De bedoeling van dit artikel is tweeledig. Sommige differentiaalvergelijken zijn dermate gecompliceerd dat het uiterst lastig, zo niet onmogelijk is deze analytisch op te lossen. Daarnaast komt het voor dat men met differentiaalvergelijkingen geconfronteerd wordt terwijl men daar niet of slecht mee bekend is. Voor bovengenoemde gevallen is dit artikel bedoeld.
In het verleden heb ik de hier onder beschreven methode enkele malen met succes toegepast, het ging in dit geval om differentiaalvergelijkingen die analytisch, althans door mij, niet op te lossen waren. De onderstaande rekenmethode is alleen van toepassing op gewone differentiaalvergelijkingen, geen partieële, en de hoogste afgeleide van de 2e graad. Aan de hand van een simpel voorbeeld zullen we de methode beschrijven.
Iedereen is bekend met de bewegingsvergelijking van Newton:

F= ma

Dit geldt alleen als de massa niet in de tijd verandert, maar we zullen ons bij aardse snelheden houden.

Voor de versnelling a kunnen we ook schrijven:

a = d2y/dt2 = F/m

Beschouwen we een vrij vallend lichaam dan is a = g en krijgen we dus

d2y/dt2 = g omdat F = mg

d2y/dt2 kunnen we ook schrijven als d(dy/dt)/dt

of

d(dy/dti)dt ≈ (dy/dti+1 - dy/dti)/Δt

Hier staat dat de linkerterm ongeveer gelijk is aan de rechterterm. de afwijking wordt geringer naarmate de tijdstap,Δt kleiner genomen wordt. De consequentie hiervan is dat de rekentijd toeneemt. De index i betekent dat de hele term waar deze bij staat op het aanvangstijdstip bepaald dient te worden en de term met index i+1 op een tijdstip,Δt later.

We hadden al eerder gezien dat d(dy/dt)/dt = g

zodat

(dy/dti+1 - dy/dti)/Δt ≈ g

Bij g hoeven we geen index te gebruiken daar deze constant is (voor korte tijd in ieder geval!)

De volgende stap levert een van de af te leiden basisvergelijkingen op.

dy/dti+1 ≈ dy/dti +gΔt                                                      (1)

De term dy/dt is per definitie

dy/dti ≈ (yi+1 - yi)/Δt                                                          (2)

De indices hebben hier dezelfde betekenis als bij de boven gegeven vergelijkingen.

Hieruit volgt

y_i+1 ≈ y_i + (dy/dt_i) Δt

Deze twee afgeleide vergelijkingen vormen de basis voor de numerieke oplossingsmethode. De berekening gaat als volgt, met vgl.(2) berekent men uit y0 en dy/dt0 de waarde voor y1. Vervolgens uit dy/dt0 en g de waarde voor dy/dt1. Daarmee is de eerste rekenlus doorlopen en kunnen we deze stappen herhalen met de nieuwe gevonden waarden als uitgangspunt. We stoppen het proces zodra de gewenste eindtijd bereikt is. Het handigst is om deze berekeningen op een PC uit te voeren, maar ook op een programmeerbare rekenmachine gaat het prima, is alleen vaak wat lastiger te programmeren omdat je al gauw in conflict komt met het aantal programmastappen.

Voor het bovenstaande voorbeeld is gekozen omdat dit ook analytisch eenvoudig op te lossen is. Bekijkt men de antwoorden in beide gevallen dan blijkt er een afwijking te zijn. Dit vindt zijn oorzaak in het feit dat we er in het begin van uitgegaan zijn dat de snelheid in het eerste tijdsinterval nul was. We veranderen daarom van strategie, we gaan nu eerst met een halve tijdstap de snelheid berekenen. d.w.z. dy/dt1/2 uit dy/dt0 en g. Vervolgens wordt uit y0 en dy/dt1/2 een waarde voor y1 berekend met een hele tijdstap. De overige berekeningen worden voor zowel y als dy/dt weer met een hele tijdstap gedaan en we krijgen enigszins gewijzigde basisvergelijkingen.

Deze worden

dy/dt1/2 ≈ dy/dto + gΔt/2                                                               (3)

dy/dti+1/2 ≈ dy/dti-1/2 + gΔt                                                           (4)

en

yi+1 ≈ yi + (dy/dti+1/2)Δt                                                                 (5)

Vergelijkt men een berekening verricht met deze 3 vergelijkingen dan blijkt de overeenkomst met de analytische oplossing zeer goed te zijn.

Tot nu toe zijn we uitgegaan van een constante, n.l. g. Het blijkt echter dat de methode even zozeer bruikbaar is voor ingewikkelder vergelijkingen. Daartoe behoeven we alleen de vergelijking in een vorm te schrijven dat hij voldoet aan:

d²y/dt² = C

Hierbij is C een willekeurige functie van alle grootheden, zoals constanten, y, t, x, dy/dt, dy/dx etc. De vergelijkingen (3) t/m (5) krijgen daardoor de vorm

dy/dt1/2 ≈ dy/dto + CΔt/2                                                              (6)

dy/dti+1/2 ≈ dy/dti-1/2 + CΔt                                                          (7)

en

yi+1 ≈ yi + (dy/dti+1/2)Δt                                                                (8)

 

