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
© Copyright 2024