(Fr.: méthode des éléments finis; Du.: Verfahren der endliche Elementen; Eng.: finite element method), werkwijze om, in het algemeen via numerieke methoden, een (benaderende) oplossing te verkrijgen voor (partiële) differentiaalvergelijkingen.
De methode is oorspronkelijk ontwikkeld voor toepassing in de mechanica, en ontstaan nadat rond 1950 computers werden ingeschakeld bij het onderzoek naar de sterkte, de stijfheid en het trillingsgedrag van constructies; tegenwoordig wordt de methode steeds meer toegepast in andere takken van wetenschap, o.a. voor de oplossing van elektromagnetische problemen. Met de elementenmethode is een instrument verkregen voor de numerieke behandeling van zeer complexe constructies of velden in zeer ingewikkelde geometrieën.
De grondgedachte is het vervangen van de conventionele analytische oplossing van (partiële) differentiaalvergelijkingen door het oplossen van stelsels algebraïsche vergelijkingen die uit deze differentiaalvergelijkingen zijn afgeleid. In de mechanica komt dit neer op een behandeling van continue constructies met behulp van methoden die voor staafconstructies al veel langer bekend zijn.
Staafconstructies zijn van nature opgebouwd uit discrete elementen: de staven. Voor continue constructies als platen, schalen enz. wordt nu ook een opdeling in elementen gemaakt. Voor de eigenschappen van de elementen is een benadering te vinden met behulp van variatieprincipes. De oplossing van het probleem is daarom uiteraard een approximatie van de exacte uitkomst.
De fundamentele onbekende in het probleem is gewoonlijk de waarde van de verplaatsingen in een discreet aantal punten (knopen) die meestal samenvallen met de hoeken van de elementen. In het geval van een statische belasting op een elastische constructie wordt de vector van alle onbekende verplaatsingen opgelost uit een stelsel algebraïsche lineaire vergelijkingen waarin de belasting voorkomt als rechterlid. Bij trillings- en stabiliteitsonderzoek moet een eigenwaardeprobleem worden opgelost. Omdat verplaatsingen de rol van onbekende spelen, spreekt men wel over de verplaatsingsmethode (displacement method). Incidenteel wordt ook met onbekende krachten gewerkt (force method).
Wiskundig gezien is de elementenmethode te rubriceren als een generalisatie van differentierekening (finite difference method). De geldende partiële differentiaalvergelijkingen worden in de elementenmethode gediscretiseerd met een differentieschema. Terwijl de constructeur in de elementen werkelijke constructiedelen ziet die worden samengevoegd tot een constructie, is een elementennet voor de wiskundige een onderverdeling van een gebied in subgebieden. De elementenmethode is daarom eveneens toepasbaar voor potentiaalproblemen bijv. bij elektromagnetische verschijnselen, warmtediffusie, consolidatie van grond enz.
Hierna wordt de uiteenzetting vooral toegespitst op het lineair-elastische gedrag van constructies onder invloed van belastingen die onafhankelijk zijn van de tijd. Volledigheidshalve wordt aan het slot kort ingegaan op trillingen, stabiliteit en niet-lineair gedrag.
Uitgangspunten.
Een constructie wordt opgevat als een lichaam van lineair-elastisch materiaal dat is onderworpen aan de invloed van uitwendige krachten. In het algemene driedimensionale geval bestaat de over het volume V verdeelde uitwendige krachtenvector p uit de drie componenten px, py en pz per volume. Elk punt (x, y, z) van het lichaam zal daardoor de corresponderende verplaatsingen ux, uy en uz ondergaan, samen de vector u vormend. De uitwendige krachten p zijn in evenwicht met de inwendige normaalspanningen σxx, σyy, σzz en schuifspanningen σxy, σyz, σzx. Deze inwendige spanningen σ gaan gepaard aan vervormingen ε, te weten de rekken εxx, εyy, εzz en de hoeken van afschuiving λxy, λyz, λzx.
De randvoorwaarden op de rand van het lichaam worden onderscheiden in twee typen: op een deel van de rand zijn de uitwendige krachten p voorgeschreven (en dus de spanningen) en op het resterende deel de verplaatsingen u. In de verplaatsingsmethode wordt gebruik gemaakt van het verband tussen ε en u en van het verband tussen σ en ε: ε = Duu; σ = Sεε.
Hierin is Du een differentiaaloperator en Sε een matrix met de materiaalconstanten uit de constitutionele vergelijkingen.
Er is nog een derde relatie die het evenwicht beschrijft tussen σ en p. Met de eerder genoemde relaties tussen enerzijds σ en ε en anderzijds ε en u worden deze evenwichtsvergelijkingen partiële differentiaalvergelijkingen in u, met differentiaties naar de coördinaten x, y, z (zie Elasticiteit: randvoorwaardeproblemen). In de elementenmethode worden deze differentiaalvergelijkingen omgevormd tot algebraïsche, meestal via variatieprincipes die equivalent zijn met de partiële differentiaalvergelijkingen. Voorbeelden daarvan zijn het principe van minimale potentiële energie en het principe van virtuele arbeid.
Hier zal het eerstgenoemde principe worden gevolgd.
De potentiële energie Epot is gedefinieerd (zie Elasticiteit: energieprincipes) als:
Epot = Eel − ∫∫∫V uTp dV
De elastische energie Eel hierin is:
Eel = ½∫∫∫ 𝜀TSεεdV
Het superscript T duidt het transponeren van een vector aan. In deze uitdrukking voor Epot wordt de uitwendige belasting p op de rand opgevat als een bijzondere vorm van de volumekrachten p in de integraal over het volume.
Verplaatsingsveld.
Om met deze uitgangspunten een benaderingsoplossing van elasticiteitsproblemen te verkrijgen, wordt een kinematisch toelaatbaar veld aangenomen voor de verplaatsingen u, waarin een eindig aantal parameters, de vrijheidsgraden v, nog moet worden opgelost. Voor het kiezen van dit verplaatsingsveld wordt het lichaam opgedeeld in elementen. Een ruimtelijke constructie kan bijv. bestaan uit tetraëderelementen en een vlakke constructie (platen, schalen) uit drie- of vierhoekige elementen. In afb. 1 is een tweetal elementen geschetst. De vrijheidsgraden worden in het algemeen gedefinieerd in knopen op de rand van een element; meestal vallen deze samen met de hoeken van de elementen maar ook knopen halverwege een zijde zijn mogelijk. De vrijheidsgraden van alle knopen die in de constructie voorkomen vormen samen de vector v. Binnen één element wordt ue\x, y, z) geïnterpoleerd op de vrijheidsgraden ve van dit element:
ue(x, y, z) = Be(x, y, z)ve
Hierin bevat ve de waarden van ue in een discreet aantal punten (de knopen) van het element. De vrijheidsgraden die in een bepaalde knoop zijn gedefinieerd gelden voor alle elementen die daar samenkomen, zodat in principe een compatibel verplaatsingsveld voor het hele lichaam is geconstrueerd uit de stuksgewijze interpolaties per element.
Berekening van verplaatsingen via het energieprincipe.
Met het verplaatsingsveld is het mogelijk de potentiële energie van de gehele constructie te schrijven als een kwadratische uitdrukking in de vrijheidsgraden v. Dit kan per element worden uitgevoerd, omdat de totale potentiële energie Epot de som is van bijdragen Eepot van de afzonderlijke elementen.
De in het verband ε = Duu voorkomende differentiaaloperator kan per element worden uitgevoerd op de interpolatiematrix Be die daarmee overgaat in de elementdifferentiematrix De. De vervormingen εe per element zijn nu uitgedrukt in de vrijheidsgraden ve:
εe = Deve
De termen in de differentiematrix zijn in beginsel weer functies van x, y en z. De elastische energie Eeel per element is dan te schrijven als:
Eeel = ½veTSeve
waarin de stijfheidsmatrix van één element
Se is:
Se = ∫∫∫Ve DeTSeεDedV
Het deel Epot,p van Epot dat het gevolg is van de uitwendige krachten wordt per element:
Eepot,p = − veTke
waarin de consistente belastingvector ke is:
ke = ∫∫∫Ve BeTpdV
Per element is de potentiële energie Eepot derhalve uitgedrukt in de vrijheidsgraden van het element:
Eepot = ½veTseve − veTke
De potentiële energie Epot van de gehele constructie kan nu worden bepaald als:
Epot = ½vTSv − vTk
Hierin heeft de vierkante matrix S de rang van v; S wordt opgebouwd als een assemblage uit alle elementmatrices Se. Evenzo wordt k samengesteld uit de vele elementvectoren ke.
De eis dat de potentiële energie Epot minimaal is, betekent dat de eerste variatie δEpot nul moet zijn voor elke willekeurige variatie δv ≠ 0 van het verplaatsingsveld:
δEpot = δvT(Sv − k) = 0
Aangezien δv ongelijk nul is moet gelden: Sv = k. Dit is het gezochte stelsel lineaire algebraïsche evenwichtsvergelijkingen waaruit de onbekende verplaatsingen v in de knopen kunnen worden opgelost.
Het stelsel vergelijkingen Sv = k is gebaseerd op een verplaatsingsveld waarvan stilzwijgend is aangenomen dat het de waarde nul aanneemt op het deel van de rand waarop de verplaatsingen zijn voorgeschreven; deze voorgeschreven waarden vu kunnen echter ook ongelijk nul zijn.
Berekening van spanningen en oplegreacties.
Als de verplaatsingen v zijn opgelost uit Sv = k kunnen de spanningen σe per element worden berekend. Het verband σe = Sεεe en εe = Deve geeft direct het recept:
σe = SeεDeve
De produktmatrix is in beginsel een functie van x, y en z. Afhankelijk van de orde van de interpolatiematrix worden de waarden van σe in een aantal punten (meestal de hoeken) van het element berekend.
Ter plaatse van de voorgeschreven verplaatsingen kunnen de oplegreacties worden berekend. Deze zijn perfect in evenwicht met de uitwendige belasting k, ongeacht de nauwkeurigheid van het gebruikte element mits is voldaan aan de consistentie-eis dat een star-lichaambeweging kan worden beschreven.
Consistentie-eisen.
De gevonden oplossing is een benadering als de gekozen interpolatie voor u de exacte oplossing niet kan beschrijven. Bij voortgaande netverfijning zal de oplossing convergeren naar de exacte oplossing als aan een aantal consistentie-eisen is voldaan. Zo moeten de verplaatsingen op de rand van een element continu zijn, moet de gekozen interpolatie alle mogelijke verplaatsingen van het element als star lichaam kunnen beschrijven en dient met de gekozen interpolatie elke homogene vervormingstoestand (constante rekken) te kunnen worden weergegeven. Wiskundig gesproken is de set benaderingsoplossingen dan compleet. Als aan de genoemde drie eisen is voldaan, zal bij toenemende netverfijning de exacte oplossing steeds beter worden benaderd in die zin dat de aanvankelijk te kleine verplaatsingen monotoon aangroeien tot de exacte waarde.
Overzicht en evaluatie.
Samengevat komt een berekening met de elementenmethode op het volgende neer. De constructie wordt in elementen opgedeeld en een verplaatsingsveld wordt geïnterpoleerd op nog nader te bepalen vrijheidsgraden v in de knopen. Per element wordt de stijfheidsmatrix Se berekend en de bijdrage ke aan het rechterlid bepaald. De belastingvector ke is de naar de knopen afgevoerde belasting. Een dergelijke belasting kan uitwendig aangrijpen of zijn opgelegd in de vorm van initiële spanningen of rekken (bijv. thermische belasting). Na assemblage van alle matrices Se tot S en alle bijdragen ke tot k wordt het rechterlid gecorrigeerd voor eventuele voorgeschreven verplaatsingen vu ongelijk nul. De verplaatsingen v worden vervolgens opgelost, waarna de spanningen σ in de elementen en de oplegreacties ku eveneens zijn te berekenen. Programma’s ter uitwerking van de elementenmethode hebben daarom gewoonlijk de structuur van afb. 2.
De uitvoer omvat in het algemeen de verplaatsingen in de knopen, de spanningen in de elementen en de oplegreacties ter plaatse van voorgeschreven verplaatsingen. Deze vormen een veld dat continu is over de elementranden; de spanningen vertonen daar in het algemeen discontinuïteiten (afb. 3).
De beschikbare programma’s zijn standaardpakketten voor een brede schakering van constructies. De ‘elementenbibliotheek’ omvat in de regel elementen voor vlakke spanningstoestanden, vlakke vormveranderingstoestanden, ruimtelijke spanningsvraagstukken, axiaalsymmetrische problemen, plaatbuiging en schalen. Ook combinaties komen voor. De vlakke elementen zijn driehoekig, of vierhoekig en de ruimtelijke elementen hebben zijvlakken van die vorm.
Elementen voor buiging zijn niet op eenvoudige wijze te maken, omdat dan moeilijk aan alle consistentie-eisen kan worden voldaan. Daarom zijn alternatieve technieken opgekomen waarbij een onderstelling voor het momentenverloop kan worden gedaan (Reissner). Ook bij buigingsproblemen in schalen hebben zulke elementen hun nut bewezen.
Trillingsonderzoek.
Het onderzoek van het trillingsgedrag vereist het medebeschouwen van de tijdafhankelijke massakrachten. Volgens het beginsel van d’Alembert kan de berekening op dezelfde wijze als voor een statische belasting worden uitgevoerd, door de dynamische sleepkrachten:
pd = − ϱü
als extra uitwendige belasting in de evenwichtsvergelijkingen te betrekken. De grootheid ϱ is de dichtheid (massa per volume) en ü bevat de drie versnellingscomponenten. In de elementenmethode worden alleen de vrijheidsgraden afhankelijk van de tijd ondersteld, zodat de interpolaties worden:
ue(x,y,z,t) = Be(x,y,z) ve(t)
Het afvoeren van de dynamische belasting naar de knopen resulteert in geconcentreerde massakrachten kd voor de constructie:
kd = − Mv ̈
waarin de massamatrix M wordt geassembleerd uit de afzonderlijke elementmatrices Me:
Me = ∫∫∫Ve BeTϱBe dVe
op dezelfde wijze als dat gebeurt voor de stijfheidsmatrix. De evenwichtsvergelijkingen voor de constructie luiden nu:
Sv + Mv ̈ = k
Het rechterlid is een opgedwongen belasting die afhankelijk is van de tijd.
De partiële differentiaalvergelijkingen voor de continue mechanica zijn aldus vervangen door een stelsel simultane gewone differentiaalvergelijkingen waarin alleen afgeleiden naar de tijd voorkomen. Door het rechterlid gelijk aan nul te nemen en v harmonisch te veronderstellen ontstaat het klassieke eigenwaardeprobleem voor de eigentrillingen. De uitvoer van programma’s kan de eigenfrequenties omvatten alsmede de eigentrillingsvormen (afb. 4). Met behulp hiervan kan vervolgens de responsie voor opgelegde periodieke (al dan niet harmonische) belastingen k worden bepaald. De oplossing wordt nauwkeuriger naarmate meer elementen worden gebruikt.
Stabiliteitsonderzoek.
Het onderzoek naar de stabiliteit van een evenwichtstoestand wordt geïllustreerd voor een in zijn vlak belaste dunne plaat. In de evenwichtstoestand heersen membraankrachten die de potentiële energie stationair maken. Er zal evenwichtsverlies door plooi optreden als de stationaire waarde niet een minimum is; daarom moet de tweede variatie van de potentiële energie worden onderzocht bij een λ-voudige belasting. Deze is kwadratisch in de plooivervormingen δv:
δE2pot = ½δvT(S − λSn)δv
De matrix S is de stijfheidsmatrix voor plaatbuiging die ook in de lineaire theorie zou optreden en Sn is een toegevoegde matrix die van de aanwezige membraankrachten afhangt en positieve termen heeft bij drukkrachten.
De aard van het evenwicht is volgens de theorie van de matrices te bestuderen aan de hand van het eigenwaardeprobleem:
(S − λSn)δv = 0
Als de kleinste eigenwaarde λ groter is dan 1, dan is δ2Epot negatief en daardoor de evenwichtstoestand instabiel. De berekende plooibelasting is nauwkeuriger naarmate meer elementen worden gebruikt. Afb. 5 geeft een indruk van de te bereiken nauwkeurigheid met een eenvoudig element.
Niet-lineair gedrag.
Met de elementenmethode kunnen niet-lineaire berekeningen worden uitgevoerd door belastingincrementen te lineariseren. Zowel niet-lineariteit van het materiaal (plasticiteit) als geometrische niet-lineariteit kan op deze wijze worden behandeld. Voor deze toepassingen zij naar de literatuur verwezen.