Drukberekening in een axisymmetrisch model van het hoofd met de eindige elementen methode

Drukberekening in een axisymmetrisch model van het hoofd

met de eindige elementen methode

Citation for published version (APA):

Günther, O. (1993). Drukberekening in een axisymmetrisch model van het hoofd met de eindige elementen methode. (DCT rapporten; Vol. 1993.059). Technische Universiteit Eindhoven.

Document status and date: Gepubliceerd: 01/01/1993

Document Version:

Uitgevers PDF, ook bekend als Version of Record

Please check the document version of this publication:

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.

• The final author version and the galley proof are versions of the publication after peer review.

• The final published version features the final layout of the paper including the volume, issue and page numbers.

Link to publication

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain

• You may freely distribute the URL identifying the publication in the public portal.

If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:

www.tue.nl/taverne

Take down policy

If you believe that this document breaches copyright please contact us at:

openaccess@tue.nl

Drukberekeningen in een

axisymmetrisch model van het

hoofd met de

eindige elementen methode

Otwin Günther

Rapport WFW: 93.059

Technische Universiteit Eindhoven Faculteit Werktuigbouwkunde

Vakgroep Fundamentele Werktuigkunde Begeleiders

Fons Sauren Maurice Claessens Juni 1993.

Samenvatting

In dit onderzoek is het gebruik van het eindige elementen pakket DIANA getoetst voor de modellering van het menselijk hoofd. Dat is gebeurd aan twee axisymmetrische modellen

w a ~ n emm wcrd uitgegaan dat het hoofd uit een schedel bestaat met daarin de hersenen. De schedel weïd als een lineair elastische vaste stof gemodelleerd. De hersenen werden in het eerste geval gemodelleerd als een lineair elastische vaste stof die door de meegegeven materiaalparameters een vloeistofkarakter kreeg. In het tweede geval is er geprobeerd werkelijk een vloeistof van te maken. De vaste stof modellen werden met verschillende materiaalparameters doorgerekend voor tevens twee verschillende soorten belastingen. Deze twee belastingen waren achtereenvolgens een blokvormige drukpuls en een driehoe- kige drukpuls.

De specifieke interesse ging uit naar de druk op verschillende plaatsen in het model als functie van de tijd. Het bij dit onderzoek gevonden verband tussen druk en tijd kon vergeleken worden met resultaten uit de literatuur. Op grond van de verkregen resultaten

kon de waarde van de modellering in DIANA vastgesteld worden. Hieruit is gebleken dat

de vloeistof modellering niet mogelijk is in de huidige vorm van het gebruikte pakket. De modellering als vaste stof, die op een vloeistof leek leverde echter hoopvolle resultaten.

Frequent gebruikte symbolen

Frequent gebruikte symbolen

P K G C E V CF

4c

E h R U F

P

Q> E

z

J B t z E -c V

Q

p'

%

dichtheid compressiemodulus glijding smodulus geluidssnelheid modulus van Young constante van Poisson Cauchy spanningstensor constante 4" orde tensor Green-Lagrange rektensor w armtegeleidingscoëfficiënt rotatietensor verlengingstensor deformatietensor 2" Piola-Kirchhoffspanningstensor

specifieke elastische energie rek

rekrichting invariant

linkse Cauchy-Green rektensor tijd tijd lineaire rektensor gradiënt operator volumebelasting randbelasting positievector

massamatrix vaste stof dempingsmatrix vaste stof stijfheidsmatrix vaste stof kolom met verplaatsingen vloeistofkrachten op vaste stof

externe krachten op vaste stof druk dichtheid vloeistof oppervlak massamatrix vloeistof dempingsmatrix vloeistof stijfheidsmatrix vloeistof puntkracht massa 1

Inhoudsopgave

Inhoudsopgave

Hoofdstuk 1

Hoofdstuk 2

Inleiding

. . .

1

Beschrijving van de modellen van Merchant en Ruan

2.1 Modelbeschrijving van Merchant

. . .

.

. .

.

. . .

.

. . .

2

2.2 Modelbeschrijving van Ruan

. . .

.

. . .

.

. . .

.

. . .

.

. . .

.

. . .

4 Modellering

3.1 Vaste stof modellering

.

. . .

.

. .

. .

. .

. .

. . .

.

. .

.

. . . .

. .

. . .

.

.

5

3.2 Vloeistofmodellering

. . .

. .

. . .

.

. . .

.

. .

. .

. .

.

. . .

.

. . .

7 Resultaten van de berekeningen

4.1 Blokvormige drukpuls

. .

.

. . .

.

. . .

.

. .

.

. . .

.

. . .

.

. .

9 4.2 Driehoekige drukpuls

.

. . .

.

. . .

.

. . .

.

. . .

.

. . .

.

. . .

12

Bespreking

5.1 Bespreking van de resultaten uit de blokvormige

drukpuls

. . .

16 5.2 Bespreking van de resultaten uit de driehoekige

17 Hoofdstuk 3 Hoofdstuk 4 Hoofdstuk 5 drukpuls

. . .

Hoofdstuk 6 Conclusies en aanbevelingen 6.1 Conclusies

. . .

. .

. .

. .

.

.

. .

. . .

.

. .

. .

. .

.

.

. . .

.

. .

. .

. .

.

.

. . . .

. .

.

. .

18 6.2 Aanbevelingen

. . .

.

. . .

.

. . .

.

. .

.

. . . .

. .

. .

.

18 Referenties

. . .

.

. .

.

. . .

.

. . .

.

. . .

.

. . .

.

. . .

.

.

20 Appendix A Materiaalmodel

. .

.

. . .

.

. . .

.

. .

.

. . .

.

. . .

.

. . .

A.l Appendix B Analysemethode

.

. .

.

. . .

.

. .

.

. . .

.

. . .

.

. .

.

. . .

.

. .

. .

. .

.

. .

. . . .

.

. .

.

B.1 Appendix C

Hoofdstuk

1

Binnen de Technische Universiteit Eindhoven is een onderzoek opgestart is dat tot doel heeft de ontwikkeling van een drie dimensionaal continuümsmodel van het hoofd waarmee het gedrag ten gevolge van stootbelastingen bestudeerd kan worden. Het uiteindelijke doel is trachten het letsel te voorspellen ten gevolge van extreme belastingen. Dit onderzoek moet de mogelijkheid bieden om concrete uitspraken te doen over de gevoeligheid van de modelresultaten voor de mate van detaillering van deelstructuren, het belang van de constitutieve eigenschappen van deze deelstructuren alsmede de interactie tussen de structuren. Ook moeten er uitspraken gedaan kunnen worden over het belang van de

randvoorwaarden ter plaatse van de overgang hoofd-nek.

Dit onderzoek heeft tot doel om inzicht te krijgen in de mogelijkheden en beperkingen van het eindige elementen pakket DIANA voor het doorrekenen van twee in de literatuur

gepubliceerde eenvoudige rotatiesymmetrische modellen van het menselijk hoofd (Mer- chant et al. [9]; Ruan et al. [12]). In beide modellen wordt de schedel voorgesteld als een gesloten dunwandige bol met homogeen, isotroop lineair elastisch materiaalgedrag.

Merchant et al. [9] modelleren de schedelinhoud als een incompressibele vloeistof. Zij gebruikten als belasting een blokvormige drukpuls. De responsie werd in hun onderzoek bepaald met behulp van de eindige differentie methode.

Ruan et al. [12] representeerden de schedelinhoud als een homogene, isotrope niet visceuze vloeistof. Zij gebruikten de eindige elementen methode om de responsie op een driehoekige drukpuls te bepalen. In beide modellen werden geen kinematische randvoor- waarden, anders dan die voortvloeien uit het axisymmetrisch karakter, opgelegd.

