Manual för ett litet FEM-program i Matlab

KTH HÅLLFASTHETSLÄRA
Manual: SMALL FE PROGRAM
Manual för ett litet FEM-program i Matlab
Programmet består av en m-fil med namn SMALL_FE_PROG.m och en hjälp-fil
för att plotta resultat som heter PLOT_DEF.m.
Input
För att köra programmet behöver man anropa det med information om
problemet man vill lösa givet på en specifik form som programmet kan
förstå. Förslagsvis skriver man ett eget Matlab-script (se bifogat exempel
RunFEprog.m) som etablerar input och anropar SMALL_FE_PROG.m. Följande input behövs: Nodes, Elements, ElemType, MtrlProp, CrossSection,
PrescribedU, PointLoads, och VolumeLoads angivna i den ordningen.
Strukturen på de olika input kommer att förklaras i relation till exemplet
som illustreras i Figur 1 nedan.
y
L
1
3!
e2
!
e6
e3
5
e5
e7
e1
2
0.134L
!
!
e4
!
4
x
0.5L
L
1.5L
2L
2.5L
Figur 1: Enkelt exempel med 5 noder och 7 element.
Nodes
Problemet har numNodes noder, då blir

x1
y1
 x2
y2

Nodes =  .
..
 ..
.
xnumN ynumN





(numN×2)
där xi , yi är koordinaterna för nod nummer i. För vårt exempel i Figur 1 så
blir


0.5 1.0
 1.0 0.134 



Nodes = 
 1.5 1.0 
 2.0 0.134 
2.5 1.0
där L = 1 [längdenhet] har antagits.
Carl Dahlberg, Feb 2015
1
KTH HÅLLFASTHETSLÄRA
Manual: SMALL FE PROGRAM
Elements
Alla elementtyper har två noder (en nod i varje ände), förutom Bar3 som
har en tredje nod i mitten. Matrisen Elements har en rad för varje element
där de noder som definierar elementet listas i ordning. För vårt exempel i
Figur 1 så blir


← e1 står på rad 1
1 2
 1 3 
← e2 står på rad 2


 2 3 
osv . . .



2
4
Elements = 


 4 3 


 3 5 
. . . till och med
4 5
← e7 på rad 7
ElemType
Det här är en text-sträng som anger vilken elementtyp som ska användas.
Det finns 4 olika element implementerade: Bar2, Bar3, Beam2 och Frame2.
Element
Bar2
Bar3
Beam2
Frame2
Antal noder
2
3
2
2
Lokala DOF
u
u
w, θ
u, w, θ
Globala DOF
Ux , Uy
Ux , Uy
Ux , Uy , Θz
Ux , Uy , Θz
Teori
Stång
Stång
Balk
Stång & balk
Tabell 1: Elementtyper.
Elementtyp Bar2 är ett linjärt stångelement där den lokala koordinaten
0 ≤ ξ ≤ 1, för övriga element gäller −1 ≤ ξ ≤ 1. Den kvadratiska stången
Bar3 är definierat så att nod 1 och nod 2 har ξ1 = −1 och ξ2 = 1 och nod
3 befinner sig mellan dessa så att ξ3 = 0. Element Bar3 bör brukas med
viss försiktighet då det kan vara lite knepigt att föreskriva tillräckligt många
randvillkor för att undvika att elementet beskriver en mekanisk led.
MtrlProp och CrossSection
Två kolumn-vektorer med längd numElem (antal element). I MtrlProp står Emodulen för materialet i varje element. I CrossSection är stången/balkens
tvärsnitt beskrivet via sidlängden h av den kvadrat som antas beskriva tvärsnittet. Antagandet om kvadratiskt tvärsnitt kan någorlunda enkelt ändras
genom att man gör några mindre ändringar i koden – mer om detta senare.
Carl Dahlberg, Feb 2015
2
KTH HÅLLFASTHETSLÄRA
Manual: SMALL FE PROGRAM
PrescribedU
Betrakta Figur 2 där exempelproblemet har beskrivits med randvillkor och
laster och kommer att användas för att exemplifiera de följande input till
programmet. Låt oss börja med föreskrivna frihetsgrader (DOF).
P
P
!
!
1
5
!
P
3
2
4
!
!
!
!
Figur 2: Exempelproblemet med laster och randvillkor.
Matrisen PrescribedU ska innehålla en rad för varje föreskriven DOF och
har tre kolumner. Den första kolumnen anger vid vilken nod man föreskriver
en DOF, den andra anger i vilken riktning1 och den sista kolumnen anger
vad det föreskrivna värdet är (ofta 0).
Från Figur 2 ser man att vid nod 2 så är både U2x = U2y = 0 och vid
nod 4 tillåts ingen förskjutning i y-led, dvs U4y = 0, vilket ger


