Mathematische modellen ter beschrijving van de dynamica van een vliesklepprothese

van een vliesklepprothese

Citation for published version (APA):

Veugelers, W. J. T. (1987). Mathematische modellen ter beschrijving van de dynamica van een vliesklepprothese. (DCT rapporten; Vol. 1987.063). Technische Universiteit Eindhoven.

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

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

providing details and we will investigate your claim.

Afdeling der Werktuigbouwkunde

Vakgroep Fundamentele Werktuigbouwkunde (WFW)

Afdeling der Technische Natuurkunde Vakgroep Transportfysica (NT)

Mathematische modellen ter beschrijving van de dynamica van een

vliesklepprothese.

W.J.T. Veugelers tussen ti j d s rapport

rapportnummer WFW: 87 .O63

Begeleiding: ir. J.B.A.M. Horsten

dr. ir. A.A. van Steenhoven

In het kader van het interunj.versitair project Hartklepprothesen wordt.

getracht een mathematisch model op t e stellen, d a t de klepvliesbeweging van en de vloeistofstromiiig rcmd een vliesklepprothese gedetailleerd beschrijft.

Dit: rapport bespreekt: enkele, u i t de literatuur bekende, quasi-één-

dimensionale analytische modellen en een numerkk model, d i e een besciiri j-

ving geven van d e dynamica van een vliesklepprothese.

De quasi-~~n-dimensionale modellen behandelen de vliesklep twee-

dimensimaal ' als een s t a r r e p l a a t rotereiid om hek aan1iecht;ingspunt. Buiging

van het v l i e s wordt dus n i e t beschreven.

Het numerieke model kan wel de buiging van het v l i e s beschrijven. Per

ti. j d s t a p wordt de vloeistofstroming en vliesbeweging ontkoppeld berekend.

H i e r b i j wordt de vloeistof stroming opgelost i n een a r b i t r a i r e Lagrange-

Euler-formulering of i.n een Euler-formulering

,

de vli.esbewegj.ng i n een

Lagrange-formulering. Met numerieke model l i j k t geschikt voor de modellering van een vìieuk1epproi:hese.

2 Natuurlijke aortaklep 2

3 Quasi-t%n-dimensionale modellen

3 . 1 Enkele yuasi-ken-dimensionale modellen

3.2 Evaluatie

4 V1oeist:of -st:ructuur interact:i.emodel

4 . 1 Probleemstelling 4.2 Impulsvergelijking 4 . 2 . 1 Vloeistof 4 . 2 . 2 Structuur 4 . 3 Meshdefinitie 4 . 4 Behoudswetten in de ALE-beschrijving 4 . 5 Eindige elementenbenadering 4 . 5 . 1 Stokes-probleem

4.5

-

2 Navier-Stokes probleem; Newton iteratie 4 5 , 3 Boetefunct.i.emethode

4.6 Di.screkisatie van hek vloeistofsysteem

in

de ALE- beschri. jviny

4.7 Discretisa.tie van de structuur 4.8 Ti jddiscretisatie

4.8

.

1 CJ-metl-iode

4 - 8 i 2 Stabiliteit van het vloeistofsysteem

4.8.3 Stabiliteit van het: structuursysteem

4.8.4 Meshpartj.tie

4.8.5 Stabiliteit van het gekoppelde systeem

4.9 Praktische problemen en onduidelijkheden

4 . 9 . 1 Automatische aanpassing van ùe v1oeistc)fmesl-i 4.9.2 Convect:j.eve term

4 . 9 . 3 Dirichlet randvoorwaarde voor de vloeistof 4.9 . 4 Neumann randvoorwaarde voor de structuur

4 . 1 0 Algoritme 5 5 14 16 17 18 18 19 20 2 1 2 4 24 28 30 33 36 38 38 40 43 43 45 47 47 49 49 50 5 1

5 Conclü:iie

6 Voortzetting