Bedrijving van de modellen van Merchant en Ruan

Hoofdstuk 2

Beschrijving van de modellen van Merchant en Ruan

2.1 Modelbeschrijving van Merchant

Het model van Merchant et al. [9] is grotendeels gebaseerd op een model van Engin et al. [4]. Deze was van mening dat het menselijke hoofd kon worden opgevat als een bol gevuld met vloeistof, welke incompressibel werd verondersteld. Daarnaast werden de bol en de vloeistof homogeen en isotroop verondersteld. De fysische aannamen die hierbij gemaakt werden zijn dat de schedel een gesloten bol en ideaal sferisch is. In werkelijkheid heeft de schedel heeft een grote opening in het achterhoofd, het foramen magnum,

waardoor de medulla oblongata samengaat met de spinale baan. Andere gaten zoals de optische foramina zijn verwaarloosbaar. Verder is het hoofd niet een vrij opgehangen object zoals in dit model, maar met spieren, pezen en banden via de nek aan de romp verbonden.

De geometrie van het hoofdmodel werd gebaseerd op het gemiddelde van een hoofd van een twintigjarige van onbekend geslacht. De diameter bedraagt in dat geval 16 cm, en de dikte van de schedel is dan 0,38 cm.

De fysische eigenschappen van de inhoud, die de hersenen moet voorstellen, lijken op een incompressibele vloeistof en benaderen de materiaaleigenschappen van water. Als

eigenschappen voor de schedel werden aangenomen: -dichtheid ~ = 2 1 4 0 kg*m"

-compressiemodulus K=9,23* lo9 Pa -glijdingsmodulus G=5,53*109 Pa -geluidssnelheid c=2785 m*s-'

Voor de vloeistof werden de volgende waarden aangenomen: -dichtheid Q=lOOO kg*m-3

-compressiemodulus K=1,96*109 Pa -glijdingsmodulus G=O Pa

-geluidssnelheid c=1400 m*s-'

Opmerkelijk hierin is het opgeven van waarden voor de compressie-en glijdingsmodulus. Meestal worden deze waarden opgegeven in hun equivalenten, namelijk de modulus van Young en de constante van Poisson. Ze zijn op de hier beschreven wijze echter aangege- ven omdat een modulus van Young met de waarde nul nietszeggend is. Het gebruikte

Beschrijving van de modellen van Merchant en Ruan

materiaalmodel voor dit onderzoek was de wet van Hooke: volkomen elastisch, homogeen en isotroop materiaalgedrag.

De analyse die Merchant et al [9] op het beschreven model uitvoerden was gebaseerd op de belasting van Engin [3], die een impulsbelasting gebruikte. Merchant benaderde deze impuls door een blokvormige drukpuls met een amplitude van 3,76 MPa te zetten op de schedel onder een tophoek van 15" (fig. 2.1).

fig. 2.1: Model van Merchant.

De belasting was gedurende 0,0614 ms aanwezig. Er waren geen randvoorwaarden zodat het model vrij kon bewegen. De analyse werd uitgevoerd met de eindige differentiemetho- de.

Het hoofd werd in elementen verdeeld met de hulp van hyperboolbeschrijvingen en ellipsen. De basislijn van het geometrische model is de symmetrieas, zodat de element doorsneden dwarsdoorsneden zijn van toroïden (fig. 2.1). De twee aan de buitenkant

gelegen elementnjen vormen de schedel. De schedel bestaat uit vierentwintig bij twee elementen. Dit is schijnbaar gedaan om de belasting op één element te laten werken. Dat volgt namelijk uit de volgende deling: 180"/24=7.5".

Er dient bij dit model echter een opmerking geplaatst te worden. Merchant geeft niet duidelijk aan hoe hij in zijn model de interactie tussen schedel en schedelinhoud heeft gemodelleerd.

Beschrijving van de modellen van Merchant en Rum

2.2 Modelbeschrijving van Ruan

Het tweede model dat gebruikt is, is dat van Ruan et al. [12]. Dit model is achteraf gekozen om het in dit onderzoek ontwikkelde model te valideren vanwege de slechte overkomsten met het model van Merchant. Dit model is ook axisymmetrisch en -bestaat uit een gesloten bol met daarin een niet visceuze vloeisto€. De wand is gemodelleerd als een dunne laag en de eigenschappen van dit materiaal benaderen die van craniaal bot. De vloeistof werd compressibel verondersteld met de dichtheid en compressiemodulus van water. Dit geeft het volgende voor de materiaaleigenschappen:

voor de schedel: -dichtheid

-modulus van Young E=6,5*109 Pa -constante van Poisson v=0,22 en voor de hersenen:

-dichtheid ~=1040,0 kg*m-3 -modulus van Young E=6,67*104 Pa -constante van Poisson ~ = 0 , 4 8

De dimensies van het model zijn gebaseerd op een 50 percentiel menselijk hoofd. De impactanalyse werd gedaan met het eindige elementenpakket ANSYS. De belasting,

onder een tophoek van 25", bestond in dit geval uit een driehoekige puls met een

maximum van 0,21 MPa en een duur van 10 ms. Omdat er geen kinematische randvoor- waarden waren opgelegd, was het model vrij om te transleren onder de gegeven belasting

e

= 14 1 2,O kg*m-3

De elementengeometrie die hierbij gebruikt werd was als volgt gecre- eerd. De schedel werd gemodelleerd met één laag elementen. De elemen- tenverdeling van de inhoud werd gemaakt door ellipsen en parabool- fig. 2.2: Model van Ruan et al.

Ruan et al. hebben nagelaten om de afmetingen mee te geven in hun modelbeschrijving. Ook bij deze modellering moet een kanttekening gemaakt worden, want

Modellering

Hoofdstuk 3

Modellering

3.1 Vaste stof modellering

Het model dat voor dit onderzoek werd gebruikt is eveneens een bol met daarin een op een vloeistof lijkende vaste stof (fig.

3.1). Er worden geen kinematische beper-

kingen aan het model opgelegd. De geometrie van dit axisymmeû-ische model is volledig gelijk aan het model van Merchant et al. [9]. De buitendiameter van

vaste stof

het hoofd is gelijk aan 16 cm en de dikte

van de schedel is gelijk aan o,38 cm.

fig. 3.1: Vereenvoudigd model van de sche-

del en de hersenen.

Voor zowel de schedel als de hersenen geldt de wet van Hooke:

U =4C:E (3.1)

Het materiaal is homogeen en isotroop. Een beschrijving van het materiaalmodel is weergegeven in appendix A.

De materiaaleigenschappen zijn voor de schedel hetzelfde als die van Merchant et al. [9]

en luiden:

-dichtheid e=2140 l~g*rn-~ -modulus van Young E=1,38283*10" Pa

-constante van Poisson v=0,25

Modellering

Voor de schedelinhoud ontstaan voor de door Merchant gebruikte materiaalpaameters enige problemen omdat de modulus van Young en de constante van Poisson als volgt afhangen van de compressiemodulus en de glijdingsmodulus:

3K-2G V =

2(3K+G) (3.3)

Dit zou voor de hersenen het volgende resultaat geven: -modulus van Young E=O Pa

-constante van Poisson ~=0,50

Deze waarden zijn onmogelijk omdat door deze getallen numeriek niet-oplosbare proble- men ontstaan. Daarom is de keuze gevallen op de volgende waarden:

-dichtheid

e=

1000 kg*m-3 -modulus van Young E=i03 Pa -constante van Poisson ~ = 0 , 4 9 -geluidssnelheid c=4,14 m*B'

Een lage modulus van Young komt overeen met een lage glijdingsmodulus wat weer representatief is voor een vloeistof als water waar de hersenen vaak mee vergeleken worden. De geometrie die gecreëerd is is te vinden in figuur 3.2.