← Nod 2, Ux , = 0
2 1 0
← Nod 2, Uy , = 0
PrescribedU =  2 2 0 
4 2 0
← Nod 4, Uy , = 0
PointLoads
Punktlaster som beskrivs enligt samma metod som för PrescribedU, dvs
först listas noden, sedan riktningen2 och sist kraft/moment-värdet. Krafterna
i Figur 2 ger


1 1 −0.707
 1 2
0.707 

PrescribedU = 
 3 2
1.000 
5 1
1.000
där P = 1 [kraftenhet] har antagits.
Riktningen anges med en siffra, 1-3, där 1 betyder Ux , 2 motsvarar Uy och 3 indikerar
rotationsfrihetsgraden, dvs Θz (om den finns med i problemet).
2
Riktningen anges med en siffra, 1-3, där 1 : Fx , 2 : Fy och 3 : Mz (om den finns med
i problemet).
1
Carl Dahlberg, Feb 2015
3
KTH HÅLLFASTHETSLÄRA
Manual: SMALL FE PROGRAM
VolumeLoads
Slutligen finns möjligheten att föreskriva konstanta volymslaster, dvs en
kraft/volymsenhet i varje element. De anges i VolumeLoads som två kolumner, där varje rad representerar [Kx Ky ] för varje element. Om vi t.ex. tänker
oss att elementen i vårt exempelproblem har densiteten ρ och påverkas av
en tyngdkraftsacceleration −g i y-led så kan vi ange detta som


0.0 −1.0
← element 1
 0.0 −1.0 
← element 2


VolumeLoads =  .

.
..
.. 
 ..
