Finite Element Formulation for One-Dimensional Problems

large strain finite element formulation for one dimensional membrane elements finite element formulation of two dimensional dynamic analysis axisymmetric
GregDeamons Profile Pic
GregDeamons,New Zealand,Professional
Published Date:03-08-2017
Your Website URL(Optional)
Comment
5 Finite Element Formulation for One-Dimensional Problems Wehavenowpreparedalloftheingredientsneededforformulatingthefiniteelementequations:(1)the weakform,whichisequivalenttothestrongformwewishtosolve,and(2)thefiniteelementweightand trialfunctions,whichwillbepluggedintotheweakform.Sowearereadytodevelopthefiniteelement equationsforthephysicalsystemswehavedescribedinChapter3:heatconduction,stressanalysisandthe advection–diffusionequation.ThisisthelaststepintheroadmapinFigure3.1.Thisstepisoftencalledthe discretization,aswenowobtainafinitenumberofdiscreteequationsfromtheweakform. The procedure is similar to the one we used in Example 3.3. We first construct admissible weight functionsand trialsolutionsintermsofarbitraryparameters.However,inthefiniteelementmethod,the parametersarethenodalvaluesofthefunctions.Fromthearbitrarinessofthenodalvaluesfortheweight function,wethendeducethefiniteelementequations,whicharelinearalgebraicequations.Weoftencall thesethediscreteequations thesystemequations;instressanalysis,theyarecalledthestiffnessequations. Thefiniteelementanalysisprocedureisoftenbrokenupintofoursteps: 1. preprocessing,inwhichthemeshisconstructed; 2. formulationofthediscretefiniteelementequations; 3. solvingthediscreteequations; 4. postprocessing,wherethesolutionisdisplayedandvariousvariablesthatdonotemanatedirectlyfrom thesolutionarecalculated. Inonedimension,preprocessingandpostprocessingarequitestraightforward,sowewillhavelittletosay about these in this chapter. However, in multidimensional problems, these are quite challenging and importantstepsforusersofsoftware. 5.1 DEVELOPMENT OF DISCRETE EQUATION: SIMPLE CASE Inordertominimizetheabstractnessofthisdescription,wefirstconsiderthespecificproblemdiscussedin Section3.2,withafiniteelementmodelconsistingoftwolinearelementsasshowninFigure5.1a.Ascanbe seen,atx¼ 0,theproblemhasatraction(natural)boundarycondition,andanessentialboundarycondition isappliedatx¼ l.NodesontheessentialboundaryarenumberedfirstasshowninFigure5.1a. TheweakformhasbeendevelopedinChapter3andisgivenasfollows. A First Course in Finite Elements J. Fish and T. Belytschko 2007 John Wiley & Sons, Ltd ISBNs: 0 470 85275 5 (cased) 0 470 85276 3 (Pbk)94 FINITE ELEMENT FORMULATION FOR ONE-DIMENSIONAL PROBLEMS Γ t Γ u (a) 3 2 1 x x = 0 3 x = l (1) (2) 1 N (x) 2 N (x) 1 N (x) 3 (b) 1 2 3 x (1) (2) u u 3 2 (c) 2 u 1 3 1 x (1) (2) Figure 5.1 (a)Two-elementmesh,(b)globalshapefunctionsand(c)anexampleofatrialsolutionthatsatisfiesan essentialboundarycondition. FinduðxÞamongthesmoothtrialsolutionsthatsatisfytheessentialboundaryconditionuðlÞ¼ u such 1 that  Z  Z T l l  dw du T T   AE dx w bdxðw tAÞ ¼ 0 8wðxÞ with wðlÞ¼ 0: ð5:1Þ  dx dx 0 0 x¼0 Intheabove,wehavetakenthetransposeoftheweightfunctions;aswðxÞisascalar,thisdoesnotchangethe valueoftheexpression,butitisnecessaryforconsistencywhenwesubstitutematrixexpressionsforwðxÞ oritsderivative. TheprocedurewewillfollowissimilartoExample3.3:Wewillevaluatetheweakformforthefinite elementtrialsolutionsandweightfunctions.Then,byinvokingthearbitrarinessoftheweightfunctions,we willdeduceasetoflinearalgebraic(discrete)equations. Thefiniteelementweightfunctionsare h wðxÞ w ðxÞ¼ NðxÞw; ð5:2Þ where  denotes an approximation and NðxÞ is the matrix of shape functions. For this mesh, wðxÞ¼ w N ðxÞþw N ðxÞþw N ðxÞ.Thefiniteelementtrialsolutionsareapproximatedbythesame 1 1 2 2 3 3 shapefunctions: h uðxÞ u ðxÞ¼ NðxÞd: ð5:3Þ Forthismesh,uðxÞ¼ u N ðxÞþu N ðxÞþu N ðxÞ.Noticethatwerefertotheweightfunctionsandtrial 1 1 2 2 3 3 solutions inthepluralcaseasthereareinfinitelymany;ourtaskwillbetofindthattrialsolutionwhich satisfiestheweakform.VariousshapefunctionsweredevelopedinChapter4,andtheprocedurewewill developwillbeapplicabletoallofthem,butfirstwewillfocusonthetwo-nodeelementwithlinearshape functions.Thesefiniteelementshapefunctions,aswelearnedinChapter4,aresufficientlysmoothtobe employedintheweakform.DEVELOPMENT OF DISCRETE EQUATION: SIMPLE CASE 95 Thetrialsolutionsmustbeconstructedsothattheysatisfytheessentialboundarycondition.Thiscanbe easilyaccomplishedbyletting u ¼u : ð5:4Þ 1 1 Theothernodaldisplacementsareunknownandwillbedeterminedbythesolutionoftheweakform.The globalshapefunctionsareshowninFigure5.1(b).Noticethattheyarethetentfunctionswehavedescribed in Chapter 4. The finite element approximation is a linear combination of these shape functions. An exampleofafiniteelementtrialsolutionisshowninFigure5.1(c).Becauseof(5.4)andthesmoothnessof thefiniteelementapproximation,allofthetrialsolutionsareadmissible. Ontheessentialboundary,theweightfunctionsmustvanish.Tomeetthisrequirement,weset w ¼ 0: ð5:5Þ 1 Theremainingnodalvalues,w andw ,arearbitrary,astheweightfunctionsmustbearbitrary. 2 3 TheelementandglobalmatricesarerelatedbygathermatricesjustasinChapter2,sowehave e e e e w ¼ L w; d ¼ L d: ð5:6Þ Thegathermatricesfollowfromtherelationbetweenlocalandglobalnodenumbers. As the finite element functions and their derivatives have kinks and jumps at the element interfaces, respectively (see Figure 3.5), efficient integration of theweak form (5.1) necessitates evaluation of the e e integralover½0;lasasumofintegralsoverindividualelementdomains½x ;x .Sowereplacetheintegral 1 2 overtheentiredomainin(5.1)bythesumoftheintegralsovertheelementdomains: ()    Z e Z e n T el x x X e e 2 2  dw du e e eT eT e   A E dx w bdxðw A tÞ ¼ 0 ð5:7Þ  e dx dx e x x x¼0 e¼1 1 1 wherewehaveplacedsuperscript‘e’ontheweightandtrialfunctionstoindicatethatthesearethepartsof thosefunctionsthatpertaintoelemente.Ineachelemente,theweightfunction(5.2)andtrialsolution(5.3) canbewrittenas e du e e e e e u ðxÞ¼ N d ; ¼ B d ; dx  ð5:8Þ T e dw eT eT eT eT eT w ¼ w N ; ¼ w B ; dx e e where d and w are given in terms of the global nodal values by (5.6). Equation (5.8) is the same approximationas (5.2)and (5.3),and thesefunctions arealsoadmissible.Theyarealocalizationof the globalapproximationstotheelements;theyfollowfromthefactthatinelemente,theglobalNandelement e shapefunctionsN areidentical(seeFigure4.8).Henceforthinthisbook,wewillwritethefiniteelement approximationsattheelementlevelintheform(5.8);theessentialboundaryconditionswillbemetonthe e e globallevelanditwillbeimplicitthatd andw aregivenintermsoftheglobalnodalvaluesby(5.6). Substituting(5.8)into(5.7)gives 8 9 e e x x 2 2 Z Z n = el X eT eT e e e e eT eT e  w B A E B dxd  N bdxðN A tÞ ¼ 0: ð5:9Þ x¼0 fflfflfflzfflfflffl e¼1 e e x x f e  1 1 fflfflfflfflfflfflfflfflfflfflfflfflfflzfflfflfflfflfflfflfflfflfflfflfflfflffl fflfflfflfflfflfflfflzfflfflfflfflfflfflffl : ; e K f e 96 FINITE ELEMENT FORMULATION FOR ONE-DIMENSIONAL PROBLEMS Intheabove,wehavedefinedtwomatricesthatwillbeveryusefulinthefiniteelementmethod(FEM): (i) theelementstiffnessmatrix e x 2 Z Z e eT e e e eT e e e K ¼ B A E B dx¼ B A E B dx; ð5:10Þ e  e x 1 (ii)theelementexternalforcematrix e x 2 Z Z   e eT eT e eT eT e    f ¼ N bdxþðN A tÞ ¼ N bdxþðN A tÞ ð5:11Þ x¼0  e e   t e fflfflfflfflfflfflfflzfflfflfflfflfflfflffl x fflfflfflfflfflfflfflzfflfflfflfflfflfflffl 1 e e f  f  e e e where istheportionoftheelementboundaryonthenaturalboundaryandf andf in(5.11)arethe t   elementexternalbodyandboundaryforcematrices,respectively.Theelementmatriceswillplaythesame key roles as in the analysis of discrete systems in Chapter 2: They are the building blocks of theglobal equations.Wewillexaminethesematricesforstressanalysisandheatconductioninmoredetaillater.In (5.10) and (5.11), the far right-hand side expressions use a notation that we will introduce in the next section. Substituting(5.10)and(5.11)into(5.9)andusing(5.6)gives n n el el X X T eT e e eT e w L K L d L f ¼ 0: ð5:12Þ e¼1 e¼1 InderivingEquation(5.12)recallthatwisnotafunctionofxandisaglobalmatrix,andhenceitcanbetaken e outsideofthesummationsymbol.Moreover,thescatteroperatorL isnotafunctionofx,butiselement dependent.Therefore,ithasbeentakenoutoftheintegral,butshouldremaininsidethesummationoverthe elements. If you compare the first sum in (5.12) to Equation (2.25), the expression can be recognized as the assembledsystem(stiffness)matrix n el X eT e e K¼ L K L : ð5:13Þ e¼1 The system matrix for the differential equation is assembled by exactly the same operations as for the discretesystems:matrixscatterandadd,whichisalsoequivalenttodirectassembly.Itshouldbestressed that we do not need to perform the large matrix multiplications indicated above to assemble the global matrices.TheassemblyprocessesareidenticaltotheassemblyprocedureswehavelearnedinChapter2. Thesecondtermin(5.12)istheassembledexternalforcematrix n el X eT e f ¼ L f : ð5:14Þ e¼1 Thisisthecolumnmatrixassemblyoperation.Itconsistsofacolumnmatrixscatterandaddandisactually easiertolearnthanmatrixassembly;itwillbeillustratedintheexamplesthatfollow. SubstitutingEquations(5.13)and(5.14)intoEquation(5.12)yields T w ðKdfÞ¼ 0 8w except w ¼ wðlÞ¼ 0; ð5:15Þ 1ELEMENT MATRICES FOR TWO-NODE ELEMENT 97 wherewehaveindicatedthearbitrarinessofthenodalvalues,w,whichemanatesfromthearbitrarinessof theweightfunctionsinthestatementoftheweakform(5.1)andtherestrictiononw,(5.5).Let r¼ Kdf; ð5:16Þ whereriscalledtheresidual.Then(5.15)becomes T w r¼ 0 8w except w ¼ 0: ð5:17Þ 1 IfwewriteEquation(5.15)forthespecificmodelinFigure5.1,wehave w r þw r ¼ 0; 2 2 3 3 wherethefirsttermhasdroppedoutbecausew ¼ 0.Astheaboveholdsforarbitraryw andw ,wecan 1 2 3 deducethatr ¼ r ¼ 0,butwecannotsayanythingaboutr ,andinfact,asitistheunbalancedforceat 2 3 1 node1,soitisthereactionforce.Ifwewritetheequations,weobtain 2 3 2 32 3 2 3 r K K K u f 1 11 12 13 1 1 4 5 4 54 5 4 5 r¼ 0 ¼ K K K u  f : ð5:18Þ 21 22 23 2 2 0 K K K u f 31 32 33 3 3 Rearrangingthetermin(5.18)gives 2 32 3 2 3 K K K u f þr 11 12 13 1 1 1 4 54 5 4 5 K K K u ¼ f : ð5:19Þ 21 22 23 2 2 K K K u f 31 32 33 3 3 Equation(5.19)isasystemofthreeequationswiththreeunknowns,u ,u andr .ItissimilartoEquation 2 3 1 (2.27)derivedinChapter2.Varioussolutionproceduressuchaspartitionandpenaltymethodshavebeen discussedinChapter2.Forinstance,usingthepartitionapproach,thenodaldisplacementsu andu are 2 3 foundfirstbysolving  K K u f K u 22 23 2 2 21 1 ¼ ; K K u f K u 32 33 3 3 31 1 followedbythecalculationofunknownreactionsatnode1: 2 3 u 1 4 5 r ¼ f ½K K K  u : 1 1 11 12 13 2 u 3 Liketheequationsfordiscretesystems,Equation(5.19)canbeviewedasequationsofdiscreteequilibrium atthenodes.Theleft-handsideisthematrixofinternalforcesandtheright-handsideisthatoftheexternal forcesandreactions.Notethatthestiffnessmatrixin(5.19)isstillsingular.However,thepartitionapproach doesnotrequireitsinversion. 5.2 ELEMENT MATRICES FOR TWO-NODE ELEMENT e e Consider a two-node linear element with constant cross-sectional area A and Young’s modulus E subjectedtolineardistributionofbodyforcesasshowninFigure5.2.Inthissection,wederivetheelement stiffnessmatrixandtheexternalforcematrix.98 FINITE ELEMENT FORMULATION FOR ONE-DIMENSIONAL PROBLEMS b 2 b 1 e e E A 1 2 x e e x x 1 2 Figure 5.2 Two-nodeelementwithlineardistributionofbodyforce. RecallthatinSection4.1weshowedthatthetwo-nodeelementshapefunctionsandtheirderivativesare givenas   e e xx xx 1 2 1 e e e N ¼ ¼ ðx xÞðxx Þ ; e e e e 2 1 e x x x x l 1 2 2 1 ð5:20Þ  d 1 1 1 e e B ¼ N¼ ¼ 11 : e e e dx l l l Theelementstiffnessmatrixisthen e e e x x " " x 2 2 2 Z Z Z e e 1 1 1 1 A E e eT e e e e e K ¼ B A E B dx¼ A E ½11 dx¼ ½11 dx e e 2 e l l 1 ðl Þ 1 fflfflfflfflfflzfflfflfflfflffl e e e x x x fflfflfflfflfflzfflfflfflfflffl 1 1 1 e B eT B 0 1 " e e 1 1 A E B C e e ¼ x x ; A 2 2 1 e fflfflfflzfflfflffl ðl Þ 11 e l  e e A E 1 1 e K ¼ : ð5:21Þ e 11 l NotethatthisresultisidenticaltothatforthebarelementderivedinChapter2basedonphysicalarguments. Inotherwords,thestiffnessmatrixofthetwo-nodeelementwithconstantcross-sectionalareaandconstant Young’smoduluswhenderivedfromtheweakformisidenticaltothatobtainedbyphysicalarguments.It then may occur to you, why go to all this trouble? The reason is that for higher order elements and in multidimensions,theproceduresdescribedinChapter2donotwork,whereastheweakformcanbeapplied tohigherorderelementsandtwoandthreedimensions. Wenowturntotheevaluationoftheexternalnodalbodyforces,thefirstterminEquation(5.11): e x 2 Z e eT f ¼ N bðxÞdx:  e x 1 Asthebodyforcedistributionislinear,itcanbeexpressedintermsoflinearshapefunctionsas  b 1 e bðxÞ¼ N b; b¼ : b 2APPLICATION TO HEAT CONDUCTION AND DIFFUSION PROBLEMS 99 Table 5.1 Terminology for finite element matrices. Matrices Elasticity Diffusion Heatconduction K Stiffness Diffusivity Conductance f Force Flux Flux d Displacement Concentration Temperature Theelementbodyforcematrixisthengivenas e e x x" 2 2 Z Z 2 e e e 1 ðx xÞ ðx xÞðxx Þ e eT e 2 2 1 f ¼ N N dxb¼ dxb  2 2 e e e e ðl Þ ðx xÞðxxÞðxx Þ e e 2 1 1 x x 1 1  e 21 b l 1 ¼ : 6 12 b 2 e Itcanbeseenthatthesumofforcesactingontheelementisl ðb þb Þ=2,whichisexactlytheintegralof 1 2 thebodyforceovertheelementdomain,i.e.thetotalforce.Asexpected,forb ¼ b ,halfoftheforcegoes 1 2 tonode1andhalftonode2. 1 5.3 APPLICATIONTOHEATCONDUCTIONANDDIFFUSIONPROBLEMS Theexpressionsforheatconductionandotherdiffusionequationscanbeobtainedbyjustreplacingthe fieldsandparametersusingtheconversiontableintroducedinChapter3(Table3.2).Theterminologyofthe matricesinthediscreteheatconductionanddiffusionequationsissummarizedinTable5.1.Theelement matricesaregivenby Z e eT e e e K ¼ B A  B dx; e  Z   e eT eT e   ð5:22Þ f ¼ N f dxþðN A Þ  e e   fflfflfflfflfflfflfflzfflfflfflfflfflfflffl  fflfflfflfflfflfflfflfflfflzfflfflfflfflfflfflfflfflffl e e f  f  withparametersdefinedbyusingtheequivalencesgiveninTable3.2. Example5.1. Heatconduction Wewillfirstuseaheatconductionproblemtoillustratehowthefiniteelementprocedureisapplied.This example will illustrate the construction and solution of the finite element equations and discuss the accuracyoffiniteelementsolutions.Mostoftheproceduresand discussion inthisexampleapplyequally tostressanalysis. 1 Considerabarwithauniformlydistributedheatsourceofs¼5Wm .Thebarhasauniformcross- 1 2  1 sectionalareaofA¼ 0:1m andthermalconductivityk¼2W C m .Thelengthofthebaris4m.  2 TheboundaryconditionsareTð0Þ¼ 0 Candqðx¼ 4Þ¼5Wm asshowninFigure5.3.Dividethe problemdomainintotwolineartemperaturetwo-nodeelementsandsolveitbytheFEM. Preprocessing Westartbynumberingthenodeson .ThefiniteelementmeshisshowninFigure5.4. T 1 Recommended for Science and Engineering Track.100 FINITE ELEMENT FORMULATION FOR ONE-DIMENSIONAL PROBLEMS −2 q (x = 4)n = q = 5 Wm T (x = 0) = T = 0°C x −1 x = 0 s = 5 W m x = 4 m Figure 5.3 ProblemdefinitionofExample5.1. Elementconductancematrix Thetwo-nodeelementshapefunctions,theirderivativesandtheresultingconductancematrix(replace e e E byk in(5.21))  Z e e A k 1 1 e eT e e e K ¼ B A k B dx¼ e 11 l e  werederivedinSection5.2.Notethatthisresultissimilartothebarelement,exceptthatYoung’smodulus isreplacedbyconductivity. Forelement1,wehave ð1Þ ð1Þ ð1Þ ð1Þ x ¼ 0; x ¼ 2; l ¼ 2; ðAkÞ ¼ 0:2; 1 2 "" 1 1 0:1 0:1 0:2 ð1Þ K ¼ ¼ ; 2 11 0:10:1 andsimilarlyforelement2:  0:1 0:1 ð2Þ K ¼ : 0:10:1 Conductancematrix The(global)conductancematrixisobtainedbythematrixassemblyoperation: n el X eT e e ð1ÞT ð1Þ ð1Þ ð2ÞT ð2Þ ð2Þ K¼ L K L ¼ L K L þL K L : ð5:23Þ e¼1 Wecanusedirectmatrixassemblytoobtainit,buttoshowthatthetwoproceduresareidenticalwewill first obtain the global conductance matrix by the above equation. We will assemble the entire con- ductancematrixwithouttakingintoaccounttheessentialboundaryconditions.Thismeansthatjustasin Chapter 2, we will obtain equations for which the right-hand side contains unknowns. However, by assembling all of the equations, we will be able to evaluate the boundary flux matrix at the essential boundaries. Γ T (1) (2) 2 1 3 x x = 0 x = 2 x = 4 1 2 3 Figure 5.4 FiniteelementmeshofExample5.1.APPLICATION TO HEAT CONDUCTION AND DIFFUSION PROBLEMS 101 Thegatheroperatorsforthetwoelementsare 2 3 2 3 "" T 1 ð1Þ T 10 0 T 1 1 6 7 ð1Þ ð1Þ 4 5 d ¼ ¼ ¼ T ¼ L d; 4 5 2 ð1Þ T 01 0 2 T 2 T 3 2 3 2 3 "" T 1 ð2Þ T T 01 0 2 1 6 7 ð2Þ ð2Þ 4 5 d ¼ ¼ ¼ 4T 5¼ L d: 2 ð2Þ T 00 1 3 T 2 T 3 Thescatteroftheconductancematricesgives 2 3 2 3  10 0:1 0:10 ð1Þ 0:1 0:1 10 0 ð1ÞT ð1Þ ð1Þ 4 5 4 5 K ¼ L K L ¼ 01 ¼ 0:10:10 0:10:1 01 0 00 000 2 3 2 3  00 000 ð2Þ 0:1 0:1 010 ð2ÞT ð2Þ ð2Þ 4 5 4 5 K ¼ L K L ¼ 10 ¼00:1 0:1 0:10:1 001 01 0 0:10:1 Thetotalstiffnessisobtainedbyaddingthescatteredelementstiffnessesgivenabove 2 3 0:1 0:10 ð1Þ ð2Þ 4 5 K¼ K þK ¼ 0:10:2 0:1 : ð5:24Þ 0 0:10:1 In practice, the above triple products are not performed, but rather a direct assembly, as previously describedinChapter2,isemployed.Thedirectassemblyfortheprocessisshownbelow   0:1 0:1 ½1 0:1 0:1 ½2 ð1Þ ð2Þ K ¼ K ¼ 0:10:1 ½2 0:10:1 ½3 ½1½2½2½3 Theresultingglobalconductancematrixis 2 3 0:1 0:10 ½1 4 5 K¼ 0:10:2 0:1 ½2: 0 0:10:1 ½3 ½1½2½3 Thismatrix,obtainedbydirectassembly,isidenticalto(5.24) Boundaryfluxmatrix  The element boundary flux matrices are calculated by (5.11) where t has been replaced by qaccordingtoTable3.2   e eT e eT eT  f ¼ðN A qÞ ¼N ðxÞ0:15¼0:5N ðx Þ: 3 3   e  q Note thatthe shapefunctionsfor element1(shownin Figure5.5)vanishon .Onlyshapefunctions q that are nonzero on the natural boundary will contribute to the nodal boundary flux. Therefore, in q computing the boundary flux matrix we need to consider only those elements that are on the natural boundary.102 FINITE ELEMENT FORMULATION FOR ONE-DIMENSIONAL PROBLEMS (1) (1) N N 1 2 1 q = 5 at x 3 A = 0.1 = constant 3 2 1 x 0 x x x (1) (2) 1 2 3 0 Figure 5.5 Shapefunctionsforelement1. Usingtheaboveequation,theelementboundaryfluxmatricesforthetwoelementsare "   ð1Þ N ðx Þ 0 0 ½1 3 ð1Þ 1 f ¼0:5 ¼0:5 ¼  ð1Þ 0 0 ½2 N ðx Þ 3 2 "  ð2Þ N ðx Þ 0 0 ½2 3 ð2Þ 1 f ¼0:5 ¼0:5 ¼  ð2Þ 1 0:5 ½3 N ðx Þ 3 2 Thescatter(ordirectassemblyprocess)thengivestheglobalboundaryfluxmatrix: 2 X eT e f ¼ L f ;   e¼1 2 3 2 3 2 3 10 00 0 ½1   0 0 6 7 6 7 6 7 f ¼4015 þ4105 ¼4 0 5½2:  0 0:5 00 01 0:5 ½3 NotethatthisresultisthesameasassigningðAqÞj tothenodewherethefluxisprescribedandzeroat  q allothernodes.Inthisway,theboundarymatrixcanbecomputeddirectly. Sourcefluxmatrix TheelementsourcefluxmatrixisderivedinSection5.2andisgivenas e x n en Z  e l 21 s 1 e eT f ¼ N sdx¼ :  6 12 s 2 e x 1 wherebin(5.11)hasbeenreplacedbysaccordingtoTable3.2.Sinces ¼ s ¼ s,theabovereducesto 1 2  e l s 1 e f ¼ :  2 1 Itcanbeseenthathalfoftheheatgoestonode1andhalftonode2.Thisalsofollowsfromthefactthatthe integraloflinearshapefunctionsovertheelementcanbecomputedasanareaofatrianglewithheight equalto1andthebaseequaltotheelementlength;thisfollowseasilyfromFigures5.5and5.6. ð1Þ ð2Þ Inthepresentexample,l ¼ l ¼ 2ands¼ 5,whichgives  5 ð1Þ ð2Þ f ¼ f ¼ :   5 Theelementsourcefluxmatrixisthenassembled: 2 3 2 3 2 3   2 10 00 5 X 5 5 eT e 4 5 4 5 4 5 f ¼ L f ¼ 01 þ 10 ¼ 10 :   5 5 e¼1 00 01 5APPLICATION TO HEAT CONDUCTION AND DIFFUSION PROBLEMS 103 (2) (2) N N 1 2 1 q = 5 at x 3 A = 0.1 = constant 1 2 3 x 0 (1) (2) x x x 1 2 3 0 Figure 5.6 Shapefunctionsforelement2. Inpractice,adirectassemblyisusedinstead:  5 ½1 2 3 ð1Þ f ¼ 5 ½1  5 ½2 4 5  ) f ¼ 5þ5 ½2:  ð2Þ 5 ½2 5 ½3 f ¼  5 ½3 Partitionandsolution Theglobalsystemofequationsisgivenby 2 32 3 2 3 2 3 2 3 2 3 0:1 0:10 0 5 0 r r þ5 1 1 4 54 5 4 5 4 5 4 5 4 5 0:10:2 0:1 T ¼ 10 þ 0 þ 0 ¼ 10 : 2 0 0:10:1 T 5 0:5 0 4:5 3 Sincenode1isontheessentialboundary,wepartitionafterthefirstrow,whichgives     0:2 0:1 T 10 T 145 2 2 ¼ ) ¼ : 0:10:1 T 4:5 T 190 3 3 Postprocessing Thetemperaturegradientisgivenas 2 3 0  ð1Þ dT 1 100 6 7 ð1Þ ð1Þ ¼ B L d¼ ½11 41455¼ 72:5; dx 2 010 190 2 3  0 ð2Þ 010 dT 1 6 7 ð2Þ ð2Þ ¼ B L d¼ ½11 145 ¼ 22:5: 4 5 dx 2 001 190 1 Notethatthetemperaturegradientispiecewiseconstantand,aswillbeseeninplottingit,aC function. Evaluationofsolutionquality The finite element solution will now be compared to the exact analytical solution. This type of comparison can be done only for some simple problems (primarily in one dimension) for which the exactsolutionisknown. WestartfromthestrongformfromChapter3:  d dT Ak þs¼ 0; 0 x l; dx dx  2 d dT d T 0:2 þ5¼ 0 ) ¼25; 2 dx dx dx dT dT 5 Tð0Þ¼ 0; qð4Þ¼k nj ¼ 5 ) ð4Þ¼ ¼2:5: x¼4 dx dx 2104 FINITE ELEMENT FORMULATION FOR ONE-DIMENSIONAL PROBLEMS Integratingthegoverningdifferentialequationgives 2 d T dT dT ¼25) ¼25xþc ) ð4Þ¼2:5¼254þc ) c ¼ 97:5: 1 1 1 2 dx dx dx Theexpressionforthetemperatureisobtainedbyintegratingthetemperaturegradient,whichgives dT 2 ¼25xþ97:5) T¼12:5x þ97:5xþc ; 2 dx 2 Tð0Þ¼ 0)12:5ð0Þ þ97:5ð0Þþc ¼ 0) c ¼ 0: 2 2 Thus,theexacttemperatureandtemperaturegradientare ex dT ex 2 T ¼12:5x þ97:5x; ¼25xþ97:5: dx Figure5.7comparestheFEMsolutionwiththeexactsolution.Itcanbeseenthatthenodaltemperatures fortheFEMsolutionareexact.Thisisanunusualanomalyoffiniteelementsolutionsinonedimension anddoesnotoccurinmultidimensionalsolutions.ItisexplainedinHughes(1987) p.25.Notethatthe essential boundary condition is satisfied exactly. This is not surprising as the trial solution was constructed so as to satisfy the essential boundary condition. In finite element solutions, essential boundaryconditionswillalwaysbesatisfiedexactly. Figure 5.8 compares the derivative of the finite element solution with the exact derivative (the derivative is proportional to the flux). As can be seen from Figure 5.8 and as mentioned before, the 1 derivative is aC function; the derivative of the temperature and hence the flux in the finiteelement solutionisdiscontinuousbetweenelements.AspointedoutinFigure5.8,thenaturalboundarycondition at x¼ 4 is not satisfied by the finiteelement solution. However, wewill see in other examples and in exercises that the natural boundary condition is met more accurately as the mesh is refined. Thus, althoughwedonothavetoconstructthefiniteelementapproximationstosatisfythenaturalboundary conditions,theyaremetapproximately. Itisalsoinformativetoseehowwelltheheatconductionequationismetbythefiniteelementsolution. Recalltheheatconductionequation(3.12)andsubstitutethefiniteelementsolutionforthetemperature: 2 d Ak ðNðxÞdÞþsðxÞ¼ errðxÞ: ð5:25Þ 2 dx ex 2 Tx =– 12.5 + 97.5x T h T x Figure 5.7 Comparisonoftheexactandfiniteelementsolutionsoftemperature.DEVELOPMENT OF DISCRETE EQUATIONS FOR ARBITRARY BOUNDARY CONDITIONS 105 Figure 5.8 Comparisonoftheexactandfiniteelementsolutionsoftemperaturegradient. In the above, we have replaced the zero on the RHS of the heat conduction equation by ‘err’ as the deviation from zero is indicative of the error in the finite element solution. The first term of Equation (5.25)willvanishinsideanelement,astheshapefunctionsarelinearinx.Therefore,insidetheelements, theerrorintheheatconductionequationwillbe errðxÞ¼ sðxÞ for x6¼ x : I Thiserroractuallyappearstobequitelargeandfurthermorewouldnotdecreasewithrefinementofthe mesh.Thebehavioratthenodesismorecomplicatedandwillnotbeconsideredhere. Thus,boththenaturalboundaryconditionandthebalanceequationsaremetapproximatelyonlyby the finiteelement solution.However, it can be shown that the finiteelement solution convergesto the exact solution as the mesh is refined, although this is not readily apparent from the weak form. ConvergenceofthefiniteelementsolutiontotheexactsolutionisdiscussedinSection5.6. 5.4 DEVELOPMENT OF DISCRETE EQUATIONS FOR ARBITRARY BOUNDARY CONDITIONS Wewill now consider the development of the finiteelement equations for theweak form with arbitrary boundaryconditions,Equation(3.49).Forconvenience,wewriteitagain: finduðxÞ2 U suchthat   Z Z T  dw du T T   AE dx w bdxðw AtÞ ¼ 0 8w2 U : ð5:26Þ 0  dx dx  t   ConsiderthefiniteelementmeshshowninFigure5.9.Theelementscanbeofanysize,andaswewillsee later, smaller elements are usually used where they are needed for accuracy. The nodes onthe essential boundaryarenumberedfirstas wewillusethepartitioningmethoddescribedinChapter2.Theactualdata106 FINITE ELEMENT FORMULATION FOR ONE-DIMENSIONAL PROBLEMS (1) (2) … n e … el x 2 n 1 …… I np x = ax = b Figure 5.9 Finiteelementmeshinonedimension. neednotbeofthatform,asthenodescanberenumberedintheprogram;mostcommercialsoftwaredonot usepartitioning.Butforthepurposeofthefollowingdevelopment,itisassumedthattheessentialboundary nodesappearfirstinallmatrices. Having selected the finite element mesh and constructed smooth approximation functions over individual element domains (5.8), we now express the integral over  in (5.26) as a sum of integrals overelementdomains: 8 9  Z Z n T = el e e X  dw du e e eT eT e   A E dx w bdxðw A tÞ ¼ 0 8w2 U ; ð5:27Þ 0  : e; dx dx  e¼1 t e e   e e where  are the element domains; integration over  is equivalent to integration over the interval e e ½x ;x . 1 nen Wewillusethesameglobalapproximationsfortheweightfunctionsandtrialsolutions,(5.2)and(5.3), respectively.Todealwitharbitraryboundaryconditions,wewillpartitiontheglobalsolutionandweight functionmatricesas     d w 0 E E d¼ ; w¼ ¼ : w w d F F F Thepartofthematrixdenotedbythesubscript‘E’containsthenodalvaluesontheessentialboundaries.As  indicated by the overbar on d , these values of the solution are set to satisfy the essential boundary E conditions,sotheycanbeconsideredasknown.Thesubmatricesdenotedbythesubscript‘F’containallthe remaining nodal values: these entries are arbitrary for the weight function and unknown for the trial solution.Theresultingweightfunctionsandtrialsolutionswillthereforebeadmissible. Substituting(5.8)into(5.27)gives 8 9  Z Z nel = X  eT eT e e e e eT eT e   w B E A B dx d  N bdxðN A tÞ ¼ 0 8w : ð5:28Þ F  : e;  e¼1 t e e   Notethat(5.28)isforarbitraryw asw isnotarbitrarybutinsteadmustvanish. F E e e e e Substituting(5.10)and(5.11)into(5.28)andusing(5.6),w ¼ L wandd ¼ L d,gives 2 3 6 7 n n el el X X 6 7 T eT e e eT e 6 7 w L K L d L f ¼ 0 8w : ð5:29Þ F 6 7 4 5 e¼1 e¼1 fflfflfflfflfflfflfflfflfflfflfflfflfflzfflfflfflfflfflfflfflfflfflfflfflfflffl fflfflfflfflfflfflfflfflfflzfflfflfflfflfflfflfflfflffl K fDEVELOPMENT OF DISCRETE EQUATIONS FOR ARBITRARY BOUNDARY CONDITIONS 107 Theabovesystemcanbewrittenas T w r¼ 0 8w ; ð5:30Þ F wherer¼ Kdf asin(5.16). PartitioningrinEquation(5.30)congruentwithwgives  r T E T T ½w w  ¼ w r þw r ¼ 0 8w : ð5:31Þ E E F F F E F r F Asw ¼ 0andw isarbitrary,itfollowsfromthescalarproducttheoremthatr ¼ 0.Equation(5.16)can E F F thenbewritteninthepartitionedformas    r K K f d E E EF E E r¼ ¼  ; T 0 K K d f F F EF F whereK ,K andK arepartitionedtobecongruentwiththepartitionsofdandf. E F EF Theaboveequationcanberewrittenas   K K d f þr E EF E E E ¼ : ð5:32Þ T K K d f F F F EF Usingthetwo-stepapproachdiscussedinSection5.1,wefirstsolvefortheunknowndiscretesolutiond by F usingthesecondrowintheabove: T  K d ¼ f K d : ð5:33Þ F F F E EF Onced isknown,theunknownreactionscanbecomputedfromthefirstrowof(5.32): F  r ¼ K d þK d f : ð5:34Þ E E E EF F E For purposes of postprocessing, the displacements and stresses are computed in each element using Equation(5.8)andthestress–strainlaw: e e e e e e e u ðxÞ¼ N ðxÞd; ðxÞ¼ E ðxÞB ðxÞd : e e e TheelementnodalvaluesareobtainedbythegatheroperatorL usingd ¼ L d. An important part of postprocessing is the visual depiction of these results. These are invaluable in interpretingtheresultsandassessingwhetherthemodelisappropriateandhasbeensolvedcorrectly.The variety and richness of visualization in one-dimensional problems is limited, but we will see that visualizationintwodimensionsisquiteimportant. Example5.2. Taperedelasticbar ConsideraproblemofanaxiallyloadedelasticbarasshowninFigure5.10.Dimensionsareinmeters. Solvefortheunknowndisplacementandstresses withafiniteelement (nel=3,nel=1)meshconsistingof asinglethree-nodeelement(n ¼3,n ¼1)asshowninFigure5.11. en el108 FINITE ELEMENT FORMULATION FOR ONE-DIMENSIONAL PROBLEMS E = 8 Pa A = 2x u(x =2) = 0 P(x = 5) = 24 N x t (x = 6) = 0 x = 2 −1 b = 8 N / m x = 6 Figure 5.10 Geometry,loadsandboundaryconditionsofExample5.2. Recallthattheelementshapefunctionsforthethree-nodequadraticelementare ð1Þ ð1Þ ðxx Þðxx Þ ðx4Þðx6Þ 1 ð1Þ 2 3 N ¼ ¼ ¼ ðx4Þðx6Þ; 1 ð1Þ ð1Þ ð1Þ ð1Þ ð2Þð4Þ 8 ðx x Þðx x Þ 1 2 1 3 ð1Þ ð1Þ ðxx Þðxx Þ ðx2Þðx6Þ 1 ð1Þ 1 3 N ¼ ¼ ¼ ðx2Þðx6Þ; 2 ð1Þ ð1Þ ð1Þ ð1Þ ð2Þð2Þ 4 ðx x Þðx x Þ 2 1 2 3 ð1Þ ð1Þ ðxx Þðxx Þ ðx2Þðx4Þ 1 ð1Þ 1 2 N ¼ ¼ ¼ ðx2Þðx4Þ; 3 ð1Þ ð1Þ ð1Þ ð1Þ ð4Þð2Þ 8 ðx x Þðx x Þ 3 1 3 2 andthecorresponding B-matrixis ð1Þ ð1Þ ð1Þ dN 1 dN 1 dN 1 ð1Þ ð1Þ ð1Þ 1 2 3 B ¼ ¼ ðx5Þ; B ¼ ¼ ð4xÞ; B ¼ ¼ ðx3Þ; 1 2 3 dx 4 dx 2 dx 4 1 ð1Þ B ¼ ½ðx5Þð82xÞðx3Þ: 4 Stiffnessmatrix Theelementstiffnessmatrixisgivenby 2 3 ðx5Þ x 3 6 Z Z 1 1 6 7 ð1Þ ð1ÞT ð1Þ ð1Þ ð1Þ K ¼ K¼ B A E B dx¼ 4ð82xÞ5ð2xÞð8Þ ½ðx5Þð82xÞðx3Þdx 4 4 x 2 1 ðx3Þ 2 3 2 6 xðx5Þ xðx5Þð82xÞ xðx5Þðx3Þ Z 6 7 2 6 7 ¼ dx: xð82xÞðx5Þ xð82xÞ xð82xÞðx3Þ 4 5 2 2 xðx3Þðx5Þ xðx3Þð82xÞ xðx3Þ 1 2 3 x (1) (1) (1) P x = 2 x = 4 x = 6 1 2 3 Figure 5.11 FiniteelementmeshofExample5.2.DEVELOPMENT OF DISCRETE EQUATIONS FOR ARBITRARY BOUNDARY CONDITIONS 109 Itcanbeseenthattheintegrandiscubicðp¼ 3Þ.Sothenumberofquadraturepointsrequiredforexact integration is 2n 1 3, i.e. n  2, that is, two-point Gauss quadrature is adequate for exact gp gp integrationoftheintegrand.TheJacobianis ba J ¼ ¼ 2: 2 Writingxintermsofandtransformingtotheparentdomain,wehave 2 3 6 1 Z Z 4 5 fðxÞdx¼ 2 fðxðÞÞd¼JW fðxð ÞÞþ W fðxð ÞÞ ¼ 2½fðxÞþfðx Þ; ð5:35Þ 1 1 2 2 1 2 z z 2 1 1 1 where  1 x ¼ xðÞ¼ 4þ2 ¼ 4þ2 pffiffiffi ¼ 2:8453; 1 1 1 3  1 x ¼ xðÞ¼ 4þ2 ¼ 4þ2 pffiffiffi ¼ 5:1547: 2 2 2 3 Using(5.35),K isgivenby 11 6 Z 2 2 2 K ¼ xðx5Þ dx¼ 2ð2:8453ð2:84535Þ þ5:1547ð5:15475ÞÞ¼ 26:667: 11 2 Thestiffnessmatrixisgivenby 2 3 2 3 26:67 32 5:33 26:67 32 5:33 4 5 4 5 K¼ 85:33 53:33 ¼ 32 85:33 53:33 : sym 48 5:33 53:33 48 Notethatthestiffnessmatrixissymmetricandthesumofthetermsineachrow(orcolumn)isequalto zero. The latter follows from the fact that under rigid body motion (for instance, when the nodal displacementsareallequalto1)theresultingnodalforcesmustbezero. Bodyforcematrix Thematrixofthenodalbodyforcesisobtainedbyaddingthecontributionsfromthedistributedloadingb (firsttermin(5.36))andthepointforceP(secondtermin(5.36)). x 3 Z ð1Þ eT eT f ¼ f ¼ N bdxþðN PÞj : ð5:36Þ  x¼5  fflfflfflfflfflfflfflzfflfflfflfflfflfflffl x 1 contribution from the point force ThederivationdetailsofthenodalbodyforcesarisingfrompointforcesaregiveninAppendixA5.Note thatthesecondtermin(5.36)consistsofaproductoftheelementshapefunctionsevaluatedatthepoint where the point force is acting and the value of the point force (positive if it acts in the positive x-direction).Forinstance,ifthepointforceisactinginthemiddleofalinearelement,thevalueofthe shapefunctioninthemiddleishalf,sohalfoftheforceflowstoeachnode.110 FINITE ELEMENT FORMULATION FOR ONE-DIMENSIONAL PROBLEMS Inthepresentexample,(5.36)gives 2 3 2 3 6 Z 0:125ðx4Þðx6Þ 0:125ðx4Þðx6Þ 4 5 4 5 f ¼ 0:25ðx2Þðx6Þ 8dxþ 0:25ðx2Þðx6Þ 24:  0:125ðx2Þðx4Þ 0:125ðx2Þðx4Þ 2 x¼5 Two-pointGaussquadratureisneededbecausethefunctionisquadratic,so 6 Z fðxÞdx¼ 2½fðxÞþfðx Þ: 1 2 2 Thus, 2hi3 ð1Þ ð1Þ 2 3 2 N ðxÞþN ðx Þ 1 2 3ð54Þð56Þ 1 1 6 7 6 ð1Þ ð1Þ 7 4 5 f ¼ 8 þ 6ð52Þð56Þ  4 2½N ðxÞþN ðx Þ5 1 2 2 2 ð1Þ ð1Þ 3ð52Þð54Þ 2½N ðxÞþN ðx Þ 1 2 3 3 2 3 2 3 2ðð2:84534Þð2:84536Þþð5:15474Þð5:15476ÞÞ 3 6 7 6 7 ¼ 4ðð2:84532Þð2:84536Þþð5:15472Þð5:15476ÞÞ þ 18 4 5 4 5 2ðð2:84532Þð2:84534Þþð5:15472Þð5:15474ÞÞ 9 fflfflfflzfflfflffl 2 3 2 3 2 3 sum¼ 24 5:33 3 2:33 4 5 4 5 4 5 ¼ 21:33 þ 18 ¼ 39:33 : 5:33 9 14:33 fflfflfflfflfflfflzfflfflfflfflfflffl fflfflfflzfflfflffl 84 24 Notethattheboundaryforcematrixvanishes,exceptforthereactionatnode1.ThustheRHSof(5.32)is: 2 3 r þ2:33 1 4 5 fþr¼ 39:33 : 14:33 Theresultingglobalsystemofequationsis wherewehavepartitionedtheequationsafterthefirstrowandcolumn.Thereducedsystemofequations are: T  K d ¼ f K d : F F F E EF fflfflfflzfflfflffl 0 Solvingtheabove   85:33 53:33 u 39:33 u 2:1193 2 2 ¼ ) ¼ : 53:33 48 u 14:33 u 2:6534 3 3TWO-POINT BOUNDARY VALUE PROBLEM WITH GENERALIZED 111 s 36– 4x 15 x 10 24– 4x x 5 x 2 4 6 Figure 5.12 Comparisonofthefiniteelement(solidline)andexactstresses(dashedline)forExample5.2. Postprocessing Once the nodal displacements have been calculated, the displacement field can be obtained by (5.3). Writingthisequationforthethree-nodeelementgives 2 3 0 ð1Þ ð1Þ ð1Þ 6 7 ð1Þ u¼ N u þN u þN u ; d¼ d ¼42:11935: 1 2 3 1 2 3 2:6534 1 1 1 uðxÞ¼ ðx4Þðx6Þð0Þþ ðx2Þðx6Þð2:1193Þþ ðx2Þðx4Þð2:6534Þ 8 4 8 2 ¼0:19815x þ2:24855x3:7045: Thestressfieldisgivenby du d ð1Þ ð1Þ ð1Þ ð1Þ ðxÞ¼ E ¼ E ðN d Þ¼ EB d dx dx 2 3 0 1 6 7 ¼ 8 ½ðx5Þð82xÞðx3Þ 2:1193 ¼3:17xþ17:99: 4 5 4 2:6534 Estimationofsolutionquality Forbrevity,onlythequalityofthestresswillbeassessed.Astheproblemisstaticallydeterminate,the exactstress fieldcan be calculated from the axial force pðxÞ by dividing it by the cross-sectional area pðxÞ ex  ¼ .Figure5.12comparestheFEsolutionofthestressfield(shownwithasolidline)withthe 2x exactstressfield(shownwithadashedline).NoticethattheFEstressfielddoesnotcapturethejumpthat occursatthelocationofthepointforce. 5.5 TWO-POINT BOUNDARY VALUE PROBLEM WITH GENERALIZED 2 BOUNDARY CONDITIONS Wewillnowconsideratwo-pointboundaryvalueproblemwithgeneralizedboundaryconditions.Wewill firstconsiderthepenaltymethod(Equation(3.62)),followedbythepartitionmethod(Equation(3.63)). 2 Recommended for Advanced Track.112 FINITE ELEMENT FORMULATION FOR ONE-DIMENSIONAL PROBLEMS Inthepenaltymethod,theessentialboundaryconditionsareconsideredasalimitingcaseofthenatural boundary conditions; thus, the natural boundary extends over the entire boundary. The weak form is repeatedhereforconvenience: 1 findðxÞ2 H suchthat Z Z   dw d 1    A dx wf dxwAðbðÞÞ ¼ 0 8w2 H ; ð5:37Þ  dx dx    wherethefieldsandparametersaredefinedinTable3.2. Inthisapproach,therearenoessentialboundaryconditions,soallofthenodalvaluesindandwarefree. Integratingtheweak(5.37)overelementdomainsandsubstitutinginterpolants(5.8)intotheweakform yields 8 9 Z  Z  n = el X   eT eT e e e e eT e e e eT eT e     w B  A B dxd þðN A bN Þ d  N f dxðN A ðþbÞÞ ¼ 0 8w:   : ; e e   e¼1 e e   ð5:38Þ e where isaportionofelementboundaryonexternalboundary.Wedefinethefiniteelementmatrices:  Z  e eT e e e eT e e  K ¼ B  A B dxþðN A bN Þ  e  e  ð5:39Þ  Z  e eT eT e    f ¼ N f dxþðN A ðþbÞÞ :  e  e  e e e e Substituting (5.39) into (5.38), using w ¼ L w, d ¼ L d and defining global matrices by (5.13) and (5.14)givesthediscreteweakform T w r¼ 0 8w; ð5:40Þ whereristheresidualmatrixdefinedin(5.16).Duetoarbitrarinessofw,itfollowsthat r¼ Kdf ¼ 0 or Kd¼ f: ð5:41Þ In (5.41), no partitioning or node renumbering is required; the essential boundary conditions are easily enforcedbyselectingbtobealargepenaltyparameter. Wenowturntothepartitionmethod,whichwasusedinSection5.1.Thegeneralweakformisstated(see Box3.6)as findðxÞ2 U suchthat Z Z   dw d    A dx wf dxwAðbðxÞððxÞðxÞÞÞ ¼ 0 8w2 U : ð5:42Þ 0  dx dx