Hierover valt op te merken dat de twee buitenste rijen de schedel voorstellen. Deze elementen indeling vertoont grote gelijkenissen met die van Merchant et al. [9]. De verschillen zitten vooral in de buurt van het middelpunt en de symmetrieas, die aange- bracht is om de rekeninspanning te verkleinen. Problemen die verwacht zouden kunnen worden zijn de twee lijnstukken binnen het vier-knoopselement die in elkaars verlengde liggen (zie pijltjes in fig. 3.2).

Bij de modellering in DIANA zijn rotatiesymmetrische, isoparametrische vier-hoops ele- menten gebruikt.

3.2 Vloeistofmodellering

Er is ook sprake geweest van een daadwerkelijke vloeistof- modellering. De geometrie en de elementen waren exact geiijk aan de vaste stof mocieiiering. Het enige verschii van dit model met het vaste stof model is dat er in DIANA een andere materiaalgegevens ingegeven moesten worden. Voor een vloeistof vereist de invoer de waarde voor de geluids- snelheid in het medium. De gekozen waarde komt overeen met de waarden van Merchant et al. [9].

-geluidssnelheid c=1400 m*s-'

Daarnaast moesten de knooppunten op de interface tussen de vaste stof en de vloeistof dubbel gemodelleerd worden. De eerste keer als onderdeel van de bol en de tweede maal fig* 3.2: Eindige

model van het hoofd. als onderdeel van de schedelinhoud. Tussen deze twee sets knooppunten werden de vloeistof-vaste stof interface ele- menten (zogenaamde FLUSTR elementen) geplaatst, die maar één dimensie hebben en wel langs dat oppervlak.

De analyse die hierop werd uit gevoerd zou worden is dezelfde als de analyse op het vaste stof model. Een blokvormige dnikpuls met een amplitude van 3760 kPa en een tijdsduur

van 0,0614 ms. De analyse kon echter niet uitgevoerd worden om de volgende redenen.

Binnen het eindige elementen pakket DIANA bestaat de mogelijkheid om aan vloeistof vaste stof structuren te rekenen. De analyse vindt plaats met de volgende methode.

Voor de vaste stof luidt de impulsvergelijking.

Msü

+CsU

+Ksu

+&=f,""'<t) (3.4)

Voor de vloeistof luidt de golfiergelijking.

V 2 p L p 1

C

De vergelijkingen (3.4) en (3.5) worden dan met de hulp van randvoorwaarden omgeschre- ven tot het volgende stelsel vergelijkingen (DIANA [2]):

(3.6)

Dit stelsel kan worden opgelost door de volgende uitdrukkingen in te vullen. Tevens vindt hier de overgang naar het frequentiedomein plaats:

Uit vergelijking (3.6) kan nu de druk p geëlimineerd worden zodat er een vergelijking

overblijft waarin alleen de verplaatsing u van de vaste stof een onbekende is. Van de vloeistof komen alleen de materiaaleigenschappen er nog in voor.

Andere redenen zijn dat de berekeningen plaatsvinden in het frequentiedomein en dat de berekeningen plaatsvinden voor kleine trillingen rond een evenwichtsstand.

Resultaten van de berekeningen

Hoofdstuk 4

4 w O 3500 3000 2500 $. 2000- 3 B 1500 1000 500 O

Resultaten van de berekeningen

- - - - - - O 0.02 0.04 0.06 0.08 0.1 0.12 O14 0.16 4.1 Blokvormige drukpuls

fig. 4.1: Blokvormige drukpuls.

De eerste impactanalyse werd uitgevoerd met een blokvormige drukpuls zoals geschetst in figuur 4.1 met het eindige elementen pakket DIANA, dat wanneer het materiaalmodel niet expliciet wordt gegeven rekent met de wet van Hooke. De simulatietijd was 0.16 ms. Nadere informatie over de analyse staat in appendix B. De impact belasting had een constante amplitude van 3760 kPa en een tijdsduur van 0.0614 ms. Merchant koos deze waarden om de impulsbelasting te benader- den die Engin [3] had gebruikt. Tijdens het

actief zijn van de drukpuls kan een drukgolf tweederde van de afstand langs de symme- trieas afleggen. Door de afwezigheid van kinematische randvoorwaarden kon het model vrij transleren langs de symmetrieas. De druk werd op verschillende plaatsen in het model berekend. Op drie plaatsen in de schedel: direct achter de plaats waar de belasting werd aangebracht, in het daar tegenoverliggende punt en tenslotte in een punt halverwege de schedelomtrek. Naast de berekeningen in de schedel werd de

druk

ook nog berekend in de hersenen. Ook hier weer op drie plaatsen: vlak achter de zone waar de belasting aangrijpt, aan de tegenoverliggende zijde en in het midden van de hersenen. De punten zijn

aangegeven in figuur 4.2, waarin ook de belasting staat aangegeven.

Door het model op te vatten als een star lichaam kan met de tweede wet van Newton een schatting worden verkregen van de maximale afstand waarover het systeem transleert tijdens de simulatietijd. De werkelijk optredende verplaatsing zal kleiner zijn omdat een gedeelte van de energie wordt opgeslagen als elastische (verv0rmings)energie.

1 1 m 2

Resultaten van de berekeningen

O

fig. 4.2: De elementen waar drukberekeningen plaatsvin- den.

weergegeven in de figuren 4.3

Het resultaat hiervan is dat de verplaatsing van de achterzij- de na 0.0614 ms kleiner moet zijn dan 0.987*10-6 m. Ook

is het mogelijk om het model te toetsen aan de voortplan- tingssnelheid van de drukpuls. Daar het geluid een voort- plantingssnelheid heeft, kan de drukpuls niet overal gelijk- tijdig geregistreerd worden. De looptijd kan eenvoudig bepaald worden door het quotiënt te nemen van de af te leggen afstand en de geluidssnelheid in het medium. Dit levert voor de punten 2, O , 1 de volgende aankomsttijden van de drukpuls. Voor punt 2 is dat O ms, voor punt O

0.045 ms en voor punt 1 0.09 ms.

Naast de gegeven getalswaarden moet ook het coup-contre- coup verschijnsel optreden. Dat houdt in dat door de druk- belasting aan de voorzijde er een positieve druk onstaat in de vloeistof aan de voorzijde en een negatieve druk aan de achterzijde. De resultaten van de numerieke berekeningen, de druk-tijd responsies in de verschillende elementen staan tot en met 4.8.

Zij geven het drukverloop weer gedurende 0,15 ms waarvan tijdens de eerste 0,0614 ms

een belasting aanwezig is.

Het dnilorerlmp in punt 2

I

tijd [ms]

fig. 4.3: Drukverloop in punt 2.

tijd [ms]

fig. 4.4: Drukverloop in punt O.

Voor het element dat direct achter het impactgebied ligt, punt 2, is het drukverloop te zien

in figuur 4.3. De drukgolf is hierin direct te zien omdat de afstand tot de impact bijna nul

is. Opmerkelijk is echter het fluctuerende karakter van de respons. Verder is te zien dat tijdens de aanwezigheid van de belasting de druk blijft dalen en dat na het wegne-

Resuïîaten van de berekenkgen 3Ooo 1m O

-

B

-ux)o -3000 -4oM) ... ... ... ... ... ... _ i + ... ... ... ... ...

;

... _ + < 4 O 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 O 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16

Het dnikveriwp in punt 3

O 0.005 -o w5 ... ... ... ... ... o1 ... i ...

-

-..o,ois f ... ... ... 2 _,,,_,,__,.,.~_,,i i... ... ... o,o25

*

, 403 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 -0.035 tijd [ms]