.
0.0 −1.0
← element 5
där ρg = 1 [kraftenhet/volym] har antagits.
Output
Programmet returnerar i sin nuvarande version fyra saker: U, F, Sigma och
numNodalDof. De tre första beskriver lösningen där U är frihetsgraderna, F
är nodkrafterna och Sigma är spänningen i varje element. Den sista output
numNodalDof är antingen 2 eller 3 och anger antalet frihetsgrader per nod
och används endast av hjälprutinen PLOT_DEF.m för att visa deformationen.
Man kan givetvis returnera fler saker ur programmet om så skulle behövas, t.ex. om man vill ha tillgång till styvhetsmatrisen.
Då spänningen i ett element inte behöver vara konstant kan det tyckas
något märkligt att programmet endast returnerar ett spänningsvärde per element. Det är designat så för att underlätta och göra programmet så enkelt
som möjligt och man kan givetvis tänka sig mer komplicerade varianter. Därför krävs dock en genomgång om vad spänningsvärdet faktiskt representerar
för de olika elementen.
Sigma för element Bar2
Eftersom detta element är linjärt i sin beskrivning av förskjutningen u så
kommer töjningen ε att vara konstant. Således blir även spänningen σ = Eε
en konstant och det finns inget problem med att återge endast ett värde för
varje element. Spänningen beräknas som σ = EBde .
Sigma för element Bar3
Eftersom detta element är kvadratiskt i sin beskrivning av förskjutningen
u så kan töjningen ε variera linjärt. Spänningen utvärderas för ett antal
punkter, ξi , längs med elementet: σi = EBi de . Det spänningsvärde som har
störst magnitud, dvs som är max |σi | väljs och returneras
sedan med sitt
√
tecken. Utvärderingspunkterna är valda till ξ1,2 = ±1/ 3, men kan givetvis
ändras.
Carl Dahlberg, Feb 2015
4
KTH HÅLLFASTHETSLÄRA
Manual: SMALL FE PROGRAM
Sigma för element Beam2
I balkteori kan töjningen givetvis variera längs med elementet (i lokal xriktning) men det även finns ett implicit linjärt beroende på den lokala zriktningen, dvs genom tjockleken. Utböjningen är en funktion av x, w =
w(x), men enligt kinematik för balkar får man att töjningen även är en
2
funktion av z enligt ε = ε(x, z) = −E ddxw2 z. I programmet har det antagits
att balken har ett symmetriskt tvärsnitt och därför kommer spänningen vara
lika stor på över och undersidan, men med olika tecken. För enkelhetens
skull presenteras det positiva värdet, men man ska alltså vara medveten om
att balken upplever en lika stor tryckspänning på andra sidan. Spänningen
beräknas som σ = E|Bde | h2
Sigma för element Frame2
Här gäller att den totala spänningen är en kombination av spänningen från
normaldeformationen och böjdeformationen, σ(x, z) = σn (x) + σb (x, z). Eftersom σn kommer vara en konstant överlagrad (superponerad) spänning
ovanpå bidraget från σb så kommer inte σ(x, z) att nödvändigtvis vara symmetrisk i z-led. Dvs, det är inte givet att det blir lika mycket tryck på ena
sidan som det blir drag på den andra. Därför beräknas bägge värdena och
den med störst magnitud väljs och returneras med tecken för att kunns skilja
på drag/tryck.
Koden
Koden består av en huvudfunktion, SMALL_FE_PROG, på rad 8 till 112. Merparten av dessa rader är kommentarer för att göra koden mer förståelig.
Några av kommentarerna saknas3 då de ingår som frågor i inlämningsuppgift 3. Huvudfunktionen beskriver flödet genom en FE-lösare och anropar
ett antal hjälpfunktioner som förhoppningsvis har fått beskrivande namn.
Nedan följer några korta kommentarer angående en del av funktionerna
ELEMENT STIFFNESS MATRIX
Funktionen på rad 147-182
R returnerar matrisen ke . Lite slarvigt kan man
säga att operationen ke = BT CBkJkdV redan är genomförd analytiskt för
varje elementtyp i förväg (detta är möjligt eftersom de inte är isoparametriska
element).
VOLUME FORCE
Funktion på rad 186-202 där volymslaster hanteras. Information från transformationsmatrisen T används för att dela upp vektorn K = [Kx Ky ]T i
3
Jag kommer lägga upp en version med alla kommentarer efter deadline för HU3.
Carl Dahlberg, Feb 2015
5
KTH HÅLLFASTHETSLÄRA
Manual: SMALL FE PROGRAM
komponenter som verkar längs och tvärs elementet. Detta därför att stänger bara kan hantera krafter längs sin utsträckning och balkar endast tvärs.
Elementet FRAME2 kan hantera bägge.
AREA och AREA MOMENT OF INERTIA
Två väldigt korta funktioner på raderna 236-247 som returnerar tvärsnittsegenskaper. Som det är kodat nu så antas ett kvadratiskt tvärsnitt, men detta
kan enkelt ändras till något annat om man vill. För spänningsutvärdering
har det antagits att tvärsnittet är symmetriskt så om man vill koda in t.ex.
ett trekantigt tvärsnitt så bör man också “ta höjd” för detta i funktionen
EVALUATE_STRESS vid elementen BEAM2 och FRAME2.
SHAPEFNC AND BMTRX
Funktion på rad 329-375 som returnerar N och B för ett givet värde på
den lokala koordinaten ξ. Eftersom de element som är implemeterade inte
behöver B för att etablera ke så används denna funktion endast vid beräkning av spänning (egentligen töjning), enligt ε = Bde i en stång och
ε = −Bde z i en balk, och sedan givetvis σ = Eε. Se kommentarer i funktionen EVALUATE_STRESS för detaljer. Så som programmet är skrivet just nu så
används inte N någonstans.
CHECK INPUT
Den sista funktionen i filen (men som faktiskt anropas först av alla). Den
returnerar antal frihetsgrader per nod och antalet noder per element. Dessutom gör den en “sanity check” på input som bör fånga upp och varna för
några av de vanligaste misstagen man kan tänkas göra i input. Det finns
dock fortfarande oändligt med möjligheter kvar att göra fel på. . .
Carl Dahlberg, Feb 2015
6