Appen(1j.x A

Appendix B

Literatuurlijst:

Ik tijrlsaffiankelijke coëff iciënten in het model van Wbppermann(l385)

Dynamica van een roterend s t a r plaatje

54

55

57

58 61

1 INLEIDING

Het doel van het project Hartklepprothesen is het opstellen van

ontwerpspecificaties waaraan een vezelversterkte vliesklepprothese moet voldoen. Binnen het deelproject Openings- en sluitingsgedrag is het

onderzoek vooral gericht op de klepvliesbeweging van en vloeistofstroming rand een vliesklepprothese. Getracht wordt een mathematisch model op te stellen, dat de klepvliesbeweging en de vloeistofstroming gedetailleerd beschrijft. Via een parameterstudie zal vervolgens worden geprobeerd tot een optimale vliesklepprothese te komen. In dit kader is, tijdens de

beginperiode van mijn afstudeeronderzoek, dit rapport geschreven, met als doel inzicht te krijgen in een mathematisch model, dat de dynamica van de vliesklepprothese kan beschrijven.

Eerst zal worden ingegaan op de natuurlijke aortaklep (hoofdstuk 2). Daarna komen enkele quasi-4dn-dimensionale modellen aan de orde, die een

analytische beschrijving geven van de vliesbeweging als functie van de tijd (hoofdstuk 3 )

.

Hoofdstuk 4 bespreekt een numeriek vloeistof-vlies

interactiemodel. De toepasbaarheid van deze mathematische modellen wordt besproken in hoofdstuk 5 . Tenslotte wordt een voorstel gedaan voor de voortgang van het onderzoek (hoofdstuk 6 ) .

2 NATUURZT JKE AORTAKLEP

De aortaklep is een ongespierde uitlaatklep van het hart en bevindt zich tussen de lichaamsslagader (aorta) en de linkerhartkanier, zie figuur 2 . 1 .

Zij bestaat in wezen uit drie delen: drie dunne vlksjes (0.6 mm), hun bevestiging aan de klepwand (aortaklepring) en achter ieder: vliesje een zakvormige uitzetting van de klepwand (sinus van Valsalva). Een schematische weergave wordt gegeven in figuur 2 . 2 .

Bij de mens is de diameter van de aorta,

1

tot 2 cm achter de aortaklep, ongeveer 30 mm. De straal van de sinus i s ongeveer gelijk aan de straal van de aorta, de lengte van de vliezen ongeveer anderhalf maal de straal van de sinus (Reid, 1 9 7 0 ) .

Ongeveer 4cin keer per secc~nde trekt. de linkerhartkamer zich samen (systole) en ontspant deze zich (diastole). De drie vliesjes sluiten tijdens de

diastole de aorta af. De vloeistofstroom door de aorta heeft een verloop zoals getekend in figuur 2 . 3 .

I I

94 o,$ '

I

Fig. 2 . 3 Het, verloap van de aortaflow a l s functie van de tijd (Milnor, 1 9 8 2 ) .

Tijdens de systole perst de linkerhartkamer het bloed door de geopende aortaklep de aorta in, en is er dus sprake van een voorwaartse stroming. Deze stroming neemt geleidelijk af als gevolg van de afnemende

samentrekkingskracht van het. hart (vloeistofvertrayingsfase). Tijdens deze fase beginnen de vliezen al geleidelijk naar elkaar toe te bewegen. Op het, einde van de systcde is daardoor het doorstrominysoppervlak van de aorta reeds voor 75% door de klepvliezen afgesloten, en is er nog slechts een

I

Q .

Fig. 2 . 1 De ligging van de aorta- klep tussen de linkerkamer en de aorta.

6.

hartkamer

Fig. 2.2 a. Schematische weergave van de aortaklep gezien vanuit de aorta.

b. Schematische weergave van de aorta, de aortaklepring en één sinusholte (van Steenhoven 1981).

kleine terugstroming nodig om de klep volledig t e doen sluiten (Bellhouse &

Talbot, 1 9 6 9 ) . D i t heeft t o t gevolg d a t zo weinig mogelijk bloed weer terugstroomt i n het hart, hetgeen uiteraard van groot belang i s voor een e f f i c i ë n t werkend hart. Het g e l e i d e l i j k sluiten van de aortaklep b l i j k t

vooral mogelijk gemaakt t e worden door de aanwezigheid van de sinus-van- Valsalva-holte achter het v l i e s . Is deze holte n i e t of nauwelijks aanweziy dan wordt het mechanisme, waartoe ook een vereffening van de druk i n de sinus-holte behoort, verstoord en z a l de klep zich abrupt onder de invloed van de terugstroming sluiten (Bellhouse & Talbot, 1969; van Steenhoven,

1 9 7 9 )

,

De karakteristieke kengetallen voor een instationaixe stroming ( f i g u u r 2.3)

z i j n tiet Keynoldsgetal en Stroutlalgetal. Deze worden hier als volgt gedefinieerd:

VR

Re =

-

_'.J

en S t = R

-

VT

met V de maximale snelheid, de pieksnelheid tijdens de systolefase

(-1.2

ms-' )

,

R de straal van de s i n u s (-1.5 cm)

,

v de kinematische v i s c o s i t e i t van het Bloed (*4*10-' m2s-') en T de t i j d s d u u r van de vertragingsfase (-0.18

5 ) . D i t levext;

3 gUASS EEN-DIMENSIONALE M0I)ELLEN

Dit hoofdstuk bespreekt enkele quasi-&&-dimensionale modellen, die een analyti.sche beschrijving geven van cle vliesbeweging. ne geometrie van de gemodelleerde aortaklep is twee-dimensionaal. Bet vlies blijft recht en roteert vrij om tiet, aanhechti.ngspunt

.

De vloeistof stroming door het model wordt &n-dj.mensionaal verondersteld

, d .

w. z . de druk en de snelheid zijn alleen afhankelijk van cle plaatscodrdinaat x en de tijd t (zie bijvoorbeeld figuur 3.1 ) I De stroming is laminair, incumpressibel, niet-visceus; en w c x d t

dus beschreven met de &n-dimensionale hxmpressibele Eulervergelijkingen,

I 2 P (behoud van impuls 1 '3U >u

-

_{a t}

_{+ U K}

_{+ y =}

_{= o}

(behoud van massa 1

met u: snelheid, p: druk, A: oppervlakte van doorstroming.

Getracht wordt een vergelijking op te stellen, d i e de vliesbeweging als

functie van d e tijd beschrijft,

3 . 1 Enkele suasi-een-dimensionale modellen

De verschillen in de theorie4n uiten zich hoofdzakelijk in de aannamen voor

de vloej.stofs.trciining aan de sinuszijde van de vlj.ezen.

in

B e l i h ~ ~ s - _ & _ T a i b ~ t ( l 9 ~ 9 ~

wordt de stromj.ng in de sinus gemodelleerd met behulp van een gemodificeerde Hill's wervel. De haufdstroming wordt hen- dimensionaal niet-visceus verondersteld, Volgens van Steenhoven & van I)onyen(1979f resulteert dit model in een drukverschil over de vliezen dat vrkj grool: is ten opzichte van de kleine massa van de vlkzen. Om deze reden

z a l het model hier verder niet worden besproken.

~ a n _ - - - i n - o ~ e n _ - _ ~ a n - ~ ~ ~ ~ e n ( l 9 ~ 9 )

bestuderen het slultingsgedrag uan een twee-dimensionaal model van de aortaklep. In fi.guur 3.1 wordt; heli

aortaklepmodel schematisch weergegeven. 0bserva.I.i es tonen aan d a t voor kleine St:rouhal-waarden ( < O . 15)

,

het v l i e s , i n (le beginfase van : j l u i t i n g ,

vrijwel recht b l i 3 f t : en roteert rondom he% aanhechtingspunt. Ook z i j n dan de snelheden i n de s i n u s veel kleiner dan i n de hoofdstroming. Op grond hiervan wordt aangenomen dat de druk i n de s i n u s ccinstant i s en g e l i j k aan de druk

aan h e t uit:ei.nde van het v l i e s . De geometrie van de s i n u s wordt; n i e t

exp1ici.e.i: meegenomen.

o

L

+ %

F i g 3.1 Schematische weergave van het 2D-aortaklepmodel volgens van Steenhoven & van Dongen( 1979) ,

Definieer een sluitingsratío

Verder wordt aangenomen dat

I

1

- A }

< <

1 ,

zodat in a l l e vergelijkingen termen

van de orde { l - A ) kunnen worden verwaarloosd. Onder deze aaname i s de

heweyingsvergeli~king van het v l i e s slechts geldig i n de beginfase van slui.ti.ng van de aortaklep. Vanwege de verwaarlooshare massa van liet v l i e s I

kan over: fret v l i e s geen drukverschil worden opgebouwd, Deze c o n d i t i e wordt n i e t Zokaal opgedrukt i maar aangenomen wordt I dat het gemiddelde

drukverschil over fret: vlj.es n u l i s .

2

U i t behoud van massa volgt

U i t behoud van impuls volgt.

met po= p(x=O,t).

2

Na substitutie en verwaarlozing van de termen û(1-A) volgt;

Door liet gemiddelde drukverschil over liet vlies nul te stellen, i.e

L

waarin de druk in de sinus gelijk is genomen aan de druk aan het uiteinde van het vlies, volgt een lineaire differentiaalvergeli7king voor A :

dh

met beginvoorwaarden %=O, A=l o p t=O

.

Bovenstaande vergelijking beschrijft de vliesbeweging als functie van de tijd, i.n de beginfase van slui.ti.ng (kleine waarden voor l-A). Bi.j een gegeven snelheid liO( ti) wordt de vergelijking numeriek opgelust.

Een vergelijking v m bovenstaande theorie met de experimentele resultaten van Bellhouse( 1969) wordt gegeven i.n figuur 3 , 2 . Belltiousef 1969 ) en

Eeflbouse & Talbot(1969) hebben experimenten uitgevoerd aan een 3D-model van de aortaklep, met, skarre wanden en vlhsjes van dun nylon (0.1 mm dik)

Bellhouse( 1969) met de behandelde theorie.

o f experiment;

,

theorie (van Cteenhoven van Donyen, 1979).

@ p ~ e ~ g m n ( 1 9 8 5 ) beschrijft zowel het openingsgedrag als het sluitingsgedrag

van de aortaklep. De aanwezigheid van de sinus wordt wel meegenomen. Het is mogelijk de afmeting van de sinus te varieren. De druk wordt aan weerszijden van het vlies berekend (dit in tegenstelling t o t van Steenlaoven & van

Dongen(1979), die de druk aan de sinuszijde van het vlies uniform en

constant veronderstelden). Fiyuur 3 . 3 geeft schematisch de geometrie van het aortaklepmodel, De gebruikte symbolen worden aangegeven in figuur 3 . 4 .

Definieer een sluitingsratio

met Ag de doorstroomopening op x=ûi en Ai(xLlt) de minimale doorctroomopening tussen de vliezen (op x=xJJ).

F i g . 3 . 3 De geometrie van het. SD-aortaklepmodel, v o l g e n s Wippermann(l985).

x

Fi.g. 3 . 4 De g e b r u i k t e synibolen i n d e modelvorming (Wippermann, 1 9 8 5 ) .

ne doorstroomopeniiig tussen de vliezen wcirdt gegeven door

De doorslroomopening aan de sinuszijde van de vliezen door

met A, ( xI t;) de daorstrr>omopenj.ng in de sinussen

De factor ni maak het mogelijk de afmeting van de sinus te variëren; m=l geeft de tna.ximaìe afmeting o p x=L. (Voor dit $kn-dimensj,onale model i s de

vorm en afmeting van de sinus vooi x>I, niet van belang.) Dus

Ex: wordt onderschei.d gemaakt tussen stroming binnen de vliezen (index i) en stroming aan de sinuszijde van de vliezen (index a ) ; op ieder tijdstip geldt. dat

De vliezen worden verondersteld yeen massa t e irehhen; voor het vlies geldt dan een niomenteneveiiwicht

1

X

met s de col2rdinaat langs het vlies, S=r,.- , xL

We hebben nu ach% niet-lineaire vergelijkingen ( 1 t/m 8 ) met acht onbekenden

Ba, Ai, uaf uif pa,

afhankelijkt de laatste twee enkel van t.

De veryeli. jkinyen worden i n de volgende drie situaties gelinearlseercl

xLf en A . ne eerste zes onbekenden zijn van x en t

A z I , A 4 . 5 , AZO

1) h z l

Na substitutie en linearisatie volgt de bewegingsvergelijking

duo

f l t f en f zijn functies van uo, en m (appendix A )

2 3

Indien de aanname wclrdt gemaakt, zoals bij van Steenhoven & van Danyen(l979) i dat

m.a.w. dat m+= (oneindig yroke si.nus) i dan reduceert de beweglnysvergelijking tot

mek voor de coëfficiënten c1 i c 2 ' c

Vergelijk deze met: de coëfficiënten van van Steenhoven & van Dongen(1979) en c4 resp. S,10/3, 5 / 2 , 10/3.

3 ? 6 / 3 , 4, 8 1 3 , 4

Hek kleine verschil is te wi3ten aan de aaname van een momentenevenwicht in plaats van een krachtenevenwicht:, zoals bij van Steenhoven SL van

Dongen( 1979).

2 ) A z 0 . S

De bewegingsvergeli3king wordt dan

g en g zijn functies van u

-

en m

E i l i 2 3 O' dt (appendix A )

3 ) A30

