Dankzij de komst van computers ligt een nieuwe, derde, weg voor ons open: numerieke simulatie (ook wel computersimulatie genaamd). We proberen hierbij om de Navier-Stokes vergelijkingen met wiskundige rekenmethoden op de computer op te lossen. Het betreffende vakgebied heet numerieke stromingsleer (in het Engels: Computational Fluid Dynamics - CFD). Eigenlijk is dit vakgebied al zo’n honderd jaar geleden begonnen, toen een zaal vol ‘rekenmeisjes’ de rol van ‘menselijke’ computer vervulde[1]. Een ‘dirigent’
zorgde voor synchronisatie van de berekeningen, en tussenresultaten werden via blaadjes papier uitgewisseld. Figuur 3 toont een resultaat van moderne CFD: een berekening van een grote golf die over de voorplecht van een schip slaat. Hoe zo’n berekening in elkaar zit, wordt hieronder uitgelegd. In dit voorbeeld zijn de scheepsbouwers geïnteresseerd in de krachten die zo’n grote golf uitoefent op de scheepsconstructie; ze willen dat bovendien in een
voor hen begrijpbare vorm. Visualisatie van de bladzijden vol getallen die uit een berekening komen, is daarmee een belangrijk onderdeel van CFD.
Modelleren
Evenals de theoretische ‘potlood-en-papier’ aanpak begint computersimulatie bij het beschrijven van de stroming die we willen analyseren. En ook nu worden de Navier-Stokes vergelijkingen nog vaak vereenvoudigd om de rekentijden niet al te groot te laten worden. Deze stap, waarbij je afweegt welke fysica voor het probleem belangrijk is en welke niet, heet modelleren. Het resultaat is een wiskundig model waarin een aantal aspecten van de stroming verwaarloosd zijn. Zo wordt bijvoorbeeld het stromende medium onsamendrukbaar
verondersteld (geen schokgolven) en/of de stroperigheid (viscositeit) van de stroming wordt verwaarloosd. In een dergelijke situatie kan de stroming worden beschreven met behulp van een (snelheids)potentiaal, analoog aan de potentiaal uit de elektriciteitsleer; de stroming wordt veroorzaakt doordat de potentiaal van plaats tot plaats
verschillend is. In termen van het weer mag je ook aan de stroming ten gevolge van drukverschillen denken, waarbij je aanneemt dat de lucht
rechtstreeks van hoge druk naar lage druk stroomt. Met deze aanname verwaarloos je de invloed van de draaiing van de aarde, die beschreven wordt door de bekende wet van Buys-Ballot: op het noordelijk halfrond krijgt de stroming een afwijking naar rechts. Dit is nou een voorbeeld van modelleren, in dit geval bedoeld om de formules hieronder niet al te moeilijk te laten worden.
Discretiseren
Vervolgens wordt het wiskundig model gediscretiseerd op een rekenrooster. Het hele stromingsgebied wordt daarbij opgedeeld in kleine blokjes; op ieder van deze blokjes worden de behoudswetten uit het wiskundig model toegepast.
voor een Fokker 50. Om de discretisatie verder uit te leggen wordt voor de eenvoud het rooster twee- dimensionaal en rechthoekig verondersteld, met gelijke maaswijdtes h in beide richtingen; zie figuur 5. Voor iedere roostercel wordt nu discreet een behoudswet opgesteld. Nb. De snelheidscomponenten in horizontale en verticale richting geven we aan met u en v, de potentiaal met P, en met ‘kompasindices’ worden de buurcellen aangeduid. Massabehoud in een roostercel houdt in dat de totale doorstroming door alle vier zijwanden op nul moet uitkomen. De doorstroming door bijvoorbeeld de rechterwand van de gearceerde cel C kan worden benaderd door de snelheid loodrecht op deze wand te vermenigvuldigen met de lengte van de wand, en wordt daarmee u hO . Sommatie over alle vier de
zijwanden levert nu de discrete behoudswet:
u h v h u h v hO + N − W − Z =0 (1)
De snelheid loodrecht op de zijwanden van de cel wordt benaderend gevonden door het potentiaal- verschil tussen de aangrenzende cellen te delen door de afstand (dus in feite door te kijken naar de helling van de potentiaal). Op de rechterzijwand leidt dit tot
u P P h
O= O− C
Voor de andere drie zijwanden gelden soortgelijke formules. Dit ingevuld in (1) geeft:
PO+PN+PW+PZ−4PC=0 (2) Dit is de klassieke ‘vijfpuntsformule’ voor de potentiaalvergelijking, ook wel vergelijking van Laplace genoemd. Op deze wijze ontstaat voor elke cel één discrete vergelijking. Na toevoeging van randvoorwaarden kan de potentiaal in de cellen worden gevonden door het volledige stelsel vergelijkingen op te lossen. Bij discretisatie zie je
FIGUUR 2 De Navier-Stokes vergelijkingen voor onsamendrukbare stroming
FIGUUR 3 Berekening van een grote golf die over een schip slaat.
2 4 0
euclides nr.4 / 2005
oplossing bevat daarmee een fout, de zogenaamde discretisatiefout. Deze fout kan worden verkleind door verfijning van het rooster. In ons voorbeeld leidt halveren van de maaswijdte tot een vier keer kleinere fout, maar ook tot een vier keer groter stelsel vergelijkingen, wat meer rekeninspanning vereist om dit stelsel op te lossen (bij de hieronder beschreven oplosmethoden zelfs 16 keer meer). Een afweging tussen nauwkeurigheid en rekeninspanning is hier vereist.
Itereren
Het aldus ontstane stelsel, het numeriek model, is vaak erg groot. Een rooster met twee miljoen cellen (zoals tegenwoordig gebruikt voor het oplossen van de Navier-Stokes vergelijkingen) leidt tot een stelsel met tien miljoen vergelijkingen. Het oplossen hiervan is erg ‘duur’, en vereist de slimheid van de numericus en de snelheid van de computer. Het oplossen gebeurt meestal iteratief, waarbij altijd eerst een schatting voor de oplossing, zeg P P= (oud), wordt gemaakt. De eenvoudigste iteratieve methode ontstaat daarna door uit (2) een verbeterde waarde voor PC te bepalen via
PC(nieuw)=1 PO(oud)+PN(oud)+PW(oud)+PZ(oud
4( ))) (3)
Dit proces, al bekend uit het begin van de 19e eeuw en genoemd naar Jacobi, kan herhaald worden tot het verschil met de gewenste oplossing voldoend klein is geworden.
Deze methode is weliswaar eenvoudig te
programmeren, maar niet erg snel. Op eenvoudige wijze is een verbetering te maken door in het iteratievoorschrift (3) meteen gebruik te maken van de nieuwe waarden zodra die er zijn. Als we bijvoorbeeld linksonder in het rooster beginnen, dan zijn de waarden voor PW en PZ al vernieuwd op het moment dat we aankomen in cel C, en kunnen we het voorschrift aanpassen tot
PC(nieuw)=1 PO(oud)+PN(oud)+PW(nieuw)+PZ(n
4( iieuw)) (4)
Met de theorie achter de numerieke wiskunde kan bewezen worden dat deze aanpassing, de methode van Gauss-Seidel, twee keeer sneller tot resultaat leidt dan de methode van Jacobi.
Valideren
De resultaten worden zichtbaar gemaakt in voor de ontwerper te interpreteren vorm, zoals de drukverdeling rond een vleugelprofiel; zie figuur 6. De druk aan boven- en onderkant van de vleugel is weergegeven; het oppervlak tussen de twee krommen is de draagkracht die de vleugel ondervindt. Voor het bepalen van de waarde van de resultaten (validatie) worden ze vaak vergeleken met experimenten. In de figuur staat een vergelijking van een berekende drukverdeling met een windtunnelmeting (linksonder in figuur 6). De krommen zijn duidelijk verschillend; welke is de beste? Van de berekening weten we al dat die slechts bij benadering goed is, omdat we de Navier-Stokes vergelijkingen vereenvoudigd hebben. Maar ook de windtunnelmeting is niet correct, omdat de wanden van de tunnel de stroming verstoord hebben en vanwege schaaleffecten (we hebben het vliegtuig kleiner gemaakt, maar de luchtmoleculen niet). Hiervan hebben we geen last als de metingen in vrije vlucht uitgevoerd worden. Nu zien we dat de berekening nauwkeuriger is dan de windtunnelmeting, en goed overeenkomt met de werkelijkheid.