waarop dit gebeurt stemt niet overeen met het tijdstip dat berekend werd voor het in de buurt gelegen schedelelement. De absolute maximale druk die hier echter berekend wordt

-

_{is erg laag namelijk nog maar zo’n 0,03}

kPa. Figuur 4.6 zou overeenkomsten moe- ten vertonen met het onderzoek van Me- chant et al. [9] (fig. 4.9). Deze overeen- komst is echter ver te zoeken. Ook in deze figuur neemt de

druk

toe gedurende de tijd dat de belasting aanwezig is. Zodra deze is verwijderd daalt de

druk.

Omdat dit ele- ment dichtbij het impactgebied ligt is er geen aanloopstuk te zien voor de drukgolf. tueerde punt. Aan deze karakteristiek valt weinig te zien omdat de drukken zeer laag fig. 4.9: Drukverloop uit onderzoek van van het in het midden van het hoofd gesi- Merchant et al.

zijn. Er is geen sprake van een drukgolf omdat de drukpuls nog niet gearriveerd kan zijn. De drukgolf kan namelijk pas na 19 ms gesignaleerd worden. Naast de berekeningen van de drukken zijn er ook berekeningen voor de verplaatsingen en rekken uitgevoerd. De grafische voorstellingen hiervan zijn terug te vinden in appendix C .

2 o

1:

n

t---

- 2 0

o 0 2 O 06 o1 O 14 In figuur 4.7 is het drukverloop te zien

Time. mcec In I O n 3 u) u) L -10 n 250 200 - - - o 2 4 6 8 10 12 14 16 4.2 Driehoekige drukpuls

met een driehoekige drukpuls. Een grafische voorstelling van deze belasting is te zien in figuur 4.10 en 4.1 1. De impact heeft een duur van 10 ms en loopt de eerste helft van de tijd op naar 218 kPa om daarna af te nemen naar nul. De analyse werd uitgevoerd met dezelfde materiaalgegevens zoals die door Ruan et al. [12] werden opgegeven. Het

er kan natuurlijk wel een soortgelijke trend in de drukvariatie verwacht worden. Ook aan dit model zijn verder geen kinematische randvoorwaarden opgelegd behalve dat op t=O de snelheid en de verplaatsing nul zijn.

De druk werd niet zoals bij de vorige analyse OP zes plaatsen berekend, maar slechts op

vier plaatsen. Bij Ruan ging het alleen om de drukberekening in de elementen in het im- pactgebied en aan de tegenoverliggende zijde. Deze drukberekening is uitgevoerd voor twee schedelelementen en twee vloeistofelementen in de nabijheid van de hierboven beschreven plaatsen.

De twee uitkomsten voor de vloeistofelementen kunnen vergeleken worden met de resultaten van Ruan et al. [12].

De resultaten van Ruan zijn te vinden in figuur 4.16. De

elementen waar de druk voor is uitgerekend staan in figuur

4.1 1. Er kunnen afschattingen gemaakt worden voor de ver-

plaatsing als star lichaam door de tweede wet van Newton toe te passen. De werkelijk optredende verplaatsing zal kleiner moeten zijn omdat een deel van de energie wordt opgeslagen als elastische (verv0rmings)energie.

F

6m

u=-t3 b' OSt55.0m~

nijde(contrecoup) Het resultaat hiervan is dat na 5 ms de verplaatsing kleiner

moet zijn dan 3.3 mm. De tijden van de aankomst van de

fig. 4.11: De elementen

waar drukberekeningen drukpuls kunnen slechts ruwweg gecontroleerd worden plaatsvinden. omdat de simulatietijd veel groter is dan in de voorgaande

analyse (0.15 ms en 15

ms).

Deze tijd is voor de achterzij- de ongeveer 0.12 ms. De resultaten voor de verschillende elementen staan weergegeven in

de figuren 4.12 tot en met 4.15.

Opmerkelijk aan figuur 4.15, welke de respons weergeeft voor het element in de schedel

aan de contrecoup zijde is het sterke fluctuerende gedrag. Vooral de amplitude van de fluctuatie neemt sterk toe nadat het punt van de maximale belasting bereikt is. Het signaal heeft een frequentie van ongeveer 4 kHz. Dit is vrij hoog. De vraag rijst dan ook of dit

werkelijk zo is of dat er sprake is van numerieke problemen binnen het algorithme. Figuur 4.14 met het drukverloop van het element in de vloeistofstructuur, punt 1, bij de

contrecoup zijde laat niet zo'n sterk fluctuerend gedrag zien. Wel is in de buurt van t=O te zien dat er een gedurende een aantal honderdsten van milliseconden nog geen drukopbouw

dat het drukverloop in het element een fase achterstand heeft op de aangebrachte belasting. Dit blijkt uit het feit dat het minimum in fig. 4.14 later komt dan het maximum in fig.

4.13.

tijd [ms]

fig. 4.12: Drukverloop in punt 3.

tijd [ms]

fig. 4.13: Drukverloop in punt 2.

De drukresponsies van zowel het schedelelement als het vloeistofelement aan de coup zijde zijn weergegeven in de figuren 4.12 en 4.13.

Het druk verlmp m punt 1 Het dnik vcrlwp in punt O 16 14 12 10 8

-

- I ,