De beweginysvergeli3king wordt dan

hl, 11 en h zijn functies van uo,

at-,

Hierbij is aangenomen dat

2 3 (appendix A )

Want voor AzO (in Ret. begin van de openingsfase c.q. op het eind van de sluitingsfase) i s de beweging van de vliezen vrijwel onafhankelijk van de afmeting van de sinus.

De zo verkregen drie lineaire differentiaalvergelijkingen, geldig op drie bepaalde tijdstippen, nl. als Azl, Az0.5 en AZO, wcmien met behulp van een

gewc,gen-gemiddelde-methode aan elkaar gefit:

,

zodanig dat twee

O . 5<h<

1

.

Door toepassing van een gewogen-gemiddelde-methode worden de

dj.f ferentiaalvergelj. jkingen niet-1i.neai.r. Hel: resultaat,

O < h < O . 5

met, de beyi.nvoorwaarden %=O, h=O op t = O .

I n f iyuur 3.5 word% het theoretiuche resultaat (doorgetrokken curve)

vergeleken met de experimentele resultaten ( o ) van Bellhouse & Talbot( 1%9)

.

De onderbroken curve is berekend met dezelfde theorie, maar met een

vertraagde s t a . r t t i j d t f t /T=0.05); de overeenkomst met de experimentele

resultaten i s i e t s beter. I n het. vroege begin van de openingsfase buigen de vlj.esjes i n de experimentele klep een k1ei.n beetje naar voren, maar r k klep

b l i j f t tijdens d i t kort t i j d s i n t e r v a l gesloten. Daarna bewegen de v l i e s j e s snel naar hun volledig geopende stand t o e . D i t is niet; ZO b i j de

gemodelleerde klep; de v l i e s j e s (die recht z i j n ) bewegen onmiddellijk a l s de stromhy begint. I n het geval van de klepopening met: een vertraagde

s t a r t t i j d , worden de beginvoorwaarden d h/d.t2= O, h=O o p t=t gebruikt.

ti Li

2

I .2 c

j

O 0.25 0.50 0.75 I .o0

Time i l T

Fig. 3 . 5 Tijdsafhankelijke beweging van de nortaklep: o experimen- tele resultaten van Bellhouse fi Talbot(1969);

-,

bewe- ging berekend met de behandelde theorie;

----

,

resultaat met hetzelfde model maar met een vertraagde starttijd (Wipper- mann, 1985).

3 . 2 Evaluatie

De besproken theorieën verschillen hoofdzakelijk in de aannamen voor de vloePstofstroming aan de sinuszijde van de vliezen.

Bellhüu:je & Talbot( 1 9 6 9 ) modelleren de stroming achter de vliezen tnet een Hill's wervel. Hun mathematische beschrijving resulteert evenwel in een fysisch niet acceptabel drukverschil over de vliezen.

Van Steenhoven & van Dongen(l979) veronderstellen verwaarloosbare snelheden in de sinus, waardoor de druk aan de sinuszijde van het vlies uniform en gelijk wordt genomen aan de druk aan het uiteinde van het vlies. Hun beweyinysvergelj.jking i s slechts geldig in de beginfase van sluj.tiny. Een bevredigende overeenstemming tussen theorie en experiment wordt gevonden. Wippermannfl985) veronderstelt geen constante druk in de sinussen; de druk aan de sinuszijde van de vliezen wordt berekend uit de verplaatsing van de

v l o e k k o f tijdens de openings- en sluitingsfase van de klep. Twee

bewegingsvergelijkingen worden verkregen, die samen het openings- en sluitingsyedrag beschrijven. In het limiet-geval yaat de

over tot de bewegingsvergeli2king van van Cteenhoven & van Donyen(1979). De overeenstemming tussen theorie en experiment is bevredigend, maar niet

4 VLOEISTOF-STRUCTUUR INTERACTTENODEI,

ïn di:t hoofdstuk wordt een numeriek vloeistof -structuur interactiemodel besproken, met als uitgangspunt het werk van Befytschko(l982)

respectievelijk Donca ( 1982 )

.

Zowel Belytschko( 1982) a l s Donea( 1982) zijn gebkeresseerd in het gedrag van onderdelen in een kernreactor bij een exp1osj.e in de reactorkern. De

beacliouwde structuur is ondergedompeld in de koelvloeistof. Door beiden wordt een arbitraire Lagrange-Euler (ALE) eindige elementenmethode

gepresenteerd, die de respons van een vloeistof-structuur interactiesysteem

moet voorspellen. 1% een ALE-formulering kunnen de knooppunten van de vloeistofmesh onafhankelijk van de vloeistofs%romj.ny worden bewogen. Dit leidt; tot een eenvoudige en nauwkeurj.ge behandeling van het contactoppervlak tussen de vloeistof en structuur. Eelytschko(l982) geeft een afleiding van de behoudswetten voor een compressi.bele, vi.skeuze vloei.stof in de ALE- beschrijving. Ook wordt een eindige elementenformulering van de