Het is vaak handig om een rekenschema te tekenen, een voorbeeld is gegeven in Fig.1 Dit schema is gemaakt om een volgend voorbeeld toe te lichten uit de mechanica en wel de sterkteleer. Gevraagd wordt de doorbuiging en de doorbuigingshoek te berekenen van een bladveer met rechthoekige doorsnede. Een probleem dat triviaal lijkt omdat dit te vinden is in vrijwel elk handboek. Alleen wordt daarbij stilzwijgend verondersteld dat het gaat om kleine doorbuigingen, meestal niet expliciet aangegeven zodat men behoorlijk de mist in kan gaan. We zullen dan ook zowel het vereenvoudigde geval als het nauwkeurige geval behandelen. Terwille van de overzichtelijkheid gaan we op een andere notatie over.

y" = d2y/dt2 en y' = dy/dt

 

 

Een bladveer van elastisch materiaal welke aan de wet van Hooke voldoet zal dusdanig deformeren, hetzij door zijn eigen gewicht, hetzij door een belasting, dat de kromtestraal van de elastische kromme evenredig is met het buigende moment M.

K = y"/[1 + (y')2]3/2 = M/EI

Hierin is

E     de Youngs modulus

I       het traagheidsmoment van de dwarsdoorsnede van de veer t.o.v. een

        horizontale lijn gaande door het zwaartepunt van de doorsnede

y      de ordinaat van de elastische kromme


Bovenstaande uitdrukking staat bekend onder de naam Bernoulle - Euler wet.

In de praktijk gaat men ervan uit dat de doorbuiging gering zal zijn en wordt de term (y')2 verwaarloosd, waardoor de uitdrukking overgaat in

y" = M/EI

Beschouwen we nu een bladveer met lengte 1, een gelijkmatige belasting g(x) N/m en een puntbelasting F N op het niet ingeklemde uiteinde dan vinden we voor het inwendige moment op een afstand x van de inklemming

 

M(x) = _x^l∫(ξ - x)g(ξ)dξ + F(l-x)

 

Uitwerking van de integraal levert op 1/2g(l-x)2 zodat de differentaalvergelijking geschreven kan worden als

y" = [1/2g(l-x)2+ F(l-x)]/EI

Integratie levert op

y' = [gl2x/2 - glx2/2 + gx3/6 + Flx - Fx2/2 + C1]/EI

Voor x = 0 is y' = 0 zodat C1 eveneens 0 is

Nogmaals integreren levert

y = [gl2x2/4 - glx3/6 + gx4/24 + Flx2 /2 -Fx3/6 +C2]/EI

Voor x = 0 is y = 0 zodat ook C2 = 0

Met behulp van de twee bovenstaande vergelijkingen is het dus mogelijk voor elk punt de hoek en de doorbuiging te berekenen. Meestal zijn we alleen geinteresseerd in de doorbuiging en de hoek op het uiteinde van de veer, d.w.z. op x = l

Substitutie hiervan in de vergelijkingen en vereenvoudiging levert voor de hoek

φ = gl3/6EI + Fl2/2EI

en

y = gl4/8EI + Fl3/3EI

De eerste termen rechts van het gelijkteken zijn de hoek en doorbuiging t.g.v. het eigen gewicht en de laatste termen de hoek en doorbuiging t.g.v. de puntbelasting.

Om de invloed van de verwaarlozing van de al eerder genoemde term te bepalen is m.b.v. bovenomschreven methode voor eenzelfde veer de hoek en doorbuiging berekend.

Voor C in de vergelijkingen (6) en (7) is gesteld

C = [1 + (y')2]3/2 . [{g(l-x)2/2 + F(l-x)}/EI]

Omdat een grafiek duidelijker is dan een reeks getallen zijn de uitkomsten verwerkt in een grafiek. Fig.2

Duidelijk is hierin te zijn dat bij kleine doorbuigingen de verwaarlozing van de bewuste term dit geen invloed heeft op de resultaten, bij grotere belastingen treedt er een duidelijk verschil op.

 

Ter complementering


E    =    21E10 N/m2 voor staal

I      =    bh3/12 voor een rechthoekige doorsnede b breed en h hoog


Deze gegevens zijn overigens ook te vinden in elk werktuighouwkundig handboek.

Het schrijven van een computerprogramma of het programmeren van een rekenmachine zal wellicht niet veel problemen opleveren, wel moet men erop bedacht zijn dat de computer de hoeken in radialen berekent en niet in graden. Indien men dat wenst zal men daar een aparte subroutine voor moeten schrijven.

Hopelijk vormt dit artikel een aanleiding om eens wat meer met een numerieke aanpak te proberen, analytisch is natuurlijk de snelste en meest zekere weg, maar ook de numerieke aanpak heeft zijn bekoring en verdient zeker ook aandacht.

 



Toelichting op de grafiek.

 

Bij de analytische berekening is de noemer in de Bernoulle - Euler uitdrukking niet meegenomen, daarentegen bij de numerieke berekening wel. Uit de grafiek is duidelijk te zien dat verwaarlozing bij kleine doorbuiging geen noemenswaardige afwijking oplevert (F = 0.1N) maar een zeer duidelijke bij een grotere belasting (F = 1.0N)

 

---