2 4 2 O -2 -4 O 2 4 6 8 10 12 14 16 O 2 4 6 8 10 12 14 16 tIJd [ml tild [ml

fig. 4.14: Drukverloop in punt 1. fig 4.15: Drukverloop in punt O. Voor het schedelelement (punt 3) ligt het absolute maximum ook in dezelfde orde van grootte als de belasting. Ook zijn hier geen grote fluctuaties te zien. De

dnik

in het schedelelement (fig. 4.12) heeft bijna hetzelfde verloop, hetzij gespiegeld, als de belas-

tingsdruk. Enige verschillen zijn het niet zo gladde verloop van de drukvariaties en de absolute maximum waarde ligt iets hoger dan bij de aangebrachte belasting. Daarentegen

Resultaten van de berekeningen

is er geen sprake van faseachterstand, want de maxima verschijnen op hetzelfde tijdstip. De grafieken met de verplaatsingen zijn opgenomen in appendix C .

PresCure De overeenkomsten tussen de figuren 4.13

en 4.14 enerzijds en figuur 4.16 anderzijds zijn het overeenkomende gedrag voigens het coup-contrecoupverschijnsel. Am de

ene zijde een positief drukverloop en aan de andere zijde een negatief drukverloop. De orde van grootte stemt ook in redelijke mate overeen. Het opvallendste verschil is het sterke fluctuerende verloop dat Ruan

( F ? L ?O00 I600 I200 800 boo 0.0 -LOO

-800 et al. gevonden hebben.

-1200 - 1600

0 1 2 3 4 5 6 1 8 9 1 0

Time(ms)

fig. 4.16 : Drukverloop uit onderzoek van Ruan et al.

Hoofdstuk

5

Bespreking

5.1 Bespreking van de resultaten aait de blokvormige drukpuls

Uit de figuren van de rekken (C.13 en C.14) blijkt dat ze voor zowel de Green-Lagrange

rek als de lineaire rek steeds in dezelfde orde van grootte zijn. De ordegrootte is ongeveer enkele duizendsten. De lineaire theorie is dus een juiste aanname geweest.

Opvallend is de enorme lage

druk

in het schedelelement (fig. 4.3) in de nabijheid van de

impact; ruim twee maal zo groot in absolute zin als de impact zelf. Dit is onverklaarbaar want juist op deze plaats wordt compressie verwacht. Verder is opvallend aan dit element het sterk fluctuerende gedrag van de

druk.

Het signaal heeft een frequentie in de buurt van

300 kHz. Dit is niet direct verklaarbaar. Er kunnen meerdere dingen een rol spelen zoals de hoge stijfheid van het schedelmateriaal en het zou veroorzaakt kunnen worden door numerieke problemen. In de andere gebieden is er nauwelijks sprake van fluctuaties. In deze elementen is bovendien de

druk

veel lager omdat zij verder van de bron verwijderd zijn.

Voor het element halverwege de schedel (fig. 4.4) zit er nog een grote uitslag in de dnik,

dit verschijnsel zou kunnen komen door terugkaatsing.

Het resultaat van punt 4 (fig 4.6) stemt niet overeen met het resultaat van Merchant et al.

Waar Merchant sterke fluctuaties vindt, komt uit deze berekening een langzaam aangroei- ende druk om na het wegnemen van de belasting weer in te zakken. Zelfs de maximale

druk stemt niet overeen, die is in het onderzoek van Merchant et al. van dezelfde orde van

grootte als de belasting. In deze analyse is de maximale

dnik

een aantal decaden lager. Omdat de geometrie hetzelfde is, net als de eigenschappen voor het schedelmateriaal, kan het alleen liggen aan de materiaalgegevens van de hersenen. En dat is natuurlijk heel goed mogelijk, want daarin zitten hele grote verschillen. De analyse methode mag geen invloed hebben. Hoewel Merchant rekende met de eindige differentie methode mag de berekening met de eindige elementen methode niet zoveel verschil opleveren ondanks het feit dat op een bepaald punt in beide analyses aannamen worden gedaan.

Van het in de buurt van de contrecoupzone gelegen element (fig. 4.8) valt nog op te

merken dat het niet op tijd begint te trillen. Er zijn al eerder dan de berekende 0.09 ms trillingen waar te nemen. Dit is zeer opmerkelijk want het is onmogelijk dat de

druk

zich sneller voortplant dan met de geluidssnelheid in het medium. Een verklaring is dan sok

5.2 Bespreking van de resultaten uit de driehoekige drukpuls

Ook hier dient eerst gekeken te worden naar de afschattingen voor de verplaatsingen. De met de analyse berekende waarden zijn fors lager dan de waarden die bij de afschatting naar voren kwamen.

De rekken in deze analyse waren ook zeer klein zodat het geoorloofd was de lineaire theorie te hanteren. De verplaatsingen loodrecht op de belastingsrichting zijn bijzonder klein en met name voor de elementen in de buurt van de contrecoupzone is deze verplaat- sing zeer sterk fluctuerend in de tij& er is sprake van een frequentie van om en nabij de 6

kHZ.

Van de drukberekeningen kunnen er twee worden vergeleken met het onderzoek van Ruan

et al,

[%a].

Dat zijn de drukvariaties in de hersenen, punt 1 en punt 2 (fig. 4.13 en 4.14)

enerzijds en figuur 4.16 anderzijds. Allereerst zijn er een aantal verschillen. Ten eerste komen er bij Ruan sterke fluctuaties naar voren, bij deze berekening niet. En er is een grotere faseverschuiving bij 5 ms. Daarnaast liggen de maximale drukken wat hoger, maar zijn wel van dezelfde orde van grootte.

Een verklaring hiervoor moet gezocht worden in de gebruikte belasting en de geometrische verschillen, want de materiaaleigenschappen zijn volledig aan elkaar geljjk. Ruan belastte zijn model onder een andere tophoek (25") en hij is vergeten de afmetingen te publiceren. Dit is verreweg het belangrijkste argument. Aan de meshverdeling zou ook nog wat te

wijten kunnen zijn. Ruan gebruikt veel minder elementen, namelijk 138, waarvan 14 voor de schedel. Terwijl in deze opdracht gebruik werd gemaakt van 288 elementen, waarvan 48 voor de schedel.

Van het drukverloop voor de twee schedelelementen (fig. 4.12 en 4.15) valt nog te zeggen dat dat bij de impact (fig.4.12) hetzelfde verloop heeft als de impact met uitzondering van de hogere maximale waarde en de spiegeling in de tijdas. Dit is zeer vreemd, temeer dat op dit punt juist de

druk

zou moeten toenemen. Het andere schedelgedeelte aan de contrecoup zijde heeft weer een sterk fluctuerend gedrag (fig. 4.15). Dit is trouwens precies omgekeerd aan de berekeningen in de vorige analyse, want daar waren juist aan de voorzijde sterke trillingen waarneembaar in tegenstelling tot de achterzijde. De trillingen hebben een frequentie van 3 kHz. Een verklaring lijkt niet zo erg eenvoudig. Dit onder- werp verdient dan ook meer aandacht.

Conclusies en aanbeveiingen

Hoofdstuk 6

Conclusies en aanbevelingen

6.1 Conclusies

Er zijn een aantal problemen gesignaleerd, waarvoor echter geen verklaring kon worden gevonden. Van het model van Merchant et al. [9] konden de resultaten niet in een

redelijke mate gereproduceerd worden omdat binnen DIANA een vloeistof niet gemodel- leerd kan worden. Het enige dat aan dit model geverifieerd kon worden was de voortplan- tingssnelheid van de onibuls. Daar kan van gezegd worden dat dat niet in svereenstem- ming gebeurt met de analytische waarde. De drukpuls loopt voornamelijk via de schedel omdat daar de geluidssnelheid vele malen groter is dan in de hersenen. De waarde van het model is hiermee niet tot nul gereduceerd omdat een validatie is uitgevoerd met een model van Rum et al. [12]. Hieruit is gebleken dat de drukkarakteristieken voor de coup zijde en de contrecoup zijde in grote mate overeenstemmen. Dit betekent dat DIANA de mogelijk- heid bezit om het menselijk hoofd te modelleren, met de beperking dat het alleen een vaste stof modellering is. De problemen die gesignaleerd zijn kunnen natuurlijk ook een numerieke oorzaak hebben.

Om in DIANA aan een vloeistofmodellering voor de hersenen te beginnen lijkt tot op heden niet haalbaar vanwege het grote aantal aanpassingen dat binnen de module FLUSTR gemaakt zou moeten worden. Deze module is momenteel alleen geschikt in dit onderzoek

om eventuele berekeningen in het frequentiegebied uit te voeren.

6.2 Aanbevelingen

Uit dit onderzoek om een model te creëren van het menselijk hoofd volgen een aantal aanbevelingen. Zowel de vaste stof modellering als de vloeistofmodellering is mogelijk. De vloeistofmodellering zal echter veel meer tijd vergen omdat de module FLUSTR

gemodificeerd zal moeten worden. Er kan ook nog gekeken worden naar de mogelijkheden

die het eindige elementenpakket MARC biedt. Dit houdt dus in dat hetzelfde model nog een keer doorgerekend zou moeten worden. Op grond van de modificeermogelijkheden binnen DIANA wordt echter de voorkeur aan DIANA gegeven. Verder kan geprobeerd worden om de interface tussen schedel en schedelinhoud te modelleren met de hulp van contactelementen.

Conclusies en aanbevelingen

Daarnaast kan er op de reeds ingeslagen weg verder gegaan worden om de mogelijkheden te verkennen die DIANA bij de modellering van het menselijk hoofd biedt.

Naast al deze mogelijkheden zou er ook gekeken kunnen worden naar een ander materiaal- model. De wet van Hooke zal ongetwijfeld niet de benadering van het materiaalgedrag beschrijven, omdat het slechts een lineair model is. Het zou ook zeer aan te raden zijn om

meerdere berekeningen uit te voeren zoals daar zijn een frequentieresponsie analyse om di; te vergelijken met de resultaten van een zelfde analyse maar dan van een lege bol. In het frequentiegebied zou er ook nog een eigenfiequentieanalyse mogelijk zijn waarvan de resultaten vergeleken moeten worden met een lege bol. Wat in dit onderzoek niet expliciet is meegenomen is een berekening van de rekken, dat zou voor de toekomst dan ook ten zeerste aan te bevelen zijn.

Het verschijnsel waaraan in deze opdracht ook gedacht is reflectie. Dit zou eventueel op kunnen treden in de analyse van het model met de blokvormige drukpuls. Dit is echter een complex probleem en het verdient dan ook de voorkeur om dit verschijnsel te onderzoeken aan een veel eenvoudiger model, wat dan naar alle waarschijnlijkheid in niets meer

overeenkomt met de modellering van het hoofd. Dit lijkt noodzakelijk gezien de weinige ervaring die er is met golfvoortplanting.

Voor het modelleren zijn er natuurlijk ook nog tal van andere mogelijkheden. Er is al

gezegd dat dit model al een complex model is en interessant genoeg binnen het vakgebied van de toegepaste mechanica. Toch beschrijft het maar heel summier het hoofd. In de toekomst moet er dan ook aan gewerkt worden om het model completer te maken. Dat kan alleen gebeuren door in het model niet zoveel te verwaarlozen. Specifiek voor dit model zou gekeken kunnen worden naar het effect van meshverfijning of een grovere mesh, De mesh lijkt in de resultaten voldoende ve6ijnd. Dit kan ook geconcludeerd worden uit de vergelijking met de literatuur (Merchant et al. [9], Ruan et al. [12]). Maar het zou

natuurlijk zo kunnen zijn dat een locale meshverfijning in de nabijheid van de impact tot betere resultaten leidt. Uit de resultaten van de berekening met de driehoekige drukpuls volgde een grote gelijkenis met de resultaten van Ruan. Het is dan misschien de moeite waard om één van zijn andere modellen door te rekenen. Hiermee wordt dan direct de stap gemaakt naar een complexere geometrie, wat de werkelijkheid beter benadert.

Als laatste lijkt het ten zeerste aan te raden dat in het belang van het onderzoek er meer aandacht geschonken wordt aan de numerieke problemen die eventueel zouden kunnen optreden. Daar is in deze opdracht niet naar gekeken, maar moet zeker een aandachtsge- bied worden binnen het onderzoek.

Referenties

Referenties

Bathe, K.J, Finite Element Procedures in Engineering Analyses, Prentice Hall, Englewood Cliffs, 1982.

DIANA, User’s Manual Draft, Volume 10: Interaction Analysis, 1991.

Engin, Ali E, The axisymmetric response of a fluid-filled spherical shell to a local radial impulse -a model for head injury, J. Biomechanics, 2: 325-341, 1969. Engin, Ali E. and King Liu, Y, Axisymmetic response of a fluid-filled spherical shell in free vibrations,

J.

Biomechanics, 3: 11-22, 1970.

Hall, D.E, Basic Acoustics, John Wiley & Sons, Inc, New York, 1987.

Hamdi, Mohamed Ali, Ousset, Y. and Verchery G, A displacement method for the analysis of vibrations of coupled fluid-structure systems, int. J. for Numerical Methods in Engineering, 13: 139-150, 1978.

Khalil, T.B. and Hubbard, R.P, Parametric study of head response by finite element modeling, J. Biomechanics, 10: 119-132, 1977.

Krawietz, A, Materialtheorie, 1986.

Merchant, H.C. and Crispino, A.J, A dynamic analysis of an elastic model of the human head, J. Biomechanics, 7: 295-301, 1974.

Petyt, M. and Lim, S.P, Finite element analysis of the noise inside a mechanically

excited cylinder, Int. J. for Numerical Methods in Engineering, 13: 109-122, 1978. Oomens, C.W.J, Constitutieve modellen, dict. nr. 4687, Technisch Universiteit Eindhoven, Faculteit Werktuigbouwkunde, januari 199 1.

Rum, J.S, Khalil, T, King, A.1, Human Head Dynamic Response to Side Impact by Finite Element Modeling, J. Biomechanical Engineering, 1 13: 276-283, augustus

1991.

Schreurs, P. J.G, Continuuumsmechanica, dict. nr. 4612, Technische Universiteit Eindhoven, Faculteit Werktuigbouwkunde, 1990.

van Veen, B, Akoestische berekeningen met de eindige elementen methoden, afstudeerverslag, Universiteit Twente, faculteit der werktuigbouwkunde, augustus

Materiaaimodel

Appendix A

Materiaalmodel

Kenmerkend voor elastisch materiaalgedrag is dat de momentane spanningstsestand volledig bepaald wordt door de momentane deformatie en dat na het wegnemen van de belasting geen blijvende vervorming valt te constateren, Uitgaande van de eis dat de spanningstoestand in een punt van het materiaal niet beïnvloed mag worden door starre translatie en/of starre rotatie, kan de onderstaande algemene constitutieve relatie voor

o

op het momentane tijdstip t worden afgeleid:

o=di.H(u).R (A. 1)

waarbij H een functie is van de verlengingstensor U en R de rotatietensor is uit de polaire decompositie van F.

Andere formulering van de constitutieve vergelijking

De algemene constitutieve vergelijking voor (T luidt:

o=R-H(u).R

Het is gebruikelijk om bovenstaande constitutieve vergelijking anders te formuleren. Hierbij wordt van de polaire decompositie van F, de deformatietensor, en E, de Green- Lagrange rektensor gebruik gemaakt.

F=RU (A.3) 1 1 2 2 E=-(F "*F

-0

=-(U2-P> waaruit volgt: RzF-U-1 (A.4)

Voor de constitutieve relatie kan hiermee geschreven worden:

o=F.U-'det(u>(det(u>H(U))U-l .F (A.7)

Het linkerlid van deze vergelijking wordt de 2" Piola-Mirchhoff spanningstensor P genoemd. De algemene constitutieve vergelijking voor elastisch materiaalgedrag wordt hiermee:

P=G(E) (A.9)

Isotroop materiaalgedrag

Ten gevolge van de elastische deformatie zal in het materiaal elastische energie worden opgeslagen. Deze inwendig opgenomen energie is te schrijven als een functie van de deformatietoestand, gekarakteriseerd door de Green-Lagrange rektensor E. Voor @(E), de specifieke elastische energie ofwel de elastische energie per eenheid van massa kan genoteerd worden:

(A. 10)

waarbij e i de hoofdrekken en 13,~ de hoofdrekrichtingen zijn.

Voor isotroop materiaal geldt dat de materiaaleigenschappen in alle richtingen hetzelfde zijn. Dat houdt in dat de opgeslagen elastische energie per massaeenheid alleen afhangt van de hoofdrekken e i.

(A. 1 i)

Het is gebruikelijk om nu <D op te vatten als een functie van de drie invarianten van E. S1(E) = a ( E ) =e +e2 + e (A. 12)

(A. 13)

Matenaaimodel

Uit themsdynmhche beschouwingen kan de volgende relatie afgeleid worden tussen de specifieke elastische energie en de 2" Piola-Kirchhoff spanningstensor.

waarin

eo

de dichtheid van het materiaal in de referentietoestand is. Dit geeft dan voor isotroop materiaal:

(A. 17)

waarbij geldt voor de drie afgeleiden van de invarianten naar de Green-Lagrange rekten- sor:

dJ1

-

=I dE -=JII-E

dJz

dE -=J&-J,E+E~ f l 3 dE Substitutie geeft: (A. 18) (A.19) (A.20) (A.21) (A.22)

Materiaalmodel

Fysisch lineair, isotroop, elastisch materiaalgedrag

Hier is sprake van als het materiaalgedrag beschreven wordt door een lineair verband tussen de spanningstensor en de rektensor. Hierbij moet altijd vermeld worden welke spannings- en rektensor gebmpki worden.

Uitgaande van het verband tussen P en E voor isotroop materiaalgedrag

P=a,,I+a1E+a2E2 (A.23)

kan de lineaire relatie afgeleid worden.

Voor een lineair verband moet gelden:

a,

=C,

waarbij co en

c,

constanten zijn. Er geldt nu:

P =c,,tr(E)I+c,E (A.24)

Gebruik makend van de Lamé-constanten h en p wordt dit vaak als volgt geschreven.

met 2 =K--G Ev

A=--

2vG - 1-2v (l+v)(l-2v) 3 p=G= E 2( 1

+v>

(A.25)

waarbij E de elasticiteitsmoduhs,

v

de dwarscontractiecoefficient, K de compressiemodu- lus en G de glijdingsmodulus is. Het lineaire verband tussen P en E kan ook geschreven

Materiaalmodel

worden met de hulp van een constante 4”orde tensor.

P=C,~T-(E)~+C,E=CJI:E

+ ~ , 4 1 : ~

P

= (cJI+c;I) :E =4C:E

(A.26)

(A.27)

Constitutieve vergelijking voor de Cauchy spanningstensor

Uitgaande van het lineaire verband tussen P en E kan onderstaande constitutieve vergelij- Icing voor de Cauchy spanningstensor worden afgeleid

1 2

C=-(Y,(tï(B) -3)B+y,(B-I).B)

Hierin is B de linkse Cauchy-Green rektensor, gedefinieerd als

B =F.F (A.28) (A.29) met 1

y

= (det(B))-”co

Deze constitutieve vergelijking is duidelijk niet-lineair in B. Indien mag worden veronder- steld dat de rekken en de afschuivingen klein zijn, geldt:

c

=R

.(~,tï(E)l +c1E) .R (A.30)

met R de rotatietensor uit F=R- U.

Als ook de rotaties klein zijn volgt:

Appendix

B

Analysemethode

Ten gevolge van van een deformatie, die veroorzaakt wordt door externe krachten treden er in het materiaal spanningen op. De spanningstoestand in een materieel punt wordt volledig beschreven door de Cauchy spanningstensor. De deformatie moet zodanig verlopen dat op elk tijdstip t,,Szlt aan een aantal fysische wetmatigheden voldaan wordt. Deze wetmatigheden zijn de balanswetten. Hier zijn er een viertal van: de wet van behoud van massa, de impuls-balanswet, de impulsmoment-balanswet en de wet van behoud van energie. Beschouw een drie-dimensionaal lichaam met een volume V(t), waarop de volumebelasting

4

(xf

t) en de randbelasting

fl

(xft) werken. Het buitenoppervlak wordt gedefinieerd door A(t). Stel:

waarin respektievelijk staan gegeven de verplaatsingsvector, de rektensor en de spannings- tensor. Dan wordt ingeval van een lineaire deformatie als ook de rotaties klein zijn de lineaire rektensor gegeven door:

1

2

E($t) =-((9

2)

+ (V

ay}

waarin 3 de gradiëntoperator is.

De impuls-balanswet levert de bewegingsvergelijking en die wordt geschreven als:

Hierin is

+e

de dichtheid.

De wet van behoud van massa levert in geval van kleine verplaatsingen op dat de dichtheid constant verondersteld kan worden. De impulsmoment-balanswet levert:

De constitutieve relatie voor elastisch materiaal luidt:

waw 4C! eer, constmts vierde OX-& t m s m 4s. Voor een isotroop en homogeen materiaal zijn

twee constanten-zoals de compressiemodulus en de glijdingsmodulus evenals de modulus van Young en de constante van Poisson onvoldoende om het materiaal te karakteriseren. Vervolgens wordt de dynamische randvoorwaarde gedefinieerd met de hulp van de stelling van Cauchy welke zegt dat in elk punt van een lichaam een tensor bestaat, de zogenaamde spanningstensor, waarmee de spanningsvector op elk vlak door dat punt bepaald kan worden.

waarin p’ de spanningsvector is,

c

de Cauchy spanningstensor en 2 de eenheidsbuitennor- maal op het betreffende vlak.

Daarnaast kunnen de volgende kinematische randvoorwaarden aangenomen worden.

waarin

$0

en

$0

de initiële verplaatsings- en snelheidsvector zijn.

Het stelsel vergelijkingen (€3.1) tot en met (€3.6) is een oplosbaar probleem en er kan de verandering in de tijd mee uitgerekend worden voor het verplaatsingsveld, het r e h e l d alsook de spanning. Oplossingen in analytische vorm zijn wel te produceren zij het voor eenvoudige geometneën en externe belastingen. Echter voor de meeste practische

problemen is er een benadering nodig. Er zijn verschillende numerieke benaderingstechnie- ken zoals bijvoorbeeld de eindige differentie methode die door Merchant et al. [9] werd gebruikt om dit zelfde probleem op te lossen. Een andere mogelijkheid is de eindige elementen methode. Er volgt nu een korte beschrijving van de methode voor dynamisch lineair elastische problemen.

Het gebied V wordt onderverdeeld in een eindig aantal kleinere gebieden, die elementen genoemd worden. De geometrie van de elementen, en dus ook het volume V, wordt

gedefinieerd door een verzameling van N punten welke de knooppunten zijn. In de verplaatsingsformulering van de eindige elementen methode is elk verplaatsingsveld gedefinieerd door een verzameling van continue functies in termen van de knooppuntsver-

plaatsingen. Voor één element i kan de verplaatsingsfunctie geschreven worden als:

X(q),t) =") -if (t) 03.7)

Eierin is

E(.?)

de r ~ a î x i i rlsi de vomfitmcties er, d ( t ) de gegemm.liseerde. verglaattings-

vector. Convergentie van de oplossing is verzekerd door de volledigheid van de vormfumc- ties wat inhoudt dat er geen rekken van elementen worden berekend ten gevolge van de verplaatsing als star lichaam.

De gediscretiseerde dynamische bewegingsvergelijkingen kunnen afgeleid worden uit het principe van virtuele arbeid of een variatie methode zoals het principe van Hamilton, wat geschreven kan worden als:

waarin

6

een variatie operator is en L de Lagrangiaan.

Voor het dynamische randvoorwaarde probleem geldt de volgende relatie:

Substitutie van de vergelijkingen (B.l) en (B.2) in (B.9), waarna vergelijking (B.7) wordt gebruikt en het uitvoeren van het minimalisatieproces zoals dat gedefinieerd is in vergelij- king (B.8) geeft een stelsel van tweede orde lineaire differentiaalvergelijkingen welke geschreven kan worden als:

(B.lO) waarin:

-

A4 de massamatrix is,

-

K

de stijfheidsmatrix is,

f

(t)

Het stelsel bewegingsvergelijkingen wordt opgelost met het integratieschema Newmark B

impliciet dat ontwikkeld is door Newmark. Integratie van het stelsel vergelijkingen (B. 10) levert de verplaatsingskolom a (t) die de totale verplaatsing beschrijft van de

knooppunten. De verplaatsingsvector T(t) die de verplaatsing als star lichaam en de deformaties van elk willekeurig punt binnen V beschrijft wordt verkregen met vergelijking

(B.7). Vervolgens worden de rekken en de spanningen berekend volgens respektievelijk (B.2) en (B.4). Nu moet er echter nog een uitdrukking gevonden worden voor de hydrosta- tische ikuk want dat is wat vergeleken kan worden met voorgaande studies. Voor een incompressibel materiaal geldt dat er geen volumeverandering optreedt. In praktijksituaties blijkt dit slechts bij benadering op te gaan. De hydrostatische druk wordt daarna op de volgende wijze berekend.

Voor een incompressibel isotroop elastisch materiaal kan de volgende relatie afgeleid worden voor de Cauchy spanningstensor (Krawietz [8] en Oomens Ell]):

o =

-pI+od=_tr(o)l+od 1 (B. 1 1)

3

hierin is p de hydrostatische druk en deze is dus gelijk aan:

1 1

(B. 12)

Verplaatsingsresultaten van de berekeningen

Appendix

C

Verplaatsingsresultaten

van de

berekeningen

C.1 Blokvormige drukpuls

Naast de drukresultaten zijn er ook resultaten met betrekking tot de verplaatsing van het model. Er is gekeken naar de twee richtingen waarin het lichaam zich kon bewegen; de x-

richting is de richting loodrecht op de belasting en de y-richting de richting parallel aan de

symmetrieas. Voor alle zes punten zijn beide verplaatsingen uitgerekend en deze zijn weergegeven in de figuren C.l tot en met C.12.

~ 1 0 - 3 x-verplaatsmg van punt 2 y-verplaatsmg van punt 2

6 O O02 0.04 0.06 O08 0.1 0.12 0.14 O 16

tijd [ms] tijd [mS]

figuur C.l: x-verplaatsing van punt 2. figuur C.2: y-verplaatsing van punt 2.

tijd [ms] tijd [ms]

Verplaatsingsresuiîaten van de berekeningen 2.5 2 1.5 I E 1 .s 0 5 E $ 9 0 . M z -0.5 -1 -1.5 ... ... ... ... _ ... +-....- _ ... ; ... ... ... ... ... ... _ j ...c... O 0.02. 0.04 0.06 0.08 0.1 0.12 0.14 0.16

x i 0 3 x."erplaaaing van punt 4

xi05 y-verpiaatring van punt 1

O 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 tijd [ms]

figuur C.6: y-verplaatsing van punt 1.

tijd [ms] tijd [ms]

figuur C.7: x-verplaatsing van punt 4. figuur C.8: y-verplaatsing van punt 4.

In deze figuren is duidelijk te zien dat d-e punten verder van de impactzone gelegen er enige tijd over doen voordat ze gaan bewegen. Opmerkelijk is verder dat de y-verplaatsing van punt nul een bijna constante snelheid heeft en dus een versnelling gelijk aan nul. Vrijwel alle verplaatsingen verlopen glad in tegenstelling tot een verloop in de druk.

Daarnaast zijn alle verplaatsing zeer klein; in ieder geval kleiner dan de analytisch berekende waarden want de verplaatsing van de contrecoup zijde (fig. C.6) is kleiner dan

de afgeschatte waarde.

Zeer vreemd is de verplaatsing in de positieve y-richting ( tegen de richting van de belasting in) van de in de hersenen gesitueerde punten (fig. C.8, C.10 en C.12). Een

Verplaatsingsresulîaten van de berekeningen

tijd [ms] tijd [ms]

figuur C.9: x-verplaatsing van punt 5. figuur C.10: y-verplaatsing van punt 5.

x10' x."erplaasingvan punt 3

tijd [m] tijd [ms]

figuur C.ll: x-verplaatsing van punt 3. figuur C.12: y-verplaatsing van punt 3. Opmerkelijk aan twee paar van de bovenstaande figuren is dat ze sterk op elkaar lijken, hetzij in spiegelbeeld ten opzichte van de tijdas (fig. C.2 en C.8 en fig. C.6 en C.12). Dit zijn punten die zeer dicht bij elkaar liggen. Het eerste paar ligt in de buurt van de impact en het tweede paar in de buurt van de contrecoupzijde. De bijbehorende x-verplaatsing (fig.

C.

1 en C.7 en fig. C.5 en C. 11) lijken daarentegen geenszins op elkaar.

Uit de waarden voor de verplaatsingen kan niet direct geconcludeerd worden dat er aan de lineaire theorie wordt voldaan met betrekking tot de rekken. Daar kan alleen wat over

gezegd worden als de rek wordt berekend volgens Green-kagrange en de lineaire rek. De figuren met de rek als functie van de tijd voor punt 4 en punt 3 zijn respektievelijk de figuren C. 13 en

C.

14.

De rek is berekend voor het lijnstukje van de twee bij elkaar in de buurt gelegen punten in de vloeistof.

De rek van Dunt 4 ~ 1 0 - 3 De rek van punt 3

tijd [mS] tiid [msl

figuur C. 13: De lineaire-en Green-Lagran-

ge rek van punt 4.

figuur C. 14: De lineaire-en Green-Lagran- ge rek van punt 3.

De doorgetrokken lijn is de lineaire rek, de gestippelde lijn de Green-Lagrange rek. Uit beide grafieken blijkt dat beide rekdefinities bijna tot hetzelfde resultaat leiden. Het was dus inderdaad gerechtvaardigd om de lineaire theorie te gebruiken.

C.2 Driehoekige drukpuls

In deze analyse zijn ook een aantal verplaatsingen in x-richting en y-richting berekend. Al deze verplaatsingen staan in de figuren C.15 tot en met C.22. Over deze karakteristieken

valt het volgende op te merken. De verplaatsingen in de y-richting zijn voor alle elemen- ten gelijk. Het is dan ook gerechtvaardigd te concluderen dat de rekken in de y-richting heel klein zullen zijn. Dit blijkt namelijk uit de y-verplaatsing voor de elementen nul en één en de elementen twee en drie. Deze verplaatsingen zijn zo op elkaar gelijkend dat er nauwelijks enige verlenging optreedt, wat op zijn beurt weer inhoudt dat de twee gebruikte rekdefinities aan elkaar gelijk zullen zijn.

Verplaaîsingsresuiten van de berekeningen

y-verplaatsing van punt 3

tijd [ms] tijd [ms]

figuur C.15: x-verplaatsing van punt 3.

~ 1 0 - 3 x-verplaatsmg van punt 2

figuur C.16: y-verplaatsing van punt 3.

y-verplaatsing van punt 2

-0.4 - 3 -0.6 - L. g -1 - -1.2 - O 2 4 6 8 10 12 14 16 O 2 4 6 8 10 12 14 16 tlJd [ms] t1ld [ml

figuur C.17: x-verplaatsing van punt 2. figuur C. 18: y-verplaatsing van punt 2.

Opmerkelijk in deze grafieken en dan met name van de verplaatsingen in de x-richting is dat de verplaatsing in de x-richting voor punt O zeer sterk fluctueert en zeer klein is en de

rekken in x-richtingen dus ook want de verplaatsing van de symmetrieas in de x-richting is gelijk aan nul. Uit de grafieken blijkt dat de werkelijk optredende verplaatsing in de y- richting naar verwachting veel kleiner is dan de afgeschatte waarde.

Opvallend is verder het gedrag van de punten in de hersenen (fig. C.17 en C.19). Ze

bewegen enigszins met een golfbeweging. De maximale waarde die bereikt wordt ligt hoger dan de waarde die de schedelpunten bereiken. Dit ligt aan het feit dat de schedelin- houd een veel lagere modulus van Young heeft dan de schedel

De verplaatsingen van de elementen in de nabijheid van de contrecoup staan weergegeven in d e nu volgende figuren. Het eerste paar beschrijft de verplaatsing van het hersenelement en het laatste paar de verplaatsing van de schedelzone.

~ 1 0 - 3 x-vemiaatsinf! van Dunt 1 y-verplaatsing van punt 1

tijd [ms]

figuur C . 19: x-verplaatsing van punt 1.

xi05 x-vepiaatsing van punt O

tijd [ms]

figuur C.20: y-verplaatsing van punt 1.

y-verplaatsing van punt O

tijd [ms] tijd [ms]