vloei.stofvergelij kingen gegeven. Donea ( 1982 1 behandelt eerst de

behoudswetten voor een compressibele, niet-viskeuze vloeistof in de ALE- beschrijving. Daarna volgt de ruj.mtelijke discretisatie o p b a s k van de

eindige elementenmethode, de numerieke tijdintegratie van de semi- gediscretj.seerde Vergelijkingen en een algoritme om de vloeistofmesh

autcmatisch aan t e passen. Reiden geven enkele numerieke resultaten, die de door hun gepresenteerde methode toelichten.

ïn het hier t e bespreken vloeistof -struc2-.uur ixiteractiemodel I wordt de

vlnel.stof i.ncc3mpressj.bel en viskeus veronders t:eM. Uitvoerig wordt; ingegaan

op het oplossen van de vloeictof stroming in de ALE-beschrijving

.

Hiervoor wordt een boetefunctiemethode gebruikt, Op de modellering van de structuur wordt minder diep ingegaan.

4 . 1 Probleemstelf inq

Getracht wordt het, in figuur 4.1 geschetste, 2D-dynamisch probleem numeriek op te lossen. Het stelt een vereenvoudigd SD-aortaklepmodel voor. Het model bestaat uit starre wanden, waarin een flexibele structuur is opgehangen (een klepvlies)

.

De vloeistofstroming is instationair, en stroomt van links naar rechts. Onder invloed van de instationaire stroming beweegt de structuur. De structuur beïnvloedt op haar beurt de vloeistofstroming.

Fig. 4.1 Het SD-aortamodel.

Het oplossen gebeurt op basis van een eindige e,ementenmethode. Hetgeen inhoudt dat liet probleem in de ruimte wordt gediscretiseerd, waarna de resulterende differentiaalvergelijkingen in de tijd geïntegreerd worden. In navolging van Belytschko(l980) wordt de vloeistof en structuur als

afzonderlijke subsystemen behandeld. Beide systemen worden per tijdstap ontkoppeld beschreven en geïntegreerd. De enige koppeling vindt plaats via

De beschrijvende vergelijking voor zowel de vloeistof als de structuur is de impulsvergelijking,

met p de dichtheid,

2

de snelheid, t de tijd, f een uitwendige kracht per massaeenheid, g

-

een spanningstensor. De vergelijking is gedefinieerd in het 1 a bora tor i uni- coö r d ina t en s te 1 s e 1 (

x

,

y

1

.

4 . 2 . 1 Vloeistof

De vloeistof wordt incompressibel en Newtons beschouwd, met een constante viscositeit. Hiervoor geldt de constitutieve relatie,

met p de druk, q dynamische viscositeit,

1

- eenheidstensor en 0

-

de def ormatiesneltieidstensor

,

Na substitutie in de impulsvergelijking

( 1 )

volyt de bekende Navier- Stokesvergelijking

Samen met de continuîteitsvergelijking

T7.g = o

4 . 2 . 2 Structuur

Voor de structuur kunnen verschillende materiaalmodellen worden gekozen.

Voor een elastisch materiaal geldt bijvoorbeeld

met

E

-

de rektensor,

9

waarbij ui de uitwijking i s ,

( 7 )

4.3 Meshdef l i l i t i e

V o w zowel de vloeistof a l s de structuur moet een meshverdeling worden

gekozen. Indien de knooppunten vast i n het laboratoriumstelse1 liggen, wordt gesproken van een Euler-mesh. R i j een Lagrange-mesh volgen de knooppunten de bewegende materhedeeltjes. Bej.de meshdef j.nitj-es hebben hun voor- en

nadelen. Voor de analyse van een flexibele strucfxur is het gebruikelijk u i t .

t e gaan van een Lagrange-beschri. jving

(convectj.eve kermen verdwijnen)

.

Bi7 grote deformaties kan een Lagrange-mesh echter t:e sterk vervormd worden. Herhaaldelijk opniew definiëren van de mesh

( rezoni.ng ) vereist: veel rekentijd en introduceert onnauwkeurigheden doordat

herhaaldelijke interpolatie naar de nieuwe knooppunten moet. plaatsvinden.

B i j grote verplaatsingen, Z(3al:i ook b i j een vloelstofst~romj.ny, i.s een Euler-

beschrilving beter geschikt.

indien b i j interactieproblemen voor de structuur een Lagrange- en voor de v1oei:itof een Euler-bescXii-ijvi.ng wordt gebruikt, beweegt: de

Lagrange/structuur-mesh v r i j door de Euler/vloeistof -mesh. Daar i n het algemeen de knooppunten van beide meshes ni.et samenvallen

,

is het opleggen

van randvoorwaarden aan het contactoppervlak moeili j k (Belytschko, 1980)

.

Om

deze moeilijkheid ke omzei-len wordt de vloeistofmesh op een door de

gebruiker voorgeschreven wijze bewogen. Gesproken wordt van een a r b i t r a i r e Lagrange-Euler (ALE) mesh (ook wel quasi-Eder me:ih genoemd) (Belytschko,

1982; Donea, 1982). I n een ALE-mesh bewegen i n ù e buurt van h e t

contactoppervlak de vloeistofknooppunten, deze beweging i s gerelateerd aan de beweging van de structuur. Gezorgd wordt dak op het contactoppervlak de vloej.stofknoopp~~nten samenvallen met de structuurknooppunten. De behandeliny

van h e t ccmtac%c)ppervlak i s nu eenvoudig.

4.4 Bekoudswett.en i n de A1,E-bel;chri.ivinq

Definieer

xi.

als de laboratorium (Euler) coördinaten, Xi als de materie (Layrange) cobrdinaten. Daarnaast wordt een referentie-codrdinaten~~steem

xi

gedef i n h e r d I d a t onafhankelijk wordt: voorgeschreven als functie van plaats

en t i j d (Belytschko, 1982)

x ₁

.

wordt o p ieder t i j d s t i p bekend verondersteld.

Een willekeurige functie f (g,t) h e e f t dan drie tijdafgeleiden (Bird, 1960)

-

de materiële ti jdafgeleide (of substantiele t i j d a f geleide)

-

de partiële t i j d a f g e l e i d e

-

de referentie ti jdafgeìeide (of t o t a l e t i j d a f g e l e i d e )

De snelheid van een materiedeeltje wordt gegeven door 3 x;

-

Y

=

a t

X;

= tons+

Definieer analoog een referentie ( g r i d ) snelheid

De r e l a t i e tussen de t i j d a f g e l e i d e n wolcdt gegeven door ( B i r d

,

1960)

Combinatie van beide levert de samenhang tussen de materiële (substantiële) tijdafgeleide en de referentie (totale) tijdafgeleide,

Beschouw een volume V(xi) met oppervlak S en normaal g, vast in de

referentie-ruimte. De integraalformuleringen van de behoudswetten worden dan gegeven door (Thompson 1972 )

-

behoud van massa

-

behoud van impuls

Definieer de Jacobiaan (functionaaldeterminant)

en de verschilsnelheid

De met vgl. ( 1 9 ) resp. (20) corresponderende partiële differentiaal- vergelijking is dan (Donea, 1982)

-

impulsvergelijking

Onder toepassing van de continufteitsvergelijking ( 2 3 ) is de impulsvergelijking ( 2 4 ) te schrijven als

En onder toepassing van de identiteit (Truesdell, 1960)

is de continu1t:eitsveryelijking (23) te schrijven als

Vergelijkingen (25) en (27) zijn de behoudswetten voor impuls

respectievelijk massa in de ALE-beschrijving. Als

xi=

X. (Lagrange-mesh) dan

-

vG=

typerend voor de behoudswetten in de Lagrange-beschrijving. Als

xi=

G

(Euler-mesh) dan

2

= Q ; verkregen worden de behoudswetten gedefinieerd in 1

; Ue cowectievr teraen ver6wijrien uit de behoudswetten. D i t i s ' i het laboratoriumstelsel.

4.5 Eindise elementenbenaderinq

Onder toepassing van de methode van Galerkin kunnen de met de behoudswetten

(251 en ( 2 7 ) in de eindige elementenbenadering corresponderende

vergelijkingen worden verkregen. Deze methode wordt hier uitgewerkt voor een incompressibele, Newtonse vloeistof beschreven in een Euler-mesh. Voor een uitvoerige bespreking wordt verwezen naar Cuvelier(l986).

Definieer een gebied R met rand

r

vast in de Euler-ruimte (xi). Binnen dit gebied stroomt de vloeistof. De Navier-Stokesveryeli~kiny (4) en de

continufteitsvergelijking ( 5 ) beschrijven de vloeistofstroming binnen R ;

4 . 5 . 1 Stokes-probleem

Een extra moeilijkheid bij het oplossen van het stelsel vergelijkingen (28) yeeft de aanwezigheid van de convectieve term (g.V)g i die niet-lineair is.

Omdat slechts op een eenvoudige wijze een lineair stelsel vergelijkingen kan worden opgelost, zal deze term gelineariseerd moeten worden (paragraaf

Beschouw eerst het Stokes-probleem; waarbij de Navier-Stokesvergelijkhg uit

(28) vervangen is door de Stokes-benadering ( 2 9 ) ,

Volgens de eindige elementenbenadering wordt het gebied R onderverdeeld in een eindig aantal elementen. Het snelheidsveld en de druk worden in elk element benaderd met basisfuncties. Dit levert de volgende benaderings- oplossingen

Hierin is vi(xfy) de basisfunctie voor de snelheid, en $ . ( x , y ) de

basisfunctie voor de druk. ïndex i telt de knooppunten af. In totaal zijn er N snelheidsknooppunten en M drukknooppunten. Voor de basisfuncties geldt-

1

Opgemerkt wordt, dat in de benaderingsoplossingen (30) het tijdsafhankelijke deel in de knooppunten wordt meegenomen. M.a.w. per tijdstap moeten de

onbekende knooppuntswaarden ui, ( t )

,

vi

( t ) en pi ( t 1 worden berekend.

De Dirichlet randvoorwaarden en de drukvoorwaarde uit stelsel ( 2 8 ) worden in

(30) verwerkt, er resulteert (de snelheidsknooppunten 1 í 2 í . . . í R liggen op I')

P

Mu,

O

o

Volgens de methode van Galerkin wordt het Stokes-probleem herschreven i

j

2 9 ,

...,

M

-,

N

j

= R t i ,

_ _

( 3 3 )

P a r t i ë l e integratie van ( 3 3 ) en invullen van de interpolatiepolynomen ( 3 1 1 ,

waarbij voor de eenvoud wordt aangenomen d a t g,fx'yIt)=gZ(xiy,t)=O en

p l ( t ) = O I l e v e r t tiet onderstaande s t e l s e l Galerkin vergelijkingen

j=

R i l , .

..,

N (341 ï n mat r ixno ta t i e O M V " O LL4 O O

met submatrices onbekenden

-

a!!

u

=

a t

I

al!

-

-

3t

v =

en rechterlidvectoren

Hogelijk suggereert matrixnotatie ( 3 5 ) ten onrecht.e, dat er vijf

onafhankelijke onbekenden (81QIu,v en p) zi-jnl die uit bovenstaand stelsel

( 3 5 ) kunnen worden opgelost. 2i en Q zijn in de tj.jd gekoppeld met u resp. v. Deze koppeling komt later aan de orde (paragraaf 4.8). Uiteindelijk wordt een sc)ortgelj.jk stelsel als ( 3 5 ) opgelost, met drie onbekenden (u,v en p).

In verkorte notatie wordt stelsel (35) als volgt geschreven

met M: massamatrix, componenten MUU, Mvv S : diffusiematrix, componenten Suu, Svv L: continufteitsmatrix, componenten L J,

-

V: snelheidskolommatrix, Componenten

u,

- P: drukkolommatrix, componenten E

-

F: rechterlidkolommatrix, componenten f,, f, PU' PV

4.5.2 Navier-Stokes Probleem: Newton iteratie

Voor fiet Navier-Stokes probleem moet de niet-lineaire convectieve term

( 1 . v ) ~

worden gelineariseerd, opdat het gediscretiseerde stelsel eenvoudig kan worden opgelost. Dit gebeurt met een iteratiemethode. Bijvoorbeeld de Picard iteratie of Newton iteratie (Cuvelier, 1986). Bij de Newton iteratie wordt er van uitgegaan dat de oplossing op de I-de iteratieslag gelijk is

aan de oplossing op de vorige slag plus een correctie,

Dan volgt voor de convectieve term

Onder verwaarlozing van de kwadratische term in

rkieerde Navier-Stokes vergeli.jking, na de f-de s l a g , er a l s volgt u i t ziet de gelinea-

Volgens de methode van Galerkin volgt dan het onderstaande stelsel lineaire differentiaalvergelijkingen. Voor de eenvoud wordt weer yenomen dat

9, ( x f y f t)=g2(x,y,t)=0 en p, (t)=O. met In

-

1-1

s v

verkorte notatie I

i

NVV, NuV, Nvu

met

N(Y)

: convectiematri.x, componenten NUU,

-

F

en M, S , L, 1, zijn reeds gedefinieerd in stelsel ( 3 6 )

: rechterlidkolommatrix, componenten f,

,

f,, f,, f4

1,

Pf

N(Y) en E zijn tijdsafhankelijk, en Pi, S, L tijdsonafhankelijk.

4.5.3 Boetefunctiemethode

Bovenstaande methode ter oplossing van de (Navier)-Stokesvergelijking wordt de directe methode genoemd; de snelheid en druk worden gekoppeld opgelost en

volgen direct uit de matrixvergelijking. Doordat het op te lossen stelsel groot is, kost dit veel rekentijd en geheugenruimte. Een efficiëntere rekenmethode geeft, de boetefunctiemethode. Het principe hierbij is om de

continufteitsvergelijking te vervangen door (Cuvelier, 1986)

- 1 - 1

met T de boetefunctieparameter (dimensie kgm s 1, een hele grote parameter

T-. Hi.erdoc)r kunnen de drukcmbekenden uit de gcdiscretiseerde Navier-

Stokesvergelijking worden geëlimineerd, zodat de snelheids- en drukbepaling ontkoppeld zijn. Dit wordt kort toegelicht (Cuvelier, 1984) aan de hand van het onderstaande boetefunctiestelsel

(43)

Aangenomen wordt dat de oplossing van dit boetefunctiestelsel voor grote T

convergeert naar die van het oorspronkelijke stelsel ( 2 8 ) . Het boetefunctiestelsel wordt gelineariseerd m.b.v. de Newton iteratie,

r i V

-

t

[ S + N [ v )

-

t

rLTD0'L]y

=

-

F

Volgens de metliode van Galerkin volgt dan het onderstaande stelsel lineaire

differentiaalvergeli.jkingen. Voor de eenvoud wordt weer genomen dat

( 4 7 )

o

M,"

Nvl4

&

L'

_{k P}

5""

+

N""

I:,

7

1,"

-

DPP In verkorte notatie

Het grote voordeel is dat de snelheids- en drukbepaling zijn ontkoppeld. Eerst wordt uit de matrixvergelijking een benaderingsoplossing voor de snelheid geconstrueerd, daarna wordt de druk bepaald. Ten opzichte van de directe methode leidt dit tot kleinere stelsels vergelijkingen, met minder onbekenden (de druk valt eruit). Hetgeen minder rekentijd en geheugenruimte vereist. Een nadeel is dat de boetefunctieparameter T goed gekozen moet

worden. Enerzijds mag T niet te klein zijn, omdat anders de

benaderingsoplossing niet naar de echte oplossing convergeert. Anderzijds mag 't niet te groot gekozen worden. Bij grote T kan de resulterende matrix

f S i

N ( 1 )

i TL D LI singulier worden (Cuvelier, 19851, hetgeen de conditie van het stelsel verslechterd. Bovendien zal door optelling van een matrix met grote getalwaarden, vanwege de machinenauwkeurigheid, cijferverlies optreden. In de praktijk wordt vaak (afhankelijk van het op d e computer beschikbare decimalen) gekozen (Cuvelier 1986)

T - 1

met Q de dichtheid, U een karakteristieke snelheid en L een karakteristieke lengte.

4 ~ 6 n i s c r e t i a a t i e van het; vloei.st;of systeem i n de ALE-beschri-ivinq

Het discretiseren van de bekoudswetten i n de ALE-beschrijving ( v g l . (25) en

( 2 7 ) ) gaat op analoge wijze als get:oond j.n pasagraaf 4 . 5 voor de

bchoudswetten i n de Euler-beschrijving

.

Beschouw hiertoe een gebied R rand

r

vast i n de referentie-ruimte

( x

.

)

.

Bimen d i % gebied st.roc)mt een

incompressibele, Newtonse vloeistof ( v g l .

(2

) )

.

De vloei :;tof si-.roming binnen R wordt beschreven door

*

met:

*

J.

*

G

O p de convectieve term (1-x

1

.Vy na is dit s t e l s e l identiek met skelsel ( 2 8 ) . De convectieve term i s t e schrijven als

G

met gridsnelheid y bekend. x.Vy wordt gelineariseesd met, de Newton

i t e r a t i e . De gelincariseerde Navi.er-Stokesvergelijking z i e t erf na de I-de

*

Het; gebi.ed R

wordt onderverdeeld in een eindig aantal elementen, Volgens de meíAiode van Galerkin volgt dan het onderstaande stelsel lineaire diff erentiaalverge- li jkinyen. V o o r de eenvoud wordt weer genomen dat y

1

(x,

y , t ) = y ₂

(x,

y f t ) =O en

dat. in de tijd bekend wordt veronderst.eld (zie vgl. ( I O ) )

,

p i (t)=O.

M,l4

o

O

u

Op en NG na is _vv dit. s t e l s e l identiek met stelsel (40)

(53)

ïn verkorte notatie

*

* *

Doordat. het integratiegebied R tijdsafhankelijk is ( R = R f t ) ) , zijn ook M,

S en L tijdsafhankelijk, dit i.t. i;. s t e l s e l ( 4 1 ) . Hierbij dient opyemerkt te worden dat in de eindige elementenbenadering, de in de matrixvergelijking voorkomende hteyralen, goed per element kunnen worden bepaald. De

*

k'

met k tellend over de elementen e

Alleen die elementen, die in de tijd van vorm veranderen, leveren een ti3dsafhankelijke elementbi2drage.

Ook hier kunnen, met behulp van de boetefunctiemethode, sneliieids- en

drukbepaling worden ontkoppeld. Voor

2

wordt dan geschreven (zie vg1.(48))

invullen in (55) levert het boetefunctiestelsel

4 . 7 Discretisatie van de structuur

Voor de analyse van de structuur is het gebruikelijk uit t e gaan van een Lagrange-beschrijving (Xi). De convectieve term uit de impulsveryelijkiny

(25) vervalt dan,

De uitwijking van het materiaal i s gedefinieerd als (vy1.(8))

* *

Beschouw een lineair elastisch materiaal lvg1.(9)) met volume R en rand

r

vast in de Lagrange-ruimte (X.). Stel dat

r

is opgedeeld in en

rl

i dus r

bewegen, en dat op

rl

als volgt geformuleerd

* *

* *

* *

**

* *

* *

* *

1.

* *

=

ro

+

r

a*

.

Neem aan dat langs ITo het materiaal niet kan

een oppervlaktekracht werkzaam is. Dit probleem wordt

* *

De Neumann randvoorwaarde 0.n kan opgesplitst worden in een component an in

de richting van de normaal en een component ot in de tangentiële richting.

-

met un=&.n en

ut=&.t

Gediscretiseerd volgens de eindige elementenmethode levert (60) het onderstaande stelsel

met M: massamatrix K: stijfheidsmatrix

-

U: riitwijkingskolommatrix

4 . 8 Ti-iddiscretisatie

In het voorafgaande is de ruimtelijke discretisatie van het beschrijvende stelsel vergelijkingen besproken. Dit gebeurde volgens een eindige

elementenmethode. Het resulterende stelsel Vergelijkingen moet nog in de tijd worden gefntegreerd. Hiervoor wordt een O-differentiemethode gebruikt.

4 . 8 . 1 @-methode Beschouw het stelsel

Gediscretiseerd met de O-differentiemethode, levert (v.d. Vosse, 1987)

Specifieke waarden voor 6 zijn:

O=O: Euler-expliciet ( rechterlid f op het 'oude' niveau)

@=I: Euler-impliciet (rechterlid f op het 'nieuwe' niveau)

8=1/2: Crank-Nicolson (ook wel trapeziumregel genoemd)

Voor Pineaj.re vergelijkingen hebben Euler-expliciet en -impliciet beide een nauwkeurigheid van O(At), Crank-Nicolson û(At 2

1 .

Toepassing van de O-methode voor de tijdintegratie van het boetefunct~iestelsel ( 5 7 ) levert

Met dit stelsel is dus uit 1 op tijdstip tn-l,

X

op tijdstip tn te bepalen. Het oplosproces moet beginnen met (de beginvoorwaarde). De bepaling van

X

G wordt besproken in paragraaf 4.9.1. De druk wordt achteraf bepaald met

Indien tijdstap At zo k1ei.n is, dat

yn-'

heel dicht bij

Yn

ligt, kan gesteld worden dat

W.a.w es wordt niet per tijdstap geftereerd, maar er wordt slechts een

Newton iteratieslag toegepast.

Toepassing van de @-methode voor de tijdintegratie van het lineaire structuursysteem (62) levert

met beginvoorwaarden

U

_-

t = o )

=

o

(A

4 . 8 . 2 Stabiliteit van het vloeistofsvsteem

Om beter te kunnen begrijpen wat de invloed is van @ bij het oplossen van het vloeistof-stelsel (651, wordt het onderstaande stelsel lineaire

differentiaalvergelijkingen van de eerste orde beschouwd (v.d. Vosse, 1987)

Hierin is A een (NxN)-matrix met retile tijdsinafhankelijke coëfficitinten, met N onafhankelijke eigenvectoren met de eigenwaarden Aii' i=1,

...,

N.

Het stelsel differentiaalvergelijkingen (68) èn de volgens de Q-methode gediscsetiseerde benadering hiervan, wordt stabiel genoemd, wanneer een eindige fout go in de beginoplossing

u.

resulteert in een eindige fout g(t)

in de oplossing Q(t) voor alle t. Deze foutvoortplanting is onderzocht door

v.d. Vosse(1987). Stelsel ( 6 8 ) is stabiel a l s (v.d. Vosse, 1987)

In figuur 4 . 2 worden, voor het interval 0<Qcl, de stabiliteitsgebieden

van

AiAt geschetst. Uit de fi.guur blijkt dat voor OLotû.5 het @-schema slechts

vvorwaardelijk stabiel is * Bijvoorbeeld als 8=0 (Euler-expliciet)

< I

Hetgeen impliceert dat At zodanig gekozen moet worden dat AiAt ligt, in het complexe vlak, binnen de cirkel met straal

1

en middelpunt (-1'01. in het geval dat de eigenwaarden van A grout zijn (maar wel negatief), zijn relatief kleine tijdstappen At vereist. Voor O. 5<0(1 is het O-schema voor alle A i A t stabiel (onvoorwaardelijke stabiliteit)

.

I m C X A t l 820.5

t4

9 1 O. 5 O

t-4

Fig. 4.2 Ctabiliteitsgebieden voor de Q-methode voor complexe resp. reële eigenwaarden (v.d. Vosse, 1 9 8 7 ) .

In figuur 4.3 is voor Q=1 resp. Q=O.5 de versterkingsfactor ci uitgezet als functie van Re[AiAt]. Hieruit blijkt dat voor @=O. 5 (Crank-Nicolson) bij gr0t.e negatieve eigenwaarden de verstesking naar -1 gaat, daarentegen gaat voor 0=1 (Euler-impliciet) de versterking naar O . Hetzelfde geldt voor complexe eigenwaarden met een overheersend groot negatief reëel deel. Bij een toenemend j.magj.nair deel van de eigenwaarden neemt: de versterkingsfactor voor de Crank-Mico1sc)nmethc)de toef voor de Euler-implicietmethode neemt deze factor af.

H.a.w. hij grote negatieve eigenwaarden zal het gebruik van Crank-Nicolson verstoringen (bijvoorbeeld ontstaan vanwege een inschakelverschijnsel) nauwelijks uitdempen en aanleiding kunnen geven tot een oscillerende oplossing. (Hierbij wordt opgemerkt, dat door v.d. Vosse(1987) alleen

oscillaties in de drukbepnling zijn gevonden, niet in hek snelheidsveld.) In tegenstelling tot de Crank-Nicolsonmet~iode is het dempend vermogen van de Eufer-Bmplicietmethode (bij grote negatieve eigenwaarden) uitstekend. Maar het; is juist; ongewenst oscillaties

,

die inherent zijn aan de oplossing

(bijvoorbeeld bij de wervelstraat van Von Karman) uit te dempen. Crank- Nicolson verdient dan de voorkeur. Be op1o;jstrategj-e voor zulk een f yslsch probleem is als volgt; begin met B=1 (uitstekende demping van inschakel- verschijnselen), verander na een of twee tijdstappen deze waarde in @=0.5,

Het gebrL1j.k van een boetefunctieparameter heeft grote negatj-eve eigenwaarden tot gevolg (van de orde TAt). Het toepassen van een expliciete Eulermethode eist dan, ten behoeve van de stabiliteit:, zeer kfehe tijdstappen (van de

orde l / r (v.d. Vosse, 1 9 8 7 ) ) . Voor het oplossen van een boetefunctiestelsel

zal dus altijd een impliciete methode worden gekozen, zoals de impliciete Eufer- en Crank-Nicolsonmethode. Bij de keuze van de ene o f de andere methode zal, m. b. t. het dempend vermogen, rekening gehouden moeten worden met het gestelde hierboven.

ReiXAtl ReikAt] IrnC x 8 4 O -4 A t 3 E d e r implicit C 1 O -1 4 O -4 Crank-Nicolson

Fig. 4 . 3 Versterkingsfactor c als functie van AAt voor de

Euler-implicietmethode (a) en de Crank-Nicolsonmethode (b),

VOQT respectievelijk complexe en reële eigenwaarden.

Voor de Mavier-Stokesvergelijking levert het boetefunctiestelsel ( 4 7 ) een

niet-symmetrische matrix met complexe eigenwaarden (v.d. Vosse, 1 9 8 7 ) .

Indien de bijdrage van de convectieve term verwaarloosbaar is (Stokes- benadering), wordt een symmetrische matrix met grote negatieve reële eigenwaarden gevonden. Indien de convectieve term sterk overheerst (zoals bij het oplossen van de wervelstraat van Von Karman)

,

worden eigenwaarden met een dominerend imaginair deel verkregen.

4.8.3 Stabiliteit van het structuursysteem

Stelsel (62) is geldig voor een lineair systeem, zcmùer demping

..

-

P

-i-

-

Toepasshg van een expliciete tijddiscretisatiemethode op stelsel (71) is slechts stabiel als (Hughes, 1978a, 1978b)

(de maximale eigenfrequentie van het systeem) de maximale

met "max

eigenwaarde van (Belytschko, 1979; Park, 1980)

Een afschatting voor wmax wordt gegeven door (Belytschko, 1978)

í 72)

(739

met c de voortpfantingssnelheid van het geluid,

1

de minimale afstand tussen de knooppunten. Dit kan zeer kleine tijdstappen opleveren. Bijvoorbeeld voor water c=1.48*10' m/s, a l s 1=3*110-3

m

dan Qt12*10-6 s .

Vaak wordt dan ook voor de tijddiscretisatie van het structuursysteem een impliciet schema toegepast (Befytschko, 1980). Hughes( 1978a, 1978b) heeft aangetoond, dat toepassing van een impliciete ti jddiscretisatiemethode op stelsel (71 ) onvoctrwaardefijk stabiel is.

3

4 . 8 . 4 Meshpartitie

In principe zijn er twee integrati.eschema's voor de semi-gediscretiseerde vergelijkingen ( (57 f en ( 6 2 ) ) mogelijk: expliciete en impliciete.

Indien de massamatrix diagonaal is, hoeft onder toepassing van een expliciet schema niet een stelsel vergelijkingen worden opgelost. De oplossing op het nieuwe tijdstip volgt via eenvoudige mathematische bewerkingen u i t de

oplossing op het vorige tijdstip. Maar numerieke stabiliteit vereist in het- algemeen zeer kleine tijdstappen. Net een impliciet schema moet wel een stelael vergelijkingen worden opgelost (ook als de massamatrix dj.agonaa1

is). Omdat in het algemeen voor lineaire problemen een impliciete methode onvoorwaardelijk stabiel is, kunnen in tegenstelling tot een expliciei; schema daarentegen veel grotere tijdstappen worden toegepast.

Voor de subsystemen kunnen verschillende tijdinteg~atieschema's wmden gekozen (meshpartitie

1

; bijvoorbeeld expliciet-expliciet (E-E) expliciet- impliciet (E-ï) en impliciet-impliciet (I-I) (Belytschko, 1980). Voor een &en-dimensionale meshverdeling wordt een en ander in figuur 4.4 grafisch voorgesteld.

t

E - E with Interpolation (e) t

Explicit

-

Implicit

n +1 n + i

n - X n -X

I

-

I with Extrapolation (@I

I

n + i

n - - r X

Fig. 4.4 De 'informatiestroom' in de expliciet-impliciet(E-I),

expliciet-expliciet(E-E~ en impliciet-impliciet(1-1~ par- titie (Eelytschkü, 1980).

Bij een 'E-E partitie wordt zowel voox de vloeisto£ als voor de structuur een expliciet tijdintegratieschema gebruikt. Ten behoeve van de stabiliteit is de tijdstap voor de structuur i.h.a. kleiner dan voor de vloeistof. Eerst wordt de integratie met de grootste tijdstap uitgevoerd, vervolgens kan de integratie met de kleinere tijdstap plaatsvinden. Dit tweede subsysteem wordt zo vaak gefntegreerd totdat het tijdstip is bereikt waarop de oplossing van het eerste subsysteem bekend is. Voor de tussenliggende

tijdstippen j.s via interpolatie eenvoudig de randvoorwaarde voor het tweede subsysteem vast te stellen.

Wanneer een E-I partitie wordt toegepast, wordt het eerste subsysteem expliciet: gefntegreerd. De oplossing op het contactoppervlak levert de

randvoorwaarde voor de impliciete integratie van het tweede subsysteem. De

tijdstap voor ieder van de subsystemen i s i n principe even groot:.

B i j een

1-1

p a r t i t i e worden bei de subsystemen impliciet gelntegreerd

.

Hierhi. j doek zich een probleem voor b i j het opleggen van de randvoorwaarde aan liei, contactoppervlak. Door middel van extrapolati e z a l de randvoorwaarde voor het eerst: t e integreren subsysteem vastges t:eld moeten worden. De i n tiet:

contactpunt geëxtrapoleerde waarde i s niet de eindwaarde, maar wordt berek.end en aangepast a l s het tweede subsysteem iii de t i j d worclt

gelntegreerd

.

4.8.5 Ctabi1i.tei.t van het qekoppelde systeem

E i 3 een X-f en E-E p a r t i t i e i s niet voldoende informatie aanwezig om de

dif~erenl:i.aalvergeli jkingen op t e kunnen lossen. R i j een

1-1

partit3.e is een extrapolatie j.n de kmoppunten op het contactoppervlak v e r e i s t , b i j een E-E

pari5.tj.e een h t e r p o l a t i e ( z i e paragraaf 4 . 8 . 4 1 a Deze extrapolaties en interpolaties befnvloeden de stabi.1itei.t van de t i j d i n t e g r a t i e (Belytschko

,

1 9 7 9 ) . Hiernaar z a l gekeken moeten worden. B i j een E-J p a r t i t i e (waarbij de

1;ijdstap voor beide partj.i:ies even groot: i s ) i s voldoende Informatie

aanwezig om h e t gekoppelde systeem op t e lossen. Door verscheidene auteurs is aangetoond d a t een E-l: p a r t i t i e voorwaardelijk :;tabiel i s (hi.jv.

Belytschko, 1978; Hughes, 1978a). De t i j d s t a p voor de t i j d i n t e g r a t i e van de t o t a l e E-1 p a r t i t i e wordt bepaald door de vereiste t i j d s t a p voor de

e x p l i c i e t e p a r t i t i e .

Vanwege de numerieke s t a b i l i t e i t wordt het, i n d i t hctofdstuk besproken,

v l ~ e i s t ~ f s y s t u e ~ , ( 5 7 ) en s t r ~ ~ t ~ u r ~ y z t e e m ( 6 2 ) i ~ p l i c i e t g e h t e g r e e r d ( z i e

paragraaf 4 . 8 . 2 en 4.8.3). Waar numerieke s t a b i l i t e i t van de afzonderli.jke

subsystemen garandeert nog niet (vanwege de beiiodiycle extrapola%ie) een s t a b i e l e ï-1: parti.t:ie

een

1-1

p a r t i t i e achaars.

Volgens Park(1980) is, voor een grote klasse vloeistof-structuur interactieproblemen (waarbij beide subsystemen lineair z i j n ) , een

ï-1

p a r t i t i e onvoorwaardelijk s t a b i e l

,

mits een geschi.kte extrapolntie voor de

randvoorwaarden wordt gekozen. Belytschko ( 1979) k e f t aangetoond d a t voor

een bepaalde klasse vloei.stof-structuur interactiemodellen (waarbij zowel

voor de vloeistof a l s voor de structuur een l i n e a i r systeem zonder demping is genomen (vgl. ( 7 1 ) 1

,

de ï-I: p a r t i t i e onder de extrapolatie

n = u n+

1

vloeistof structuur

U

onvoorwaardelijk instabiel i s . Park( 1980) heeft voor het j.nteracI;i.emodel van

Eelytschko ( 1979) een extrapolatie gefndentificeerd, die een onvoorwaardelijk

stabiele

1-1:

p a r t i t i e oplevert;

n-

1

= 2un - u

n+

1

vloeistof structuur structuur

U

Voor zover o p ( l i b moment t e beoordelen i s , v a l t het i n d i t hoofdstuk

behandelde vloeistof-structuur interactiemodel niet onder de docx Relyt:jchko ( 1979 ) gefdentif iceerde klasse

mogelijk i s een u i t s l u i t s e l t e geven c i f h e t behandelde vloeistof-structuur

interactiemodel numeriek s t a b i e l is en zo j a , onder welke extrapolatie, Dil:

vereisi; nog nader onderzoek. Bewi-jsvc)ering van Park(l98O) i s n i e t begrepen. Behandeling conf:actoppervlak door Belytschko ( 1 9 7 3 ) is nog onduj.delijk.

4 9 PrakZ.isclie problemen en onduideli-ik.heden

:In de A&E-f ormulering beweegt de vloeistofmesh, deze beweging is gerelateerd aan de beweging van de structuur. Per t i j d s t a p wordt de beweging van de mesh berek.end u i t de beweging van de structuur. Gezorgd wordt dat op he%

contactoppervlak de vloeistofknoopp~~nten samenvallen met de s.tructuurknoop- punten, daarnaast worden interne vloeistofknooppunten zodanig bewogen dat lokale vecstorj.nyen van de mesli zo klein mogelijk blijven.

De bedoeling i s zowel het vloeistof systeem als het structuursysteem implM.et t e integxeren a Hierbij wordt; het vloeistof sy:iteem het eerst;

gehtegreerd, daarna volgt h e t structuursysteem. Aan het vloekkofsysteem worden kinematische randvoorwaarden voorgeschreven, aan het structuursysteem dynamische.

4 . 9 . 1 Automatj sche aanpassing van de vloeistofmesh

Voor praktische implementatie van de ALE-beschrijving

,

is een automatische

aanpassing van de vloeistofmesh v e r e i s t . U i t de oplossing van het;

structuursysteem (op t i j d s t i p t ) is de gridsneliieid van de vlnej.:itofme:jh

f o p t i j d s t i p tn) af t e leiden. Beschouw

,

i n onderstaande f i g u u r

,

vloei s t o f -

knooppunt i:

.

n

A

F i g . 4 . 5 Aanpassing van de vloeistofmesh.

Belytschko ( 1982) s c h r i j f t de yridsne1hej.d van de Interne knooppunten al:;

mei;

(de u i t w i j k i n g van kiicrcrppunt B op het contactoppervlak) volgt u i t de

oplossing van het structuursysteem.

Donea(1382) s t e l t een veel ingewikkelder schema voor,

met N het aantal knooppunten die met knooppunt

1,

v i a de elementzijden en

dj.agonalen, z i j n verbonden. L i s de afstand tussen knooppunt

r

en het hieraan verbonden kncmppunt J,

fi

i s de verplaatsing van de knocrppuiiten. De

eer:;te term geeft aan, dat de gridsnelheid van knooppunt T o p ti.jtl:jt;ip tn

g e l j . j k genomen wordt aan de gemiddelde gridsnelheid van de naburige knooppunten o p liet vorige t i j d s t i p tn-'. De tweede term is een ccrrrectie hierop. De correctj.eterm verhoogt de gridsnelheid indien de af stand tt.i:i>ien

twee natruriye knooppunten t e klein dreigt t e worden.

De aanpassing van de knooppuntscoBrdinaten van de voeistofmesfi volgt u i t i

I J

Een andexe mogelijkheid i s om per t i j d s t a p een nieuwe vloeistofmesh t e

maken (rekening houdende met. de u i t w i j k i n g van üe structuur), waarbij alleen

de c::dxdi.r,aten van be knooppunten veranderen, Uit het verschj.1 van de oude

(3") en nieuwe

(zn"

1 meshknooppunten, volgt de gridsnelheid

xGf

vergelijking ( 7 7 ) ) ,

( z i e

Voor de modellering van een hartklepprothese i s een zuivere ALE-mesh

waarschijnlijk niet toepasbaar. Dit doordat tijden25 het openen en sluiten van de klep, het vlcreistof domein een t e grote vormverandering ondergaat. Deze i s met: een ALE-mesh moeilijk t e behandelen (Horsten, 1 9 8 7 h ) , Beschouw bijvoorbeeld i n f i g u u r 4.5 de vloeistof elementen i n de doorstrocrmopening

onehdig klein. Om de vc~ordelen van de ALE-mesh t o c h t.e kunnen beì~c~ucten, i s

een combi.na.l;ie van ALE-mesh met rezoning de meest voor de hand liggende k.euze. ne AIIE-%echiii.ek li2k2, wel geschikt voor de uitwerking van

testproblemen, waarbij het v l i e s r e l a t i e f kleine Verplaatsingen ondergaat ~

4 . 9 . 2 Convectieve kerm

G

De convectieve term (y-y 1 . V x wordt, z o a l s aangegeven i n paragraaf 4 . 6 ,

gelineariseerd. Hieruit volgt de convectiematrix

N(1-Y

1

,

waarin j! bekend

wordt verondersteld. Na ti jddiscretisatie wordt, de term N(x -y

verkregen ( z i e v g l . (6'5)

.

Doordat het vloeistofsysteem tiet eerst wordt gefntegreerd, i s

YG,"

onbekend. Echter indien t i j d s t a p A t zo klein is, dat

1

G G

n G,n _)yn

n-

1

_{heel d i c h t b i j}

_xn

_{l i g t , kan gesteld worden d a t}

Hetgeen inhoudt: d a t er slechts &n Newi:on j.teratieslag wordt toegepasi;

,

c1aarnaa:jt: i:; yri.dsne1liej.d

YGf

bekend,

0 P Om s t e l s e l ( 6 5 ) op t e kunnen lossen, moet het oude snelheidsveld ( j ! n-1)

de nieuwe knooppunten ten t i j d e tn bekend z i j n . Deze knooppunteii bewegen met. een snelheid yGín. Doordat verondersteld wordt dat Gfn-' heel dicht b i j YGf

nieuwe knoi)ppunten

xn

g e l i j k i s aan het, op de c)ude knooppunten &"-li

berekende snelheidsveld Yn-' d a t mek een snelheid

Y

l i g t , kan worden aangenomen d a t het oude snelheidsveld heersend op de i s verplaatst..

4 . 9 . 3 Dirichlet. randvoorwaarde voor de vloeistof

Vanwege de 'no-slip' conditie voor viskeuze vloeistcjf fen geldt, d a t aan het contactoppervlak (i' ) de snelheid van de vloeistof (y ) g e l i j k is aan de

snelheid van tie structuur

(x')

(kinematische randvoorwaarde)

,

F

conk

E S

-

Doordat hef: vloeistof systeem het eerst wordt yefnteyreerd

extrapolatie i n de t i

$1

van de kinematische randvoorwaarde voor het

vloeistof systeem v e r e i s t . De eenvoudigste keuze i s

i s een

Vanwege de s t a b i l i t e i t van het t o t a l e oplosproces keuze noodzakelijk ( z i e paragraaf 4 . R . 5) ; bijvoorbeeld

i s mogelijk een andere

I I

I’cont “cont 1 rcont

4.5.4 Neumann randvocrrwaarde voor de structuur

De integratie van het structuursysteem volgt na de integratie van het

vloeistofsysteem. De dynamische randvoorwaarde

-

0.n {vergelijking { 6

1

) 1 wordt

bepaald uit de oplossi.ng (y) van het vloeistofsysteem,

met v _n =x.g en vt=x.t;

Het is wog onduidelijk of b i j de bepaling van de spannl.i~gen

(e’!)

rekmin.;

gehouden moet worden met de beweging van de vloeistofmesh (E G

1 .

Dit: vereist:

4. IO A l qo r i tm e

Een mogelijk algoritme voor het oplossen van het vloeistof -structuur

interact:ieprohleem i.n de ALE-formuleri.ng i zou het onderstaande schema kunnen

z i j n ( t i j d s t a p A t i s voor beide subsystemen even g r o o t ) .

1 )

Per t i j d s t a p n wosdt het v1oeistofsysl;eem ( 6 5 ) opgelost, verkregen wordt: de sne1hej.d yn op de knooppunten.

2)

Bereken u i t di.1; snelheidsveld de normaal- en schuif spanni-ngen, die op de structuur werken (paragraaf 4.9.4).

3 1 Het: stsuctuursysteem ( 6 7 ) wordt opgelost:,

verkregen wordt de u i t w i j k i n g

an

o p de knooppunten. 4

1

Bepaal de gridsnelheid

2"

van de vloeistof kmoppunten

,

(paragraaf 4 . 9 s 1 1 .

5) Pas de cobrdinaten van de vloeistofmesh aan;

6 ) Extrapoleer i n de tAjd de kinematische randvoorwaarde voos het:

vloeistof systeem aan het, contactoppervlak. u i t de u i t w i j k i n g

van de structuur (paragraaf 4 . 3 . 3

.

7 ) Ga verder met de volgende t i j d s t a p n = n t l .

Voor stap

1

) zal de boetefunctiemetlic?de uit SEFRAN aangepast: moeten worden

aan de ALE-formulering

.

2)

i s binnen CEPRAN uitvoerbaar. Vc>c)r een f l e x i b e l v l i e s i.s st:ap 3 ) standaard met SEPKAN niet mogelijk. Voor 4) en 5 ) zal

mogelijk per geval een geschikte oplossing gekozen moeten worilen. Over h e t j u i s t e schema voor 6 1 bestaan nog fundamentele onzekerheden.

<

4 . I 1 Updated Euler me%kode

Een andere mogelijk.heid om de klepvliesbeweging numeriek op t e lossen, is de

updated Euler methode. Hierbij wordt het vloeiskof systeem opgelost i n een Euler-mesh. De stxuctuur bi Ijvoorbeeld i n een Tagrange-mesh

.

Na af loop van een ti.jdstap wordt de vloeist(>fme>jh aangepast aan de beweging van de structuur. Hierbij mogen nieuwe elementen worden yecreeerd

.

Per kiljdstap wordt gezorgd dat op het contactoppervlak de vloeistofknooppunten samen- vallen met de structuurknooppunten. B i j deze methode skaat de vloeistofmesh

per t i j d s t a p dus s t i l (dit i n tegenstelling t o t de AtE-mesh)

.

Ue

wisselwerking tussen vloeistof en structuur geschiedt op geheel analcryc wijze als beschreven i n paragraaf 4.9.3 en 4.9,4. Waar voor het nplossen van

de Navier-Stokes vergelijking i n de nj.euwe Euler-mesh i s echter wel een

interpolatie van het; discrete snelkiiìsveld van de oude naar de nieuwe

kncloppunten v e r e i s t .

Beschouw het i n de ruimte en i n de t i j d gediscretiseerde s t e l s e l ( 6 5 )

(waarbii gridanelheid

'

'

Y

g e l i j k is aan n u l ) . De kiiooppuiiten waarop snelheid

p-

I

worden bepaald. Om s t e l s e l ( 6 5 ) op t e kunnen lossen, moet liet oude

snelheldsveltl

(Yn-')

o p de nieuwe knooppunten ten t i j d e tn bekend z i j n . D i t

v e r e i s t ini:erpolatie van h e t discrete snelheidsveld yn-' , Hiervoor kunnen de

interpolatiepolynt~men ( 30) worden gebruikt:.

Aan een interpolatie z i j n onnauwkeurigheden inherent. Doordat per ti jdctap een interpolatie wordt toegepast, zal d i t tot een cumulatie van

oniiauwk.eurigheden leiden f d i t i n tegenstelling t o t de ALE-f ormulering )

.

i s bepaald vallen niet samen met de knooppunten waarop snelheid

Xn

moet,

-

Een alyoritnie wor het oplosser, van het v?veistof-stïuctUur fnte*âctiepnj-

bleem i.n de updated Euler methcxie, zou er a l s volgt u i t kunnen z i e n ( t i j d s t a p A t i s voor beide subsystemen even grcjot) ;

i 1 Per t i j d s t a p n wordt iiet vïoeistofsysteem ( 6 s ) opgelost (met EG=o) i verkregen wordt de snelheid

xn

op de knooppunten.

2)

Bereken u i t d i t snelheidsveld de normaal- en schuifspanningen, d i e