1. No category

Det Teknisk Naturvidenskabelige Fakultet
Aalborg Universitet
Titel: Stivhedsanalyse af aluminium – Virkelighedens teori eller teoriens virkelighed?
Tema: Analyse og design af bærende konstruktioner
Projektperiode: 7. semester, 2. september 2005 - 23. december 2005
Projektgruppe: C 128
__________________________
Mads-Peter Hansen
__________________________
Søren Mosegaard Hansen
__________________________
Per Højbjerg Kristensen
__________________________
Ivar Chr. Bjerg Pedersen
__________________________
Albert Pætursson
__________________________
Lars Saaugaard
Vejleder: Rune Brincker
Oplagstal: 8
Sideantal hovedrapport: 116
Sideantal appendiks: 37
Synopsis
I dette projekt undersøges spændingerne og
flytningerne for en række cirkulære homogene og porøse forsøgsemner udsat for to
modstående jævnt fordelte belastninger.
Forsøgsemnerne er udført i aluminium.
Idet geometrien af forsøgsemnerne er af en
sådan art, at simple beregningsmetoder ikke
er anvendelige, opstilles der både analytiske
og numeriske modeller til bestemmelse af
spændinger og flytninger.
For at vurdere anvendeligheden af de analytiske og numeriske modeller sammenlignes
resultaterne med spændinger og flytninger
bestemt eksperimentelt.
Til bestemmelse af spændingerne og flytningerne i de homogene og porøse forsøgsemner indgår hhv. homogene og porøse
materialeparametre. De homogene er bestemt eksperimentelt ved måling på en homogen klods, hvorimod de porøse er bestemt ved at sammenligne forsøg på en porøs skive med analytiske og numeriske estimater af materialeparametrene.
1
2
Forord
Rapporten er udarbejdet på baggrund af temaet Analyse og design af bærende konstruktioner, og som et led i uddannelsesforløbet for 7. semesters studerende ved det Teknisk- Naturvidenskabelige Fakultet på Aalborg Universitet under Instituttet for Bygningsteknik.
Målgruppen for rapporten er læsere med et grundlæggende kendskab til overbygningen
indenfor ingeniørfaget konstruktion.
Rapporten omfatter de væsentligste forudsætninger, der ligger til grund for de anvendte
teorier samt resultater, der enten er bestemt eksperimentelt eller ved programmering i Matlab. Efter rapporten følger et appendiks, hvori forsøgsbeskrivelser og resultater for de udførte forsøg er nærmere beskrevet. Henvisninger til appendiks er indført i rapporten, hvor
det findes nødvendigt med uddybende information. Endvidere er der vedlagt en CD-rom
med ”CD-bilag” indeholdende bl.a. programmeringer i Matlab, beregninger i Excel regneark og rapporten i pdf-format. Ønskes der kendskab til de udførte beregninger, henvises
der i de enkelte afsnit til CD-bilag.
Figurer og tabeller er i rapporten nummereret fortløbende i hvert afsnit og ledsaget af forklarende tekst og evt. kildehenvisninger. Eksempelvis er den tredje figur i kapitel 6 i rapporten angivet som figur 6.3. Tilsvarende er figur 6 i appendiks A angivet som figur A.6.
Formler anvendt i rapporten er ligesom figurer og tabeller fortløbende nummeret. Således
angives formel nummer 6 i afsnit 4 som (4.6).
Ved henvisning til ”CD-bilag” angives henvisningen ved matlab-fil og excel-fil for hhv.
Matlab programmeringer og Excel regneark efterfulgt af et nummer.
Kildeangivelsen af rapporten er opdelt i faglitteratur og øvrige. Bagerst i rapporten forefindes en litteraturliste med tilhørende kildehenvisninger. Placeres kildehenvisningen før et
punktum henviser den til den foregående sætning, og placeres den efter et punktum henviser den til det foregående afsnit. Såfremt der henvises direkte til en kilde i teksten skrives
eksempelvis ” …ifølge Concepts and application of finite element analysis…”, og henvisninger forekommer derfor i to former. Enten som direkte refererede inde i teksten eller som
kildereference efter en sætning eller et afsnit. Kildehenvisningerne og deres referencer er
angivet på følgende måde:
Faglitteratur
Henvisning:
Litteraturliste:
Øvrige
Henvisning:
Litteraturliste:
[Forfatter efternavn et al, Udgivelsesår]
Forfatter efternavn, Fornavn
Titel, Forlag, udgivelsesår, ISBN:
[Efternavn, år]
Efternavn, fornavn
Udgiver, dato/år
I rapporten er der til beregningerne anvendt beregningsprogrammerne Microsoft Excel,
Matlab samt Calfem, der er en toolbox til Matlab.
3
Indhold
Den vedlagte CD-rom indeholder følgende:
• Rapport og appendiks i pdf-format
• Microsoft Excel regneark anvendt til beregningerne
• Matlab programmering til modellering af modellerne
I rapporten er anvendt følgende notation:
Symbol
[]
{}
[ ]T
[ ]-1
,/
∂
∂u
; eksempel u , x =
∂x
∂x
Forklaring
Rektangulær eller kvadratiske matricer
Søjlematrice
Transponeret matrice
Inverteret matrice
Partiel differentiation
Udover ovenfor nævnte notation er der i afsnit 2.2.4 og 2.2.5 anvendt indeksnotation, hvor
latinske indeks i antager værdierne 1, 2 og 3, og græsk indeks α antager værdierne 1 og 2.
4
Indhold
1 Indledning .......................................................................................................... 9
1.1
Beskrivelse af forsøgsemner.................................................................. 10
2 Estimering af materialeparametre.................................................................. 13
2.1
Beskrivelse af enhedscelle..................................................................... 13
2.2
Analytisk estimering af materialeparametre ........................................... 15
2.2.1
Ashby og Gibsons cellulære metode .............................................. 15
2.2.2
Voigt estimat ................................................................................... 17
2.2.3
Reuss estimat ................................................................................. 19
2.2.4
Dilute estimat .................................................................................. 20
2.2.5
Den selvkonsistente metode ........................................................... 25
2.3
Numerisk estimering af materialeparametre .......................................... 26
2.4
Sammenligning af materialeparametre .................................................. 29
3 Analytiske modeller ........................................................................................ 33
3.1
Analytisk løsning af cirkelskive............................................................... 34
3.1.1
Airy’s spændingsfunktion ................................................................ 36
3.1.2
Centralt belastet cirkelskive ............................................................ 37
3.1.3
Jævnt belastet cirkelskive ............................................................... 51
3.1.4
Sammenligning af enkelt last og jævnt fordelt last .......................... 57
3.1.5
Opsamling....................................................................................... 59
3.2
Analytisk løsning af cirkelring ................................................................. 59
3.2.1
Flytningsfeltet.................................................................................. 59
3.2.2
Beregning af konstanter.................................................................. 62
3.2.3
Resultater ....................................................................................... 64
3.2.4
Opsamling....................................................................................... 64
3.3
Delkonklusion......................................................................................... 65
5
Indhold
4 Teori for elementmetoden............................................................................... 67
4.1
Indledende betingelser for FEM ............................................................. 67
4.1.1
Konstitutiv betingelse ...................................................................... 67
4.1.2
Kinematisk betingelse ..................................................................... 68
4.1.3
Statisk betingelse ............................................................................ 69
4.1.4
Kompatibilitetsbetingelse ................................................................ 70
4.2
Interpolation og formfunktioner............................................................... 70
4.2.1
4.3
Potentiel energi i elastisk legeme ........................................................... 72
4.3.1
4.4
Pascals trekant................................................................................ 72
Spændingsberegning ...................................................................... 74
Elimination af frihedsgrader.................................................................... 74
5 Elementtyper.................................................................................................... 79
5.1
Isoparametrisk CST-element Q3............................................................. 79
5.2
Isoparametrisk Q8-element..................................................................... 83
5.2.1
Jacobi Transformation..................................................................... 85
5.2.2
Numerisk integration ....................................................................... 86
5.3
Isoparametrisk Q4-element..................................................................... 88
5.4
Elementsvagheder ................................................................................. 89
5.4.1
Spændingsberegning ...................................................................... 90
6 Numeriske modeller ........................................................................................ 91
6.1
Konvergensanalyse................................................................................ 92
6.2
Numeriske resultater .............................................................................. 95
6.2.1
Cirkelskive....................................................................................... 95
6.2.2
Cirkelring ......................................................................................... 96
6.3
Porøs cirkelskive opbygget af superelementer....................................... 98
6.3.1
6
Opbygning af superelement ............................................................ 99
Indhold
6.3.2
Konvergensanalyse .......................................................................103
6.3.3
Bestemmelse af flytninger..............................................................105
6.4
Delkonklusion........................................................................................106
7 Konklusion......................................................................................................107
7.1
Materialekonstanter...............................................................................107
7.2
Flytninger ..............................................................................................108
7.3
Spændinger ..........................................................................................110
8 Litteraturliste ..................................................................................................115
7
8
1
Indledning
Undertiden kan det blive nødvendigt at anvende konstruktionselementer af en sådan udformning, at almindeligt kendte bjælketeorier ikke giver en brugbar løsning ved dimensioneringen af disse. Endvidere kan konstruktionselementerne være opbygget af et materiale
med en kompleks struktur så som beton, som også kaldes et kompositmateriale. Disse materialer indeholder som bekendt flere forskellige materialer med hver deres egenskaber, og
idet det er yderst besværligt at modellere disse ned i mindste detalje, er det fordelagtigt at
kunne bestemme en gennemsnitlig materialeegenskab for f.eks. beton. Til dette findes der
flere forskellige estimeringsmetoder, som bl.a. vil blive gennemgået i denne rapport.
I dette projekt betragtes hhv. en homogen cirkelskive og cirkelring samt en porøs cirkelskive og cirkelring, der alle er udført i aluminium. Følgende vil disse under et omtales som
forsøgsemnerne.
Normalt forbindes et porøst materiale med et materiale, hvori der findes et stort antal porer,
der enten indeholder luft eller vand, hvis dimensioner er så små, at de ikke kan ses med det
blotte øje. I denne rapport har ordet porøsitet en lidt anden mening, idet de i rapporten behandlede aluminiumsemner, hvori der er boret 8 mm huller, opfattes som emner af et porøst materiale.
For forsøgsemnerne ønskes det at bestemme spændingerne og flytningerne forårsaget af to
modstående lastpåvirkninger. Dette lader sig ikke umiddelbart gøre ved simple beregningsmetoder, hvorfor der opstilles både analytiske og numeriske modeller. Formålet er at
sammenligne de analytiske og numeriske modeller med eksperimentelt bestemte flytninger
og spændinger, og ud fra sammenligningen vurdere nøjagtigheden af modellerne, hvor
usikkerheden på eksperimenterne tages i betragtning.
Der er ikke bestemt spændinger i samtlige forsøgsemner, hvorfor der i tabel 1.1 skabes
overblik over hvilke metoder, der anvendes til at bestemme spændinger hhv. flytninger i de
forskellige forsøgsemner. Af tabel 1.1 ses det f.eks., at det er muligt at sammenligne spændingerne i den homogene cirkelskive, da disse er bestemt ved de to modeller samt eksperimentelt. Derimod er det ikke muligt at sammenligne spændingerne i den homogene cirkelring, da det kun er for den numeriske model disse kan bestemmes. Med udgangspunkt i
tabel 1.1 kan en afsluttende sammenligning af flytnings- og spændingsresultaterne foretages.
9
1.1 Beskrivelse af forsøgsemner
Tabel 1.1. Angiver hvilke metoder, der kan anvendes til at bestemme flytninger hhv. spændinger i forsøgsemnerne. Ja hhv. nej angiver om der er bestemt / ikke bestemt flytningerne og spændingerne.
Analytisk
Flytninger
Spændinger
Numerisk
Flytninger
Spændinger
Forsøg
Flytninger
Spændinger
Homogen
cirkelskive
Ja
Ja
Ja
Ja
Ja
Ja
Homogen
cirkelring
Porøs
cirkelskive
Porøs
cirkelring
Ja
Nej
Ja
Ja
Ja
Nej
Ja
Nej
Ja
Nej
Ja
Nej
Ja
Nej
Ja
Nej
Ja
Nej
I de efterfølgende modelberegninger er kendskabet til aluminiumets materialeparametre
nødvendig. For det homogene materiale er materialeparametrene bestemt eksperimentelt,
og forsøgsbeskrivelse samt resultater fremgår af appendiks B.1. Det porøse materiales materialeparametre er derimod bestemt ved en række analytiske estimater samt numerisk og
eksperimentelt. Det er valgt indledningsvis i rapporten at præsentere de analytiske og numeriske estimater, der afslutningsvis sammenlignes med de eksperimentelt bestemte porøse materialeparametre, for hvilke forsøgsbeskrivelse og resultater er at finde i appendiks
B.2. I sammenligningen vælges de porøse materialeparametre, der arbejdes videre med i
modellerne for de porøse forsøgsemner.
1.1 Beskrivelse af forsøgsemner
Som beskrevet i afsnit 1 er formålet at bestemme flytningerne og spændinger i de fire forsøgsemner. Følgende beskrives forsøgsemnernes geometri, og hvorledes flytningerne bestemmes eksperimentelt. Endvidere beskrives forsøgsemnerne anvendt til bestemmelse af
de homogene og porøse materialeparametre, hvilket gøres indledende.
Figur 1.1. Billeder af hhv. den homogene klods (tv.) og den porøse kvadratiske skive (th.), der er anvendt til
bestemmelse af materialeparametrene.
10
1 Indledning
De homogene og porøse materialeparametre bestemmes ud fra forsøgsemnerne vist på
figur 1.1, hvor flytningerne af klodsen og den porøse kvadratiske skive hhv. bestemmes
vha. straingages og flytningsmålere. For en nærmere beskrivelse henvises til appendiks B.1
og B.2.
De fire forsøgsemner, for hvilke de analytiske, numeriske og eksperimentelle flytninger
sammenlignes, fremgår af figur 1.2. Det bemærkes, at porøsiteten både for cirkelskiven og
cirkelringen er indlagt som et kvadratisk net placeret ens på de to forsøgsemner.
ry
ry
b
L
b
L
c
c
a
ry
ry
b
ri
ri
L
b
L
L
L
a
Figur 1.2. Homogene og porøse forsøgsemner med angivelse af længder og radier.
11
1.1 Beskrivelse af forsøgsemner
Af tabel 1.2 fremgår de fire forsøgsemners geometriske størrelser.
Tabel 1.2. Forsøgsemnernes geometriske størrelser.
Beskrivelse
ry
ri
a
b
c
t
L
12
[mm]
[mm]
[mm]
[mm]
[mm]
[mm]
[mm]
Ydre radius
Indre radius
Radius af hullerne
Indbyrdes vandrette og lodrette afstand mellem centrum af hullerne
Vandrette og lodrette afstand fra centrum af emnet til et hul
Tykkelsen af forsøgsemnerne
Afstanden mellem målepunkterne for flytningsmålerne
Størrelse
79
24
4
16
8
20
96
2
Estimering af materialeparametre
I dette afsnit estimeres materialeparametrene for det porøse materiale ved en række analytiske modeller samt én numerisk model. Resultaterne for disse indgår slutteligt i en sammenligning med materialeparametrene bestemt eksperimentelt, og ud fra sammenligningen
vurderes hvilke af de bestemte materialeparametre, der anvendes efterfølgende.
Den analytiske estimering af materialeparametrene gøres iht. følgende fem metoder
• Ashby og Gibsons cellulære metode
• Voigts estimat
• Reuss’ estimat
• Dilute estimatet
• Den selvkonsistente metode
Estimeringen af materialeparametrene på et numerisk grundlag foretages ved anvendelse af
Finite Element Method (FEM), som er en approksimativ løsning, der for en finere og finere
diskretisering af elementer konvergerer mod det eksakte. Beregningerne foretages i computerprogrammet Matlab.
2.1 Beskrivelse af enhedscelle
Estimeringen af materialeparametrene bygger generelt set på et repræsentativt volumen
element (RVE), hvad enten det gælder de analytiske eller numerisk estimater, dog med
undtagelse af Ashby og Gibsons cellulære metode. Derfor opstilles der i det følgende et
RVE, som er gældende for alle estimaterne. RVE anvendes oftest i forbindelse med kompositmaterialer til beskrivelse af hvilke delmaterialer de består af samt i hvor stor en
mængde de findes i kompositmaterialet. Det materiale, der er det fremherskende, kaldes
for matrixmaterialet.
RVE udgør pr. definition en så lille del af materialets samlede struktur, som det realistisk
er muligt, hvilket vil sige, at hvis der foretages en detaljeret analyse af RVE, så skal det
beskrive materialets egenskaber på makroskopisk niveau.
RVE opstilles for det tidligere omtalte porøse materiale, som er en aluminiumsskive med
huller, hvorfor det må antages, at kompositteorien kan anvendes tilnærmet. Det fremgår af
figur 2.1, at den porøse aluminiumsskive er opbygget med centrum af hullerne i et kvadratisk net. RVE bestemmes derfor som et udsnit af skiven, der repræsenterer et hul placeret
13
2.1 Beskrivelse af enhedscelle
dobbelt symmetrisk i et kvadratisk udsnit. Det vælges i det følgende at betegne RVE for
enhedscellen.
Enhedscellen er ud fra hullernes indbyrdes afstand i prøvelegemet bestemt til at have en
sidelængde på 16 mm, og hullerne en diameter på 8 mm, figur 2.1. Tykkelsen ud af planet
er 20 mm.
b = 16
y
x
n
64
b = 16
M
a=4
64
Figur 2.1. Enhedscelle anvendt til estimering af materialeparametrene for et porøst materiale med angivelse
af volumenerne ΩM og Ωn. Alle mål i [mm].
For enkelte af estimaterne anvendes porøsiteten f af enhedscellen, hvorfor denne bestemmes som porevolumenet i forhold til det totale volumen, formel (2.1).
f =
Vp
V
=
π ⋅ a2 ⋅ t
(2.1)
b2 ⋅ t
hvor
Vp er porevolumenet [mm3]
V er det totale volumen [mm3]
a er radius af hullet [mm]
b er enhedscellens sidelængde [mm]
t er tykkelsen af enhedscellen [mm]
Porøsiteten bestemmes til
π ⋅ ( 4mm ) ⋅ 20mm
2
f =
14
(16mm )
2
⋅ 20mm
=
π
16
0,1964
2 Estimering af materialeparametre
2.2 Analytisk estimering af materialeparametre
I de følgende afsnit foretages en estimering af de effektive materialeparametre for det porøse materiale. Det gøres ud fra før nævnte analytiske metoder, som bygger på betragtning
af enhedscellen, afsnit 2.1. Afsnittene er alle skrevet ud fra noterne 2, 5, 6, 8 og 9 af Henrik
Myhre Jensen, Professor ved Aalborg Universitet, hvis ikke andet er angivet. [Myhre,
2005]
2.2.1 Ashby og Gibsons cellulære metode
Ved et cellulært materiale forstås et materiale med en høj porøsitet, hvilket vil sige omkring 90 – 95 % [Myhre 1, 2005]. Ashby og Gibsons metode er en simpel model til estimering af elasticitetsmodulet E for et cellulært materiale. Metoden er udledt på baggrund af
materialer, hvis mikrostruktur minder om et bjælke-søjlesystem, hvor der er en stor porøsitet, figur 2.2.
Figur 2.2. Cellulært materiale med stor porøsitet.
I det efterfølgende er det antaget, at Ashby og Gibsons metode kan anvendes til estimering
af elasticitetsmodulet for enhedscellen, som er vist på figur 2.1. Det bør nævnes, at metoden kun giver et overslag på elasticitetsmodulet.
Metoden giver et estimat, der er en funktion af hhv. det porøse og det homogene materiales
densitet. Estimatet beregnes af formel (2.2).
E
porøs
=E
homogen
⎛ ρ porøs ⎞
⎜ homogen ⎟
⎝ρ
⎠
2
(2.2)
hvor
Eporøs er elasticitetsmodulet for det porøse materiale [MPa]
Ehomogen er elasticitetsmodulet for det homogene materiale [MPa]
ρporøs er densiteten af det porøse materiale [g/mm3]
ρhomogen er densiteten af det homogene materiale [g/mm3]
15
2.2 Analytisk estimering af materialeparametre
Densiteten af det homogene materiale er i laboratoriet bestemt til 2.703 kg/m3. Til bestemmelse af det porøse materiales densitet beregnes af formel (2.3) volumenet af det homogene materiale i enhedscellen.
V homogen = b 2 ⋅ t − a 2 ⋅ π ⋅ t
(2.3)
Ved at anvende målene for enhedscellen angivet i afsnit 2.1 bestemmes Vhomogen til
V homogen = (16 ⋅16 ⋅ 20 ) mm3 − ( 4mm ) ⋅ π ⋅ 20mm = 4.115mm3
2
Derefter bestemmes massen af det homogene materiale mhomogen af formel (2.4).
m homogen = ρ homogen ⋅V homogen
m homogen = 2, 7 ⋅10−3
g
mm3
(2.4)
⋅ 4.115mm3 = 11,12 g
Endeligt beregnes ρporøs af formel (2.5).
ρ porøs =
ρ
porøs
m homogen
V porøs
(2.5)
11,12 g ⋅106
=
= 2.172 mkg3
16mm ⋅16mm ⋅ 20mm
Heraf estimeres elasticitetsmodulet E som et forhold mellem det homogene materiale og
det porøse materiales densitet iht. formel (2.6) til
2
E
porøs
=E
homogen
⎛ 2,17 ⎞
homogen
⋅⎜
⎟ = 0, 65 ⋅ E
2,
70
⎝
⎠
(2.6)
Elasticitetsmodulet blev i laboratoriet bestemt til 79 GPa, og derfor estimeres det porøse
elasticitetsmodul iht. Ashby og Gibsons metode til 51,4 GPa.
16
2 Estimering af materialeparametre
2.2.2 Voigt estimat
Voigt-estimatet er et estimat opstillet af Voigt i 1889 til beregning af stivhedstensoren for
et kompositmateriale. I metoden antages en homogen tøjningstilstand i alle bestanddele af
kompositmaterialet. I det efterfølgende tages udgangspunkt i enhedscellen vist på figur 2.1.
Den gennemsnitlige spænding σij i materialet kan beregnes af formel (2.7).
σ ij =
1
σ ijlokal dV
∫
VV
(2.7)
hvor
V er det samlede volumen af enhedscellen [mm3]
σijlokal er den lokale spænding i de enkelte bestanddele af materialet [N/mm2]
Erstattes σijlokal i formel (2.7) med den konstitutive relation givet ved formel (2.8)
lokal
σ ijlokal = Eijkl
⋅ ε kllokal
(2.8)
fås formel (2.9)
σ ij =
⎞
1⎛
M
M
c
c
⎜ ∫ Eijkl ⋅ ε kl dV + ∫ Eijkl ⋅ ε kl dV ⎟
V ⎝ V −Ω
Ω
⎠
(2.9)
hvor
V-Ω er volumenet af matrixmaterialet [mm3]
Ω er volumenet af de enkelte bestanddele [mm3]
EM er stivhedstensoren for matrixmaterialet
Ec er stivhedstensoren af de øvrige bestanddele
Med et matrixmateriale forstås det overordnede materiale, der binder de øvrige bestanddele
sammen. I dette tilfælde består matrixmaterialet af aluminium. Da der jf. Voigts metode
antages homogen tøjningstilstand i materialet gælder følgende relation
ε klM = ε kl1 = ε kl2 = ε kl
Derfor kan formel (2.9) omskrives til formel (2.10)
σ ij =
n
1⎛
M
c ⎞
V
E
−
Ω
+
Ωc Eijkl
(
)
∑
ijkl
⎜
⎟ ⋅ ε kl
V⎝
c =1
⎠
(2.10)
Udtrykket i parentesen er et udtryk for stivhedstensoren for kompositmaterialet som helhed
angivet af Voigt. Da de enkelte bestanddeles lokale stivhedstensorer ikke afhænger af volumenet, skrives udtrykket i parentesen i formel (2.10) om til formel (2.11).
17
2.2 Analytisk estimering af materialeparametre
Voigt
=
Eijkl
Ω n
V − Ω M Ω1 1
Eijkl +
Eijkl + ... + n Eijkl
V
V
V
(2.11)
Forsøgsemnerne af aluminium er lidt specielle, da materialet kun består af aluminium og
luft. Formel (2.11) reduceres derfor til formel (2.12)
Voigt
Eijkl
=
V − Ωcirkulært hul
V
alu
Eijkl
+
Ωcirkulært hul
V
luft
Eijkl
(2.12)
For enhedscellen med huldiameteren 8 mm estimeres stivhedstensoren til følgende, da sidste led af formel (2.12) falder bort grundet Eijklluft = 0.
Voigt
ijkl
E
16mm ⋅16mm ⋅ 20mm − 42 mm 2 ⋅ π ⋅ 20mm alu
alu
=
Eijkl = 0,80 ⋅ Eijkl
16mm ⋅16mm ⋅ 20mm
Stivhedstensoren Eijkl kan for et homogent og isotropt materiale med plan spændingstilstand skrives som følgende [Cook, et. al., 2002]
⎡
⎤
⎢1 ν
0 ⎥
⎥
E ⎢
ν 1
0 ⎥
[ E] =
2 ⎢
1 −ν ⎢
1 −ν ⎥
⎢0 0
⎥
2 ⎦
⎣
(2.13)
Heraf følger iht. ovenstående beregninger
E Voigt
(
1 − ν Voigt
)
2
⎡
⎢ 1
ν Voigt
⎢ Voigt
1
⎢ν
⎢
⎢ 0
0
⎣
⎤
⎥
⎥
E alu
0 ⎥ = 0,8
1 − ν alu
1 −ν Voigt ⎥⎥
2 ⎦
0
( )
2
⎡
⎢ 1 ν alu
⎢ alu
1
⎢ν
⎢
⎢ 0
0
⎣
⎤
⎥
⎥
0 ⎥
1 −ν alu ⎥⎥
2 ⎦
0
Af ovenstående kan opstilles følgende to ligninger til bestemmelse af EVoigt og νVoigt
E Voigt
(
1 − ν Voigt
)
2
2
⋅ν Voigt = 0,8 ⋅
E Voigt
(
1− ν
Voigt
)
E alu
= 0,8 ⋅
( )
1 − ν alu
(2.14)
2
E alu
( )
1− ν
alu
2
⋅ν alu
(2.15)
Indsættes de i laboratoriet fundne værdier af elasticitetsmodulet og Poisson’s forhold på
hhv. 79 GPa og 0,3 fås følgende værdier for EVoigt og νVoigt
EVoigt = 63, 2GPa
ν Voigt = ν alu = 0,3
18
2 Estimering af materialeparametre
2.2.3 Reuss estimat
Reuss-estimatet er analogt til Voigt-estimatet med den forskel, at der antages en homogen
spændingstilstand i stedet for en homogen tøjningstilstand. Fremgangsmåden er den samme som for Voigts estimat, hvorfor middeltøjningen i materialet opskrives som
ε ij =
⎞
1 lokal
1⎛
M
c
ε ij dV = ⎜ ∫ Cijkl
⋅ σ klM dV + ∫ Cijkl
⋅ σ klc dV ⎟
∫
VV
V ⎝ V −Ω
Ω
⎠
(2.16)
hvor
V-Ω er volumenet af matrixmaterialet [mm3]
Ω er volumenet af de enkelte bestanddele [mm3]
CM er fleksibilitetstensoren for matrixmaterialet
Cc er fleksibilitetstensoren af de øvrige bestanddele
Da der antages homogen spændingstilstand i alle punkter, omskives formel (2.16) til formel (2.17), da følgende gælder
σ klM = σ kl1 = σ kl2 = σ kl ⇒
n
1⎛
M
h ⎞
ε ij = ⎜ (V − Ω ) Cijkl + ∑ Ωc Cijkl
⎟ ⋅ σ kl ⇒
V⎝
c =1
⎠
(2.17)
Reuss
⋅ σ kl ⇒
ε ij = Cijkl
Reuss
=
Cijkl
Ω
V − Ω M Ω1 1
n
Cijkl + ⋅ Cijkl + ... + n ⋅ Cijkl
V
V
V
Af formel (2.17) følger det, at Reuss’ estimatet er et estimat af fleksibilitetstensoren, der er
lig den inverse af stivhedstensoren, hvilket vil sige at følgende relation gælder
(
Reuss
Reuss
Cijkl
= Eijkl
)
−1
(2.18)
Af formel (2.17) bestemmes et estimat for fleksibilitetstensoren CReuss for enhedscellen
(16 ⋅16 ⋅ 20 ) mm3 − ( 4mm ) ⋅ π ⋅ 20mm ⋅ C alu + ( 4mm ) ⋅ π ⋅ 20mm ⋅ C luft
=
ijkl
(16 ⋅16 ⋅ 20 ) mm3
(16 ⋅16 ⋅ 20 ) mm3 ijkl
2
Reuss
ijkl
C
2
(2.19)
Fleksibiliteten Cluft går imod uendelig, hvilket svarer til, at stivheden går imod nul. Derfor
giver Reuss-estimatet ikke noget brugbart resultat for porøse materialer.
19
2.2 Analytisk estimering af materialeparametre
2.2.4 Dilute estimat
Følgende udledes Dilute estimatet for et generelt tilfælde mht. RVE. Tøjningerne bestemt
for det generelle tilfælde anvendes efterfølgende til at bestemme Dilute estimatet for enhedscellen.
Dilute estimatet forudsætter, at porerne er spændingsfrie, hvilket betyder, at porerne forstyrrer spændingsbilledet i RVE, da spændingerne i porerne er lig med nul. Spændingerne
kan derfor opskrives som
σ ij =
1 S
1
ti x j dS = ∫ σ ijlokal dV
∫
V S
VM
(2.20)
hvor
S er den ydre overflade af RVE
tiS er kræfter, der virker på S
M er volumenet af matrixmaterialet
Formel (2.20) medfører, at de gennemsnitlige lokale spændinger i RVE er lig med de påførte globale spændinger.
Det forudsættes endvidere, at matrixmaterialet er lineærelastisk og isotropisk samt, at der
kun forekommer små flytninger.
S
S1
V
ni
S1
ni
M
Figur 2.3. RVE med det totale volumen V og angivelse af matrixmaterialets volumen M, ydre overfladerand
S, indre overfladerande S1 og enhedsnormalen ni.
For et homogent materiale kan tøjnings-flytningsrelationen beskrives ved formel (2.21).
ε ij = 12 ( ui , j + u j ,i )
hvor
ui,j og uj,i er flytningsgradienter
20
(2.21)
2 Estimering af materialeparametre
En flytningsgradient i et materiale med mikrostruktur regnes ækvivalent med den gennemsnitlige flytningsgradient over volumenet af RVE. Den gennemsnitlige flytningsgradient
bestemmes som
ui , j =
1 lokal
u
dV
V V∫ i , j
(2.22)
Ved at benytte Gauss divergenssætning omskrives volumenintegralet, formel (2.22), til et
fladeintegral, formel (2.23).
ui , j =
1 lokal
u n j dS
V ∫S i
(2.23)
hvor
nj er den udadrettede enhedsnormalen til den ydre overflade S
Ved indsættelse af formel (2.23) i (2.21) fås følgende udtryk for tøjningen i et homogent
materiale
ε ij =
1
2V
∫ (u
lokal
i
n j + u lokal
ni ) dS
j
(2.24)
S
I modsætning til spændingerne kan den gennemsnitlige lokale tøjning ikke siges at være lig
med den globale tøjning for et homogent RVE. Det skyldes de spændings frie porer, der
betyder, at tøjningerne i disse ikke kan bestemmes af formel (2.24), hvorfor der til tøjningerne bestemt af formel (2.24) skal adderes tillægstøjninger. Den globale tøjning εij for et
porøst RVE kan derfor udtrykkes som
ε ij =
1
ε ijlokal dV = ε ij + ε ijc
∫
VM
(2.25)
hvor
ε ij er de globale tøjningerne for et homogent RVE
εijc er tillægstøjningerne pga. porerne i et porøst RVE
De globale tøjninger for det homogene RVE udtrykkes ved formel (2.26) og tillægstøjningerne ved formel (2.27).
ε ij = Cijklσ kl
ε ijc = −
1
2V
(2.26)
∫ (u
lokal
i
n j + u lokal
ni ) dS
j
(2.27)
S1
hvor
Cijkl er den globale fleksibilitetstensor
21
2.2 Analytisk estimering af materialeparametre
Dilute-estimatet for enhedscellen
I det følgende udledes Dilute-estimaterne vha. spændings-tøjningsrelationen for et plant
spændingstilfælde samt tillægstøjningerne.
For et homogent materiale udtrykkes spændings-tøjningsrelationen i det plane spændingstilfælde ved Hooke’s udvidede lov, formel (2.28).
ε αβ =
1
EM
( (1 +ν )σ
M
αβ
−ν M δαβ σ γγ )
(2.28)
hvor
EM er matrixmaterialets elasticitetsmodul [MPa]
νM er Poisson’s forhold af matrixmaterialet [-]
εαβ er tøjningen i retningen xx, yy eller xy
σαβ er spændingen i retningen xx, yy eller xy
σγγ er summen af spændingerne både i retningen xx og yy
δαβ er Kroneckers delta
Kroneckers delta og σγγ bestemmes af (2.29).
σ γγ = σ xx + σ yy
⎧1 hvis α = β
⎩0 hvis α ≠ β
δαβ = ⎨
(2.29)
Normal- og forskydningstøjningerne for det homogene materiale bestemmes følgende til
ε xx =
σ xx −νσ yy
E
∧ ε yy =
σ yy −νσ xx
E
∧ ε xy =
1 +ν
σ xy
E
Spændings-tøjningsrelationen for det porøse materiale fås ved superposition af formel
(2.28) og tillægstøjningerne.
I enhedscellen bestemmes tillægstøjningerne ved at antage énakset spændingstilstand i enhedscellen, figur 2.4. For énakset spændingstilstand gælder der, at enhedscellen er belastet
på to modstående sider langt fra det cirkulære hul, hvilket ikke helt er tilfældet for den anvendte enhedscelle. Tillægstøjningerne i enhedscellen med det cirkulære hul udledes ud fra
formel (2.27), men for nemheds skyld omskrives flytningerne u først til polære koordinater, da det er flytningerne på randen af hullet der bestemmes. Formel (2.27) omskrives til
(2.30).
ε ijc = −
hvor
1
2V
∫ (u
lokal
i
n j + u lokal
ni ) dS = −
j
S1
1
2b 2
2π
∫ (u
lokal
i
n j + u lokal
ni ) a dθ
j
0
n1 = − cos (θ ) ∧ n2 = − sin (θ )
u1lokal = ur cos (θ ) − uθ sin (θ ) ∧ u2lokal = ur sin (θ ) − uθ cos (θ )
22
(2.30)
2 Estimering af materialeparametre
u1 og u2 er flytningerne i et kartesisk koordinatsystem
ur og uθ er flytningerne i et polært koordinatsystem
y
σ xx
σ xx
x
a=4
Figur 2.4. Enhedscellen udsat for énakset spænding.
Flytningerne i polære koordinater er givet ved formel (2.31), som er udledt ud fra de totale
spændinger for et cirkulært hul i en skive udsat for énakset spænding.
σ xx ⎛
⎞
⎛
a4 ⎞
4a 2
a2
ur =
cos ( 2θ ) + (1 + ν ) + (1 −ν ) r ⎟
⎜ (1 + ν ) ⎜ r − 3 ⎟ cos ( 2θ ) +
2E ⎝
r ⎠
r
r
⎝
⎠
σ xx ⎛
⎞
⎛
a4 ⎞
2a 2
uθ = −
(1 −ν ) ⎟ sin ( 2θ )
⎜ (1 + ν ) ⎜ r − 3 ⎟ r +
2E ⎝
r ⎠
r
⎝
⎠
(2.31)
Formlerne i (2.31) er generelle udtryk for spændingerne i et hvilket som helst sted i skiven,
men da det kun er spændingerne på randen af hullet, der er interessante, sættes r = a. Det
medfører, at flytningerne kan skrives som
σ xx a
(1 + 2 cos ( 2θ ) )
E
2σ a
uθ = − xx sin ( 2θ )
E
ur =
(2.32)
Tillægstøjningerne bestemmes nu ved integration af formel (2.30) over hele hullets rand
[0,2π]. Integrationen er foretaget i Matlab, matlab-fil a11, og deraf fås følgende tøjninger
for en énakset belastning i x1-retningen
ε xxc =
3σ 11 f
;
E
ε yyc = −
σ 11 f
E
;
ε xyc = 0
Tilsvarende findes tillægstøjningerne for en énakset belastning i x2-retningen til
ε xxc = −
σ 22 f
E
;
ε yyc =
3σ 22 f
;
E
ε xyc = 0
23
2.2 Analytisk estimering af materialeparametre
Slutteligt bestemmes forskydnings-tillægstøjningerne ved at rotere koordinatsystemet 45º i
forhold til det originale og påføre spændingerne -σ12 og σ12 i hhv. den nye x2-retning og x1retning. Det medfører følgende tillægstøjninger
ε xxc = 0;
ε yyc = 0;
ε xyc =
4σ 12 f
E
Heraf fås følgende relationer for spændings-tøjningsrelationen for det porøse materiale ved
superposition af tøjningerne for det homogene materiale og tillægstøjningerne.
σ xx −νσ yy
3σ xx f σ yy f
−
E
E
E
σ −νσ xx 3σ yy f σ xx f
ε yy = yy
+
−
E
E
E
4σ f
1 +ν
ε xy =
σ xy + xy
E
E
ε xx =
+
Spændings-tøjningsrelationerne for det porøse materiale kan skrives tilsvarende formel
(2.28) ved formel (2.33).
ε αβ =
1
Ed
( (1 +ν ) σ
d
αβ
−ν d δαβ σ γγ )
(2.33)
hvor
Ed er det estimerede elasticitetsmodul efter Dilute [MPa]
νd er det estimerede Poissons forhold efter Dilute [-]
Det estimerede elasticitetsmodul og Poissons forhold bestemmes efter formlerne i (2.34).
Ed =
E
1+ 3 f
νd =
ν+f
1+ 3 f
(2.34)
hvor
f er porøsiteten af enhedscellen
Endelig bestemmes de estimerede materialeparametre for det porøse materiale, hvor materialeparametrene for det homogene materiale er E = 79.000 MPa og ν = 0,3 bestemt ved
forsøg, afsnit B.1.
24
Ed =
79.000 MPa
= 49.700 MPa
1 + 3 ⋅ 0,1964
νd =
0,3 + 0,1964
= 0,31
1 + 3 ⋅ 0,1964
2 Estimering af materialeparametre
2.2.5 Den selvkonsistente metode
I modsætning til Dilute estimatet tager den selvkonsistente metode hensyn til den vekselvirkning, der er mellem hullerne. Den selvkonsistente metode er en approksimativ metode,
hvor de estimerede materialeparametre bestemmes ud fra de effektive materialeparametre
af kompositmaterialet og ikke materialeparametrene for matrix materialet, som det gælder
for Dilute estimatet.
Selve udledningen af de estimerede materialeparametre Ek og νk for den selvkonsistente
metode er tilsvarende udledningen for Dilute estimatet. Derfor kan spændingstøjningsrelationen for det porøse materiale ved den selvkonsistente metode opskrives direkte, hvor Ek og νk er de ubekendte, der skal bestemmes. Spændings-tøjningsrelationen er
givet ved
ε xx =
ε yy =
ε xy =
σ xx −νσ yy
+
E
σ yy −νσ xx
+
E
3σ xx f σ yy f
−
Ek
Ek
3σ yy f
Ek
−
σ xx f
Ek
4σ f
1 +ν
σ xy + xy
E
Ek
som igen kan udtrykkes ved Hooke’s udvidede lov
ε αβ =
1
Ek
( (1 +ν )σ
k
αβ
−ν k δ αβ σ γγ )
(2.35)
For den selvkonsistente metode bestemmes det estimerede elasticitetsmodul og Poisson’s
forhold efter formlerne i (2.36).
Ed = E (1 − 3 f )
ν d = ν (1 − 3 f ) + f
(2.36)
De estimerede materialeparametre for det porøse materiale, hvor E, ν og f er de samme
som for Dilute estimatet, er bestemt til
Ed = 79.000 MPa (1 − 3 ⋅ 0,1964 ) = 32.450 MPa
ν d = 0,3 (1 − 3 ⋅ 0,1964 ) + 0,1964 = 0,32
Det skal bemærkes, at formlerne i (2.36) kun er gyldige for porøsiteter af enhedscellen,
som er under 0,33, da større porøsiteter giver negative værdier af Ek og νk.
25
2.3 Numerisk estimering af materialeparametre
2.3 Numerisk estimering af materialeparametre
I det følgende estimeres elasticitetsmodulen og Poisson’s forhold for det porøse materiale
vha. en numerisk model, som, de analytiske estimater, bygger på enhedscellen, afsnit 2.1.
Den numeriske model for enhedscellen opbygges i programmet MatLab. Da enhedscellen
er symmetrisk omkring begge akser, er det valgt kun at modellere en fjerdedel af denne.
Det kan lade sig gøre, fordi der ikke forekommer flytninger vinkelret på symmetri akserne
[Cook, et. al., 2002]. Ved at begrænse modelleringen til en fjerdedel af enhedscellen er det
muligt at lave en finere opdeling og få et nøjagtigere resultat på samme beregningstid.
Derved udnyttes computerkraften bedre til at opnå en mere præcis estimering af E og ν. I
figur 2.5 ses den del af enhedscellen, som modelleres i MatLab.
y
8
4
x
Figur 2.5. En fjerdedel af enhedscellen, som modelleres i Matlab. Alle mål i mm.
Ved modelleringen vælges det at regne med plan spændingstilstand, hvilket vil sige, at
spændingerne ud af planen antages at være lig nul.
Ved beregning af E og ν for det porøse materiale er der i modellen valgt at give den lodrette rand til højre en tvungen flytning, i x-aksens positive retning, på 1 mm. Randbetingelserne for en fjerdedel af enhedscellen er vist på figur 2.6, hvor de simple understøtninger på
hhv. x-aksen og y-aksen skyldes symmetri. Understøtningsforholdende på den øverste
vandrette rand er indført således, at elasticitetsmodulet og Poisson’s forhold kan bestemmes ved to ligninger med to ubekendte. Det skyldes, at den øverste vandrette rand er fastholdt i lodret retning, hvorved tøjningen i y-retningen bliver lig nul. Dvs. εyy = 0.
u=1
y,ε yy
x,ε xx
Figur 2.6. Enhedscellens randbetingelser og den påtvungne flytning.
26
2 Estimering af materialeparametre
Til estimering af E og ν benyttes Hooke’s udvidede lov, der beskriver sammenhængen
mellem spændinger og tøjninger, og som gælder for lineære, homogene og isotrope materialer.
Idet spændingen i enhedscellen som før nævnt regnes plan, kan Hooke’s udvidede lov beskrives med formel (2.37).
⎡ ε xx ⎤
⎡1
⎢
⎥ 1⎢
⎢ ε yy ⎥ = E ⎢ −ν
⎢ 2ε xy ⎥
⎢⎣ 0
⎣
⎦
−ν
1
0
⎤ ⎡σ xx ⎤
⎢ ⎥
0 ⎥⎥ ⎢σ yy ⎥
2(1 +ν ) ⎥⎦ ⎢⎣σ xy ⎥⎦
0
(2.37)
hvor
εαβ er tøjningen [-]
E er elasticitetsmodulen [MPa]
ν er Poisson’s forhold [-]
σγδ er spændingen [MPa]
Da enhedscellen er fastholdt i y-retningen bliver εyy lig 0. Af formel (2.37) fås εxx til
ε xx =
1
1
⋅ σ xx − ⋅ν ⋅ σ yy
E
E
Tøjningen εyy bliver ligeledes af formel (2.37) til
ε yy =
1
1
⋅ −ν ⋅ σ xx + ⋅ σ yy = 0
E
E
Derved kan elasticitetsmodulen og Poisson’s forhold beregnes ud fra de to ligninger med
to ubekendte.
E=
σ xx 2 − σ yy 2
ε xx ⋅ σ xx
(2.38)
ν=
σ yy
σ xx
(2.39)
Tøjningen εxx er uden videre kendt grundet den tvungne flytning i x-retningen. Dermed er
det kun σxx og σyy, der skal findes til bestemmelse af E og ν.
ε xx =
u xx 1
= = 0,125
8
l
For at beregne de globale spændinger σxx og σyy i enhedscellen benyttes formel (2.40).
Denne formel er opstillet i det numeriske program, og beregningen kan ses i den vedlagte
matlab-fil n9.
27
2.3 Numerisk estimering af materialeparametre
De globale spændinger i enhedscellen beregnes som areal-vægtede spændinger i hvert element, formel (2.40)
σ αβ =
1 n i
∑ σ αβ ⋅ Ai
A i =1
(2.40)
hvor
A er enhedscellens samlede areal inklusiv hul [mm2]
σiαβ er spændingen i det i´te element [MPa]
Ai er arealet af det i´te element [mm2]
Der er valgt at inddele modellen i Q4-elementer, for hvilket der ses et eksempel på figur
2.7. Den stiplede linie viser det deformerede legeme.
Figur 2.7. Eksempel på typografien af enhedscellen opbygget i Matlab samt den deformerede enhedscelle
(stiplet linier).
Ud fra beregningerne i Matlab fås spændingerne σxx og σyy ved en diskretisering på 1.250
elementer til hhv. 6.629 MPa og 1.765 MPa forårsaget af en vandret flytning på 1 mm.
Først bestemmes Poisson’s forhold ud fra formel (2.39) til
ν=
1.765 MPa
= 0, 2663
6.629 MPa
Elasticitetsmodulen bestemmes af formel (2.38) til
( 6.629 MPa ) − (1.765 MPa )
E=
2
0,125 ⋅ 6.629 MPa
28
2
= 49.272 MPa
2 Estimering af materialeparametre
Ved den numeriske modellering af spændingerne, og derved også E og ν, er enhedscellen i
figur 2.6 delt op i 1.250 elementer, idet resultaterne, som det ses i tabel 2.1, konvergerer
ved denne inddeling. Det ses endvidere, at ν konvergerer væsentligt hurtigere end E.
Tabel 2.1. Konvergensanalyse af E og ν. Det fremgår, at både E og ν er konvergeret ved 1.250 elementer.
Elementer
E
[MPa]
ν
[-]
2
8
32
128
288
512
800
1.152
56.231
51.641
49.945
49.437
49.337
49.301
49.285
49.276
0,252
0,256
0,262
0,265
0,266
0,266
0,266
0,266
1.250
1.568
2.048
2.592
3.200
49.272
49.270
49.267
49.264
49.263
0,266
0,266
0,266
0,266
0,266
2.4 Sammenligning af materialeparametre
Sammenligningen af materialeparametrene foretages mellem de estimerede og de eksperimentelt bestemte. På baggrund af sammenligningen vurderes, hvorvidt de enkelte metoder
er anvendelige. De estimerede materialeparametre tillades en sammenligning med de eksperimentelle, da det er samme porøse prøveemne, der ligger til grund for resultaterne. De
porøse materialeparametre bestemt af forsøget fremgår af tabel 2.2, men for en nærmere
beskrivelse af forsøget henvises til appendiks B.2.
Tabel 2.2. De porøse materialeparametre bestemt eksperimentelt. Indeks f angiver forsøg.
Materialeparametre
Resultater
Ef [MPa]
ν [-]
45.000
-
Det er ikke muligt at bestemme Poisson’s forhold eksperimentelt for det porøse materiale,
da flytningerne i x-retningen ikke er målt. Derfor anvendes et af de estimerede Poisson’s
forhold. Det vælges at benytte estimatet bestemt ved den numeriske model i de videre beregninger, da det er vurderet til at give det bedste billede af Poisson’s forhold, tabel 2.3.
Det begrundes ved, at et meget porøst materiale såsom kork har et Poisson’s forhold som
er tilnærmelsesvis nul, da tøjningerne i x-retningen er forsvindende ved belastning i yretningen. I dette tilfælde for skiven af aluminium er porøsiteten på 20 %, hvorfor der bør
forekomme et Poisson’s forhold, der er mindre end for det homogene materiale, hvilket
ikke er tilfældet for hverken Voigt, Dilute eller den selv konsistente metode. For det homogene materiale er Poisson’s forhold bestemt til 0,3, appendiks B.1.
29
2.4 Sammenligning af materialeparametre
Resultaterne af de estimerede materialeparametre er gengivet i tabel 2.3 sammen med
sammenligningen af elasticitetsmodulerne udtrykt ved differensen i procent.
Tabel 2.3. Sammenligning af eksperimentelle og estimerede materialeparametre. Indeks e angiver estimeret.
Beregningsmetode
Ashby & Gibsons cellulære metode
Voigt estimat
Reuss estimat
Dilute estimat
Den selvkonsistente metode
FEM-model
Estimeret
Ee
[MPa]
51.400
63.200
49.700
32.450
49.300
[-]
Differens
Ee - Ef
[%]
0,30
0,31
0,32
0,27
14,2
40,4
10,4
27,9
9,6
νe
Ashby og Gibson’s metode bygger på cellulære materialer med høj porøsitet, hvorfor denne metode ikke vurderes særlige anvendelig for dette materiale på trods af den relative lille
differens på ca. 14 %. Det skyldes, at materialets mikrostruktur ikke minder om et cellulært
materiale, samt at porøsiteten ikke er stor nok.
Voigt og Reuss estimaterne beskriver hhv. den øvre og nedre værdi af elasticitetsmodulen,
og af tabel 2.3 fremgår det, at Voigt er en øvre værdi. Derimod kan Reuss-estimatet ikke
bestemmes, da partikelmaterialet, som i dette tilfælde er en porøsitet, består af luft, for
hvilket elasticitetsmodulet er nul. Følgende kommenteres Reuss-estimatet ikke nærmere.
Voigts estimat bygger som nævnt på en antagelse om homogen tøjningstilstand i hele cellen, hvilket må synes som en rimelig antagelse for et kompositmateriale som helhed, men
som det fremgår af figur 2.1, består enhedscellen af matrixmaterialet aluminium og et cirkulært luftfyldt hul. At den ene bestanddel af kompositmaterialet blot er luft, kan være årsagen til, at der forekommer en differens på ca. 40 % i forhold til det eksperimentelt bestemte elasticitetsmodul, tabel 2.3. Det stemmer overens med, at en eksakt værdi af Voigtestimatet kun kan bestemmes, såfremt et vilkårligt snit har samme fordeling af matrix- og
partikelmateriale, hvilket ikke er tilfældet. Endvidere fremgår det, at Poisson’s forhold ikke
ændrer sig i forhold til det homogene materiale, og som tidligere beskrevet bør Poisson’s
forhold mindskes ved porøsiteter.
Dilute-estimatet og den selvkonsistente metode bygger på forudsætningen om, at det cirkulære hul udsættes for énakset spændingstilstand i en uendelig stor skive. Ved at studere enhedscellen nærmere ses det, at den foregående forudsætning ikke er opfyldt, da afstanden
fra den indre rand til den ydre kun er af forholdet 1:1 i forhold til radien af hullet. Til trods
for dette stemmer elasticitetsmodulet for Dilute-estimatet pænt overens med det eksperimentelle. Det gælder ikke for den selvkonsistente metode, hvor der forekommer en differens på ca. 28 %. At det lige er den ene metode frem for den anden, der stemmer nogenlunde med eksperimentet kan være tilfældigt, men da differensen af det numerisk estimerede elasticitetsmodul er af samme størrelses orden som Dilute, antages disse at være de
mest korrekte værdier af de estimerede elasticitetsmoduler.
30
2 Estimering af materialeparametre
For Poisson’s forhold bliver det porøse større end det homogene både for Dilute og den
selv konsistente metode, hvilket ikke stemmer overens med, hvad der tidligere er beskrevet. Det er erfaret for Dilute-estimatet, at formlen for Poisson’s forhold, formel (2.34), kun
er gyldig, såfremt det homogene porøse Poisson’s forhold er større end én tredjedel. For
værdier under denne værdi bliver Poisson’s forhold større end det homogene Poisson’s forhold, og netop i dette tilfælde er det homogene Poisson’s forhold bestemt eksperimentelt til
0,3.
De estimerede materialeparametre bestemt ved den numeriske FEM-model vurderes at være af en størrelsesorden, der er acceptabel i forhold til de virkelige materialeparametre. Det
skyldes, at elasticitetsmodulet har en differens på mindre end 10 %, og Poisson’s forhold
antager en værdi, der er mindre end for det homogene mareriale. I modellen er der anvendt
Q4-elementer, der ved deformation bibeholder rette sidder. Det kan give anledning til en
fejlkilde, da enhedscellen indeholder et cirkulært hul, der ikke kan modelleres perfekt. Af
konvergensanalysen i afsnit 2.3 ses det dog, at fejlkilden negligeres ved den valgte diskretisering på 1250 elementer. En anden mulig fejlkilde ved FEM-modellering er locking af
elementerne, som er nærmere beskrevet i afsnit 5.4. Locking betyder, at elementerne opfører sig stivere, end de i virkeligheden er, men igen kan fejlkilden negligeres ved en fin diskretisering. Hvor stor indflydelse nævnte fejlkilder har på materialeparametrene er uvis,
men vurderes til ikke at have den store betydning.
Det eksperimentelt bestemte porøse elasticitetsmodul er behæftet med en vis usikkerhed,
som er nærmere beskrevet i appendiks B.2. Grundet måleusikkerheden på eksperimentet
kan det tænkes, at differensen mellem elasticitetsmodulet bestemt eksperimentelt og elasticitetsmodulet bestemt ved Dilute hhv. numerisk mindskes.
I de følgende beregninger af flytningerne for de porøse forsøgsemner vælges det at anvende elasticitetsmodulet bestemt eksperimentelt og Poisson’s forhold bestemt ved den numeriske estimering, tabel 2.4.
Tabel 2.4. Materialeparametrene, der anvendes ved beregning af flytningerne for de porøse forsøgsemner.
Elasticitetsmodul [MPa]
Poisson’s forhold [-]
45.000
0,27
31
32
3
Analytiske modeller
Formålet med dette afsnit er at bestemme analytiske løsninger for flytninger i cirkelskiven
og cirkelringen, der belastes centralt af to modsatrettede belastninger. Der er opstillet følgende analytiske modeller, figur 3.1.
• Homogen cirkelskive centralt belastet (a)
• Porøs cirkelskive centralt belastet (b)
• Homogen cirkelskive jævnt belastet (c)
• Porøs cirkelskive jævnt belastet (d)
• Homogen cirkelring centralt belastet (e)
• Porøs cirkelring centralt belastet (f)
(a)
(c)
(e)
(b)
(d)
(f)
Figur 3.1. Skitse af de beregnede analytiske modeller.
33
3.1 Analytisk løsning af cirkelskive
Den analytiske løsning til cirkelskiven, både centralt- og jævnt belastet, udarbejdes vha.
Airy’s spændingsfunktion og en Fourier rækkeudvikling, der tager udgangspunkt i randspændingerne. Cirkelringens analytiske løsning udarbejdes vha. et gæt på et flytningsfelt,
hvorpå princippet for stationær potentiel energi anvendes til bestemmelse af ukendte koefficienter. Ved beregningerne for de porøse emner genanvendes de homogene beregninger
med elasticitetsmodulet og Poissons forhold ændret fra de homogene materialeparametre
til de porøse Ep og νp.
3.1 Analytisk løsning af cirkelskive
I dette afsnit bestemmes spændinger, tøjninger og flytninger i cirkelskiven, der belastes
hhv. centralt af to modsatrettede enkeltkræfter og to modsatrettede jævnt fordelte belastninger på et stykke af randen, figur 3.2. Formålet med afsnittet er, vha. Airy’s spændingsfunktion og en Fourier rækkeudvikling, at opstille en analytisk løsning for den definerede
problemstilling. Af de bestemte spændinger beregnes flytningen af to punkter på skivens
lodrette diameter. Hermed kan resultaterne for flytningerne fra hhv. den analytiske og numeriske model samt de ved laboratorieforsøg bestemte flytninger sammenlignes.
Grunden til, at der opstilles beregninger for både enkelt- og fordelte belastninger er, at den
eksakte belastningsform, der finder sted ved forsøget i laboratoriet, ikke kendes. Derfor må
det forventes, at størrelsen af den eksperimentelt bestemte sammentrykning af skiven befinder sig tæt på de to beregnede sammentrykninger. På figur 3.2 er de to belastningssituationer skitseret, hvor kraften P regnes negativ ved tryk. Målepunkterne, mellem hvilke
flytningen måles er ligeledes skitseret på figur 3.2.
2ε
P
P
2ε
L
P
P
2ε
Figur 3.2. De to belastningsformer på den cirkulære skive, samt angivelse af de to målepunkter, mellem hvilke flytningen u måles.
Den fordelte belastning, der er skitseret midt på figur 3.2, er antaget som en radiær belastning, hvilket vil sige, at der ingen forskydningsspændinger påføres randen.
34
3 Analytiske modeller
Da der er tale om en cirkulær skive, beskrives spændingsbilledet nemmest ved hjælp af
polære koordinater. På figur 3.3 er de tre spændingskomponenter σrr, σθθ og σrθ derfor angivet.
σθθ
y
r=
σr θ
σrr
1
x
σrr
σrθ
σθθ
Figur 3.3. Skitse med enhedscirkel og orientering af de tre spændingskomponenter.
Funktionerne, der opstilles for spændingerne, formuleres i første omgang som funktioner
af kraften P, radius r og vinklen θ. De polære spændingskomponenter transformeres senere
om til kartesiske spændingskomponenter, idet det hermed bliver muligt at sammenligne
resultaterne med de numeriske beregninger. Der tages udgangspunkt i enhedscirklen, figur
3.3, hvorfor de udledte spændinger i første omgang kun gælder for denne. Det betyder, at r
kan varieres mellem 0 og 1, hvor r = 0 angiver spændingen i centrum af skiven, mens r =1
angiver spændinger på randen. Udtrykkene omformuleres senere, så de kan benyttes for
cirkulære skiver med vilkårlig radius R, der er belastet af en vilkårlig størrelse af kraften P.
I det følgende gennemføres følgende skridt for slutteligt at kunne bestemme den lodrette
flytning af de to punkter i skiven:
• Opstilling af spændingsudtryk hhv. i skiven og på randen af skiven vha. Airy’s spændingsfunktion
• Opstilling af Fourier rækker for spændinger på randen
• Bestemmelse af konstanter i de to funktioner til opstilling af endelige udtryk
• Konvergensanalyse
• Transformation af spændinger fra polære til kartesiske koordinater
• Bestemmelse af tøjninger i et repræsentativt antal punkter
• Bestemmelse af sammentrykningen u
35
3.1 Analytisk løsning af cirkelskive
3.1.1 Airy’s spændingsfunktion
I dette afsnit opstilles en spændingsfunktion F, der gælder generelt for de cirkulære tilfælde beskrevet i dette afsnit [Kohl, 1930]. En spændingsfunktion er defineret som ”en funktion, der baseret på koordinater er valgt således, at spændingskomponenterne bestemt ved
differentiation opfylder ligevægtsligningerne” [Myhre, 2005]. I polære koordinater er de
koordinater, der er omfattet af definitionen, radius r og vinklen θ. Ydermere gælder, at
”hvis spændingerne er afledt af en spændingsfunktion, er ligevægtsligningerne automatisk
opfyldt [Myhre, 2005].
Airy’s spændingsfunktion for de cirkulære tilfælde kan skrives på følgende form [Kohl,
1930]
F = b0 r 2 + b1r 3 cos (θ ) + d1r 3 sin (θ ) +
∞
∑(
n=2
)
∞
(
)
an r n + bn r n + 2 cos ( nθ ) + ∑ cn r n + d n r n + 2 sin ( nθ )
n=2
(3.1)
hvor
an, bo, b1, bn, cn, d1 og dn er konstanter, der ønskes bestemt
r er radius i enhedscirklen [mm]
θ er vinklen, der måles fra enkeltkraften, figur 3.4 [rad]
Foruden ligevægtsligningerne skal kompatibilitetsligningen være opfyldt, hvilket er tilfældet, hvis Airy’s spændingsfunktion opfylder den biharmoniske ligning [Myhre, 2005], beskrevet ved formel (3.2) [Fenster, et. al.,2003].
⎛ ∂2 1 ∂ 1 ∂2 ⎞ 2
∇4 F = ⎜ 2 +
+
⎟∇ F = 0
r ∂r r 2 ∂θ 2 ⎠
⎝ ∂r
(3.2)
hvor
V2F er Laplace operatoren
V2F beregnes af formel (3.3) [Fenster, et. al., 2003].
∂ 2 F 1 ∂F 1 ∂ 2 F
∇ F= 2 +
+
∂r
r ∂r r 2 ∂θ 2
2
(3.3)
I matlab-fil a1 er det eftervist, at Airy’s spændingsfunktion i form af formel (3.1) opfylder
den biharmoniske ligning, hvoraf det sluttes, at både ligevægts- og kompatibilitetsligningerne er opfyldt.
Spændingskomponenterne σrr, σθθ og σrθ beregnes af formlerne (3.4), (3.5) og (3.6)
[Fenster, et. al., 2003].
36
3 Analytiske modeller
1 ∂ 2 F 1 ∂F
σ rr = 2
+
⇒
r ∂θ 2 r ∂r
∞
σ rr = 2b0 + 2b1r cos (θ ) + ∑ (n(1 − n)an r n − 2 + (n + 2 − n 2 )bn r n ) cos ( nθ ) +
(3.4)
n=2
∞
2d1r sin (θ ) + ∑ (n(1 − n)cn r n − 2 + (n + 2 − n 2 )d n r n ) sin ( nθ )
n=2
σ θθ =
∂2F
⇒
∂r 2
∞
σ θθ = 2b0 + 6b1r cos (θ ) + ∑ (n(n − 1)an r n − 2 + (n + 2)(n + 1)bn r n ) cos ( nθ ) +
(3.5)
n=2
∞
6d1r sin (θ ) + ∑ (n(n − 1)cn r n −2 + (n + 2)(n + 1)d n r n ) sin ( nθ )
n=2
σ rθ =
1 ∂F 1 ∂ 2 F
−
⇒
r 2 ∂θ r ∂r ∂θ
∞
σ rθ = 2b1r sin (θ ) + ∑ (n(n − 1)an r n − 2 + n(n + 1)bn r n ) sin ( nθ ) −
(3.6)
n=2
∞
2d1r cos (θ ) − ∑ (n(n − 1)cn r n − 2 + n(n + 1)d n r n ) cos ( nθ )
n=2
For begge belastningstilfælde kan udtrykkene (3.4), (3.5) og (3.6) reduceres væsentligt, da
σrr og σθθ, som er symmetriske, pr. definition kan kaldes lige funktioner, mens σrθ kan betegnes som en ulige funktion [Cullen, et. al., 2001]. Hermed kan sinus-leddene i (3.4) og
(3.5), samt cosinus-leddene i (3.6) udelades. σrr, σθθ og σrθ kan derfor skrives som [Kohl,
1930]
∞
σ rr = 2b0 + 2b1r cos (θ ) + ∑ (n(1 − n)an r n − 2 + (n + 2 − n 2 )bn r n ) cos ( nθ )
(3.7)
n=2
∞
σ θθ = 2b0 + 6b1r cos (θ ) + ∑ (n(n − 1)an r n − 2 + (n + 2)(n + 1)bn r n ) cos ( nθ )
(3.8)
n=2
∞
σ rθ = 2b1r sin (θ ) + ∑ (n(n − 1)an r n − 2 + n(n + 1)bn r n ) sin ( nθ )
(3.9)
n=2
3.1.2 Centralt belastet cirkelskive
På figur 3.4 (tv.) er den centralt belastede cirkelskive og vinklen θ, der indgår i de følgende
beregninger, vist. Til højre på figur 3.4 er spændingsfordelingen σrr på randen, som funktion af θ, afbildet, som den intuitivt må se ud. Det fremgår, at spændingen på randen går
mod minus uendeligt for θ = 0, π, 2π, …, nπ, mens den for alle andre vinkler er lig nul.
Med udgangspunkt i dette opstilles i det følgende en Fourier rækkeudvikling for randspændingen, der kan beskrive denne spændingsfordeling som funktion af θ. Desuden opstilles et lignende udtryk for σrθ.
37
3.1 Analytisk løsning af cirkelskive
σrr
P
0
π
2π
θ
θ
P
Figur 3.4. Principiel spændingsfordeling af σrr på randen som funktion af θ.
En Fourier-serie for en lige funktion, også kaldet en cosinusrække, kan skrives på formen
[Cullen, et. al., 2001]
f ( x) =
⎛ nπ
A0 ∞
+ ∑ An cos ⎜
2 n =1
⎝ p
⎞
x⎟
⎠
(3.10)
hvor
p er den øvre grænse i det interval, hvor funktionen gentager sig
A0 og An er koefficienter, der kan bestemmes af formlerne (3.11) og (3.12).
A0 =
2 p
f ( x ) dx
p ∫0
An =
⎛ nπ
2 p
f ( x ) cos ⎜
∫
p 0
⎝ p
(3.11)
⎞
x ⎟ dx
⎠
(3.12)
For en ulige funktion, kaldet en sinusrække, kan Fourierrækken skrives som [Cullen, et.
al., 2001]
∞
⎛ nπ
f ( x ) = ∑ Bn sin ⎜
n =1
⎝ p
⎞
x⎟
⎠
(3.13)
hvor Bn bestemmes af formel (3.14).
Bn =
⎛ nπ
2 p
f ( x ) sin ⎜
∫
p 0
⎝ p
⎞
x ⎟ dx
⎠
(3.14)
Da funktionerne, der ønskes fundet, ikke er funktioner af x, men derimod af θ, erstattes alle
f(x)’er og x’er i udtrykkene (3.10) - (3.14) med hhv. σ(θ) og θ.
38
3 Analytiske modeller
Til bestemmelse af koefficienterne A0 og An betragtes tilfældet på figur 3.5, hvor kraften P
er delt ud på et lille stykke 2εr, der på randen svarer til 2ε.
P
P
2ε
θ
2ε
P
P
Figur 3.5. Skitse til bestemmelse af σrr.
σrr kan skrives som følgende i intervallet 0 < θ < π, hvorefter funktionen gentager sig for
hvert π.
⎧P
for 0 < θ < ε
⎪
⎪⎪ 2ε
σ rr (θ ) = ⎨ 0 for ε < θ < π − ε
⎪ P
⎪
for π − ε < θ < π
⎩⎪ 2ε
(3.15)
Det ses af (3.15), at P har dimensionen kraft pr. længdeenhed ind i planet. P angives derfor
efterfølgende med enheden [N/mm].
Koefficienten A0 bestemmes af formel (3.11) til
A0 =
2 p
σ rr (θ ) dθ ⇒
p ∫0
π −ε
π
P
2⎛ ε P
⎞
dθ + ∫ 0dθ + ∫
dθ ⎟ ⇒
⎜ ∫0
ε
π −ε 2ε
π ⎝ 2ε
⎠
2 P
A0 =
(ε + π − (π − ε ) ) ⇒
π 2ε
2P
A0 =
A0 =
(3.16)
π
39
3.1 Analytisk løsning af cirkelskive
Herefter bestemmes An af formel (3.12) til
An =
⎛ nπ ⎞
2 p
σ rr (θ ) cos ⎜ θ ⎟ dθ ⇒
∫
0
p
⎝ p ⎠
π −ε
π
P
2⎛ ε P
⎞
n
d
d
cos
θ
θ
0
θ
cos ( nθ ) dθ ⎟ ⇒
+
+
(
)
⎜ ∫0
∫
∫
ε
π
−
ε
π ⎝ 2ε
2ε
⎠
ε
π
2 P ⎛ ⎡ sin ( nθ ) ⎤ ⎡ sin ( nθ ) ⎤ ⎞
⎜⎢
An =
⎥ +⎢
⎥ ⎟⇒
π 2ε ⎜ ⎣ n ⎦ 0 ⎣ n ⎦π −ε ⎟
⎝
⎠
2 P
An =
sin ( nε ) − sin ( 0 ) + sin ( nπ ) − sin ( n (π − ε ) ) ⇒
π 2ε n
1 P
An =
( sin ( nε ) − sin ( nπ − nε ) )
π εn
An =
(
(3.17)
)
Da det sidste sinusled i (3.17) ønskes omskrevet, introduceres følgende to generelle trigonometriske formler [Edwards, et. al., 2002]
sin (α + β ) = sin (α ) cos ( β ) + cos (α ) sin ( β )
(3.18)
sin ( −α ) = − sin (α )
(3.19)
Ved hjælp af (3.18) og (3.19) kan sinusleddet i An omskrives til
sin ( nπ − nε ) = sin ( nπ ) cos ( − nε ) + cos ( nπ ) sin ( − nε ) ⇒
sin ( nπ − nε ) = 0 + cos ( nπ ) sin ( − nε )
sin ( −nε ) = − sin ( nε ) ⇒
cos ( nπ ) sin ( −nε ) = − cos ( nπ ) sin ( nε )
(3.20)
(3.21)
Hvis disse to omskrivninger indføres i (3.17), kan An skrives som
An =
1 P
( sin ( nε ) + cos ( nπ ) sin ( nε ) )
π εn
(3.22)
For ulige værdier af n bliver An lig 0, da cos(nπ) for n=1, 3, 5, … er lig -1. Derfor tages i
det efterfølgende kun lige værdier af n med i beregningen, og med den begrundelse erstattes n med 2v, hvor v = 1, 2, 3, 4, …,∞. Formel (3.22) omskrives derfor til (3.23).
An =
An =
1 P
( sin ( 2vε ) + cos ( 2vπ ) sin ( 2vε ) ) ⇒
π ε 2v
2 P sin ( 2vε )
π
(3.23)
ε 2v
Hvis ε går mod 0, hvilket svarer til at cirkelskiven reelt belastes af to enkeltkræfter, går sinus-leddet i formel (3.23) mod 1, og derfor bliver An tilnærmet lig med
40
3 Analytiske modeller
An =
2P
π
Cosinusrækken for σrr kan, med baggrund i ovenstående beregninger, omskivninger og iht.
formel (3.10), endeligt skrives som
σ rr =
∞
A0
+ An ∑ cos ( nθ ) ⇒
2
n=2
(3.24)
2P ⎛ 1 ∞
⎞
cos
2
v
σ rr = +
θ
=
+ ∑ cos ( 2vθ ) ⎟
(
)
∑
⎜
π π v =1
π ⎝ 2 v =1
⎠
P
2P
∞
Beregningen af sinusrækken for σrθ går væsentligt hurtigere, da σrθ på randen er lig 0 i hele
intervallet, formel (3.25). Heraf følger, iht. formel (3.14), at Bn er lig 0.
σ rθ = 0 for 0 < θ < π
(3.25)
Bestemmelse af koefficienter
I det foregående blev der opstillet to forskellige udtryk for hhv. σrr og σrθ, samt et udtryk
for σθθ på baggrund af Airy’s spændingsfunktion og Fourier rækkeudviklingen. Af de opstillede ligninger kan konstanterne 2b0, 2b1, a2v og b2v bestemmes.
På randen af cirkelskiven er radius r lig 1, da formlerne er udledt på baggrund af enhedscirklen. Derfor kan begge udtryk for σrr opskrives på følgende form, der dog kun gælder
for spændinger på randen. Bemærk at n er erstattet med 2v, afsnit 3.1.2.
⎧ Airy ' s funktion for randspændinger
⇒
⎩ Fourier funktion for randspændinger
σ rr = ⎨
∞
⎧
+
+
2
b
2
b
cos
θ
( ) ∑ (2v(1 − 2v)a2v + (2v + 2 − 4v 2 )b2v ) cos ( 2vθ )
1
⎪ 0
v =1
⎪
⇒
σ rr = ⎨
∞
⎪ 2 P ⎛ 1 + cos ( 2vθ ) ⎞
⎟
⎪⎩ π ⎜⎝ 2 ∑
v =1
⎠
∞
(
)
2b0 + 2b1 cos (θ ) + ∑ 2v(1 − 2v)a2 v + (2v + 2 − 4v 2 )b2 v cos ( 2vθ ) =
v =1
(3.26)
2P ⎛ 1 ∞
⎞
+ ∑ cos ( 2vθ ) ⎟
⎜
π ⎝ 2 v =1
⎠
Af (3.26) ses det umiddelbart, at 2b0 er lig P/π, da disse to led er de eneste konstante i udtrykkene. Samtidigt ses det, at 2b1 må være lig 0, da der på højre side af lighedstegnet ikke
er et led med en struktur, der afhænger af cos(θ).
Herefter er begge udtryk for σrθ opskrevet, idet udtrykkene igen kun er gældende for spændinger på randen.
σ rθ
∞
⎧
+
2
b
sin
θ
( ) ∑ (2v(2v − 1)a2v + 2v(2v + 1)b2v ) sin ( 2vθ )
⎪
=⎨ 1
⇒
v =1
⎪0
⎩
(3.27)
41
3.1 Analytisk løsning af cirkelskive
∞
2b1 sin (θ ) + ∑ (2v(2v − 1)a2 v + 2v(2v + 1)b2v ) sin ( 2vθ ) = 0
v =1
Med 2b0 og 2b1 indført i udtrykkene for σrr og σrθ haves følgende to ligninger til rådighed
til bestemmelse af a2v og b2v.
Fra σrr haves
P
π
∞
(
)
+ ∑ 2v(1 − 2v)a2 v + (2v + 2 − 4v 2 )b2v cos ( 2vθ ) =
v =1
2v(1 − 2v)a2v + (2v + 2 − 4v )b2v =
2
P
π
+
2P
π
∞
∑ cos ( 2vθ ) ⇒
v =1
2P
(3.28)
π
Fra σrθ haves
∞
∑ ( 2v(2v − 1)a
2v
v =1
+ 2v(2v + 1)b2 v ) sin ( 2vθ ) = 0 ⇒
(3.29)
2v(2v − 1)a2v + 2v(2v + 1)b2v = 0
Hermed haves to ligninger med to ubekendte, og derfor kan a2v og b2v bestemmes til
a2v = −
b2 v =
2P
1
π 2(2v − 1)
(3.30)
2P
1
π 2(2v + 1)
(3.31)
Endeligt kan de tre spændingskomponenter opskrives, idet 2b0, 2b1, a2v og b2v indsættes i
formlerne (3.7), (3.8) og (3.9). Først omskrives σrr
∞
σ rr = 2b0 + 2b1 cos (θ ) + ∑ ( 2v(1 − 2v)a2v r 2 v −2 + (2v + 2 − 4v 2 )b2v r 2 v ) cos ( 2vθ ) ⇒
v =1
⎛ ⎛ 2 P 2v(1 − 2v) ⎞ 2v − 2 ⎛ 2 P (2v + 2 − 4v 2 ) ⎞ 2v ⎞
σ rr = + ∑ ⎜⎜ ⎜ −
r
+⎜
⎟ r ⎟ cos ( 2vθ ) ⇒
π v =1 ⎝ ⎝ π 2(2v − 1) ⎟⎠
2(2v + 1) ⎠ ⎟⎠
⎝ π
P
∞
⎛ ⎛ 2 P (−4v 2 + 2v) ⎞ 2v − 2 ⎛ 2 P (−4v 2 + 2v + 2) ⎞ 2 v ⎞
σ rr = + ∑ ⎜⎜ ⎜ −
+⎜
⎟ r ⎟⎟ cos ( 2vθ ) ⇒
⎟r
π v =1 ⎝ ⎝ π (4v − 2) ⎠
(4v + 2)
⎝ π
⎠ ⎠
P
∞
∞
⎛ ⎛ 2P
⎞
2P
+ ∑⎜⎜−
( −v ) ⎞⎟ r 2v −2 + ⎛⎜ ( −v + 1) ⎞⎟ r 2v ⎟ cos ( 2vθ ) ⇒
π v =1 ⎝ ⎝ π
⎠
⎝ π
⎠ ⎠
∞
P 2P
σ rr = +
∑ vr 2v−2 + ( −v + 1) r 2v cos ( 2vθ ) ⇒
σ rr =
P
π
42
π
v =1
(
)
σ rr =
2P ⎛ 1 ∞
⎞
+ ∑ vr 2 v − 2 + ( −v + 1) r 2v cos ( 2vθ ) ⎟ ⇒
⎜
π ⎝ 2 v =1
⎠
σ rr =
2P ⎛ 1 ∞
⎞
+ ∑ vr 2 v − 2 − vr 2 v + r 2 v cos ( 2vθ ) ⎟
⎜
π ⎝ 2 v =1
⎠
(
(
)
)
3 Analytiske modeller
Dette forkortes til formel (3.32).
⎛
⎞
2 P ⎜ 1 ∞ 2v −2
2v
+∑r
σ rr =
v − ( v − 1) r cos ( 2vθ ) ⎟
⎜
⎟
π ⎜ 2 v =1
⎟
⎝
⎠
(
)
(3.32)
Efter samme fremgangsmåde kan σθθ og σrθ bestemmes, men her afgrænses til blot at vise
resultatet.
σ θθ =
σ rθ =
2 P ⎛ 1 ∞ 2v−2
⎞
+∑r
−v + (v + 1)r 2 cos ( 2vθ ) ⎟
⎜
π ⎝ 2 v =1
⎠
(
2P
π
∞
∑r
v =1
2v −2
)
( −v + vr ) sin ( 2vθ )
2
(3.33)
(3.34)
Hermed er de tre udtryk for spændingerne opstillet som funktion af P, r og θ. Indtil videre
gælder udtrykkene dog kun for enhedscirklen for det plane tilfælde. Ved at dele udtrykkene
for σrr, σθθ og σrθ med R, som er en vilkårlig radius, kan udtrykkene benyttes til at beskrive
spændingerne i en skive med vilkårlig radius. Desuden kan skiven gives en tykkelse ved at
dele udtrykket med t. Udtrykket for eksempelvis σrr kommer til at se således ud [Kohl,
1930]
⎛
⎞
2 ⋅ P ⎜ 1 ∞ 2v −2
2v
σ rr =
v − ( v − 1) r cos ( 2vθ ) ⎟
+∑r
⎜
⎟
π ⋅ R ⋅ t ⎜ 2 v =1
⎟
⎝
⎠
(
)
hvor det samme gør sig gældende for σθθ og σrθ.. Dette benyttes i det efterfølgende til at bestemme spændingerne i den homogene cirkulære skive.
y
r=
r=3
1/2
/4
R=7
9mm
x
Figur 3.6. Definition af radier. R er den virkelige radius af skiven, mens r er et forhold, der varierer mellem
0 og 1.
43
3.1 Analytisk løsning af cirkelskive
Spændinger i centralt belastet skive
De opstillede udtryk for de tre spændingskomponenter benyttes. Spændingsudtrykkene for
σrr, σθθ og σrθ er opstillet i et MatLab program, og spændingerne kan derfor beregnes i et
vilkårligt punkt i skiven som funktion af en vilkårlig størrelse af kraften P.
For at eftervise randbetingelserne er σrr afbildet for belastningen P = -1 N/mm, radius r = 1,
R = 79 mm og t = 20 mm, hvilket svarer til spændingerne på randen af den homogene skive, figur 3.7. Spændingerne er afbildet i intervallet 0 < θ < 2π. Størrelserne er indsat i
(3.32), hvorved følgende udtryk fremkommer
σ rr =
2 ⋅ (−1) ⎛ 1 400 2v − 2
⎞
v − ( v − 1)12 v cos ( 2vθ ) ⎟
+ ∑1
⎜
π ⋅ 79 ⋅ 20 ⎝ 2 v =1
⎠
(
)
(3.35)
Der er summeret 400 gange grundet konvergensanalysen.
Figur 3.7. σrr på randen afbildet som funktion af θ i intervallet 0 < θ < 2π (tryk regnes negativt).
Af figur 3.7 fremgår det, at der er en god overensstemmelse mellem formen af den beregnede værdi af σrr og den på figur 3.4 afbildede σrr. Desuden ses det, at σrr er symmetrisk,
samt at spændingen går mod minus uendeligt for θ = 0, π, …, nπ, sammenlignet med
spændingerne afbildet på figur 3.20 for r = 1/2.
44
3 Analytiske modeller
På figur 3.8 og figur 3.9 er σθθ og σrθ afbildet for samme belastning og i samme interval.
Figur 3.8. σθθ på randen afbildet som funktion af θ i intervallet 0 < θ < 2π (tryk regnes negativt).
Figur 3.9. σrθ på randen afbildet som funktion af θ i intervallet 0 < θ < 2π (tryk regnes negativt).
På figur 3.8 og figur 3.9 ses det, at σθθ er symmetrisk, mens σrθ er konstant 0 på randen,
hvilket stemmer overens med den afgrænsning, der tidligere er gjort, hvor σrr og σθθ er antaget som symmetriske lige funktioner, mens σrθ skal være lig 0 på hele randen, afsnit
3.1.2.
45
3.1 Analytisk løsning af cirkelskive
Konvergensanalyse
Som beskrevet i det foregående består spændingsudtrykkene for σrr, σθθ og σrθ af led, der
består af en uendelig sum af hhv. cosinus- og sinusled. Da det ikke er hensigtsmæssigt at
regne med en uendelig sum, er der vha. figur 3.10, figur 3.11 og figur 3.12 lavet en grafisk
konvergensanalyse. Dette er gjort ved at afbilde spændingerne σrr, σθθ og σrθ som funktion
af antal summationsled v. Radius er valgt til r = 78/79, da gennemregninger viser, at spændingerne konvergerer mod en fast værdi ved flere summationer, jo større radius r er. θ er
valgt til π/8, men beregninger viser, at det stort set er uden betydning for antal summationer, hvor stor θ vælges.
Figur 3.10. Grafisk konvergensanalyse for σrr (tryk regnes negativt).
Figur 3.11. Grafisk konvergensanalyse for σθθ (tryk regnes negativt).
46
3 Analytiske modeller
Figur 3.12. Grafisk konvergensanalyse for σrθ (tryk regnes negativt).
Figurerne viser, at σrr, σθθ og σrθ alle konvergerer mod en fast værdi ved ca. 400 summationer. Med den begrundelse udføres de efterfølgende beregninger af spændinger og tøjninger
med v = 400.
Flytninger i centralt belastet cirkelskive
Da spændingerne fra de numeriske beregninger for den homogene skive er givet i kartesiske koordinater, transformeres de fundne spændinger om fra polære til kartesiske koordinater. Dette gøres ved at benytte transformationsformlerne givet ved formlerne (3.36),
(3.37) og (3.38) [Fenster, et. al., 2003].
σ xx = σ rr sin 2 (θ ) + σ θθ cos 2 (θ ) + 2σ rθ sin (θ ) cos (θ )
(3.36)
σ yy = σ rr cos 2 (θ ) + σ θθ sin 2 (θ ) − 2σ rθ sin (θ ) cos (θ )
(3.37)
σ yx = σ xy = (σ rr − σ θθ ) sin (θ ) cos (θ ) + σ rθ ( cos 2 (θ ) − sin 2 (θ ) )
(3.38)
σθθ
σyy
σr θ
σrr
σrr
σθθ
y
σxy
σxx
σr θ
σyx
σxx
σxy
θ
x
σyx
σyy
Figur 3.13. Angivelse af spændingsretninger.
47
3.1 Analytisk løsning af cirkelskive
Det overordnede formål med dette afsnit er, som beskrevet i indledningen, at bestemme
flytningen af to punkter på skivens lodrette diameter. Dette gøres ved først at beregne
spændingen i 12 punkter langs diameteren, hvorefter tøjningen i punkterne bestemmes.
Flytningen, der måles eksperimentelt i laboratoriet, bestemmes over 96 mm symmetrisk
omkring den vandrette diameter i skiven, og derfor beregnes i det efterfølgende kun flytningen for den øvre halvdel af skiven. Stykket på 48 mm deles op i 12 stykker af 4 mm,
hvor tøjningen bestemmes midt i hvert interval. Flytningen bestemmes derefter ved en numerisk integration. De 12 punkter fremgår af figur 3.14.
P
48 mm
Punkt 12
y
Punkt 1
Figur 3.14. Angivelse af punkterne, hvor tøjningerne bestemmes.
Spændingerne i punkterne, som er nummereret fra 1 til 12, er angivet i tabel 3.1 for en størrelse af kraften P på -1 N/mm. Ved brug af formel (3.39) [Foley, 2004] beregnes tøjningen
i hvert enkelt punkt vha. Hookes udvidede lov, da der er antaget plan spændingstilstand i
skiven.
ε yy =
1
( −νσ xx + σ yy )
E
hvor
E er elasticitetsmodulet bestemt eksperimentelt i laboratoriet til 79 [GPa]
ν er Poissons forhold bestemt eksperimentelt i laboratoriet til 0,3 [-]
48
(3.39)
3 Analytiske modeller
Tabel 3.1. Spændinger og tøjninger til beregning af flytninger i den homogene cirkelskive beregnet for en
enkeltlast på -1 N/mm (tryk regnes negativt).
Punkt
y-koordinat
[mm]
σxx
[N/mm2]
σyy
[N/mm2]
εyy
[mm/mm]
2
6
10
14
18
22
26
30
34
38
42
46
0,00020
0,00020
0,00020
0,00020
0,00020
0,00020
0,00020
0,00020
0,00020
0,00020
0,00020
0,00020
-0,00060
-0,00061
-0,00062
-0,00063
-0,00065
-0,00067
-0,00070
-0,00074
-0,00079
-0,00085
-0,00092
-0,00102
-8,4.10-9
-8,5.10-9
-8,6.10-9
-8,8.10-9
-9,0.10-9
-9,3.10-9
-9,7.10-9
-10,1.10-9
-10,7.10-9
-11,5.10-9
-12,4.10-9
-13,6.10-9
1
2
3
4
5
6
7
8
9
10
11
12
Den numeriske integration er skitseret principielt på figur 3.15. Flytningen i intervallet
[0;48] på den lodrette akse svarer til arealet af de 12 søjler.
εyy
2
6
10
14
18
22
26
30
34
38
42
46
y
Figur 3.15. Skitse af tøjningerne beregnet i de 12 punkter.
Den samlede flytning u af de to punkter på den lodrette diameter af skiven beregnes af formel (3.40).
12
u = 2∑ ε yy ,i ⋅ yi
(3.40)
i =1
hvor
εyy,i er tøjningen beregnet i hvert af de 12 punkter [mm/mm]
yi er bredden af hvert interval på y-aksen, hvilket svarer til 4 mm
((
uhomogen cirkelskive, enkeltlast = 2 ⋅ −8, 42 ⋅10−9
mm
mm
)
(
⋅ 4mm + ... + −13, 6 ⋅10−9
mm
mm
))
⋅ 4mm P ⇒
uhomogen cirkelskive, enkeltlast = 9, 65 ⋅10−7 mm
49
3.1 Analytisk løsning af cirkelskive
Hermed er flytningen af de to punkter bestemt til ca. 9,65.10-7mm, idet flytningen er beregnet for en last på -1 N/mm. For en vilkårlig størrelse af kraften P med enheden [N/mm] kan
formel (3.41) anvendes til at beregne flytningen
uhomogen cirkelskive, enkeltlast = 9, 65 ⋅10−7
mm 2
N
⋅P
(3.41)
Flytninger i porøs centralt belastet cirkelskive
For at bestemme flytningerne i den porøse centralt belastede cirkelskive anvendes samme
fremgangsmåde som for det homogene materiale, idet det homogene elasticitetsmodul erstattes med det porøse elasticitetsmodul Ep bestemt fra laboratorieforsøg. Poissons forhold
for det porøse materiale er ikke fundet ved laboratorieforsøg, men er derimod bestemt numerisk ud fra enhedscellen, afsnit 2.3. Der anvendes derfor et Poissons forhold for det porøse materiale νporøs på 0,27.
Ved bearbejdning af resultater fra belastningsforsøg af det porøse materiale i laboratoriet er
det porøse elasticitetsmodul bestemt til 45 GPa for aluminiumsskiverne med 8 mm huller,
appendiks B.2.
Spændingerne i de føromtalte punkter bliver således de samme som i de tidligere beregninger for enkeltkraften. Det er kun tøjningerne, der antager varierende værdier. Dette giver naturligvis ikke et præcist resultat, men derimod et estimat på flytningen, idet der i
princippet regnes spændinger og tøjninger i punkter på diagonalen, hvor der er huller. De
beregnede tøjninger i de 12 punkter for den porøse cirkelskive fremgår af tabel 3.2.
Tabel 3.2. Spændinger og tøjninger til beregning af flytninger i porøs skive med enkeltlast (tryk regnes negativt).
Punkt
1
2
3
4
5
6
7
8
9
10
11
12
y-koordinat
[mm]
σxx
[N/mm2]
σyy
[N/mm2]
εyy
[mm/mm]
2
6
10
14
18
22
26
30
34
38
42
46
0,00020
0,00020
0,00020
0,00020
0,00020
0,00020
0,00020
0,00020
0,00020
0,00020
0,00020
0,00020
-0,00060
-0,00061
-0,00062
-0,00063
-0,00065
-0,00067
-0,00070
-0,00074
-0,00079
-0,00085
-0,00092
-0,00102
-14,7.10-9
-14,7.10-9
-14,9.10-9
-15,2.10-9
-15,6.10-9
-16,1.10-9
-16,8.10-9
-17,7.10-9
-18,7.10-9
-20,0.10-9
-21,7.10-9
-23,8.10-9
Flytningen er beregnet efter samme fremgangsmåde som angivet ved formel (3.40).
u porøs cirkelskive, enkeltkraft = 16,80 ⋅10−7
50
mm 2
N
⋅P
3 Analytiske modeller
3.1.3 Jævnt belastet cirkelskive
I det følgende bestemmes et analytisk udtryk for en jævnt belastet cirkelskive, figur 3.1, for
derved at beregne flytningen mellem to punkter på den lodrette diagonal.
Det antages, at den jævne belastning af cirkelskiven er radiær, hvilket medfører, at der ikke
forekommer forskydningsspændinger på randen. Denne antagelse stemmer ikke overens
med forsøget, hvor der vil opstå forskydningsspændinger på randen ved den lodrette jævne
belastning. Dette vurderes ikke at have stor betydning for en lille udbredelse af den jævne
belastning. Det medfører, at der gælder samme randbetingelser, formel (3.15) og (3.25), og
udledninger som dem beskrevet for den centralt belastede cirkelskive, afsnit 3.1.2, til og
med formel (3.23).
I stedet for at lade ε gå mod 0, så forbliver ε nu en variabel, der beskriver størrelsen af den
buelængde, hvor kraften P fordeles over, figur 3.2. Konstanten An i Fourier rækkeudviklingen er som beskrevet i formel (3.23), hvilket for rækkeudviklingen giver følgende udtryk
for σrr, formel (3.42)
σ rr =
P
π
2 P sin ( 2vε )
⋅ cos ( 2vθ ) ⇒
ε 2v
v =1 π
∞
+∑
∞
sin ( 2vε )
⎞
P⎛
⋅ cos ( 2vθ ) ⎟
σ rr = ⎜1 + ∑
π ⎝ v =1 ε v
⎠
(3.42)
Konstanten Bn forbliver 0, da det antages, at der ikke forekommer forskydningsspændinger, og dermed er σrθ lig 0, som vist i afsnit 3.1.2.
Bestemmelse af koefficienter
For at bestemme koefficienterne sættes randspændingerne bestemt ved Airy´s spændingsfunktion og ved Fourier rækkeudviklingen lig hinanden, formel (3.43).
∞
⎧
+
+
2
b
2
b
cos
θ
(2v(1 − 2v)a2 v + (2v + 2 − 4v 2 )b2 v ) cos ( 2vθ )
(
)
∑
1
⎪ 0
v =1
⎪
σ rr = ⎨
∞
⎪ P ⎛ 1 + sin ( 2vε ) ⋅ cos ( 2vθ ) ⎞
⎟
⎪ π ⎜⎝ ∑
εv
v =1
⎠
⎩
(3.43)
Det ses af formel (3.43), at 2b0 er lig P/π, da ingen af disse indeholder variable, og at 2b1
er lig 0, da der ikke er led, der kun afhænger af θ.
Forskydningsspændingerne σrθ ændrer sig ikke fra den centralt belastede cirkelskive, hvilket giver udtrykket beskrevet ved formel (3.27).
Som for den centralt belastede cirkelskive indføres konstanterne 2b0 og 2b1 i udtrykkene
for σrr og σrθ, hvilket giver to ligninger til bestemmelse af a2v og b2v, formel (3.44) og
(3.45).
51
3.1 Analytisk løsning af cirkelskive
Fra σrr haves
2v (1 − 2v ) a2 v + ( 2v + 2 − 4v 2 ) b2 v =
P sin ( 2vε )
π
εv
(3.44)
Fra σrθ haves
2v ( 2v − 1) a2 v + 2v ( 2v + 1) b2 v = 0
(3.45)
Ud fra disse to ligninger bestemmes a2v og b2v til
a2 v = −
b2 v =
P sin ( 2vε )
π 2 ( 2v − 1) vε
(3.46)
P sin ( 2vε )
π 2 ( 2v + 1) vε
(3.47)
Indsættes formel (3.46) og (3.47) samt konstanterne 2b0 og 2b1 i udtrykkene for spændingerne bestemt ved Airy´s spændingsfunktion, formel (3.7), (3.8) og (3.9), fås følgende udtryk for spændingen σrr i den jævnt belastede cirkelskive, formel (3.48).
σ rr =
∞ ⎛
⎛ P sin ( 2vε )
+ ∑ ⎜ 2v (1 − 2v ) ⋅ ⎜⎜ −
π v =1 ⎜⎝
⎝ π 2 ( 2v − 1) vε
P
⎞ 2v−2 ⎞
⎟ cos ( 2vθ )
⎟⎟ ⋅ r
⎟
⎠
⎠
∞ ⎛
⎛ P sin ( 2vε ) ⎞ 2 v ⎞
+ ∑ ⎜ ( 2v + 2 − 4v 2 ) ⋅ ⎜⎜
⎟⎟ ⋅ r ⎟⎟ cos ( 2vθ ) ⇒
⎜
2
2
1
π
v
v
ε
+
(
)
v =1 ⎝
⎝
⎠
⎠
⎞
⎛ sin ( 2vε ) ⎞ 2 v − 2 ⎞
P⎛ ∞ ⎛
⎟
cos
2
σ rr = ⎜ ∑ ⎜ ( 2v − 1) ⎜⎜
r
v
θ
⋅
⎟
(
)
⎟
⎟
⎟
⎟
2
1
π ⎝⎜ v =1 ⎝⎜
v
ε
−
(
)
⎝
⎠
⎠
⎠
⎞
⎛ sin ( 2vε ) ⎞ 2 v ⎞
P⎛ ∞ ⎛
⎟
r
cos
2
v
θ
⋅
+ ⎜ ∑ ⎜ 2 (1 − v )( 2v + 1) ⋅ ⎜⎜
⎟
(
)
⎟
⎟
⎟
⎟
v
v
ε
π ⎜⎝ v =1 ⎜⎝
2
2
1
+
(
)
⎝
⎠
⎠
⎠
Dette forkortes til følgende
σ rr =
∞ ⎛
⎞
⎛ sin ( 2vε ) ⎞ 2 v − 2
⎛ sin ( 2vε ) ⎞ 2 v ⎞
P ⎛
1
cos
2
r
v
r
v
θ
⋅ ⎜1 + ∑ ⎜ ⎜
⋅
+
−
⋅
⋅
⎟
(
)
(
)
⎟
⎟
⎜
⎟
⎟
⎟
π ⎝⎜ v =1 ⎝⎜ ⎝
ε
⎠
⎝ vε
⎠
⎠
⎠
(3.48)
Udtrykkene for σθθ og σrθ er bestemt til følgende, formel (3.49) og (3.50).
52
σ θθ =
∞ ⎛
⎞
⎛ sin ( 2vε ) ⎞ 2 v − 2
⎛ sin ( 2vε ) ⎞ 2 v ⎞
P ⎛
⋅ ⎜1 + ∑ ⎜ ⎜ −
⋅r
+ ( v + 1) ⋅ ⎜
⋅ r ⎟ cos ( 2vθ ) ⎟
⎟
⎟
⎟
⎟
π ⎝⎜ v =1 ⎝⎜ ⎝
ε
⎠
⎝ vε
⎠
⎠
⎠
(3.49)
σ rθ =
∞ ⎛
⎞
⎛ sin ( 2vε ) ⎞ 2 v − 2 ⎛ sin ( 2vε ) ⎞ 2 v ⎞
P ⎛
⋅ ⎜1 + ∑ ⎜ ⎜ −
⋅r
+⎜
⋅ r ⎟ cos ( 2vθ ) ⎟
⎟
⎟
⎟
⎟
π ⎝⎜ v =1 ⎝⎜ ⎝
ε
ε
⎠
⎝
⎠
⎠
⎠
(3.50)
3 Analytiske modeller
Disse udtryk for spændingerne er som for den centralt belastet cirkelskive udledt ud fra
enhedscirkelen og for at bestemme spændingerne i en cirkel med vilkårlig radius divideres
der med R og en tykkelse, figur 3.6.
Spændinger i jævnt belastet cirkelskive
Spændingerne i den jævnt belastede cirkelskive er bestemt vha. Matlab, og da der i forsøget er anvendt en ε på ca. 0,2 rad, appendiks B.3, bestemmes spændingerne for denne jævne belastning.
Der er fortaget en konvergensanalyse, matlab-fil a6, hvor antallet af summationsled er bestemt til 400, figur 3.16. Konvergensanalysen er lavet tæt på randen, da det er der, hvor
funktionerne kræver flest summationsled for at konvergere.
Figur 3.16. Grafisk konvergensanalyse for σrr (rød), σθθ (blå) og σrθ (sort), (tryk regnes negativt).
For at eftervise at spændingsfunktionerne er korrekte, kontrolleres det, om de overholder
randbetingelserne. Dette gøres ved at plotte spændingsfunktionerne som funktion af θ og
fastholde r lig 1. På figur 3.17 ses det, at randbetingelsen, formel (3.15), er overholdt, da
der er en jævn fordelt belastning i intervallet 0 < θ < ε og igen fra π-ε < θ < π. I intervallerne er den jævne belastning ikke helt konstant, figur 3.17, hvilket skyldes overshooting
og ε´s lille størrelse. Randbetingelsen om, at σrr er lig nul på den resterende del, er ligeledes opfyldt.
53
3.1 Analytisk løsning af cirkelskive
Figur 3.17. σrr på randen afbildet som funktion af θ i intervallet 0 < θ < 2π (tryk regnes negativt).
Det er tidligere antaget, at σθθ er en symmetrisk funktion, hvilket også ses af figur 3.18.
Figur 3.18. σθθ på randen afbildet som funktion af θ i intervallet 0 < θ < 2π (tryk regnes negativt).
I hele intervallet er det antaget, at forskydningsspændingen er nul på randen, hvilket den
fundne spændingsfunktion overholder, figur 3.19.
54
3 Analytiske modeller
Figur 3.19. σrθ på randen afbildet som funktion af θ i intervallet 0 < θ < 2π (tryk regnes negativt).
Det er valgt at vise spændingerne halvt inde i cirklen, hvilket er for r = 0,5, for at illustrere,
at σrθ er en antimetrisk funktion, og spændingernes forløb gennem cirkelskiven, figur 3.20.
Figur 3.20. Alle spændingerne, σrr (rød), σθθ (blå) og σrθ (sort), afbildet som funktion af θ i intervallet 0 < θ
< 2π (tryk regnes negativt).
55
3.1 Analytisk løsning af cirkelskive
Flytninger i jævnt belastet cirkelskive
Flytningerne i den jævnt belastede cirkelskive bestemmes ud fra E og ν bestemt ved forsøget med den homogene klods, appendiks B.1.
Spændingerne er som for enkeltlasten transformeret fra polære til kartesiske koordinater
vha. formel (3.36), (3.37) og (3.38). Da det er formålet at sammenligne flytningerne fra
forsøget med de analytiske løsninger, er tøjningerne bestemt, formel (3.39), i punkterne
beskrevet i figur 3.15.
Tabel 3.3. Spændinger og tøjninger til beregning af flytninger i homogen cirkelskive med jævn belastning
(tryk regnes negativt).
y-koordinat
[mm]
σxx
[N/mm2]
σyy
[N/mm2]
εyy
[mm/mm]
1
2
0,00019
-0,00059
-8,2.10-9
2
3
4
5
6
7
8
9
10
11
12
6
10
14
18
22
26
30
34
38
42
46
0,00019
0,00019
0,00019
0,00019
0,00018
0,00018
0,00018
0,00017
0,00016
0,00015
0,00013
-0,00060
-0,00060
-0,00062
-0,00063
-0,00065
-0,00068
-0,00071
-0,00075
-0,00080
-0,00085
-0,00092
-8,3.10-9
-8,4.10-9
-8,5.10-9
-8,7.10-9
-9,0.10-9
-9,3.10-9
-9,7.10-9
-10,1.10-9
-10,7.10-9
-11,3.10-9
-12,1.10-9
Punkt
Flytningen af de to punkter på den lodrette diameter i den homogene cirkelskive for den
jævnt fordelte belastning er ud fra tøjningerne i tabel 3.3 jf. formel (3.40) beregnet til
uhomogen skive , jævnt fordelt belastning = 9,14 ⋅10−7
mm 2
N
⋅P
Flytninger i porøs jævnt belastet cirkelskive
Ved bestemmelse af flytningerne i den porøse cirkelskive med den jævne belastning avendes der ligeledes, Ep fundet ved forsøg og νporøs bestemt ud fra enhedscellen. Tøjningerne i
de 12 punkter langs den lodrette diagonal, figur 3.15, er bestemt ved formel (3.39) og vist i
tabel 3.4.
56
3 Analytiske modeller
Tabel 3.4. Spændinger og tøjninger til beregning af flytninger i den porøse cirkelskive med jævn belastning
(tryk regnes negativt).
Punkt
1
2
3
4
5
6
7
8
9
10
11
12
y-koordinat
[mm]
σxx
[N/mm2]
σyy
[N/mm2]
εyy
[mm/mm]
2
6
10
14
18
22
26
30
34
38
42
46
0,00019
0,00019
0,00019
0,00019
0,00019
0,00018
0,00018
0,00018
0,00017
0,00016
0,00015
0,00013
-0,00059
-0,00060
-0,00060
-0,00062
-0,00063
-0,00065
-0,00068
-0,00071
-0,00075
-0,00080
-0,00085
-0,00092
-14,3.10-9
-14,4.10-9
-14,6.10-9
-14,8.10-9
-15,2.10-9
-15,6.10-9
-16,2.10-9
-16,8.10-9
-17,7.10-9
-18,6.10-9
-19,8.10-9
-21,2.10-9
Flytningerne er beregnet ved formel (3.40) til
u porøs skive , jævnt fordelt belastning = 15,94 ⋅10−7
mm 2
N
⋅P
3.1.4 Sammenligning af enkelt last og jævnt fordelt last
I det følgende sammenlignes spændinger, tøjninger og flytningerne hhv. fra enkellasten og
fra den jævnt fordelte belastning.
Af normalspændingerne i x-retningen på den lodrette diagonal, som er vist på figur 3.21, er
det tydeligt, at differensen mellem spændingerne fra enkellasten og fra den jævne belastning øges betydeligt, når radius nærmer sig randen, r = 1. Bestemmelsen af flytningerne
strækker sig dog kun fra centrum og 48 mm ud, figur 3.15, hvilket medfører, at differensen
mindskes.
57
3.1 Analytisk løsning af cirkelskive
Figur 3.21. Spændingerne σxx for enkelt last (rød) og jævn belastning (blå) på den lodrette diagonal (tryk
regnes negativt).
På figur 3.22 ses normalspændingerne σyy for en enkeltlast og ved jævn belastning på den
lodrette diagonal. Figur 3.22 viser tydeligt at differensen mellem normalspændingerne fra
enkeltlasten og fra den jævne belastning er lille, når der ses på intervallet fra centrum og til
en radius på 48 mm. Ud fra formel (3.39) ses det, at det er normalspændingerne i yretningen, der er de dominerende mht. tøjningen, da E er konstant og ν er 0,3, hvilket medfører, at normalspændingerne i x-retningen kun regnes med ca. en tredjedel af værdien. Det
vurderes derfor, at forskellen på flytningerne mellem enkeltlasten og den jævnt fordelte
last i dette tilfælde ikke afviger meget.
Figur 3.22. Spændingerne σyy for enkelt last (rød) og jævn belastning (blå) på den lodrette diagonal (tryk
regnes negativt).
58
3 Analytiske modeller
Sammenlignes flytningerne i den homogene skive beregnet af modellen med enkeltkraften
med resultatet fra beregningen for den jævnt fordelte belastning ses det af følgende beregning, at resultatet afviger med
uenkeltlast − u jævnt fordelt belastning
u jævnt fordelt belastning
9, 65 ⋅10−7 − 9,14 ⋅10−7
⋅100 =
⋅100 ≈ 5, 6%
9,14 ⋅10−7
3.1.5 Opsamling
Der er gennem analytiske løsninger bestemt spændinger, tøjninger og flytninger i en homogen cirkelskive udsat for en central enkeltlast og en jævnt udbredt belastning. Ligeledes
er der bestemt spændinger i en porøs cirkelskive udsat samme to belastninger. Af de bestemte spændinger er flytningen af to punkter på cirkelskivernes lodrette diameter beregnet, tabel 3.5.
Tabel 3.5. Flytninger mellem to punkter på den lodrette diagonal.
Cirkelskive
Homogen
Porøs
|uenkeltlast|
10-7[mm]
|ujævnt fordelt belastning|
10-7[mm]
9,65P
9,14P
16,80P
15,94P
Det kan konkluderes, at der i dette tilfælde ikke er betydelig forskel på om flytningerne
beregnes ud fra en enkeltlast eller den jævnt fordelte belastning. Det er derfor vurderet, at
det kun er nødvendigt at bestemme en analytisk løsning for cirkelringen belastet med en
enkeltlast.
3.2 Analytisk løsning af cirkelring
I dette afsnit beregnes en analytisk løsning til cirkelringen til bestemmelse af flytningerne.
Flytningerne bestemmes ved at gætte på et flytningsfelt for cirkelringen. Flytningsfeltets
konstanter bestemmes derefter vha. princippet om stationær potentiel energi, hvorved denne minimeres mht. hver af konstanterne.
3.2.1 Flytningsfeltet
Som flytningsfelt vælges det at tage udgangspunkt i parameterfremstillingen for en ellipse,
formel (3.51)
⎧ x ⎫ ⎡ a ⋅ cos θ ⎤ ⎡ r ⋅ cos θ ⎤ ⎡δ a ⋅ cos θ ⎤
⎨ ⎬=⎢
⎥=⎢
⎥+⎢
⎥
⎩ y ⎭ ⎣b ⋅ cos θ ⎦ ⎣ r ⋅ sin θ ⎦ ⎣ −δ b ⋅ sin θ ⎦
ellipse
cirkel
(3.51)
δu
59
3.2 Analytisk løsning af cirkelring
y
y
δb
r
b
θ
x
a
r
δa
x
Figur 3.23. Viser ellipse med de forskellige symboler. Den stiplede linie er en cirkel med radius r.
En illustration af formel (3.51) med de forskellige symboler ses på figur 3.23. Det sidste
led, δu, er det led, der omformer den udeformerede cirkel til en deformerede cirkel, som
antages at have form som en ellipse.
En skitse af cirkelringen er vist på figur 3.24, og med stiplede linier er vist, hvordan cirkelringen er deformeret fra en cirkelring til en ellipseformet ring.
y
y
P
δu P
x
x
δu P
P
Figur 3.24. Den udeformerede cirkelring med enkeltkraften P (tv.) og den deformerede cirkelring markeret
med stiplet linie (th.).
Det antages, at flytningerne både er symmetrisk omkring den lodrette og vandrette akse
gennem origo, hvorved origo bevarer den samme placering efter en deformation. Grundet
symmetri vælges der kun at se på den øvre del af cirkelringen, og der fås følgende model
for flytningsfeltet, figur 3.25.
60
3 Analytiske modeller
y
δb2
P
δb1
r2
r
r1
θ
x
δa1
δa2
Figur 3.25. Model for den analytiske beregning af cirkelskiven. Bemærk understøtningsforholdene.
Modellen overholder kompatibilitetsbetingelsen, da deformationen ikke river cirkelringen
fra hinanden, og de statiske betingelser er overholdt, da der er ligevægt mellem den påførte
kraft og reaktionerne i understøtningerne. De kinematiske randbetingelser, er overholdt ved
de givne understøtningsforhold, hvilket beskriver, at flytningen i den lodrette symmetriakse kun kan foregå i y-retningen, mens flytningen for den vandrette kun kan foregå i xretningen.
For at kunne beskrive flytningen af alle punkterne mellem de to rande gættes der på to
formfunktioner, der interpolerer flytningerne mellem randene. Der gættes på en lineær
sammenhæng mellem flytningerne, figur 3.26, som er udtrykt ved formel (3.52).
W1
W2
1
1
r1
r2
r
r1
r2
r
Figur 3.26. Viser formfunktionerne w1 og w2, der vægter henholdsvist den indre og ydre rand.
w1 =
(r2 − r )
(r2 − r1 )
(r − r1 )
w2 =
(r2 − r1 )
(3.52)
Derved kan flytningsfeltet i x- og y-retningen udtrykkes ved formel (3.53).
u x (r , θ ) = δ a1 ⋅ cos θ ⋅ w1 + δ a2 ⋅ cos θ ⋅ w2
u y (r , θ ) = −δ b1 ⋅ sin θ ⋅ w1 − δ b2 ⋅ sin θ ⋅ w2
(3.53)
61
3.2 Analytisk løsning af cirkelring
hvor
ux er flytningen i x-aksens retning afhængende af r og θ [mm]
uy er flytningen i y-aksens retning afhængende af r og θ [mm]
w1 er en formfunktion, der vægter den indre rand, se figur 3.26 [-]
w2 er en formfunktion, der vægter den ydre rand, se figur 3.26 [-]
δa1 og δb1 er flytningen af den indre rand i x- og y- retningen [mm]
δa2 og δb2 er flytningen af den ydre rand i x- og y- retningen [mm]
Flytningsfeltet transformeres over til polære flytninger, bestående af en radiær og tangentiel flytning. Dette gøres ved transformationsformlen (3.54).
⎡ur ⎤ ⎡ cos(θ ) sin(θ ) ⎤ ⎡u x ⎤
⎢u ⎥ = ⎢ − sin(θ ) cos(θ ) ⎥ ⎢u ⎥
⎦⎣ y⎦
⎣ θ⎦ ⎣
(3.54)
Det endelige gæt på flytningsfeltet, omskrevet til polære koordinater, ses af formel (3.55).
ur = δ a1 ⋅ cos 2 θ ⋅ w1 + δ a2 ⋅ cos 2 θ ⋅ w2 + δ b1 ⋅ sin 2 θ ⋅ w1 − δ b2 ⋅ sin 2 θ ⋅ w2
uθ = (δ a1 ⋅ cos θ ⋅ w1 + δ a2 ⋅ cos θ ⋅ w2 ) ⋅ − sin θ + (δ b1 ⋅ sin θ ⋅ w1 − δ b2 ⋅ sin θ ⋅ w2 ) ⋅ cos θ
(3.55)
3.2.2 Beregning af konstanter
De fire ubekendte, konstanterne δa1, δa2, δb1 og δb2, bestemmes vha. den potentielle energi
for elementet.
Den potentielle energi kan beskrives ved formel (3.56).
Πp =U −Ω
(3.56)
hvor
U er tøjningsenergien
Ω er energien forårsaget af den påførte last
Når der ikke er en initial spænding eller tøjning i legemet kan tøjningsenergien U beskrives
ud fra formel (3.57).
U =∫
V
1
T
{ε } [ E ]{ε } dV
2
(3.57)
hvor
{ε} er tøjningstensoren
[E] er den konstitutive matrice
For plan spændingstilstand er den konstitutive matrice udtrykt ved formel (3.58), [Foley,
2004].
62
3 Analytiske modeller
⎡1 ν
E ⎢
ν 1
[E] =
1 −ν 2 ⎢
⎢⎣ 0 0
0⎤
0 ⎥⎥
1−ν ⎥
2 ⎦
(3.58)
Den potentielle energi, der stammer fra lasten kan beskrives ved formel (3.59).
Ω = {U }
T
{P}
(3.59)
hvor
{U} er flytningstensoren
{P} er krafttensoren
Derved kan formel (3.56) omskrives til formel (3.60).
Πp = ∫
V
1
T
T
{ε } [ E ]{ε } dV − {U } {P}
2
(3.60)
Det vælges at regne i polære koordinater, idet det er nemmere at beskrive en cirkel i polære
koordinater frem for i kartesiske. Idet at cirkelringen har en konstant tykkelse på 20 mm og
antagelsen om, at P kun forårsager en radiær udbøjning, kan formel (3.60) i polære koordinater omskrives til formel (3.61).
θ 2 r2
1
T
{ε } [ E ]{ε } ⋅ r dr dθ − {U }T ⋅{P}
θ1 r1 2
Π p = 20 ⋅ ∫ ∫
(3.61)
De polære tøjningerne bestemmes ud fra formel (3.62), [Myhre, 2005].
⎧ ∂u
⎫
⎪ r
⎪
r
∂
⎪
⎪
⎧ε rr ⎫
⎪ ⎪ ⎪⎪ ur 1 ∂uθ
⎪⎪
{ε } = ⎨εθθ ⎬ = ⎨ + ⋅
⎬
⎪ε ⎪ ⎪ r r ∂θ
⎪
⎩ rθ ⎭ ⎪ 1 ⎛ 1 ∂u u ∂u ⎞ ⎪
r
− θ + θ ⎟⎪
⎪ ⎜ ⋅
∂
∂ r ⎠ ⎭⎪
2
r
r
θ
⎩⎪ ⎝
(3.62)
hvor
ur er den radiære flytning [mm]
uθ er den tangentielle flytning [mm]
r er radius [mm]
θ er vinklen [rad]
Udtrykket for den potentielle energi er nu opstillet, og det ses, at tøjningstensoren og dermed også den potentielle energi afhænger af de fire ubekendte konstanter. Ved at anvende
princippet om stationær potentiel energi kan udtrykkende vist i formel (3.63) opskrives.
63
3.2 Analytisk løsning af cirkelring
∂Π P
∂Π P
∂Π P
∂Π P
=0;
=0;
=0;
=0
∂δ a1
∂δ a2
∂δ b1
∂δ b2
(3.63)
Derved fås fire ligninger med fire ubekendte, hvori δa1, δa2, δb1 og δb2 kan isoleres. Disse
værdier kan så sættes ind i ligningen for flytningerne i polære koordinater, formel (3.55).
Alle disse udregninger er foretaget i Matlab, og udregningerne kan ses i matlab-fil a10.
3.2.3 Resultater
Flytningen, der måles eksperimentelt i laboratoriet, bestemmes over 96 mm, symmetrisk
omkring den vandrette diameter i skiven. For at have et sammenligningsgrundlag skal de
resultater, der bestemmes i den analytiske model derfor multipliceres med to, idet der kun
er taget udgangspunkt i den øverste halvdel af cirkelskiven. Der vælges at belaste cirkelringen med 25 kN.
Resultaterne fra de analytiske beregninger er skematiseret i tabel 3.6.
Tabel 3.6. Flytningsresultaterne for den homogene og porøse cirkelring.
Cirkelring
Homogen
Porøs
P
[kN]
r1
[mm]
r2
[mm]
r
[mm]
θ
[rad]
E
[Mpa]
ν
[-]
|2u|
[mm]
25
25
24
24
79
79
48
48
π/2
π/2
79.000
45.000
0,30
0,27
0,0238
0,0416
3.2.4 Opsamling
Det ses ud fra beregningerne, at nedbøjningen bliver næsten dobbelt så stor i den porøse
cirkelring i forhold til den homogene. Dette er også forventet, da det er stivheden, som er
halveret gennem en ændring af E og ν til Ep og νp.
Skaleres flytningen i forhold til kraften P fås følgende udtryk for flytningerne, tabel 3.7.
Tabel 3.7. Flytningerne som funktion af kraften for den homogene og porøse cirkelring.
Cirkelring
Homogen
Porøs
64
|uenkeltlast|
10-7[mm]
9,52P
16,64P
3 Analytiske modeller
3.3 Delkonklusion
Det kan konkluderes ud fra de analytiske beregninger af hhv. den centralt- og jævnt belastet cirkelskive, tabel 3.8, at der ikke er ret stor forskel i flytningerne ved de to beregninger. Det vælges derfor ikke at beregne på en jævn belastning for andre emner.
Det ses ydermere af tabel 3.8, at flytningen for cirkelringen ved en enkellast er mindre end
for cirkelskiven, hvilket ikke er forventet, da stivheden for cirkelskiven burde være større
end for cirkelringen. Det tyder på, at det gættede ellipseformede flytningsfelt ikke passer
ret godt med det virkelige flytningsfelt, hvilket i dette tilfælde giver en for stiv løsning af
cirkelringen.
Tabel 3.8. Resultater for de analytiske beregninger.
Emne
Homogen cirkelskive
Porøs cirkelskive
Homogen cirkelring
Porøs cirkelring
|uenkeltlast|
10-7[mm]
|ujævnt fordelt belastning|
10-7[mm]
9,65P
16,80P
9,52P
16,64P
9,14P
15,94P
-
65
66
4
Teori for elementmetoden
I dette afsnit udledes den generelle stivhedsmatrice anvendt i de numeriske beregninger af
cirkelskiven og cirkelringen. Inden udledningen foretages, opstilles de grundlæggende betingelser for teorien bag finite element method (FEM). Det gøres for et 2D kontinuum, idet
skiven og ringen undersøges i det plane tilfælde. Efterfølgende udledes flytningsinterpolationsmatricen, og stivhedsmatricen kan derefter opstilles ud fra betragtning af potentiel energi.
Efter udledningen af den generelle stivhedsmatrice udledes teorien for eliminering af indre
frihedsgrader i et legeme. I den forbindelse opstilles et udtryk for den reducerede stivhedsmatrice, der anvendes under opbygningen af et superelement, afsnit 6.3.1. Superelementet indgår i opbygningen af den porøse cirkelskive, for hvilken der foretages en numerisk beregning af flytningerne.
Følgende er skrevet på baggrund af Concepts and applications of finite element analysis
[Cook, et. al., 2002], hvor andet ikke er nævnt.
4.1 Indledende betingelser for FEM
De indledende betingelser, der skal være på plads, før stivhedsmatricen kan udledes i et
kartesisk koordinatsystem, er følgende
• Konstitutive betingelser
Krav til relation mellem spændinger og tøjninger
• Kinematiske betingelser
Krav til elementets flytningsfelt
• Statiske betingelser
Krav til elementets ligevægt
• Kompatibilitetsbetingelser
Krav til elementets sammenhæng
For den eksakte løsning vil ethvert lille udsnit af materialet være i ligevægt, kompabilitet
forefindes overalt, og alle randbetingelser for spændinger og flytninger vil være opfyldt.
For løsningerne bestemt efter FEM er det kun kompatibilitetsbetingelsen, der nødvendigvis
er overholdt, idet FEM er en tilnærmet metode.
4.1.1 Konstitutiv betingelse
Spændings-tøjningsrelationen kan opskrives ved formel (4.1) ud fra antagelsen om, at materialet er lineært elastisk. I dette tilfælde tages der ikke højde for initiale spændinger. Relationen gælder både for en, to og tre dimensioner.
{σ} = [E]{ε}
(4.1)
67
4.1 Indledende betingelser for FEM
hvor
{σ} er spændingstensoren udtrykt som en søjlevektor
[E] er den konstitutive matrice
{ε} er tøjningstensoren udtrykt som en søjlevektor
Relationen udskrevet for to dimensioner fremgår af formel (4.2).
⎧σ xx ⎫ ⎡ E11
⎪ ⎪ ⎢
⎨σ yy ⎬ = ⎢ E21
⎪τ ⎪ ⎢ E
⎩ xy ⎭ ⎣ 31
E12
E22
E32
E13 ⎤ ⎧ε xx ⎫
⎪ ⎪
E23 ⎥⎥ ⎨ε yy ⎬
E33 ⎥⎦ ⎪⎩γ xy ⎪⎭
(4.2)
Den konstitutive matrice, som indeholder elastiske konstanter, er symmetrisk således Eij =
Eji er gældende. For et isotropt materiale, og hvor plan spændingstilstand er gældende, dvs.
σzz = τyz = τzx = 0, kan den konstitutive matrice [E] og dens inverse [E-1] udtrykkes ved formel (4.3) og (4.4).
0
⎡1 ν
⎤
E ⎢
⎥
ν 1
0
[E] =
2 ⎢
⎥
1 −ν
⎢⎣ 0 0 (1 −ν ) / 2 ⎥⎦
[ E]
−1
(4.3)
0 ⎤
⎡ 1/ E −ν / E
⎢
= ⎢ −ν / E 1/ E
0 ⎥⎥ ;
⎢⎣ 0
0
1/ G ⎥⎦
G=
E
2(1 + ν )
(4.4)
hvor
G er forskydningsmodulet [MPa]
ν er Poisson’s forhold [-]
Ud fra omskrivning af formel (4.1) kan tøjningen bestemmes, formel (4.5).
{ε} = [ E] {σ}
−1
(4.5)
4.1.2 Kinematisk betingelse
Flytningerne i FEM beskrives ved et flytningsfelt, der beskriver deformationerne af et element såvel som flytningen af det. Ud fra flytningsfeltet kan tøjningsfeltet udledes vha. tøjnings-flytningsrelationen. For en todimensional tøjnings-flytningsrelation kan tøjningerne i
et xy-plan udtrykkes efter formel (4.6) ved at anvende partiel differentiation og lade ∆x og
∆y gå mod nul, figur 4.1.
ε xx =
(a)
68
∂u
∂x
ε yy =
(b )
∂v
∂y
γ xy = 2 ε xy =
(c)
∂u ∂v
+
∂y ∂x
(4.6)
4 Teori for elementmetoden
Normaltøjningen er defineret som længdeændringen i enten x- eller y-retningen divideret
med den oprindelige længde. Forskydningstøjningen er defineret som halvdelen af den
samlede vinkeldrejning elementet får ved samtidig længdeændring i begge retninger, figur
4.1.
∆u
∆u
∆x
∆x
∆y
∆y
∆y
∆x
∆v
∆u
∆y
∆v
∆x
(a)
(b)
∆v
(c)
Figur 4.1. Et infinitesimalt element placeret i et kartesisk koordinatsystem. (a), (b) og (c) angiver hhv. normaltøjningen i x-retningen, normaltøjningen i y-retningen og forskydningstøjningen.
Ved at anvende en matrixoperator kan tøjnings-flytningsrelationen for et 2D kontinuum
opskrives som formel (4.7)
⎡∂
⎢
⎧ε xx ⎫ ⎢ ∂x
⎪ ⎪ ⎢
⎨ε yy ⎬ = ⎢ 0
⎪γ ⎪ ⎢
⎩ xy ⎭ ⎢ ∂
⎢
⎣ ∂y
⎤
0⎥
⎥
∂ ⎥ ⎧u ⎫
⎨ ⎬
∂y ⎥⎥ ⎩ v ⎭
∂⎥
⎥
∂x ⎦
(4.7)
4.1.3 Statisk betingelse
Følgende opstilles ligevægtsligninger for et plant element, placeret i et kartesisk koordinatsystem, påvirket af kræfter, figur 4.2. Disse stammer fra spændinger på siderne samt elementkræfter. Elementkræfterne Fx og Fy er defineret som kræfter per enhedsvolumen, og
kan frembringes af f.eks. tyngdekraften og acceleration. For et element med en ensformig
tykkelse opstilles ligevægtsligningen i x-retningen, hvor kræfterne regnes positive i xretningen.
−σ xx t dx − τ xy t dx + (σ xx + σ xx , x dx) t dy + (τ xy + τ xy , y dy ) t dx + Fx t dx dy = 0
(4.8)
Tilsvarende findes en ligevægtsligning i y-retningen. Ved at reducere ligningerne fås for et
plant 2D kontinuum følgende ligevægtsligninger, formel (4.9)
x-retningen:
σ xx , x + τ xy , y + Fx = 0
y -retningen:
τ xy , x + σ yy , y + Fy = 0
(4.9)
69
4.2 Interpolation og formfunktioner
σ yy +σ yy,y dy
τ xy +τ xy,y dy
Fy
dy
τ xy +τ xy,x dx
σ xx
(y,v)
Fx
τ xy
(x,u)
τ xy
σ xx +σ xx,x dx
σ yy
dx
Figur 4.2. Spændings- og elementkræfter på et plant element med ensformig tykkelse.
FEM-elementer baseret på flytningsfelter vil ikke opfylde statiske betingelser overalt i materialet, men betingelserne vil gennemsnitligt være opfyldt således, at statiske betingelser
f.eks. er opfyldt for en model som helhed. Med andre ord vil de ydre statiske betingelser
være opfyldt, hvorimod snitkræfter ikke nødvendigvis er i ligevægt.
4.1.4 Kompatibilitetsbetingelse
Generelt udledes tøjningerne ud fra kendskab til flytningerne, men i visse tilfælde er det
muligt at bestemme flytningerne ud fra kendskab til tøjningerne. Det kræver, at tøjningskomponenterne overholder kompatibilitetsbetingelsen, som for et plant tilfælde udtrykkes ved formel (4.10). Er kompatibilitetsbetingelsen overholdt er tøjningerne integrerbare, og flytningerne kan derfor bestemmes vha. af randbetingelserne. For et lineært elastisk element kan det forklares ved, at der ikke må forekomme brud ved deformation af
elementet. [Foley, 2004]
2
2
∂ 2ε x ∂ ε y ∂ γ xy
+
=
∂y 2
∂x 2
∂y∂x
(4.10)
Kompatibilitetsbetingelsen sikrer, at deformationerne over elementgrænserne er kontinuerte. I dette projekt benyttes elementer, hvor flytningsfeltet tilnærmes vha. polynomier, hvilke altid vil opfylde kompatibilitetsbetingelserne.
4.2 Interpolation og formfunktioner
Den generelle flytningsinterpolationsmatrice [N] udledes følgende, idet de indgående
formfunktioner, partielt differentieret, beskriver tøjnings-flytningsmatricen [B], der indgår
i stivhedsmatricen [k].
Formfunktionen beskriver, hvorledes sammenhængen er mellem knuder i elementet. Formfunktionen bestemmes ved at foretage en interpolation af flytningerne mellem knuderne i
70
4 Teori for elementmetoden
elementet. Indenfor FEM anvendes der hovedsageligt polynomier som interpolationsfunktioner, hvilket skyldes, at de er kontinuerte, og at der for en enkelt værdi af funktionen kun
findes et koordinatsæt i et kartesisk koordinatsystem, der tilfredsstiller denne værdi. Et interpolerende polynomium med den afhængige variabel φ og den uafhængige variabel x kan
udtrykkes ved formel (4.11).
n
φ = ∑ ai xi
(4.11)
i =0
hvor
φ er en skalar, der beskriver flytningsfeltet
ai er frihedsgraderne i knude i
xi er en uafhængig variabel
n angiver graden af polynomiet
Formel (4.11) kan skrives på matrixform ved formel (4.12).
φ = [ X ]{a} ;
[ X] = ⎡⎣1
x
x 2 … x n ⎤⎦ ;
{a} = [ a0
a1
a2
an ]
T
(4.12)
Følgende bevises det, at flytningsfeltet φ kan beskrives ved formel (4.13).
φ = [ N ]{φe }
(4.13)
hvor
[N] er flytningsinterpolationsmatricen
{φe} er knudeflytningerne
Relationen mellem knudeflytningerne og frihedsgraderne ai kan udtrykkes ved formel
(4.14) og flytningsinterpolationsmatricen ved formel (4.15). Indsættes disse i formel (4.13)
ses det, at formel (4.12) fremtræder.
{φe } = [ A ]{a}
[ N ] = [ X][ A ]
−1
(4.14)
= [ N1
N 2 …]
(4.15)
hvor
[A] er koordinat-matricen, der indeholder [X] i hver række
Ni er formfunktionen i i´te knude
Det er vha. formfunktionerne, at flytningerne i elementet kan interpoleres. [A] matricen
opstilles for de enkelte elementtyper under udledningen af tøjnings-flytningsmatricen [B].
Til bestemmelse af [X] anvendes Pascals trekant.
71
4.3 Potentiel energi i elastisk legeme
4.2.1 Pascals trekant
Ud fra Pascals trekant kan flytningerne i elementerne beskrives ved polynomier af forskellig grad, tabel 4.1. Hvilken grad polynomiet skal antage for det enkelte element afhænger
af antallet af knuder i elementet således, at flytningerne i et CST-element med tre knuder
(Q3) beskrives ved tre led, altså et første gradspolynomium. Q4-elementet, der er opbygget
af fire knuder, beskrives ved fire led, og bliver derfor et andengradspolynomium i henhold
til Pascals trekant, osv.
For at elementerne er i stand til at håndtere stivlegemebevægelser og translation, skal polynomiet altid vælges af laveste orden, hvilket betyder, at polynomierne altid starter med 1.
Endvidere må polynomiet ikke favorisere x eller y, idet det ville betyde at flytningsfeltet
ikke længere ville være i ligevægt.
Tabel 4.1. Pascals trekant til bestemmelse af polynomiets grad i et plant tilfælde.
Pascals trekant
Grad af polynomiet
1
Konstant
x
x
x
3
y
2
2
xy
2
y
2
xy
y3
xy
1. gradspolynomium
2. gradspolynomium
3. gradspolynomium
Polynomierne kan i henhold til formel (4.12) opstilles på matrixform, og vil således for
elementer med tre, fire og otte knuder blive som følgende
Q3 :
Q4 :
Q8 :
[ X ] = [1
[ X ] = [1
[ X ] = ⎡⎣1
x
y]
x
y
xy ]
x
y
xy
x2
y2
x2 y
xy 2 ⎤⎦
4.3 Potentiel energi i elastisk legeme
For senere at kunne opstille stivhedsmatricen [k] for et legeme opstilles et udtryk for legemets potentielle energi. Elementet forudsættes lineær elastisk, og den potentielle energi Пp
opstilles mht. virtuelt arbejde, formel (4.16).
Πp =U −Ω
(4.16)
hvor
U er systemets tøjningsenergi
Ω er arbejdet udført af ydre laster
Til opstilling af stivhedsmatricen udnyttes udtrykket for sammenhængen mellem elementog knudeflytninger, formel (4.13), hvor flytningerne u interpoleres mht. knudeflytningerne
d.
{u} = [ N ]{d}
72
(4.17)
4 Teori for elementmetoden
hvor
{u} = [u
v ] er elementflytningen af et 2D element med to frihedsgrader
T
Tøjningerne bestemmes nu ved differentiation af flytninger som vist i formel (4.7). Således
kan sammenhængen mellem tøjninger og flytninger opstilles.
{ε} = [ B]{d}
(4.18)
hvor
[B] er tøjnings-flytningsmatricen defineret som [B]=[∂][N]
Ved opstilling af udtrykket for virtuelt arbejde antages små virtuelle flytninger {δu}, der
overholder kompabilitet og grænseværdier for flytningerne, formel (4.19) og figur 4.3.
∫ {δ ε} {σ} dV = ∫ {δ u} {F} dV + ∫ {δ u} {Φ} dV
T
T
T
(4.19)
hvor
{F} er kraftpåvirkning i volumenet V [N]
{Φ} er kræfter på randen [N]
Figur 4.3. Laster på element og rand anvendt ved VAP.
Formel (4.17) og (4.18) omskrives og indsættes i formel (4.19)
{δ u}
T
= {δ d} [ N ]
{δ ε}
T
T
δ Π p ({d} ) = {δ d}
T
T
= {δ d} [ B ]
T
T
(4.20)
( ∫ [B] [E][B] dV {d} −∫ [N] [F] dV − ∫ [N] [Φ] dS )
T
T
T
(4.21)
Nu kan udtrykket for stivhedsmatricen [k] opskrives og samtidig ses, at {d} er uafhængig
af koordinater, og derfor ikke optræder i integralet. Samtidig er d1, d2,…, dn uafhængige af
hinanden og løsningen findes derfor ved at kræve, at der ikke sker nogen variation i potentiel energi når {d} varieres, formel (4.22).
73
4.4 Elimination af frihedsgrader
δ Π p ({d} ) = 0 ∀ {δ d}
(4.22)
[ k ]{d} = {re }
(4.23)
hvor
T
[ k ] = ∫ [B] [E][ B] dV
(4.24)
{re } = ∫ [ N ] [F ] dV + ∫ [ N ] [Φ] dS
T
T
(4.25)
Lasttensoren {re} beskriver alle laster påført elementet udefra. Udtrykket for stivhedsmatricen (4.24) kan bruges generelt, og vil blive benyttet i afsnit 5 .
4.3.1 Spændingsberegning
Alle elementer benyttet i dette projekt defineres ud fra et gæt på flytningsfeltet. Spændinger og tøjninger bestemmes derfor ud fra de flytninger, der er bestemt ved formel (4.26).
[ k ]{d} = {r}
(4.26)
hvor
{r} er laster påført i knuderne
Derpå bestemmes tøjningerne af formel (4.18) og spændingerne af (4.5). Her skal der gøres
opmærksom på, at selv hvis flytningsfeltet bestemmes tilnærmelsesvist korrekt, vil [B] afdække og forstørre eventuelle fejl i flytningsfeltet. Det skyldes, at [B] bestemmes ved differentiation af flytningsfeltet.
4.4 Elimination af frihedsgrader
I dette afsnit beskrives teorien for eliminering af indre frihedsgrader i det generelle tilfælde, hvor legemer kan opbygges af mange forskellige elementtyper. Eliminationen af de indre frihedsgrader er mulig med forskellige former for teorier, men der anvendes her en teori, der tager udgangspunkt i potentiel energi. Afsnit 4.4 er skrevet på baggrund af [Byskov,
2002].
Den potentielle energi for et legeme opskrives efter formel (4.27).
Π p ({d} ) =
1
T
T
{d} [k ]{d} − {d} {r}
2
hvor
{d} er kundeflytningerne for et legeme
[k] er stivhedsmatricen for et legeme
{r} er lasttensoren for et legeme
74
(4.27)
4 Teori for elementmetoden
Knudeflytningerne deles op i flytninger af de indre frihedsgrader {di} og flytninger af de
ydre frihedsgrader {dy}, formel (4.28) og (4.29).
{di } = [ Ti ]{d}
(4.28)
{d } = ⎡⎣T ⎤⎦ {d}
(4.29)
y
y
hvor
[Ti] er en placeringsmatrice for de indre frihedsgrader
[Ty] er en placeringsmatrice for de ydre frihedsgrader
Placeringsmatricerne angiver, hvor hhv. de indre og ydre frihedsgrader er placeret i knudeflytningstensoren, som derpå kan udtrykkes ved addition af formel (4.28) og (4.29), formel (4.30).
{d} = [Ti ] {di } + ⎡⎣Ty ⎤⎦ {d y }
T
T
(4.30)
Den potentielle energi, formel (4.31), opskrives herefter ud fra formel (4.30).
Πp
({d } ,{d }) = 12 ⎣⎢⎡{d }
T
i
y
i
[Ti ] + {d y }
T
T
T
⎡⎣Ty ⎤⎦ ⎤ [k ] ⎡[ Ti ] {d i } + ⎡⎣Ty ⎤⎦ {d y }⎤
⎣⎢
⎦⎥
⎦⎥
T
T
− ⎡⎢{d i } [ Ti ] + {d y } ⎡⎣Ty ⎤⎦ ⎤⎥ {r}
⎣
⎦
(4.31)
Stivhedsmatricen deles op vha. placeringsmatricerne, hvilket giver fire delstivhedsmatricer, formel (4.32), (4.33), (4.34) og (4.35).
[k ii ] = [Ti ][k ][Ti ]
T
⎡⎣k iy ⎤⎦ = [ Ti ][k ] ⎡⎣ Ty ⎤⎦
(4.32)
T
(4.33)
⎡⎣k yi ⎤⎦ = ⎡⎣Ty ⎤⎦ [k ][ Ti ]
T
⎡⎣k yy ⎤⎦ = ⎡⎣ Ty ⎤⎦ [k ] ⎡⎣Ty ⎤⎦
(4.34)
T
(4.35)
hvor
[kii] er stivhedsmatricen udelukkende for de indre frihedsgrader
[kiy] er stivhedsmatricen for både indre og ydre frihedsgrader
[kyi] er stivhedsmatricen for både ydre og indre frihedsgrader
[kyy] er stivhedsmatricen udelukkende for de ydre frihedsgrader
75
4.4 Elimination af frihedsgrader
Lasttensoren deles ligeledes op vha. placeringsmatricerne, formel (4.36) og (4.37).
{ri } = [Ti ]{r}
(4.36)
{r } = ⎡⎣T ⎤⎦ {r}
(4.37)
y
y
hvor
{ri} er lasttensoren for de indre frihedsgrader
{ry} er lasttensoren for de ydre frihedsgrader
Anvendes formlerne (4.32) til (4.37) i formel (4.31) fås følgende udtryk for den potentielle
energi, formel (4.38).
Πp
⎡ [k ii ]
⎣⎡k iy ⎦⎤ ⎥⎤ ⎡ {d i } ⎤
⎢
⎥
⎥ ⎣⎢{d y }⎦⎥
⎤
⎡
⎤
k
yi
yy
⎣⎣ ⎦ ⎣ ⎦⎦
({d } ,{d }) = 12 ⎡⎢⎣{d } {d } ⎤⎥⎦ ⎢⎢ ⎡k
T
T
i
y
i
y
{d }
− ⎡⎢{d i }
⎣
T
T
y
⎡ {r } ⎤
⎤⎢ i ⎥
⎦⎥ ⎢{ry }⎥
⎣
⎦
(4.38)
Udskrives formel (4.38) fås formel (4.39).
Πp
({d } ,{d }) = + 12 {d }
T
i
y
i
+
1
2
T
[k ii ]{di } + {di } ⎡⎣k iy ⎤⎦ {d y }
T
T
1
1
d y } ⎣⎡k yi ⎦⎤ {d i } + {d y } ⎣⎡k yy ⎦⎤ {d y }
{
2
2
(4.39)
− {d i } {ri } − {d y } {ry }
T
T
Det er muligt at variere den potentielle energi med hensyn til de indre frihedsgrader, da
disse ikke er sammensat med andre elementer. Det er ikke tilladt at variere den potentielle
energi med hensyn til de ydre frihedsgrader uden at tage højde for sammenhængen med de
omkringliggende legemer, og der varieres derfor kun med hensyn til de indre frihedsgrader, formel (4.40).
δ Πp
{di } = 0
∀ {δd i }
(4.40)
Formel (4.40) udskrives til formel (4.41).
0 = {δd i } [k ii ]{d i } + {δd i } ⎡⎣k iy ⎤⎦ {d y } − {δd i } {ri } ∀ {δd i }
⇓
T
T
T
(4.41)
{0} = [k ii ]{d i } + ⎡⎣k iy ⎤⎦ {d y } − {ri }
Den potentielle energi er altid positiv, hvilket medfører, at det[kii] er forskellig fra nul, og
den kan derfor inverteres, hvilket giver følgende udtryk for flytningen af de indre frihedsgrader, formel (4.42).
76
4 Teori for elementmetoden
{di } = [k ii ] ({ri } − ⎡⎣k iy ⎤⎦ {d y })
−1
(4.42)
Denne indsættes i udtrykket for den potentielle energi, formel (4.39).
Πp
({d }) = + 12 ([k ] ({r } − ⎡⎣k ⎤⎦ {d }))
T
−1
y
ii
i
(
iy
y
[k ii ] ([k ii ] ({ri } − ⎡⎣k iy ⎤⎦ {d y }) )
−1
))
(
T
1
−1
[k ii ] {ri } − ⎡⎣k iy ⎤⎦ {d y } ⎡⎣k iy ⎤⎦ {d y }
2
T
1
−1
+ {d y } ⎡⎣k yi ⎤⎦ [k ii ] {ri } − ⎡⎣k iy ⎤⎦ {d y }
2
T
1
+ {d y } ⎡⎣k yy ⎤⎦ {d y }
2
+
(
(
− [k ii ]
−1
))
(
({r } − ⎡⎣k ⎤⎦ {d }))
T
i
iy
y
(4.43)
{ri } − {d y } {ry }
T
Stivhedsmatricen er symmetrisk, og det fremgår deraf, at formel (4.44) er gældende.
⎡⎣k iy ⎤⎦ = ⎡⎣k yi ⎤⎦
T
(4.44)
Anvendes formel (4.44) og matrixalgebra kan formel (4.43) omskrives til følgende, formel
(4.45).
Πp
({d }) = − 12 {r }
[k ii ] {ri } + {d y }
i
T
1
dy }
{
2
T
1
+ {d y }
2
−
⎡⎣k iy ⎤⎦ [k ii ]
T
−1
T
y
T
−1
{ri }
⎣⎡k iy ⎦⎤ [k ii ] ⎡⎣k iy ⎤⎦ {d y }
T
−1
(4.45)
⎡⎣k yy ⎤⎦ {d y } − {d y } {ry }
T
Hvis der varieres med hensyn til flytningen af de ydre frihedsgrader, så forbliver alle de led
der ikke afhænger af flytningen af de ydre frihedsgrader konstant, og kan derfor fjernes i
udtrykket. Det medfører, at formel (4.45) kan omskrives til formel (4.46).
Πp
({d }) = + 12 {d }
T
y
y
− {d y }
T
( ⎡⎣k
({r } − ⎡⎣k
y
)
⎤⎦ − ⎡⎣k iy ⎤⎦ [k ii ] ⎡⎣k iy ⎤⎦ {d y }
T
yy
⎤⎦ [k ii ]
T
iy
−1
−1
{ri }
)
(4.46)
77
4.4 Elimination af frihedsgrader
Formel (4.47) og (4.48) indføres.
⎡⎣k *yy ⎤⎦ = ⎡⎣k yy ⎤⎦ − ⎡⎣k iy ⎤⎦ [k ii ] ⎡⎣k iy ⎤⎦
(4.47)
{r } = {r } − ⎡⎣k
(4.48)
T
y
y
⎤⎦ [k ii ]
T
*
iy
−1
{ri }
hvor
[k*yy] er den reducerede stivhedsmatrice
{r*y} er den reducerede lasttensor
Den potentielle energi, formel (4.46), kan derpå opskrives som formel (4.49).
Πp
({d }) = + 12 {d }
T
y
y
⎡⎣k *yy ⎤⎦ {d y } − {d y } {r*y }
T
(4.49)
Sammenlignes formel (4.49) med formel (4.27), ses det, at den reducerede stivhedsmatrice,
[k*yy], svarer til stivhedsmatricen for legemet og den reducerede lasttensor, {r*y}, svarer til
lasttensoren for legemet. [k*yy] kan nu benyttes som stivhedsmatrice for hele legemet.
78
5
Elementtyper
I det følgende udledes tøjnings-flytningsmatricen [B] for tre forskellige elementer, som er
CST-elementer (constant-strain triangle), 4-knudet elementer og 8-knudet isoparametriske
elementer med fire Gausspunkter. I det følgende benævnes de tre elementer hhv. Q3, Q4 og
Q8, hvor indekset referer til det enkelte elements antal knuder. De benyttede elementer er
defineret isoparametrisk, således vil både flytninger og knudekoordinater blive bestemt af
de samme flytningsfunktioner [N]. Således kan elementerne placeres vilkårligt i et kartesisk koordinatsystem, og formen kan varieres så længe geometrien ikke udsættes for voldsomme ændringer.
5.1 Isoparametrisk CST-element Q3
For en lineær, plan trekant varierer elementets flytningsfelt lineært i et kartesisk koordinatsystem, hvorfor der ved en spændingsanalyse opstår et konstant tøjningsfelt. Heraf navnet
CST-element (constant strain triangle).
En optegning af et vilkårligt CST-elementet i et kartesisk koordinatsystem fremgår af figur
5.1, hvor der samtidig er påført knudeflytninger. En knudeflytning beskrives ved en
vandret og lodret flytning af knuden benævnet hhv. u og v. Det resulterer i, at elementet
har seks frihedsgrader.
I den følgende udledning af tøjnings-flytningsmatricen transformeres elementet over i et
referencekoordinatsystem (ξ,η), således knuderne, uanset deres tidligere placering, får koordinaterne som på figur 5.1.
Figur 5.1. En vilkårlig plan trekant placeret i et kartesisk koordinatsystem (TV) og transformationen af elementet til referencekoordinatsystemet (TH).
79
5.1 Isoparametrisk CST-element Q3
Som skrevet er den primære opgave at udlede tøjnings-flytningsmatricen [B] for et vilkårligt placeret trekantelement. Men idet tøjnings-flytningsmatricen er relationen mellem tøjningerne og knudeflytninger, udtrykt ved formel (5.1), bestemmes endvidere et udtryk for
tøjningerne, der fremgår af formel (5.10).
{ε} = [ B][d]
(5.1)
For CST-elementet bestemmes indledningsvis flytningsinterpolationsmatricen i referencekoordinatsystemet [N(ξ,η)], der bestemmes som.
⎡N
⎡⎣ N (ξ ,η ) ⎤⎦ = ⎢ 1
⎣0
0
N2
0
N3
N1
0
N2
0
0⎤
N 3 ⎥⎦
(5.2)
hvor den første række beskriver flytningerne i ξ–retningen og anden række flytningerne i
η–retningen.
Med udgangspunkt i formel (4.15) kan formfunktionerne i formel (5.2) bestemmes efter
formel (5.3).
⎡⎣ N (ξ ,η ) ⎤⎦ = [ X ] ⎡⎣ A (ξ ,η ) ⎤⎦
−1
(5.3)
Ud fra Pascals trekant opskrives [X] og referencekoordinatmatricen [A(ξ,η)]. Referencekoordinatmatricen for et CST-element opskrives som en 3×3 matrice indeholdende Xmatricen, og formfunktionerne kan ud fra koordinaterne i referencekoordinatsystemet bestemmes som
−1
[ N1
N2
⎡1 0 0 ⎤
N 3 ] = [1 ξ η ] ⎢⎢1 1 0 ⎥⎥ = [1 − ξ − η ξ η ]
⎢⎣1 0 1 ⎥⎦
⇓
0
ξ 0 η 0⎤
⎡1 − ξ − η
⎡⎣ N (ξ ,η ) ⎤⎦ = ⎢
1 − ξ − η 0 ξ 0 η ⎥⎦
⎣ 0
Da tøjnings-flytningsmatricen først bestemmes for referencekoordinatsystemet nødvendiggør det en transformationsrelation, således tøjnings-flytningsmatricen kan bestemmes for et
kartesisk koordinatsystem. Transformationsrelationen benævnes Jacobimatricen [J]. Forbindelsen mellem flytningsfeltet i referencekoordinatsystemet og i det kartesisk koordinatsystem er følgende opstillet.
⎧φ , x ⎫
⎧φ ,ξ ⎫
⎨ ⎬ = [J] ⎨ ⎬
⎩φ ,η ⎭
⎩φ , y ⎭
80
(5.4)
5 Elementtyper
Det generelle udtryk for Jacobimatricen fremgår af formel (5.5).
⎡ x,
[ J ] = ⎢ x,ξ
⎣
η
y,ξ ⎤ ⎡ ∑ N i ,ξ xi
=
y,η ⎥⎦ ⎢⎣ ∑ N i ,η xi
∑ N i ,ξ yi ⎤
∑ N i ,η yi ⎥⎦
(5.5)
Jacobimatricen for et trekantelement kan efter formel (5.5) skrives som
⎡ x1
⎡ −1 1 0 ⎤ ⎢
[J ] = ⎢
⎥ ⎢ x2
⎣ −1 0 1 ⎦ ⎢
⎣ x3
y1 ⎤
⎡− x + x
y2 ⎥⎥ = ⎢ 1 2
−x + x
y3 ⎥⎦ ⎣ 1 3
− y1 + y2 ⎤
− y1 + y3 ⎥⎦
For at bestemme tøjningerne i det kartesiske koordinatsystem er det efter formel (5.4) nødvendigt at invertere Jacobimatricen, hvilket er gjort følgende.
[J ]
−1
=
1 ⎡ − y1 + y3
J ⎢⎣ x1 − x3
y1 − y2 ⎤
; J = (− x1 + x2 )(− y1 + y3 ) − (− x1 + x3 )(− y1 + y2 )
− x1 + x2 ⎥⎦
(5.6)
hvor
|J| er lig med det[J]
Hvis |J| udregnes med hensyn til koordinaterne i figur 5.1 ses det, at den er lig med 2A,
hvor A er arealet af trekanten i et kartesisk koordinatsystem.
Følgende udledes formel (5.1) ved at betragte de kinematiske betingelser opstillet i formel
(4.6).
⎧u , x ⎫
⎧ε xx ⎫ ⎡1 0 0 0 ⎤ ⎪ ⎪
u,
⎪ ⎪
{ε} = ⎨ε yy ⎬ = ⎢⎢0 0 0 1 ⎥⎥ ⎨⎪ v, y ⎬⎪
⎪γ ⎪ ⎢ 0 1 1 0 ⎥ ⎪ x ⎪
⎦ ⎪ v, ⎪
⎩ xy ⎭ ⎣
⎩ y⎭
(5.7)
De partielle afledte af u og v mht. x og y kan også skrives som
⎧u , x ⎫
⎡ − y1 + y3
⎪u , ⎪
⎢
⎪ y ⎪ 1 ⎢ x1 − x3
=
⎨ ⎬
⎪ v, x ⎪ 2 A ⎢ 0
⎢
⎪⎩ v, y ⎪⎭
⎣ 0
y1 − y2
0
− x1 + x2
0
0
0
− y1 + y3
x1 − x3
⎤ ⎧u ,ξ ⎫
⎪ ⎪
0 ⎥⎥ ⎪u ,η ⎪
⎨ ⎬
y1 − y2 ⎥ ⎪ v,ξ ⎪
⎥
− x1 + x2 ⎦ ⎪⎩ v,η ⎪⎭
0
(5.8)
I formel (5.8) transformeres flytninger fra referencekoordinatsystemet til det kartesiske koordinatsystem vha. den udvidede inverse Jacobimatrice.
Flytningerne i referencekoordinatsystemet er bestemt ud fra formel (4.17), hvor formfunktionerne i hhv. u-retningen og v-retningen både differentieres mht. ξ og η. Flytningerne i
referencekoordinatsystemet kan bestemmes efter formel (5.9).
81
5.1 Isoparametrisk CST-element Q3
⎧u,ξ ⎫ ⎡ N1,ξ
⎪u, ⎪ ⎢ N
⎪ η ⎪ ⎢ 1,η
⎨ ⎬=⎢
⎪ v,ξ ⎪ ⎢ 0
⎪⎩ v,η ⎪⎭ ⎣⎢ 0
0
0
N1,ξ
N1,η
⎡ −1 0
⎢ −1 0
⎢
=
⎢ 0 −1
⎢
⎣ 0 −1
1
0
0
0
0
0
N 2,ξ
N 2,η
N 2,ξ
N 2,η
0
0
0
0
1
0
0
1
0
0
N 3,ξ
N 3,η
0
0
0 ⎤
0 ⎥⎥
{d}
N 3,ξ ⎥
⎥ 6×1
N 3,η ⎦⎥
⎡ u1 ⎤
0 ⎤ ⎢⎢ v1 ⎥⎥
0 ⎥⎥ ⎢u2 ⎥
⎢ ⎥
0 ⎥ ⎢ v2 ⎥
⎥
1 ⎦ ⎢ u3 ⎥
⎢ ⎥
⎢⎣ v3 ⎥⎦
(5.9)
Ved indsættelse af formel (5.8) og (5.9) i formel (5.7) fås følgende udtryk for tøjningerne.
⎡ − y1 + y3
⎧ε xx ⎫ ⎡1 0 0 0 ⎤
⎢
1 ⎢ x1 − x3
⎪ ⎪ ⎢
⎥
⎨ε yy ⎬ = ⎢0 0 0 1 ⎥
⎪γ ⎪ ⎢0 1 1 0 ⎥ 2 A ⎢ 0
⎢
⎦
⎩ xy ⎭ ⎣
⎣ 0
y1 − y2
0
− x1 + x2
0
0
− y1 + y3
0
x1 − x3
⎤ ⎡ −1 0
0 ⎥⎥ ⎢ −1 0
⎢
y1 − y2 ⎥ ⎢ 0 −1
⎥⎢
− x1 + x2 ⎦ ⎣ 0 −1
0
1 0 0
0 0 1
0 1 0
0 0 0
[B]
⎡ u1 ⎤
0 ⎤ ⎢⎢ v1 ⎥⎥
0 ⎥⎥ ⎢u2 ⎥
⎢ ⎥
0 ⎥ ⎢ v2 ⎥
⎥
1 ⎦ ⎢ u3 ⎥
⎢ ⎥
⎢⎣ v3 ⎥⎦
Udregnes de tre matricer, der udgør tøjnings-flytningsmatricen [B], kan tøjningerne i et
trekantelement bestemmes af formel (5.10).
⎡ε xx ⎤
⎡ − y3 + y2
⎢ ⎥ 1 ⎢
⎢ε yy ⎥ = 2 A ⎢ 0
⎢γ xy ⎥
⎢⎣ x3 − x2
⎣ ⎦
0
− y1 + y3
0
y1 − y2
x3 − x2
− y3 + y2
0
x1 − x3
x1 − x3
− y1 + y3
0
− x1 + x2
⎡ u1 ⎤
⎢v ⎥
0 ⎤⎢ 1⎥
⎢u ⎥
− x1 + x2 ⎥⎥ ⎢ 2 ⎥ (5.10)
v
y1 − y2 ⎥⎦ ⎢ 2 ⎥
⎢ u3 ⎥
⎢ ⎥
⎢⎣ v3 ⎥⎦
Tøjnings-flytningsmatricen [B] i det kartesiske koordinatsystem er derfor givet som
⎡ − y3 + y 2
1 ⎢
[ B] = ⎢ 0
2A
⎣⎢ x3 − x2
0
− y1 + y3
0
y1 − y2
x3 − x2
0
x1 − x3
0
− y3 + y 2
x1 − x3
− y1 + y3
− x1 + x2
⎤
− x1 + x2 ⎥
⎥
y1 − y2 ⎦⎥
0
Udledningen af [B] fremgår af matlab-fil n12.
Elementets stivhedsmatrice [k] kan nu bestemmes vha. formel (4.24) og den konstitutive
matrice vist i formel (4.3). I praksis benyttes Calfem funktionen plante.
82
5 Elementtyper
5.2 Isoparametrisk Q8-element
I et isoparametrisk element defineres knudekoordinaterne [x y] og knudeflytningerne [u v]
ud fra den samme formfunktion [N]. Isoparametriske elementer overføres derfor til et referencekoordinatsystem (ξ,η). Elementet med otte knuder gennemgået her kan have buede
sider, beskrevet ved parabler, som følge af den ekstra knude placeret på hver rand. Prisen
er en koordinattransformation vha. en Jacobimatrice [J], og at numerisk integration er påkrævet ved bestemmelse af stivhedsmatricen [k]. [Cook, et. al., 2002]
Figur 5.2. 8-knudet element i xy-system (tv), samme element i ξη-system (th).
Nummereringen af knuderne på figur 5.2 er ikke den eneste mulighed, men i følgende afsnit er udtrykket skrevet specifikt til denne nummerering.
Overført til ξη-systemet vil et element uanset formen være et kvadrat med sidelængden to,
figur 5.2. Idet elementet er isoparametrisk benyttes samme formfunktion til bestemmelse af
koordinater og flytninger, formel (5.11).
⎧⎪ x ⎫⎪ ⎪⎧ ∑ N i xi ⎪⎫
⎨ ⎬=⎨
⎬ = [ N ]{c}
⎩⎪ y ⎭⎪ ⎪⎩∑ N i yi ⎪⎭
⎪⎧u ⎪⎫ ⎪⎧∑ N i ui ⎪⎫
⎨ ⎬=⎨
⎬ = [ N ]{d}
⎩⎪ v ⎭⎪ ⎪⎩ ∑ N i vi ⎪⎭
{c} = [ x1
y1
x2
y2
…
{d} = [u1
v1
u2
v2
…
16×1
16×1
⎡ N1
N
=
[ ] ⎢
16×2
⎣⎢ 0
y7
v7
x8
y8
u8
v8
0
N2
0
N3
…
N8
N1
0
N2
…
N7
0
(5.11)
]T
(5.12)
]T
0⎤
⎥
N 8 ⎦⎥
(5.13)
Flytningsfeltet, valgt ud fra Pascals trekant, tabel 4.1, opskrives nu som matricen [X], og
heraf bestemmes formfunktioner vha. formel (5.14).
83
5.2 Isoparametrisk Q8-element
[ N ] = [ X ] ⎡⎣ A −1 ⎤⎦
[ X ] = ⎡⎣1
ξ
(5.14)
η
ξ2
ξη
η2
ξ 2η
8×1
ξη 2 ⎤⎦
(5.15)
Matricen [A] opskrives vha. [X] og koordinaterne indsættes i henhold til figur 5.2.
⎡1
⎢
⎢1
⎢
⎢1
⎢
⎢1
A
=
[ ] ⎢
⎢1
⎢
⎢1
⎢
⎢1
⎢
⎢
⎣1
−1
−1
1
1
1
−1
1
−1
1
−1
1
−1
1
1
1
1
1
1
−1
1
1
−1
1
1
0
−1
0
0
1
0
1
0
1
0
0
0
0
1
0
0
1
0
−1
0
1
0
0
0
−1⎤
⎥
1⎥
⎥
1⎥
⎥
−1⎥
⎥
0⎥
⎥
0⎥
⎥
0⎥
⎥
⎥
0⎦
(5.16)
Matricen inverteres således at formfunktionerne derefter kan bestemmes.
⎡ −0, 25
⎢
⎢ 0
⎢
⎢ 0
⎢
⎢ 0, 25
−1
[A] = ⎢
⎢ 0, 25
⎢
⎢ 0, 25
⎢
⎢ −0, 25
⎢
⎢
⎣ −0, 25
−0, 25
−0, 25
−0, 25
0,5
0,5
0,5
0
0
0
0
0,5
0
0
0
0
−0,5
0
0,5
0, 25
0, 25
0, 25
−0,5
0
−0,5
−0, 25
0, 25
−0, 25
0
0
0
0, 25
0, 25
0, 25
0
−0,5
0
−0, 25
0, 25
0, 25
0,5
0
−0,5
0, 25
0, 25
−0, 25
0
−0,5
0
0,5 ⎤
⎥
−0,5⎥
⎥
0 ⎥
⎥
0 ⎥
⎥ (5.17)
0 ⎥
⎥
−0,5⎥
⎥
0 ⎥
⎥
⎥
0,5 ⎦
Følgende formfunktioner kan nu bestemmes for Q8-element vha. formel (5.14).
N1 = −
1
(1 − ξ )(1 − η )(1 + ξ + η )
4
N5 =
1
1 − ξ 2 ) (1 − η )
(
2
N2 = −
1
(1 + ξ )(1 − η )(1 − ξ + η )
4
N6 =
1
(1 + ξ ) (1 − η 2 )
2
1
N 3 = − (1 + ξ )(1 + η )(1 − ξ − η )
4
N4 = −
84
1
(1 − ξ )(1 + η )(1 + ξ − η )
4
1
N 7 = (1 − ξ 2 ) (1 + η )
2
N8 =
1
(1 − ξ ) (1 − η 2 )
2
(5.18)
5 Elementtyper
5.2.1 Jacobi Transformation
For at kunne opskrive udtrykket [k(x,y)] skal der transformeres mellem afledte mht. xy og
ξη, hvor sammenhængen, beskrevet i formel (5.4), er gældende.
Ud fra formel (5.4), (5.5) og formfunktionerne i (5.18) kan [J] opstilles for et Q8-element.
Som eksempel opstilles her kun for formfunktionen N1.
⎡ 14 (ξ − ξ 2 ) + 12 (η − ηξ )
[J] = ⎢ 1
⎢ 4 (η − η 2 ) + 12 (η − ηξ )
⎣
⎡ x1
⎢
…⎤ ⎢ x2
⎥⎢
…⎥⎦ ⎢
⎢
⎢⎣ x8
y1 ⎤
⎥
y2 ⎥ ⎡ J 11
⎥=⎢
⎥ ⎢⎣ J 21
⎥
y8 ⎥⎦
J 12 ⎤
⎥
J 22 ⎥⎦
(5.19)
Sammenhængen mellem tøjninger og flytninger er tidligere blevet opstillet generelt, formel
(5.20). I formel (5.21) beskrives sammenhængen mellem tøjninger og flytninger vha. en
matrice ud fra de kinematiske betingelser beskrevet i afsnit 4.1.2.
{ε} = [ B]{d}
3×1
(5.20)
3×16 16×1
⎧ε xx ⎫ ⎡1
⎪ ⎪ ⎢
{ε} = ⎪⎨ε yy ⎪⎬ = ⎢0
⎪ ⎪ ⎢
⎩⎪γ xy ⎭⎪ ⎢⎣ 0
0
0
0
0
1
1
⎧u , x ⎫
0⎤ ⎪ ⎪
⎥ ⎪⎪u , y ⎪⎪
1⎥ ⎨ ⎬
⎥ ⎪ v, x ⎪
0 ⎥⎦ ⎪ ⎪
⎪⎩ v, y ⎭⎪
(5.21)
Formel (5.20) omskrives nu mht. [J] således, at tøjningerne i xy-systemet kan bestemmes
vha. formfunktionerne opstillet i ξη-systemet, hvor [0] er en 4×4 nulmatrice.
⎧u , x ⎫
⎪ ⎪
⎪⎪u , y ⎪⎪ ⎡ J -1
⎨ ⎬=⎢
⎪ v, x ⎪ ⎢⎣ 0
⎪ ⎪
⎪⎩ v, y ⎪⎭
⎧u,ξ ⎫ ⎡ N1,ξ
⎪ ⎪ ⎢
⎪⎪u,η ⎪⎪ ⎢ N1,η
⎨ ⎬=⎢
⎪ v,ξ ⎪ ⎢ 0
⎪ ⎪ ⎢
⎪⎩ v,η ⎪⎭ ⎢⎣ 0
⎧u ,ξ ⎫
⎪ ⎪
0 ⎤ ⎪⎪u ,η ⎪⎪
⎥⎨ ⎬
J -1 ⎥⎦ ⎪ v,ξ ⎪
⎪ ⎪
⎪⎩ v,η ⎪⎭
(5.22)
0
N 2,ξ
0
N 3,ξ
N 8,ξ
0
N1,η
0
N1,η
N 8,η
N1,ξ
0
N 2,ξ
0
0
N1,η
0
N 2,η
0
0
0 ⎤
⎥
0 ⎥
⎥ {d}
N 8,ξ ⎥ 16×1
⎥
N 8,η ⎥⎦
(5.23)
Formel (5.21), (5.22) og (5.23) sammenskrives nu således, at resultatet bliver udtrykt som i
formel (5.20). For udregning af [B] henvises til CD-bilag matlab-fil n13
85
5.2 Isoparametrisk Q8-element
⎡1
⎢
{ε} = ⎢0
⎢
⎢⎣0
0
0
0
0
1
1
⎡ J -1
0⎤ ⎢
⎥⎢ 0
1⎥ ⎢
⎥⎢ 0
⎥
0⎦ ⎢
⎢0
⎣
0
0
J -1
0
0
J -1
0
0 ⎤ ⎡ N1,ξ
⎥⎢
0 ⎥ ⎢ N1,η
⎥⎢
⎥
0 ⎢ 0
⎥⎢
-1 ⎥ ⎢
J ⎦⎣ 0
0
8×8
0
0
N1,ξ
N1,η
0 ⎤
⎥
0 ⎥
⎥ {d} (5.24)
⎥
N8,ξ
⎥
⎥
N8,η ⎦
⎡ B( x , y ) ⎤
⎣
⎦
hvor [J]-1 er bestemt efter formel (5.6)
I den videre notation sættes det[J] lig J. Formel (4.24) opstillet for stivhedsmatricen kan nu
benyttes, og et udtryk for stivhedsmatricen opstilles mht. det isoparametriske element.
T
⎡ k ( x , y ) ⎤ = ∫ ∫ ⎡ B( x , y ) ⎤ [ E] ⎡ B( x , y ) ⎤ t dx dy =
⎣
⎦
⎣
⎦
⎣
⎦
1 1
∫
T
∫ ⎡⎣ B(ξ ,η ) ⎤⎦ [E] ⎡⎣B(ξ ,η ) ⎤⎦ t J dξ dη
(5.25)
−1 −1
I praksis benyttes Calfem funktionen plani8e til opstilling og løsning af stivhedsmatricen
[k].
5.2.2 Numerisk integration
Prisen, for at kunne opskrive isoparametrisk, er en koordinattransformation vha. Jacobimatricen [J], og at numerisk integration er påkrævet ved bestemmelse af stivhedsmatricen
[k], udtrykket er for komplekst til at kunne løses analytisk. Her benyttes Gaussintegration
og en omskrivning af formel (4.24). Ved Gaussintegration samles værdier ved en række
valgte punkter og summeres for at bestemme det endelige integrationsresultat I. Nøjagtigheden af summationen afhænger af antallet af opsamlingspunkter og kompleksiteten af det
udtryk, der integreres.
1
I = ∫ φ d ξ ≈ W1φ1 + W2φ2 + … + Wnφn
(5.26)
−1
hvor
W er udtrykket for ”vægten” af de enkelte Gausspunkter
φ er funktionen, der integreres, transformeret til ξη-system
Vægtningsfaktoren W indsættes for at bestemme et nøjagtigt resultat med færre målinger,
og varierer alt efter, hvor mange Gausspunkter der benyttes. Således er vægten for eksempel lig to for 1. ordens og en for 2. ordens integration.
Det 8 knudet isoparametriske element integreres i to dimensioner, og derfor benyttes numerisk integration ad to omgange som vist i formel (5.27).
86
5 Elementtyper
1 1
I=
1
⎤
i
−1 −1
≈
⎡
∫ ∫ φ (ξ ,η ) dξ dη ≈ ∫ ⎢⎣∑W φ (ξ ,η )⎥⎦ dη
i
i
−1
⎡
⎤
Wi W j φ (ξi ,η j )
∑j W j ⎢⎣∑i Wi φ (ξi ,η j )⎥⎦ = ∑∑
i
j
(5.27)
Fuld integration af et Q8-element kræver fire Gausspunkter som vist på figur 5.3, og er derfor af 2. orden, fordi der er 2×2 punkter placeret symmetrisk. Hvert enkelt punkt vægtes
således, at WiWj er lig én. Det skyldes, at der arbejdes med 2. orden. Her indsættes udtrykket, der ønskes integreret. Således [k] bliver udtrykt ved en omskrivning af formel (5.25)
til formel (5.28).
2
2
T
⎡ k ( x , y ) ⎤ = J ⋅ t ⋅ ∑∑ ⎡ B(ξ ,η ) ⎤ [ E] ⎡ B ⎤
⎣
⎦
⎣
⎦
⎣ (ξ ,η ) ⎦
i =1 j =1
(5.28)
Figur 5.3. Gausspunkter for numerisk integration af Q4 og Q8.
Stivhedsmatricen kan nu opstilles for et Q8-element. Placering af Gausspunkter er valgt
således, at der for 1. orden evalueres i (0,0), 2. orden i (ξ,η) = ±1/√3 og 3. orden i (ξ,η) =
±√(0,6). Det bemærkes ved valg af punkter, at nøjagtigheden af Gaussintegrationen ikke
stiger såfremt, der benyttes flere punkter end krævet for fuld integration.
87
5.3 Isoparametrisk Q4-element
5.3 Isoparametrisk Q4-element
Isoparametriske elementer med 4 knuder udledes som et Q8 element med den undtagelse
af antallet af formfunktioner er færre og flytningerne interpoleres lineært. Det er nødvendigt at benytte en numerisk integration med 2×2 Gausspunkter. Udledningen af [B], der er
en 3×8 matrice, gennemføres her kort med de vigtigste resultater. Derefter kan formel
(5.28) benyttes til at bestemme [k].
⎧⎪ x ⎫⎪ ⎪⎧ ∑ N i xi ⎪⎫
⎨ ⎬=⎨
⎬ = [ N ]{c}
⎪⎩ y ⎪⎭ ⎪⎩ ∑ N i yi ⎪⎭
⎪⎧u ⎪⎫ ⎪⎧ ∑ N i ui ⎪⎫
⎨ ⎬=⎨
⎬ = [ N ]{d}
⎪⎩u ⎪⎭ ⎪⎩ ∑ N i vi ⎪⎭
{c} = [ x1
y1
x2
y2
x3
{d} = [u1
v1
u2
v2
u3
y3
y4 ]
T
x4
v3
u4
v4 ]
(5.30)
T
⎡ N1
0
N2
0
N3
0
N4
⎣⎢ 0
N1
0
N2
0
N3
0
[N] = ⎢
(5.29)
0⎤
⎥
N 4 ⎦⎥
(5.31)
Formfunktionerne kan opskrives individuelt, formel (5.32).
N1 =
1
(1 − ξ )(1 − η )
4
1
N 3 = (1 + ξ )(1 − η )
4
N2 =
1
(1 + ξ )(1 + η )
4
(5.32)
1
N 4 = (1 − ξ )(1 + η )
4
Jacobimatricen kan nu bestemmes vha. formel (5.5).
[J] =
1 ⎡ − (1 − η )
⎢
4 ⎢ − (1 − ξ )
⎣
(1 − η )
− (1 + ξ )
(1 + η )
(1 + ξ )
⎡ x1
⎢
− (1 + η ) ⎤ ⎢ x2
⎥⎢
(1 − ξ )⎥⎦ ⎢
⎢
⎢⎣ x4
y1 ⎤
⎥
y2 ⎥ ⎡ J 11
⎥=⎢
⎥ ⎢⎣ J 21
⎥
y4 ⎥⎦
J 12 ⎤
⎥
J 22 ⎥⎦
(5.33)
Herefter benyttes fremgangmåden, hvor formel (5.21), (5.22) og (5.23) kombineres for, at
kunne opskrive et udtryk for [B]. Stivhedsmatricen kan nu beregnes efter formel (5.28).
⎡1
⎢
{ε} = ⎢0
⎢
⎢⎣0
0
0
0
0
1
1
⎡ N1,ξ
⎢
0⎤
⎥ ⎡ J -1
0 ⎤ ⎢ N1,η
⎥⎢
1⎥ ⎢
-1 ⎢
⎥ ⎢⎣ 0
J ⎥⎦ 0
⎢
⎥
×
4
4
0⎦
⎢ 0
⎣
⎡B( x , y ) ⎤
⎣
⎦
88
0
0
N1,ξ
N1,η
0 ⎤
⎥
0 ⎥
⎥ {d}
N 4,ξ ⎥
⎥
N 4,η ⎥⎦
(5.34)
5 Elementtyper
5.4 Elementsvagheder
De benyttede elementer lider alle under en række mangler, som en bruger af elementet bør
være opmærksom på. Usikkerheden af spændingen bestemt med FEM er generelt beskrevet i afsnit 4.3.1.
FEM benytter ofte elementer, som er stivere end virkeligheden. Dette er især udtalt for Q3
og Q4-elementer. For eksempel er Q4 ude af stand til at simulere ren bøjning. Bøjning af
Q4-elementet medfører både forskydning og normalspændinger og en deformation som beskrevet i figur 5.4 (th.). Denne falske forskydning forbruger energi, og gør elementet stivere end forudsat. Denne proces kaldes på engelsk for locking.
M
M
M
M
Figur 5.4. Deformation af en blok udsat for moment (tv), og deformation af et Q4 element udsat for tilsvarende belastning (th.).
Udover locking, der opstår pga. en falsk forskydningsspænding, findes der andre måder et
net af elementer kan låse, men disse vil ikke blive diskuteret i dette projekt. Generelt kan
locking beskrives som en overdreven stivhed i en eller flere deformationsformer. Elementer vil sjældent låse således, at deformation er umulig.
Konvergering af et system af elementer er ikke umuligt, selv om elementerne låser undervejs. I stedet vil låsen betyde, at elementerne ikke kan benyttes til simulering, inden diskretiseringen er tilstrækkelig fin.
I forbindelse med bøjningsproblemer kan Q8 elementer benyttes til at undgå locking. Det
skyldes, at [N] interpoleres vha. 2. gradspolynomier således, at Q8 ikke ligesom f.eks.
CST-elementer vil medføre ”falske” forskydningsspændinger, når elementet udsættes for
bøjning. Dette skyldes, at 2. gradspolynomier i [N] bedre kan tilnærmes den bøjede deformation.
89
5.4 Elementsvagheder
5.4.1 Spændingsberegning
De generelle usikkerheder ved beregning af spændinger er tidligere blevet omtalt i afsnit
4.3.1. I det følgende nævnes kort usikkerheden ved interpolering af spændinger. Spændingerne kan ligesom flytningerne interpoleres over elementet vha. formfunktionerne under
forudsætning af, at der tages hensyn til elementdefekter. I det følgende beskrives, hvor
spændingerne generelt bør beregnes i Q4 og Q8.
• Q4 vil ved bøjning være defekt som gennemgået tidligere. Derfor bør spændingen
først beregnes i punktet ξ = η = 0, figur 5.3, og derefter interpoleres ud over elementet.
• For Q8 vil spændingerne, beregnet i de fire Gausspunkter, være de mest korrekte, og
derefter kan spændingerne interpoleres over hele elementet. Spændingsbilledet bestemt ved denne metode er jævnet, og der er ingen garanti for, at alle spændinger kan
interpoleres korrekt.
Som udgangspunkt vælges det at gå ud fra en konstant værdi af spændingerne over de enkelte elementer, med det forbehold, at spændingsbilledet først kan regnes tilnærmelsesvis
korrekt med en fin diskretisering.
90
6
Numeriske modeller
I det følgende foretages numeriske beregninger af flytningerne for følgende forsøgsemner
• Homogen og porøs cirkelskive
• Homogen og porøs cirkelring
• Porøs cirkelskive opbygget af superelementer
De numeriske modeller af forsøgsemnerne opbygges af elementer, der bygger på interpolation af flytningsfeltet iht. elementmetoden beskrevet i afsnit 4 . Modellerne løses numerisk
med MatLab og den dertil hørende toolbox Calfem. Programmeringen af de enkelte modeller findes på CD-bilag matlab-fil n2 og n4, og bygger på funktioner defineret i MatLab og
Calfem [Calfem, 1999]. Grundlæggende programmeringskommandoer anvendt i Calfem
fremgår af appendiks A.1. Principielt opbygges modellerne af cirkelringen og -skiven som
vist på figur 6.1. Opbygning og brug af superelementet beskrives i afsnit 6.3.
F
F
N
R
F
R
F
Figur 6.1. Eksempler på numeriske modeller af cirkelskiven (th.) og –ringen (tv.) med angivelse af randbetingelser, belastning og princip for opdeling i elementer (N×R).
Det vælges at udnytte symmetrien omkring både den vandrette og lodrette diagonal i cirkelskiven og –ringen, hvorfor der kun modelleres en kvart cirkelskive / -ring. Det har den
fordel, at antallet af frihedsgrader reduceres, hvorved beregningstiden mindskes. Ved kun
91
6.1 Konvergensanalyse
at modellere et kvart forsøgsemne skal det sikres, at randbetingelserne svarer til, at hele
forsøgsemnet modelleres. Det betyder, at elementerne i symmetriakserne understøttes simpelt og således, at flytninger vinkelret på symmetriakserne er lig nul, hvilket er gældende,
så længe deformationen er lineær [Cook, et. al., 2002]. Endvidere skal kun den halve last F
påføres i modellen.
6.1 Konvergensanalyse
Under opbygningen af de numeriske modeller af cirkelringen og -skiven er der benyttet tre
forskellige elementtyper. Der er i det følgende gennemført en konvergensanalyse for at udpege hvilken af disse elementer, der er mest anvendelig i det videre arbejde med de numeriske modeller.
Det er valgt at afgrænse konvergensanalysen til kun at omhandle cirkelringen, idet det vurderes, at det valgte element ligeledes er det mest anvendelige for cirkelskiven. Det skyldes,
at den eneste forskel mellem cirkelskiven og -ringen er indsættelsen af elementer i centrum.
De numeriske modeller er uanset elementtypen opbygget vha. en drejning af en række elementer omkring centrum. Derfor vil der i konvergensanalysen ikke blive skelnet mellem
forskellige opbygninger af modellerne. Konvergensanalysen foretages med tre forskellige
elementtyper, som er Q3, Q4 og Q8. Se også teorien bag elementerne i afsnit 5 .
Til bestemmelse af hvilket element, der er mest anvendeligt til modellering af forsøgsemnerne, gennemføres konvergensanalysen for nedbøjningen som funktion af antal knuder,
og som funktion af beregningstiden. Dvs. det element, der giver en konstant nedbøjning
ved færrest elementer og kortest beregningstid, bør vælges. Nedbøjningen måles i punktet
vist på figur 6.2.
Målepunkt
F
N=4
CL
R=4
Figur 6.2. Skitse af målepunktets placering og inddelingen (tv.) samt nedbøjningen δ for cirkelringen (th.).
92
6 Numeriske modeller
I konvergensanalysen inddeles cirkelringen i lige mange inddelinger over radien R som
over vinklen N, hvilket er illustreret på figur 6.2. For hver beregningsgennemgang øges
inddelingerne ved at addere med den samme konstant. Nedbøjningen bestemmes for en last
F på 10 kN, der påføres lodret i den øverste knude, figur 6.2.
Af tabel 6.1 fremgår hardware og software specifikationerne anvendt ved beregningerne.
Tabel 6.1. Hardware og software specifikationer anvendt til konvergensanalysen.
Processor
RAM
MatLab
CalFem
Pentium 4; 2,80 GHz
2,0 Gb
Version 7.0.1
Version 3.4
Resultaterne af konvergensanalysen angives i tabel 6.2 og grafisk på figur 6.3 og figur 6.4.
Tabel 6.2. Angivelse af konvergeringen for de tre elementer.
Inddelinger
N
R
[-]
[-]
0
5
10
15
20
25
30
0
5
10
15
20
25
30
Q3
Knuder / Tid
Elementer
[s]
0/0
36 / 50
121 / 200
256 / 450
441 / 800
676 / 1250
961 / 1800
0
0,1
0,3
0,6
1,5
3,2
6,6
δ
[mm]
0
0,0075
0,0083
0,0085
0,0085
0,0086
0,0086
Q4
Knuder / Tid
Elementer
[s]
0/0
121 / 100
441 / 400
961 / 900
0
0,4
1,8
5,8
δ
[mm]
0
0,0084
0,0086
0,0086
Q8
Knuder / Tid
Elementer
[s]
0/0
341 / 100
1281 / 400
2821 / 900
0
0,5
2,8
12,4
δ
[mm]
0
0,0086
0,0086
0,0086
0,0088
0,0086
Nedbøjning [mm]
0,0084
0,0082
Q3
0,008
Q4
Q8
0,0078
0,0076
0,0074
0,0072
0
500
1000
1500
2000
2500
3000
Antal knuder
Figur 6.3. Konvergenskurve for Q3, Q4 og Q8 med nedbøjningen som funktion af antal knuder.
93
6.1 Konvergensanalyse
0,009
Nedbøjning [mm]
0,0085
0,008
Q3
Q4
Q8
0,0075
0,007
0,0065
0
1
2
3
4
5
6
7
8
Tid [s]
Figur 6.4. Konvergenskurve for Q3, Q4 og Q8 med nedbøjningen som funktion af beregningstiden.
Ved at sammenholde resultaterne af nedbøjningen for de tre elementtyper fremgår det af
tabel 6.2 og figur 6.3, at Q8-elementet konvergerer ved ca. 350 knuder, mens Q3 og Q4elementerne konvergerer ved hhv. 450 og 700 knuder. At Q8-elementet konvergerer ved en
lavere diskretisering end Q3 og Q4-elementerne stemmer overens med tidligere beskrevne
elementsvagheder, afsnit 5.4, idet locking er mere udpræget for disse end for Q8-elementet.
Problemet med locking er i denne konvergensanalyse ikke af særlig stor betydning, da diskretiseringen af Q3 og Q4-elementerne ikke behøver at være væsentlig finere end Q8elementet før end problemet med locking er elimineret. Det vil sige, at der er opnået konvergens. Udover antal knuder som funktion af nedbøjningen er beregningstiden som funktion af nedbøjningen bestemt. Af denne ses det, at Q8-elementet konvergerer efter et tidsforløb på 0,5 s, hvorimod det tager hhv. 3,2 og 1,8 s før end Q3 og Q4-elementerne konvergerer. Resultatet af tidsforløbet skal dog anvendes med forbehold, da beregningerne eventuelt kan effektiviseres ved ændring af programmeringen.
Da beregningstiden for konvergeringen af de enkelte elementer ikke er væsentlig, er det
valgt at anvende Q4-elementet i de endelige modelleringer af cirkelskiven og –ringen. Det
begrundes med, at modelleringen af forsøgsemnerne med Q4-elementet er lettere end Q8elementet samtidig med, at det konvergerer før Q3-elementet.
94
6 Numeriske modeller
6.2 Numeriske resultater
Cirkelskiven og -ringen opbygges iht. konvergensanalysen af standard isoparametriske Q4elementer, hvor centrum af cirkelskiven dog opbygges af Q3-elementer, da det ikke kan
lade sig gøre med Q4. Kompatibiliteten af skiven bringes ikke i fare ved at bruge to forskellige elementer, idet formfunktionerne er lineære for begge elementer. Herudover modelleres superelementet til en fuldstændig numerisk modellering af den porøse cirkelskive.
Ved opbygning af den porøse skive vha. superelementet vælges det at anvende Q8elementer, idet disse er bedre til at modellere den indre rand i superelementet, figur 6.11.
For at sikre at modellerne konvergerer mod den eksakte løsning, er det på baggrund af
konvergensanalysen valgt at foretage beregningerne med en diskretisering på 60×60. Ved
denne diskretisering elimineres endvidere problemet med locking.
6.2.1 Cirkelskive
Som tidligere skrevet er cirkelskiven opbygget af Q3 og Q4-elementer med en opdeling på
60×60 elementer, svarende til 2.000 knuder. Opbygningen af skiven ses bedst i matlab-fil
n2 på CD-bilag, men et overblik fås i figur 6.5, hvor nedbøjningen samt spændingsbilledet
af en kvart skive er illustreret.
Figur 6.5. Skitse af cirkelskiven opdelt i elementer samt nedbøjningen for en 20×20 diskretisering og en last
på 50 kN (tv.). Bemærk, at nedbøjningen er skaleret med en faktor 100, hvor sort = udeformeret og rød =
deformeret. Plot af normalspændingen σyy [MPa] (th.) ved en last på 50 kN og en opdeling på 60×60. Bemærk, at der ikke er brugt samme diskretisering i de to figurer.
Svagheden ved den numeriske model er den store deformation omkring det punkt, hvor
kraften påføres. Denne usikkerhed kommer ligeledes til udtryk i spændingsbilledet, og derfor må aflæsninger i dette område benyttes med forsigtighed.
Til senere sammenligning af de lodrette flytningerne bestemt analytisk og ved forsøg er
den lodrette flytning både for den homogene og porøse cirkelskive som funktion af kraften
opstillet grafisk ved figur 6.6 og som et retliniet polynomium, formel (6.1). Den lodrette
flytning er bestemt i punktet (0; 48), hvilket svarer til placeringen af målepunkterne i for95
6.2 Numeriske resultater
søgsemnerne. Derfor får den lokale deformationspåvirkning, som ovenfor er beskrevet, ikke nogen betydning.
Ved bestemmelse af flytningerne i den porøse cirkelskive er materialeparametrene i den
numeriske beregning erstattet med de porøse materialeparametre, som er bestemt i afsnit
2.4. Dette er også gældende for cirkelringen.
0,08
0,07
y = 1,652E-06x - 6,167E-08
R2 = 1,000E+00
Flytning [mm]
0,06
0,05
0,04
0,03
y = 9,480E-07x + 5,143E-08
R2 = 1,000E+00
0,02
0,01
0
0
5000
10000
15000
20000
25000
30000
35000
40000
45000
F [N]
Homogen skive, lodret flytning (0;48)
Porøs skive, lodret flytning (0; 48)
Lineær (Porøs skive, lodret flytning (0; 48))
Lineær (Homogen skive, lodret flytning (0;48))
Figur 6.6. Lodrette flytninger i punktet (0; 48) som funktion af kraften for både den homogene og porøse
cirkelskive.
Som det fremgår af udtrykkene opskrevet i (6.1) er der set bort fra konstanterne, der er angivet i udtrykkene for tendenslinien i figur 6.6. Det skyldes, at det kun er for meget små
belastninger, at konstanten får betydning. Endvidere ses det af korrelationskoefficienten,
som er lig med én (R2 = 1) for tendenslinien, at punkterne ligger på en perfekt ret linie.
uhomogen,skive = 9, 48 ⋅10−7
uporøs,skive = 16,52 ⋅10−7
mm
N
mm
N
⋅F
⋅F
(6.1)
hvor
F er trykkraften [N]
u er nedbøjningen i y-retningen i målepunktet (0; 48) [mm]
6.2.2 Cirkelring
Cirkelringen er opbygget af Q4-elementerne, der i modsætning til cirkelskiven vil være udsat for større bøjning, hvorfor der er mulighed for locking, afsnit 5.4. Opdeling med 60×60
elementer, svarende til 2.000 knuder, viser iht. konvergensanalysen, at konvergeringen er
gennemført, og locking derfor er elimineret. Opbygningen af skiven ses bedst i matlab-fil
n4 på CD-bilag, men et overblik fås i figur 6.7, hvor nedbøjning og spændingsbillede af en
kvart ring illustreres.
96
6 Numeriske modeller
Figur 6.7. Skitse af cirkelringen opdelt i elementer samt nedbøjningen for en 20×20 diskretisering og en last
på 50 kN (tv.). Bemærk, at nedbøjningen er skaleret med en faktor 100, hvor sort = udeformeret og rød =
deformeret. Plot af normalspændingen σyy [MPa] (th.) ved en last på 50 kN og en opdeling på 60×60. Bemærk, at der ikke er brugt samme diskretisering i de to figurer.
Som for den numeriske model af cirkelskiven opstår der lokalt en stor deformation det
punkt, hvor kraften påføres, men igen får denne ikke stor betydning, da den lodrette flytning er bestemt i punktet (0; 48).
Til senere sammenligning af de lodrette flytningerne bestemt analytisk og eksperimentelt
er den lodrette flytning både for den homogene og porøse cirkelring som funktion af kraften opstillet grafisk ved figur 6.8 og som et retliniet polynomium, formel (6.2).
0,18
0,16
Flytning [mm]
0,14
y = 3,477E-06x + 5,000E-08
R2 = 1,000E+00
0,12
0,1
0,08
0,06
y = 1,992E-06x + 3,509E-08
R2 = 1,000E+00
0,04
0,02
0
0
5000
10000
15000
20000
25000
30000
35000
40000
45000
50000
F [N]
Homogen ring, lodret flytning (0; 48)
Porøs ring, lodret flytning (0; 48)
Lineær (Porøs ring, lodret flytning (0; 48))
Lineær (Homogen ring, lodret flytning (0; 48))
Figur 6.8. Lodrette flytninger i punktet (0; 48) som funktion af kraften for både den homogene og porøse
cirkelring.
97
6.3 Porøs cirkelskive opbygget af superelementer
Igen ses der bort fra konstanterne i udtrykkene for den lodrette flytning af cirkelringen,
formel (6.2), med samme begrundelse som for cirkelskiven, og korrelationskoefficienten
viser, at punkterne ligger på en perfekt ret linie, figur 6.8.
uhomogen,ring = 19,92 ⋅10−7
uporøs,ring = 34, 77 ⋅10−7
mm
N
mm
N
⋅F
⋅F
(6.2)
hvor
F er trykkraften [N]
u er nedbøjningen i y-retningen i målepunktet (0; 48) [mm]
6.3 Porøs cirkelskive opbygget af superelementer
I dette afsnit bestemmes flytningen af de to punkter på den lodrette diameter i den porøse
cirkelskive vha. en numerisk beregning. Den numeriske model opbygges af et antal superelementer, der hver især er opbygget af standardelementer. Superelementerne måler 16×16
mm, og skiven modelleres som angivet på figur 6.9 (th.). Det fremgår, at modellen er meget grov på randen af cirkelskiven, hvilket vil indgå i den endelige vurdering af resultatet.
Beregningen skal opfattes som et alternativ til at benytte modellen for den homogene skive, hvor elasticitetsmodulet E og Poisson’s forhold v for det homogene materiale erstattes
med E og v for det porøse materiale. Porøsiteten skabes af et kvadratisk net af borede huller med en diameter på 8 mm, hvilket fremgår af figur 6.9.
L
Figur 6.9. Den porøse cirkelskive (tv.) og den tilhørende numeriske model heraf (th.).
Som for de øvrige numeriske modeller tages der udgangspunkt i en fjerdedel af skiven, da
skiven er symmetrisk om både den vandrette og lodrette akse. Beregningsmodellen er skitseret på figur 6.10. Det fremgår, at der indgår 19 superelementer i modellen.
98
6 Numeriske modeller
y
x
Figur 6.10. Numerisk beregningsmodel.
6.3.1 Opbygning af superelement
Superelementet opbygges på baggrund af teorien beskrevet i afsnit 4.4, idet superelementet
sammensættes af et antal standardelementer, hvorefter de indre frihedsgrader elimineres.
Herved kan den globale stivhedsmatrice for superelementet opstilles, og derefter benyttes
som lokal stivhedsmatrice i den globale model vist på figur 6.10.
Det er valgt at opbygge superelementet af Q8-elementer, da det herved forventes, at modellen konvergerer tilfredsstillende ved færre antal elementer end ved brug af Q3- og Q4elementer.
7
3
Q8
6
4
8
1
N=1
N=2
5
2
N=4
Figur 6.11. Opbygning af superelement, hvor N angiver inddelinger radiært og på randen.
Superelementet er opbygget som en funktion dels af antal inddelinger på hver rand, dels af
antal radiære inddelinger, hvilket på figur 6.11 er angivet med et N. For N = 1 haves dermed fire elementer, for N = 2 haves 16 elementer, mens der for N = 4 haves 64 elementer.
Tilsvarende gælder, at der for N = 1 er 20 knuder i superelementet, mens der for hhv. N =
2 og N = 4 er 64 og 224 knuder. Relationen mellem N og antallet af elementer hhv. knuder
kan beskrives ved formlerne (6.3) og (6.4)
Antal elementer = 4 N 2
(6.3)
Antal knuder = 12 N 2 + 8 N
(6.4)
Antallet af frihedsgrader er det dobbelte af antallet af knuder.
Ideen med at opbygge et superelement er, at opstille en stivhedsmatrice for superelementet
som helhed, idet superelementet herefter kan benyttes som et nyt element til opbygning af
99
6.3 Porøs cirkelskive opbygget af superelementer
skiven. Da antallet af frihedsgrader som funktion af N hurtigt går imod et meget stort tal,
elimineres de indre knuder således, at det kun er knuderne på randen, der er beskrevet i superelementets stivhedsmatrice. Hermed er kompatibiliteten mellem de enkelte elementer
superelementer sikret. Antallet af knuder på superelementets rand kan beskrives med relationen (6.5)
Antal knuder på rand = 8 N
(6.5)
hvilket betyder, at besparelsen af frihedsgrader svarer til 24N2. For N = 4 er besparelsen
derfor lig med 384 frihedsgrader.
Elimination af frihedsgrader
I det efterfølgende opstilles den reducerede stivhedsmatrice for et superelement. Som eksempel vælges et element med N = 2. På figur 6.12 er superelementet optegnet med nummerering af frihedsgrader. Det er valgt at nummerere hhv. de ydre og indre frihedsgrader
fortløbende efter hinanden, idet det letter arbejdet med at opstille de efterfølgende matricer.
Generelt benyttes index y for frihedsgrader, der hører til ydre knuder, mens index i benyttes for indre frihedsgrader og knuder. De ydre knuder er de knuder, der findes på den ydre
rand af superelementet.
13,14
11,12
9,10
7,8
5,6
3,4
15,16
2
1
17,18
1,2
127,128
19,20
21,22
31,32
23,24
27,28
25,26
29,30
Ydre frihedsgrader : 32
Indre frihedsgrader : 96
Figur 6.12. Nummerering af frihedsgrader.
Den reducerede stivhedsmatrice [k*yy], der indeholder informationer for frihedsgraderne på
den ydre rand af superelementet, har formen
⎡⎣k *yy ⎤⎦ = ⎡⎣k yy ⎤⎦ − ⎡⎣k iy ⎤⎦ [k ii ] ⎡⎣k iy ⎤⎦
T
100
−1
(6.6)
6 Numeriske modeller
Matricerne [kyy], [kiy] og [kii], der er delmatricer fra den globale stivhedsmatrice for superelementet, ønskes fundet. Derfor findes først matricerne [Ty] og [Ti] jf. afsnit 4.4.
{d } = ⎡⎣T ⎤⎦ { d }
(6.7)
{di } = [Ti ] { d }
(6.8)
y
32 x1
96 x1
y
32 x128 128 x1
96 x128 128 x1
Da nummereringen af frihedsgraderne er foretaget som vist på figur 6.12, betyder det, at
flytningstensoren {d} kan skrives som (6.9)
⎡{d y }⎤
⎢
⎥
{ d } = ⎢ 32d×1 ⎥
⎢{ i } ⎥
128 x1
⎢⎣ 96×1 ⎥⎦
(6.9)
Heraf følger det, at [Ty] og [Ti] har følgende form for et superelement med N = 2
⎡⎣ Ty ⎤⎦ = [ I | 0
]
[ Ti ] = [ 0
]
32 x128
96 x128
(6.10)
32 x 32 32 x 96
| I
(6.11)
96 x 32 96 x 96
hvor
[I] er identitetsmatricen med 1-taller på diagonalen
[0] er en nulmatrice
[kyy], [kiy] og [kii] kan nu bestemmes.
⎡⎣k yy ⎤⎦ = ⎡⎣Ty ⎤⎦ ⎡⎣ k ⎤⎦ ⎡⎣ Ty ⎤⎦
32 x 32
T
⎡⎣k iy ⎤⎦ = ⎡⎣Ti ⎤⎦ ⎡⎣ k ⎤⎦ ⎡⎣Ty ⎤⎦
96 x 32
T
(6.13)
96 x128 128 x128 128 x 32
⎡⎣k ii ⎤⎦ = ⎡⎣ Ti ⎤⎦ ⎡⎣ k ⎤⎦ ⎡⎣ Ti ⎤⎦
96 x 96
(6.12)
32 x128 128 x128 128 x 32
T
(6.14)
96 x128 128 x128 128 x 96
Hermed er den globale stivhedsmatrice [k] for superelementet inddelt i delmatricer, hvilket
kan illustreres således
⎡ ⎡k yy ⎤ ⎡k yi ⎤ ⎤
⎢ ⎣32 x 32⎦ ⎣32 x 96⎦ ⎥
[k ] = ⎢
⎥
128 x128
⎢ ⎡⎣k iy ⎤⎦ [k ii ] ⎥
96 x 96 ⎦
⎣ 96 x 32
(6.15)
101
6.3 Porøs cirkelskive opbygget af superelementer
Det er [kyy] der er interessant at finde, men med informationer fra de indre frihedsgrader
inkluderet. Det er her, at [k*yy] beskrevet ved formel (6.6) kommer ind i billedet. Det er
valgt ikke at udskrive hhv. den globale stivhedsmatrice [k] for superelementet og den reducerede stivhedsmatrice [k*yy], da de fylder meget. I stedet opskrives relationerne (6.16)
og (6.17) mellem knudekræfter og knudeflytninger til at eftervise, at brug af den reducerede stivhedsmatrice giver identiske resultater med den globale stivhedsmatrice.
{r}
128 x1
= [k ] {d} ⇒ {d} = [k ]
−1
128 x128 128 x1
{r } = ⎡⎣k
*
*
y
32 x1
yy
32 x 32
128 x1
{r}
⎤⎦ {d y } ⇒ {d y } = ⎡⎣k * yy ⎤⎦
32 x1
(6.16)
128 x128 128 x1
32 x1
32 x 32
−1
{r }
*
y
32 x1
(6.17)
hvor
{r} er lastvektoren tilhørende den globale stivhedsmatrice
{r*y} er lastvektoren tilhørende den reducerede stivhedsmatrice
Eftervisning af reduceret stivhedsmatrice
Til eftervisning af, at brug af den reducerede stivhedsmatrice giver korrekte resultater, benyttes opstillingen på figur 6.13. Kraften F påføres i en tilfældig knude, som i dette tilfælde
er knuden med frihedsgraderne 9 og 10. Størrelsen af F er tilfældig valgt til 10 kN. Understøtningsforholdene i opstillingen fremgår ligeledes af figuren.
F
13,14
11,12
7,8
5,6
9,10
15,16
3,4
17,18
1,2
19,20
31,32
21,22
23,24
25,26
27,28
29,30
Figur 6.13. Opstilling af model til sammenligning af resultater.
102
6 Numeriske modeller
Da kraften påføres i frihedsgrad 10, kommer de to lastvektorer {r} og {r*y} til at tage følgende form
{r}
T
= [ 0 0 0 0 0 0 0 0 0 F 0 0 0 0 0 0 0 .... 0]
1 x128
{r } = [0 0 0 0 0 0 0 0 0 F 0 .... 0]
*
T
y
1 x 32
Efterfølgende benyttes formlerne (6.16) og (6.17) til at beregne flytningerne af de ydre
knuder med den angivne størrelse af lasten F. Resultaterne er angivet i tabel 6.3.
Tabel 6.3. Flytninger i frihedsgrader for hhv. global og reduceret stivhedsmatrice for en belastning på 10 kN.
Frihedsgrad
nr.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
{d}
[mm]
{dy}
[mm]
0,01104
-0,01450
0,00930
-0,01354
0,00255
-0,01419
0,00151
-0,02276
0,00572
-0,03848
0,00992
-0,02276
0,00888
-0,01419
0,00213
-0,01354
0,01104
-0,01450
0,00930
-0,01354
0,00255
-0,01419
0,00151
-0,02276
0,00572
-0,03848
0,00992
-0,02276
0,00888
-0,01419
0,00213
-0,01354
Frihedsgrad
nr.
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
{d}
[mm]
{dy}
[mm]
0,00039
-0,01450
0,00461
-0,01160
0
0
0,00323
-0,01314
0,00572
-0,01014
0,00820
-0,01314
0,01143
0
0,00683
-0,01160
0,00039
-0,01450
0,00461
-0,01160
0
0
0,00323
-0,01314
0,00572
-0,01014
0,00820
-0,01314
0,01143
0
0,00683
-0,01160
Af tabel 6.3 fremgår det, at der fås helt identiske resultater for begge beregninger. Desuden
ses det, at flytningerne er lig nul i de frihedsgrader, hvor superelementet er fastholdt.
6.3.2 Konvergensanalyse
Der er foretaget en konvergensanalyse af flytningen af frihedsgrad 1 for en kraft F påført i
frihedsgrad 10, hvilket er illustreret på figur 6.14. I analysen varieres størrelsen N med 2.
Derved kan en passende størrelse af N bestemmes. Resultatet er angivet i tabel 6.4 og på
figur 6.15.
103
6.3 Porøs cirkelskive opbygget af superelementer
F
u1
Figur 6.14. Opstilling til konvergensanalyse.
0,025
Udbøjning [mm]
0,02
0,015
0,01
0,005
0
2
4
6
8
10
12
14
16
18
20
22
24
26
28
30
Antal inddelinger N
Figur 6.15. Konvergensanalyse for flytningen af frihedsgrad 1 med en kraft på 10 kN, som er påført i frihedsgrad 10.
Tabel 6.4 viser resultaterne fra konvergensberegningen. Den akkumulerede tid henviser til
regnetiden, det har taget at udføre beregningerne i MatLab. Det fremgår, at ændringen i
flytningen u1 ikke når at gå imod en fast værdi med N = 30.
Taget regnetiden, som for N = 30 er ca. 3 minutter, og den relativt grove opbygning af cirkelskiven i betragtning, er en opdeling af superelementet på N = 12 fundet acceptabel, og
derfor vælges det at opbygge cirkelskiven af superelementet med dette antal opdelinger,
hvilket giver 3648 frihedsgrader pr. superelement. Desuden haves ikke computerkraft til at
opstille den reducerede stivhedsmatrice for superelementet med flere inddelinger end N =
12.
104
6 Numeriske modeller
Tabel 6.4. Resultater fra konvergensberegning.
N
Antal frihedsgrader
u1
[mm]
Akkumuleret tid
[sek]
∆ tid
[sek]
2
4
6
8
10
12
14
16
18
20
22
24
26
28
30
128
448
960
1.664
2.560
3.648
4.928
6.400
8.064
9.920
11.968
14.208
16.640
19.264
22.080
0,0110
0,0145
0,0163
0,0176
0,0186
0,0194
0,0201
0,0207
0,0212
0,0216
0,0220
0,0224
0,0228
0,0231
0,0234
0,14
0,29
0,75
1,93
4,28
9,26
18,48
34,19
57,99
93,92
146,10
219,11
315,50
441,73
618,38
0,15
0,46
1,18
2,34
4,99
9,21
15,71
23,80
35,92
52,19
73,01
96,40
126,22
176,65
6.3.3 Bestemmelse af flytninger
I Matlab er cirkelskiven, optegnet på figur 6.16, modelleret, og flytningen over de 48 mm
beregnes. Beregningsgangen gennemgås ikke her, men fremgår af matlab-fil n6. Da cirkelskiven er symmetrisk om både den lodrette og vandrette diameter påføres den halve størrelse af kraften F, mens det er den halve størrelse af flytningen L over 96 mm, der beregnes. Resultatet kan sammenlignes med resultaterne fra laboratorieforsøgene og de øvrige
beregninger for de homogene modeller, hvor de homogene materialekonstanter er erstattet
med porøse.
1
2
1
2
F
u
Figur 6.16. Opstilling til beregning af flytning.
Af beregningen i Matlab fås følgende relation mellem kraften F og den samlede flytning
over 96 mm
u = 17,50 ⋅10−7
mm
N
⋅F
hvor F skal indsættes med enheden [N], mens flytningen fås med enheden [mm]
105
6.4 Delkonklusion
6.4 Delkonklusion
Det kan konkluderes ud fra de numeriske beregninger, at de lodrette flytninger stemmer
overens med det forventede, idet den homogene cirkelskive har en mindre flytning og derved en større stivhed end den porøse cirkelskive. Det samme er gældende for cirkelringen,
hvor det endvidere ses, at flytningerne bliver større end for cirkelskiven, hvilket skyldes
hullet i midten.
Tabel 6.5. Resultater for de numeriske beregninger.
Emne
Homogen cirkelskive
Porøs cirkelskive
Homogen cirkelring
Porøs cirkelring
Porøs cirkelskive opbygget af superelementer
106
Lodret flytning u
10-7[mm]
9,48 F
16,52 F
19,92 F
34,77 F
17,50 F
7
Konklusion
I denne rapport er der ud fra analytiske, numeriske samt eksperimentelle metoder blevet
foretaget beregninger og estimater af elasticitetsmodul og Poissons forhold for aluminium.
Med udgangspunkt i de bestemte materialekonstanter er gennemført flytnings- og spændingsberegninger i en cirkelskive og cirkelring.
Det følgende afsnit, der samler op rapportens resultater og konklusioner, sammenligner
hvorvidt det er lykkedes, at nå frem til de samme resultater eksperimentelt, numerisk og
analytisk.
7.1 Materialekonstanter
Materialekonstanterne for aluminium er indledningsvis bestemt ved forsøg, idet der ikke
haves oplysninger om det homogene materiale. Resultaterne ses i tabel 7.1, og er af samme
størrelsesorden som opslagsværdierne for aluminium, hvor elasticitetsmodulet E er fundet
til 72.000 MPa [Teknisk Ståbi, 2003]. Med udgangspunkt i de bestemte værdier for det
homogene materiale er det estimeret, hvorledes en porøsitet på 20 % reducerer materialekonstanterne. Resultatet fremgår ligeledes af tabel 7.1.
Tabel 7.1. Materialeparametre benyttet i rapportens modeller.
ν
E
[MPa]
[-]
Usikkerhed
[%]
Homogen
79.000
0,30
-
Porøs
45.000
0,27
9,4
De porøse materialekonstanter opgivet i tabellen er bestemt eksperimentelt og numerisk.
Poissons forhold v er således fastlagt vha. den numerisk model af enhedscellen og ud fra
den betragtning, at v falder i værdi sammen med porøsiteten. Den numeriske model er den
eneste metode, der giver et fornuftigt estimat af v. Det vælges at kombinere v med et eksperimentelt bestemt E. Dette begrundes med, at ingen af de analytiske og numeriske estimater gav samstemmende resultater for både v og E. Det dertilhørende E svarer i størrelsesorden til den cellulære metode og Dilutes estimatet af E, bestemt til ca. 50.000 MPa.
Pga. uoverensstemmelsen vælges det at arbejde med det eksperimentelle resultat på E lig
45.000 MPa, med en afvigelse på ± 9,4 %, der er bestemt iht. appendiks B.2.
107
7.2 Flytninger
7.2 Flytninger
Flytningerne er bestemt mellem de to punkter både i cirkelskiven og -ringen, hvilket fremgår af figur 7.1. Flytningerne bestemt vha. analytiske og numeriske modeller sammenlignes med flytningerne målt eksperimentelt. Her benyttes de materialekonstanter som fremgår af tabel 7.1.
Stålklods
Forsøgsemne
96
96
32
32
Stålklods
Figur 7.1. Viser placering af de to flytningsmålepunkter. Alle mål i [mm].
Sammenligningen er lavet for at undersøge, hvor meget resultaterne, fundet ved de analytiske og numeriske beregninger, afviger fra flytningerne målt eksperimentelt i laboratoriet.
Alle målte og beregnede flytninger, samt afvigelser herimellem, er opstillet i tabel 7.2.
Tabel 7.2. Viser resultaterne af flytningerne mellem de to punkter langs den lodrette diagonal. Derudover er
usikkerheden på de eksperimentelle forsøg opstillet, samt afvigelserne mellem forsøg og beregningsmodellerne.
Forsøgsemne
Fejl
[%]
Flytningsudtryk
10-7 [mm]
Afvigelse
[%]
Homogen cirkelskive
Eksperimentelt
Analytisk
Numerisk
3,1
-
9,88 F
9,65 F
9,48 F
-2,3
-4,0
Porøs cirkelskive
Eksperimentelt
Analytisk
Numerisk
Superelement
3,1
-
16,30 F
16,80 F
16,52 F
17,50 F
3,1
1,3
7,4
Eksperimentelt
3,1
20,70 F
-
-
9,52 F
19,92 F
-54,0
-3,8
3,1
-
42,40 F
16,64 F
34,77 F
-60,1
-18,0
Homogen cirkelring
Porøs cirkelring
108
Analytisk
Numerisk
Eksperimentelt
Analytisk
Numerisk
7 Konklusion
Homogen cirkelskive
Ved sammenligningen ses, at flytningerne, beregnet analytisk og numerisk i den homogene
cirkelskive, er henholdsvis 2,3 og 4,0 % mindre end flytningen målt eksperimentelt. Afvigelsen på 4 % ligger uden for de beregnede fejlkilder på 3,1 %, men grundet den lille difference er det vurderet, at beregningerne er acceptable.
Porøs cirkelskive
De analytiske og numeriske beregninger af den porøse cirkelskive viser en afvigelse til den
modsatte side i forhold til den homogene cirkelskive. Afvigelsen ligger inden for det acceptable, dog ikke for superelementet, som afviger fra modellen med 7,4 %.
Homogen cirkelring
Ved den homogene cirkelring ses det, at den analytiske afvigelse er 54 % i forhold til den
eksperimentelle flytning, hvilket ikke er acceptabelt. På grund af den lille forskel mellem
eksperiment og numerisk model må dette skyldes, at den analytiske model er dårlig.
Porøs cirkelring
Ud fra konklusionerne for den homogene cirkelring er det forventet, at den analytiske model heller ikke med Eporøs kan bruges, hvilket indses med en afvigelse på 60 %. Den numeriske beregning af cirkelringen giver en nedbøjning, der afviger 18 % i forhold til de eksperimentelle beregninger. Dette er ikke acceptabelt, sammenholdt med de forudsatte fejl,
tabel 7.2.
Opsamling
Samlet kan følgende konkluderes ud fra en sammenligning af de forskellige forsøg, analytiske og numeriske modeller.
I sammenligningerne er der set bort fra den analytiske model af cirkelringen, idet at fejlen
afviger for meget, set i forhold til de resterende resultater. Dette skyldes, at der er gættet på
et flytningsfelt, der beskriver det virkelige flytningsfelt dårligt, og som derfor skal forbedres.
Sammenholdes afvigelserne fra den porøse cirkelskive med den homogene cirkelskive ses
det, at den porøse cirkelskives nedbøjning er den eneste, der er større end den eksperimentelle nedbøjning. Det stemmer ikke overens med den numeriske model, da det vides fra
elementmetoden, at denne giver resultater, der er stivere end virkeligheden, hvilket svarer
til en mindre nedbøjning. Derfor er der noget der tyder på, at det valgte porøse elasticitetsmodul er valgt for lavt således, at et af de estimerede porøse elasticitetsmoduler på omkring 50.000 MPa, der i stedet skulle være valgt.
Superelementet er benyttet til en grov modellering af cirkelskiven, og det ses, at nedbøjningen er uden for det acceptable og forudsætter et E mindre end Eporøs. Her skal tages forbehold for, at opbygningen med superelementer er grov i forhold til skivens form.
Foretages der en sammenholdelse af de numeriske afvigelser mellem den homogene og
porøse cirkelring, som gjort for cirkelskiven, ses det, at nedbøjningen for den porøse cirkelring i dette tilfælde ikke bliver større end den eksperimentelle, som nedbøjningen for
den porøse cirkelskive. Ud fra dette kan sammensluttes, at den porøse cirkelring ikke sva109
7.3 Spændinger
rer til det billede, der er givet af Eporøs, der som før nævnt er valgt for lavt. Numeriske og
målte værdier er langt fra hinanden og følger ikke tendensen mellem den homogene og porøse cirkelskive, hvor flytningen vokser. Derfor konkluderes det, at der har været nogle
grove fejl under udførelse af forsøget, og at der ses bort fra dette resultat.
7.3 Spændinger
Formålet med spændingssammenligningen er at finde ud af, hvordan spændingsbilledet for
de analytiske og numeriske beregninger i cirkelskiven er i forhold til hinanden. Disse
spændinger er derefter sammenlignet med de spændinger, der er beregnet ud fra forsøgsemnerne.
Grunden til, at det kun er valgt at sammenligne spændinger i cirkelskiven er, at dette forsøgsemne er det eneste, hvor der vha. rosettegages er bestemt tøjninger eksperimentelt, således spændingerne kan bestemmes. Gagenes placering fremgår af figur 7.2.(tv.)
y
1,4,7
ε yy
(C)
39,5
(B) 4 6
5
39,5
1 3
2
(A)
79
8
39,5
ε xx
α
x
3,6,9
ε xx
2,5,8
Figur 7.2. Placeringen af de tre rosettegages A, B og C med tøjningsmålerne 1-9 (tv.). De tre længdetøjninger, der bliver målt i hver rosettegage (th.).
Først er spændingsbillederne fra de numeriske og analytiske beregninger af spændingerne
optegnet for en enkelkraft på 25 kN. Spændingsbillederne er medtaget for at illustrere
spændingernes variation overalt i cirkelskiven. Spændingsbillederne er plottet i Matlab, og
det er valgt kun at medtage en fjerdedel af skiven, grundet spændingssymmetri omkring
den vandrette og lodrette diagonal.
110
7 Konklusion
Figur 7.3. Viser de analytiske og numeriske spændingsplot, σxx, σyy og σxy forårsaget af en enkelkraft på 25
kN for cirkelskiven. Til venstre ses analytiske spændingsplot og til højre de numeriske. Spændingerne er angivet i [MPa].
Spændingsbillederne viser et ensartet billede af spændingerne i skiven, beregnet med de to
metoder. Det ses på figur 7.3, at de største spændinger samlet set er i det område, hvor enkelkraften påføres og mindst i områderne langs den ydre rand. Derudover kan det iagttages, at σxx er konstant ned langs den lodrette diagonal, og σyy aftager langs diagonalen ind
mod centrum. Det ses endvidere, at σxy er nul langs den lodrette diagonal udover området
ved enkeltkraften. Det konkluderes, at de analytiske og numeriske spændingsbilleder stemmer overens.
For endeligt at kunne konkludere på rigtigheden af den analytiske og numeriske model,
sammenlignes spændinger i tre punkter, figur 7.2. Længdetøjningerne εxx og εyy kan bestemmes direkte ud fra gagene 1, 2, 4, 5, 7 og 8, og derefter omregnes til spændinger.
111
7.3 Spændinger
Tværtøjningen εxy kan ikke bestemmes eksperimentelt, men bestemmes ud fra transformationsformler [Foley, 2004].
Spændingerne σxx, σyy og σxy i A, B og C bestemmes ud fra spændings-tøjningsrelationen
beskrevet ved Hookes udvidede lov. Spændingernes variation med kraften er optegnet på
figur 7.4.
Figur 7.4. Viser kraften som funktion af henholdsvis σxx, σyy og σxy i de tre punkter A B og C, på cirkeskivenl
der ses på figur 7.1. Forsøg (sort), numerisk bestemt (blå) og analytisk bestemt (rød).
Det ses generelt, at spændingerne, udregnet ved de analytiske og numeriske metoder, stemmer godt overens med forsøgsresultaterne. Spændingen σxx har tilnærmelsesmæssigt samme værdi i punkt A som i punkt B, mens σyy er større i punkt A end i punkt B, hvilket også
stemmer overens med den forventede spændingsfordeling.
Det ses endvidere, at punkt C er det punkt, hvor de analytiske og numeriske metoder passer
dårligst med forsøget. Dette gælder især for σxy, men idet σxy er meget lille, set i forhold til
den påførte last, vil denne fejl ikke være afgørende for modellen.
Langs den lodrette rand af cirkelskiven er forskydningsspændingerne teoretisk set lig nul,
hvilket også passer overens med de analytiske beregninger. Ved den lodrette rand afviger
forsøgsresultatet fra nul, hvilket kan grundes måleusikkerheder. Den numeriske beregning
af forskydningsspændingen i punkt A viser sig også at være forskellig fra nul på den lod112
7 Konklusion
rette rand. Dette skyldes, at de beregnede spændinger i Matlab er midlet over det elementareal, der tangerer punkt A. Det samme er gældende i punkt B, men her er elementarealet
så lille, at det ikke giver udslag.
Det kan konkluderes, at beregningsmodellerne giver et godt billede af spændingerne, idet
de ligger meget tæt på outputtet i forsøget. Derudover fremgår det, ligesom ved flytningssammenligningerne, at de analytiske og numeriske spændinger ligger meget tæt op ad hinanden, men afviger fra forsøgsresultaterne, hvilket stemmer overens med de forventede
usikkerheder på de eksperimentelle resultater.
113
114
8
Litteraturliste
[Byskov, 2002]
Byskov, Esben; Elementary continuum mechanics for everyone – and some more, Vol. 2; Aalborg University; 2002;
ISBN: 1395-8232-U0208
[Calfem, 1999]
Department of Mechanics and Materials; CALFEM, a finite
element toolbox to MATLAB version 3.3; Lund University;
Februar 1999; ISSN: 0281-6679
[Cook, et. al., 2002]
Cook, Robert D.; Malkus, David S.; Plasha, Michael E.;
Witt, Robert J.; Concepts and applications of finite element
analysis, 4. edition; John Wiley & Sons. Inc.; 2002; ISBN:
0-471-35605-0
[Cullen, et. al., 2001]
Cullen, Michael R.; Zill, Dennis G.; Differential equations
with boundary-value problems, 5. edition; Brooks/Cole;
2001; ISBN: 0-534-38002-6
[Edwards, et. al., 2002]
Edwards, C. Henry; Penny, David E.; Calculus, 6. edition;
Prentice Hall; 2002; ISBN: 0-13-095006-8
[Fenster, et. al., 2003]
Fenster, Saul K.; Ugural, Ansel C.; Advanced strength and
applied elasticity, 4. edition; Prentice Hall; 2003; ISBN: 013-047392-8
[Foley, 2004]
Foley, Christina; Noter til indledende kontinuummekanik;
Aalborg Universitet; september 2004; ISBN 1395-8232
U0304
[Kohl, 1930]
Kohl, E.; Berechnung von kreisrunden Scheiben unter der
Wirkung von Einzelkräften in ihrer Ebene; Artikel fra Ing.Arch I; 1930
[Myhre, 2005]
Jensen, Henrik Myhre; Professor ved Aalborg Universitet;
Note 2: Simple models of cellular materials and honeycomb
materials, Note 3: Plane conditions / The airy stress function,
Note 5: Stress concentration at circular holes, Note 6: Representative volume elements, unit cells and definitions of mac115
roscopic stresses and strains, Note 8: Voigt and Reuss models, bounds on moduli and compliances, Note 9: Dilute and
self-consistent calculations of elastic moduli of materials
with cavities; 2005
[Myhre 1, 2005]
Jensen, Henrik Myhre; Professor ved Aalborg Universitet;
Forelæsningslektion 2; 2005
[Pilegaard, 2005]
Hansen, Lars Pilegaard; Noter til eksperimentelle metoder;
Aalborg Universitet, Institut for bygningsteknik; 2005
[Teknisk Ståbi, 2003]
C. G. Jensen og K. Olsen; Teknisk Ståbi, 18. udgave, 4. oplag; Ingeniøren bøger; 2003; ISBN: 87-571-2134-6
116