Biomedical Engineering - Computational Biomechanics Laboratory

Appunti lezioni completi FEM

Divided by topic

Fundamental concepts of Finite Elements for Elastic SolidsDarioGastaldiComputationalBiomechanicsLaboratory1.WhatisaFEA?Basicsteps1.Idealization2.Discretization3.Solution2.Typesofproblems3.Examples4.Theelasticproblem.StrongForm5.TheintegralformoftheproblemContent2 The Finite Element Method (FEM) allows finding an approximate numerical solution to mathematical model of a physical problem What is a FEA?3IDEALIZATIONDISCRETIZATIONSOLUTIONPhysicalSystemMathematicalModelDiscreteModelDiscreteSolutionSolution errorSolution + discretization errorSolution + discretization + modeling errorPhysical modelLinear Elasticity•Governing equations•Boundary conditionsFE ResultsWhat is a FEA? Basic Steps4FE Model Mathematical ModelPhysical ModelIdealizationModelingframework•Governing equations (set of partial differential equations)•Boundary conditionsWhat is a FEA? Basic Steps5For each element we calculate:Library of elementsAssembly[K] {u} = {F}[Ke] {ue} = {Fe}[Ke] {ue} = {Fe}[Ke] {ue} = {Fe}Analytic solutionF(x,y, z)Numerical solution discretized in nodesWhat is a FEA? Basic Steps6Discretization [K] {u} = {F}Commercial softwares:•Abaqus•Ansys•Adina•Nastran•LS-Dyna ...{u} = [K]-1{F}Solution at nodesResolutionStiffness MatrixUnknown DisplacementsLoadsKPreprocessingPostprocessingWhat is a FEA? Basic Steps7SolutionDiscretization:Divisionin elementsCompute elementalstiffnessmatrices Assembly:Global stiffnessmatrixand Load vectorImpositionofboundaryconditionsResolution of thesystem of equations:Displacements Displacementsin Global orLocalcoordinatesDeterminationof stresses,and deformations** Numerical processPreprocessing(user)Postprocess(user)•Governing equations•Boundary conditionsFEM: Steps in the solution8 Deformable SolidsBeamsTr u s s e sCablesPipesPlane StressPlane strainAxisymmetryHexahedral ElementsTe t r a h e d r a l E l e m e n t sOthersStaticDynamicTr u s s e s1D ElementsPlane Elasticity2D ElementsThree dimensional solids3D ElementsLinearNon linearModelingGeometric classification of the solidsFractureMaterialGeometricTypes of problems9Automotive (static)Examples11 Automotive (dynamic)Examples12https://www.youtube.com/watch?v=hrfcROMz2II13Biomechanics (static)Examplesmechanical compatibility assessment of knee prosthesisanalysis of contact pressure distribution for wear risk assessment 14Fatiguestrengthtest simulation(standard ISO 14879-1)Vo n Misesstress distributionFBiomechanics (static)ExamplesBiomechanics (static)Examples15 ( ) ( ) ( ) ( ) 0 1 3 1 2 6 2 12 2 23 5 2 3 1 4 2 2 3 3 1 2 2 1 < - s + s s + s + s + s + s + s + s + s + s = s a a a a a a ) ( f ij Tsai-Wostress criterion(anisotropicnon symmetricand brittlematerials) Biomechanics (dynamic)Examples16Curtesy of Prof. J.F. Rodriguez Matas (LaBS, Politecnico di Milano)Biomechanics (quasi-static)Examples17Morlacchi, Chiastra et al. EuroIntervention, 2014 Biomechanics (quasi-static)Examples18Curtesy of G. Luraghi (LaBS, Politecnico di Milano)Upper gripFemoral componentUHMWPE insertTibial Component Lower gripKNEE PROSTESYS EXPERIMENTAL SET-UPBiomechanics(dynamics)Examples 20KNEE PROSTESYS EXPERIMENTAL SET-UP FEM modelLinear ElasticityPhysicalmodelFE ModelWhat is a FEA? Basic Steps21IDEALIZATIONDISCRETIZATIONSOLUTIONPhysicalSystemMathematicalModelDiscreteModelDiscreteSolutionSolution errorSolution + discretization errorSolution + discretization + modeling error The linear elasticity frameworkStress: ijnj= ti ij= ji ·n= t = T Kinematic equations: ✏ij= 1 2 ✓@ui @xj+ @uj @xi ◆ ✏= Hu 8x2⌦ Constitutive equations: ij= Cij kl ✏kl = D✏= DHu = Su Hand Qare differential operatorsBoundary conditions: ⇢ ui=¯ ui 8x2u ijnj= ¯ti 8x2t Strong FormulationEquilibrium equations: b= Q= QSu = Au 8x2⌦ ¯ttu⌦ ⌘ @⌦ ¯u=0b The elastic problem22Mathematical modelThe Voigt NotationTensorial notationVo i g t n o t a t i o n ✏= 8>>>>>>< >>>>>>: ✏xx ✏yy ✏zz xy xz yz 9>>>>>>= >>>>>>; ij= 0 @ xx xy xz xy yy yz xz yz zz 1 A ✏ij= 0 @ ✏xx xy xz xy ✏yy yz xz yz ✏zz 1 A = 8>>>>>>< >>>>>>: xx yy zz xy xz yz 9>>>>>>= >>>>>>; The elastic problem23Mathematical model Differential operatorsDeformation ✏= Hu 3D problems2D problems1D problems ✏xx = @ @xu Kinematic Operator 0 BBBBBB@ ✏xx ✏yy ✏zz xy xz yz 1 CCCCCCA = 2 6666666664 @@x 00 0 @@y 0 00 @@z @@y @@x 0 @@z 0 @@x 0 @@z @@y 3 7777777775 0 @ u v w 1 A 0 @ ✏xx ✏yy xy 1 A = 2 64 @@x 0 0 @@y @@y @@x 3 75 ✓ u v ◆ The elastic problem24Mathematical modelConstitutive Equations 0 BBBBBB@ xx yy zz xy xz yz 1 CCCCCCA = 2 6666664 2G + 000 2G + 000 2G + 000 000 G 00 0000 G 0 00000 G 3 7777775 0 BBBBBB@ ✏xx ✏yy ✏zz xy xz yz 1 CCCCCCA = D✏ = DHu = Su Stress Operator 0 BBBBBB@ xx yy zz xy xz yz 1 CCCCCCA = 2 6666664 2G + 000 2G + 000 2G + 000 000 G 00 0000 G 0 00000 G 3 7777775 2 6666666664 @@x 00 0 @@y 0 00 @@z @@y @@x 0 @@z 0 @@x 0 @@z @@y 3 7777777775 0 @ u v w 1 A The elastic problem25Mathematical model Constitutive Equations. Special CasesPlane StressPlane StrainOne Dimension xx = E✏xx 0 @ xx yy xy 1 A = E 1⌫2 2 4 1 ⌫ 0 ⌫ 10 00 1⌫2 3 5 2 64 @@x 0 0 @@y @@y @@x 3 75 ✓ u v ◆ 0 @ xx yy xy 1 A = E (1 + ⌫)(1 2⌫) 2 4 1⌫⌫ 0 ⌫ 1⌫ 0 00 12⌫ 3 5 2 64 @@x 0 0 @@y @@y @@x 3 75 ✓ u v ◆ The elastic problem26Mathematical modelDifferential operatorsEquilibrium equations 2 64 @@x 00 @@y @@z0 0 @@y 0 @@x 0 @@z 00 @@z 0 @@x @@y 3 75 0 BBBBBB@ xx yy zz xy xz yz 1 CCCCCCA = 0 @ bx by bz 1 A 0 @ bx by bz 1 A b= Q= QDHu = QSu = Au 2 64 @@x 00 @@y @@z0 0 @@y 0 @@x 0 @@z 00 @@z 0 @@x @@y 3 75 2 6666664 2G + 000 2G + 000 2G + 000 000 G 00 0000 G 0 00000 G 3 7777775 2 6666666664 @@x 00 0 @@y 0 00 @@z @@y @@x 0 @@z 0 @@x 0 @@z @@y 3 7777777775 0 @ u v w 1 A = 0 @ bx by bz 1 A The elastic problem27Mathematical model D = 2 6666664 2G + 000 2G + 000 2G + 000 000 G 00 0000 G 0 00000 G 3 7777775 Kinematic OperatorMaterial MatrixStress Operator S= DH = 2 66666666666666666664 (2G +)@ @x @ @y @ @z @ @x (2G +)@ @y @ @z @ @x @ @y (2G +)@ @z G @ @y G @ @x 0 G @ @z 0 G @ @x 0 G @ @z G @ @y 3 77777777777777777775 H = 2 6666666664 @@x 00 0 @@y 0 00 @@z @@y @@x 0 @@z 0 @@x 0 @@z @@y 3 7777777775 The elastic problem28Mathematical modelDifferential operatorsThe strong formulation (or differential form) can be written as: Au = b 8x2⌦ ⇢ ui=¯ ui 8x2u ijnj= ¯ti 8x2t The Finite Element Method (FEM) allows to find an approximate solution to this problemIn order to use FEM we first need to find an Integral form(or weak form) of the problemTwo choices:1.Potential Energy (only for conservative systems)2.Galerkin’s Method ¯ttu⌦ ⌘ @⌦ ¯u=0b The elastic problem29Mathematical model For an elastic body, the total potential energy P is defined as U = 1 2 Z ⌦T✏dV WP = Z ⌦uTbdV Z uTtdS X i uTiPi Π= U+ WP ⇧ = 1 2 Z ⌦T✏dV Z ⌦uTbdV Z uTtdS X i uTiPi Principle of Minimum Potential EnergyFor conservative systems, of all the kinematicallyadmissible displacement fields, those corresponding to equilibrium extremizethe total potential energy. If the extremum condition is a minimum, the equilibrium state is stableKinematically admissible displacements are those that satisfy compatibility and the boundary conditionsThe integral form of the problem30Potential energy31The weak form of the problemTheminimumof potentialenergycorrespondsto thePrincipleof virtual Works: 3D solidsplanar beamsplates 32The integral form of the problemlinear elastic constitutive lawFrom differential problem to algebraic problem?!?3D Differential operator33The finite element equations 34The finite element equationsFrom global structure (U, F) to single finite element e (ue, Fe) ?!?Assembly procedure by means of booleanmatrix Ce35The finite element equationsNumerical integration…Geometrical mappingGauss integration formula… Isoparametricelements Tw o d i m e n s i o n a l p r o b l e m sDarioGastaldiComputationalBiomechanicsLaboratory1.GoverningEquation2.Discretization3.Fournodesquadrilateralelement1.Shapefunctions2.ElementStiffness3.ElementLoadVector4.ConstantStrainTriangle1.Isoparametricrepresentation2.ElementStiffness3.ForcetermContents2 Let us consider a two-dimensional solid:Displacement vector: u=[ uv ]T ¯ttu⌦ ⌘ @⌦ ¯u=0b Pxy 0 @ ✏xx ✏yy xy 1 A = 2 64 @@x 0 0 @@y @@x @@y 3 75 ✓ u v ◆ Strain-displacement relations (Kinematic operator):Constitutive equation: = D✏ (plane stress or plane strain)Equilibrium equation: Z ⌦T✏(⌘)dV = Z ⌦⌘TbdV + Z ⌘TtdS +X i ⌘TPi Governing Equations3Stresses and strains: =[ xx yy xy] ✏ =[ ✏xx ✏yy ✏xy] 0 @ ✏xx ✏yy xy 1 A = 2 64 @@x 0 0 @@y @@x @@y 3 75 ✓ u v ◆ Finite Element discretization with Bi-Linear Quadrilateral ElementsP ¯t ⌘ @⌦ u ¯u=0 4Discretization tu ⌘ @⌦ ¯u=0 PFinite Element discretization with Linear Triangular Elements ¯t 1.Each node has two degrees of freedom i.e., horizontal and vertical displacementsThe standard convention establishes the local numbering by going around the element in counterclockwise direction to avoid calculating a negative element area ⇢ 2i1 xdisp . 2iy disp .,i =1 ,2,3 xy23Bi-linear quadrilateral Element local numbering41 u1⌘ q1 v1⌘ q2 u2⌘ q3 v2⌘ q4 u3⌘ q5 v3⌘ q6 v4⌘ q8 u4⌘ q7 u(e)= 2 66666666664 u1 v1 u2 v2 u3 v3 u4 v4 3 77777777775 ⌘ 2 66666666664 q1 q2 q3 q4 q5 q6 q7 q8 3 77777777775 5Four node Quadrilateral Element i= 1, 2, 3, 42.The local displacement vector is ordered as:rs123(1,-1)(-1,-1)(-1,1)4(1,1)Cartesian systemNatural systemxy2341 u1⌘ q1 v1⌘ q2 u2⌘ q3 v2⌘ q4 u3⌘ q5 v3⌘ q6 v4⌘ q8 u4⌘ q7 6Four node Quadrilateral Element IsoparametricRepresentationNatural Element System rs123(1,-1)(-1,-1)(-1,1)4(1,1) N1= 1 4(1 r)(1 s) N2= 1 4(1 + r)(1 s) N3= 1 4(1 + r)(1 + s) N4= 1 4(1 r)(1 + s) 1234123412341234N1N2N3N47Four node Quadrilateral Element IsoparametricRepresentationShape functions u= Nu (e) x= Nx (e) u = N1u(e)1 +N2u(e)2 +N3u(e)3 +N4u(e)4 v = N1v(e)1 +N2v(e)2 +N3v(e)3 +N4v(e)4 x = N1x(e)1 +N2x(e)2 +N3x(e)3 +N4x(e)4 y = N1y(e)1 +N2y(e)2 +N3y(e)3 +N4y(e)4 N =  N1 0 N2 0 N3 0 N4 0 0 N1 0 N2 0 N3 0 N4 Geometryand displacementsare approximated using the same interpolation functions8Four node Quadrilateral Element IsoparametricRepresentationDisplacement and Geometry interpolation To calculate the strains we need to derive the displacement field with respect to the Cartesian Coordinates. Let us consider the displacement along the x-direction, u. @u @r= @u @x @x @r+ @u @y @y @r @u @s= @u @x @x @s+ @u @y @y @s @u@r@u@s ! = " @x@r @y@r @x@s @y@s # @u@x@u@y ! @u@x@u@y ! = J1 @u@r@u@s ! J= " @x@r @y@r @x@s @y@s # 9Four node Quadrilateral Element IsoparametricRepresentation. Element strainFor the quadrilateral element, the Jacobian matrix is not a constant matrixThe derivative of the displacement in the x-direction is found to be @u@x@u@y ! = 1 det J  J22 J12 J21 J11 @u@r@u@s ! Similarly for the derivative of the displacement in the y-direction @v@x@v@y ! = 1 det J  J22 J12 J21 J11 @v@r@v@s ! 10Four node Quadrilateral Element IsoparametricRepresentation. Element strain The strain field is given as 0 B@ ✏xx ✏yy xy 1 CA = 0 B@ @u@x@v@y @u@y+ @v@x 1 CA = A 0 BBB@ @u@r@u@s@v@r@v@s 1 CCCA A = 1 det J 2 4 J22 J12 00 00 J21 J11 J21 J11 J22 J12 3 5 0 BBB@ @u@r@u@s@v@r@v@s 1 CCCA = Gu (e) 11Four node Quadrilateral Element IsoparametricRepresentation. Element strain G = 1 4 2 664 (1 r)0(1 r)0(1+ r)0 (1 + r)0 (1 s)0(1+ s)0(1+ s)0 (1 s)0 0 (1 r)0(1 r)0(1+ r)0 (1 + r) 0 (1 s)0(1+ s)0(1+ s)0 (1 s) 3 775 ✏= AGu (e)= Bu (e) Gis the matrix corresponding to the derivative of the shape functionsThe strain is not any longer constant within the element 12Four node Quadrilateral Element IsoparametricRepresentation. Element strain The element stress is found as = D✏= DBu (e) where Dis D = E (1 + ⌫)(1 2⌫) 2 4 1⌫⌫ 0 ⌫ 1⌫ 0 00 12⌫ 3 5 D = E 1⌫2 2 4 1 ⌫ 0 ⌫ 10 00 1⌫2 3 5 Plane strain:Plane stress:13Four node Quadrilateral Element IsoparametricRepresentation. Element stressThe contribution to the internal work of element eis IW e= Z e✏(⌘)TdV Substituting the expressions for and , ✏ IW e = Z e(⌘(e))TBTDBu (e)tedA =( ⌘(e))TZ1 1 Z1 1teBTDB det Jdrds u(e) where we have used dA= detJdrdsto change the integration from the Cartesian system to the natural element system K(e)= Z1 1 Z1 1teBTDB det Jdrds Note that this integral cannot be solved analytically. Therefore numerical integration is required8x8 matrix14Four node Quadrilateral Element Element stiffness The contribution to body force load vector of element eis: Z e⌘TbdV =( ⌘(e))TZ eNTbtedA =( ⌘(e))TZ1 1 Z1 1teNTbdet Jdrds f(e)If the body force is constant in the element:8x1 matrix f(e)= Z1 1 Z1 1teNTdet Jdrds ✓ bx by ◆ If the body force is not constant in the element:8x1 matrix f(e)= Z1 1 Z1 1teNTN det Jdrds 0 BBBBBB@ (bx)(e)1 (by)(e)1... (bx)(e)8 (by)(e)8 1 CCCCCCA 15Four node Quadrilateral Element Force term. Body force Z e ⌘TTtedl Consider a constant traction vector T=(Tx, Ty)along the edge l23xy23Tx4Ty1rs123(1,-1)(-1,-1)(-1,1)4(1,1)TxTyNote that on edge 23, r=1and therefore N1=N4=0. In addition N2= 0.5(1-s), N3= 0.5(1+s)16Four node Quadrilateral Element Force term. Traction force The contribution of constant tractions action on edge 23 are then given as: T(e)= tel23 2 [00 Tx Ty Tx Ty 00 ]T They are directly accounted for in the global Assembly processForce term. Point forcesFor the remaining edges, similar expressions are obtained T(e)= tel12 2 [Tx Ty Tx Ty 0000 ]T T(e)= tel34 2 [0000 Tx Ty Tx Ty]T T(e)= tel14 2 [Tx Ty 0000 Tx Ty]T Edge 12Edge 34Edge 1417Four node Quadrilateral Element Force term. Traction force1.GoverningEquation2.Discretization3.Fournodesquadrilateralelement1.Shapefunctions2.ElementStiffness3.ElementLoadVector4.ConstantStrainTriangle1.Isoparametricrepresentation2.ElementStiffness3.ForcetermContents18 The contribution of constant tractions action on edge 23 are then given as: T(e)= tel23 2 [00 Tx Ty Tx Ty 00 ]T They are directly accounted for in the global Assembly processForce term. Point forcesFor the remaining edges, similar expressions are obtained T(e)= tel12 2 [Tx Ty Tx Ty 0000 ]T T(e)= tel34 2 [0000 Tx Ty Tx Ty]T T(e)= tel14 2 [Tx Ty 0000 Tx Ty]T Edge 12Edge 34Edge 1417Four node Quadrilateral Element Force term. Traction force1.GoverningEquation2.Discretization3.Fournodesquadrilateralelement1.Shapefunctions2.ElementStiffness3.ElementLoadVector4.ConstantStrainTriangle1.Isoparametricrepresentation2.ElementStiffness3.ForcetermContents18 1.Each node has two degrees of freedom i.e., horizontal and vertical displacements2.The local displacement vector is ordered as:The standard convention establishes the local numbering by going around the element in counterclockwise direction to avoid calculating a negative element area ⇢ 2i1 xdisp . 2iy disp .,i =1 ,2,3 u(e)= 2 6666664 u1 v1 u2 v2 u3 v3 3 7777775 ⌘ 2 6666664 q1 q2 q3 q4 q5 q6 3 7777775 xy123Linear Triangular Element (local numbering) u1⌘ q1 v1⌘ q2 u2⌘ q3 v2⌘ q4 u3⌘ q5 v3⌘ q6 19Constant Strain Trianglewherexy123(x1,y1)(x2,y2)(x3,y3)P(x,y)A1A2A3 N1= A1 A ,N 2= A2 A ,N 3= A3 A N1+N2+N3=1 Cartesian System Ai= 1 2A(ai+bix+ciy) 2A = x13y23 x23y13 a1= x23y3+y32x3 b1= y23 c1= x32 a2= x31y1+y13x1 b2= y31 c2= x13 a3= x12y2+y21x2 b3= y12 c3= x21 Isoparametricrepresentation: u= P 3i=1 Niu(e)i ,x = P 3i=1 Nix(e)i v= P 3i=1 Niv(e)i ,y = P 3i=1 Niy(e)i xij= x(e)i x(e)j ,y ij= y(e)i y(e)j 20IsoparametricrepresentationConstant Strain TriangleThe interpolation functions for a linear triangular element can be defined in terms of the cartesiansystem or a natural element system Natural Element Systemrs123(1,0)(0,0)(0,1) N1= r, N 2= s, N3=1 rs u= Nu (e) ✓ u v ◆ =  N1 0 N2 0 N3 0 0 N1 0 N2 0 N3 0 BBBBBBBB@ u(e)1 v(e)1 u(e)2 v(e)2 u(e)3 v(e)3 1 CCCCCCCCA ✓ x y ◆ =  N1 0 N2 0 N3 0 0 N1 0 N2 0 N3 0 BBBBBBBB@ x(e)1 y(e)1 x(e)2 y(e)2 x(e)3 y(e)3 1 CCCCCCCCA The interpolation functions can also be defined in terms of a local system r-s. Defining: r⌘ A1/A, s ⌘ A2/A Isoparametricrepresentation: u= P 3i=1 Niu(e)i ,x = P 3i=1 Nix(e)i v= P 3i=1 Niv(e)i ,y = P 3i=1 Niy(e)i 21IsoparametricrepresentationConstant Strain Triangle123N3123N1N2123Representation of Shape Functions22IsoparametricrepresentationConstant Strain Triangle To calculate the strains we need to derive the displacement field with respect to the Cartesian Coordinates. Let us consider the displacement along the x-direction, u. @u @r= @u @x @x @r+ @u @y @y @r @u @s= @u @x @x @s+ @u @y @y @s @u@r@u@s ! = " @x@r @y@r @x@s @y@s # @u@x@u@y ! @u@x@u@y ! = J1 @u@r@u@s ! J= " @x@r @y@r @x@s @y@s # 23IsoparametricRepresentation. Element strainConstant Strain TriangleFor the triangular element, the Jacobian matrix is a constant matrix J1= 1 2A  y23 y13 x23 x13 2A =det J J=  x13 y13 x23 y23 24IsoparametricRepresentation. Element strainConstant Strain Triangle @u@x@u@y ! = 1 2A " y23@u@r y13@u@s x23@u@r x13@u@s # = 1 2A " y23(u1u3) y13(u2u3) x23(u1u3) x13(u1u3) # The derivative of the displacement in the x-direction is found to be @u@x@u@y ! = 1 2A " y23@u@r y13@u@s x23@u@r x13@u@s # = 1 2A " y23(u1u3) y13(u2u3) x23(u1u3) x13(u1u3) # Similarly for the derivative of the displacement in the y-direction @v@x@v@y ! = 1 2A " y23@v@r y13@v@s x23@v@r x13@v@s # = 1 2A " y23(v1v3) y13(v2v3) x23(v1v3) x13(v1v3) # @u@x@u@y ! = 1 2A " y23@u@r y13@u@s x23@u@r x13@u@s # = 1 2A " y23(u1u3) y13(u2u3) x23(u1u3) x13(u1u3) # The strain field for a triangular element is found to be ✏= 0 B@ ✏xx ✏yy xy 1 CA = 1 2A 2 4 y23 0 y31 0 y12 0 0 x32 0 x13 0 x21 x32 y23 x13 y31 x21 y12 3 5 0 BBBBBBBB@ u(e)1 v(e)1 u(e)2 v(e)2 u(e)3 v(e)3 1 CCCCCCCCA ✏= Bu (e) The strain field is constant within the element25IsoparametricRepresentation. Element strainConstant Strain TriangleElement strain-displacement matrixThe element stress is found as = D✏= DBu (e) where Dis D = E (1 + ⌫)(1 2⌫) 2 4 1⌫⌫ 0 ⌫ 1⌫ 0 00 12⌫ 3 5 D = E 1⌫2 2 4 1 ⌫ 0 ⌫ 10 00 1⌫2 3 5 Plane strain:Plane stress:26IsoparametricRepresentation. Element stressConstant Strain Triangle The contribution to the internal work of element eis K(e)= teAeBTDB K(e)is a symmetric matrix IW e= Z e✏(⌘)TdV Substituting the expressions for and, ✏ IW e = Z e(⌘(e))TBTDBu (e)tedA =( ⌘(e))TBTDBu (e)te Z edA =( ⌘(e))TteAeBTDBu (e) =( ⌘(e))TK(e)u(e) 27Element stiffnessConstant Strain TriangleThe force term in the principle of virtual work is associated with the external workBody ForcesConcentrate ForcesTraction Forces EW = Z ⌦⌘TbdV + Z ⌘TtdS +X i ⌘TPi The contribution of element eto the force term is: EW e= Z e⌘TbdV + Z e ⌘TtdS Concentrate Loads are accounted for during the assembly process28Force termConstant Strain Triangle Considering a constant body force on the element, and using Z e⌘TbdV = Z e(⌘xbx+⌘yby)tedA ⌘x= N1(⌘x)(e)1 +N2(⌘x)(e)2 +N3(⌘x)(e)3 ⌘y= N1(⌘y)(e)1 +N2(⌘y)(e)2 +N3(⌘y)(e)3 Z e⌘TbdV =( ⌘x)(e)1 tebx Z eN1dA +( ⌘x)(e)2 tebx Z eN2dA +( ⌘x)(e)3 tebx Z eN3dA +( ⌘y)(e)1 teby Z eN1dA +( ⌘y)(e)2 teby Z eN2dA +( ⌘y)(e)3 teby Z eN3dA 29Force term. Body forceConstant Strain Triangle1N123h=1 Z eN1dA represents the volume of the tetrahedron Z eN1dA = 1 3Ae·h= 1 3Ae similarly Z eN2dA = Z eN3dA = 1 3Ae 30Force term. Body forceConstant Strain TriangleThe vector of elemental body forces is: b(e)= teAe 3 [bx by bx by bx by]T Z e⌘TbdV =( ⌘(e))Tb(e)=( ⌘(e))TteAe 3 [bx by bx by by by]T Z e⌘TbdV =( ⌘(e))Tb(e)=( ⌘(e))TteAe 3 [bx by bx by by by]T Traction forces act on the edges of the element. We could have the following traction forces acting on the edge of the element xy123p2p3Normal PressureComponent Distributionxy123Ty3Ty2xy123Tx3Tx231Force term. Traction forceConstant Strain Triangle l23 = p(x3x2)2+( y3y2)2 cx= y3y2 l23 cy= x3x2 l23 Tx1= cxp1 Tx2= cxp2 Ty1= cyp1 Ty2= cyp2 l23 = p(x3x2)2+( y3y2)2 cx= y3y2 l23 cy= x3x2 l23 Tx1= cxp1 Tx2= cxp2 Ty1= cyp1 Ty2= cyp2 l23 = p(x3x2)2+( y3y2)2 cx= y3y2 l23 cy= x3x2 l23 Tx1= cxp1 Tx2= cxp2 Ty1= cyp1 Ty2= cyp2 l23 = p(x3x2)2+( y3y2)2 cx= y3y2 l23 cy= x3x2 l23 Tx1= cxp1 Tx2= cxp2 Ty1= cyp1 Ty2= cyp2 l23 = p(x3x2)2+( y3y2)2 cx= y3y2 l23 cy= x3x2 l23 Tx1= cxp1 Tx2= cxp2 Ty1= cyp1 Ty2= cyp2 l23 = p(x3x2)2+( y3y2)2 cx= y3y2 l23 cy= x3x2 l23 Tx1= cxp1 Tx2= cxp2 Ty1= cyp1 Ty2= cyp2 l23 = p(x3x2)2+( y3y2)2 cx= y3y2 l23 cy= x3x2 l23 Tx1= cxp1 Tx2= cxp2 Ty1= cyp1 Ty2= cyp2 l23 = p(x3x2)2+( y3y2)2 cx= y3y2 l23 cy= x3x2 l23 Tx1= cxp1 Tx2= cxp2 Ty1= cyp1 Ty2= cyp2 l23 = p(x3x2)2+( y3y2)2 cx= y3y2 l23 cy= x3x2 l23 Tx1= cxp1 Tx2= cxp2 Ty1= cyp1 Ty2= cyp2 Z e ⌘TTtdl = Z l23 (⌘xTx+⌘yTy)tdl Noting that along the edge l23, hx, hy, Tx, andTycan be expressed as: ⌘x =( ⌘x)(e)2 N2+( ⌘x)(e)3 N3 ⌘y =( ⌘y)(e)2 N2+( ⌘y)(e)3 N3 Tx = Tx2N2+Tx3N3 Ty = Ty2N2+Ty3N3 For a traction force applied on edge l23: Z l23 (⌘xTx+⌘yTy)tdl =( ⌘x)(e)2 te ✓ Tx2 Z eN22dl +Tx3 Z eN2N3dA ◆ +( ⌘x)(e)3 te ✓ Tx2 Z eN2N3dl +Tx3 Z eN23dA ◆ +( ⌘y)(e)2 te ✓ Ty2 Z eN22dl +Ty3 Z eN2N3dA ◆ +( ⌘y)(e)3 te ✓ Ty2 Z eN2N3dl +Ty3 Z eN23dA ◆ 32Force term. Traction forceConstant Strain Triangle Noting: Z eN22dl = l23 3, Z eN23dl = l23 3, Z eN2N3dl = l23 6 T(e)= tel23 6 [002 Tx2+Tx3 2Ty2+Ty3 Tx2+2 Tx3 Ty2+2 Ty3]T l23 = p(x3x2)2+( y3y2)2 Z l23 (⌘xTx+⌘yTy)tdl =( ⌘(e))TT(e) (⌘(e))Ttel23 6 [002 Tx2+Tx3 2Ty2+Ty3 Tx2+2 Tx3 Ty2+2 Ty3]T Substituting:33Force term. Traction forceConstant Strain Triangle u= 2 6666666664 U1 V1 U2 V2... UN VN 3 7777777775 ⌘ 2 6666666664 Q1 Q2 Q3 Q4... Q2N1 Q2N 3 7777777775 The global assembly process is the same as described for 1D problems with the only consideration that now each node will have two degrees of freedom instead of one.The local displacement vector (at element level) and the corresponding global displacement vector for a mesh with Nnodes are To each local and global node are associated two dof ⇢ 2i1 xdisp . 2iy disp .,i =1 ,2,3 ⇢ 2i1 xdisp . 2iy disp .,i =1 ,··· ,N The correspondence between the local and global dofis established through the element connectivity array u(e)= 2 6666664 u1 v1 u2 v2 u3 v3 3 7777775 ⌘ 2 6666664 q1 q2 q3 q4 q5 q6 3 7777775 34Global AssemblyConstant Strain Triangle From the example mesh given at the beginning of this presentation, the connectivity data for the first two elements is as followse1231172618292019The global dofassociated with each element is given byeu11v12u23v24u35v3613334515235362171839403738Local numberingGlobal numberingLocal dofGlobal dofThe 6x6 stiffness matrix of element 1 will occupy the rows and columns indicated in the previous table 35Global AssemblyConstant Strain Triangle Three dimensional problemsDarioGastaldiComputationalBiomechanicsLaboratory1.GoverningEquation2.Lineartetrahedral1.Shapefunctions2.ElementStiffness3.Forceterm3.Trilinearhexahedralelement1.Shapefunctions2.ElementStiffness4.3DelementsinAbaqus5.ShearlockingandhourglasseffectsContent2 Constitutive equation: = D✏ Equilibrium equation: Z ⌦T✏(⌘)dV = Z ⌦⌘TbdV + Z ⌘TtdS +X i ⌘TPi Displacement vector: u=[ uvw ]T Strain-displacement relations (Kinematic operator): 0 BBBBBB@ ✏xx ✏yy ✏zz xy xz yz 1 CCCCCCA = 2 6666666664 @@x 00 0 @@y 0 00 @@z @@x @@y 0 @@x 0 @@z 0 @@y @@z 3 7777777775 0 @ u v w 1 A Stresses and strains: =[ xx yy zz xy xz yz ]T ✏ =[ ✏xx ✏yy ✏zz xy xz yz ]T ¯ttu⌦ ⌘ @⌦ ¯u=0b PxyzGoverning EquationsLet us consider a three-dimensional solid:31.Each node has three degrees of freedomi.e., displacements along the coordinate axis2.The local displacement vector is ordered as:xy123Linear Tetrahedron local numbering v1⌘ q2 z u1⌘ q1 4 w1⌘ q3 u2⌘ q4 v2⌘ q5 w2⌘ q6 u3⌘ q7 v3⌘ q8 w3⌘ q9 u4⌘ q10 v4⌘ q11 w4⌘ q12 u(e)= 2 6666666664 u1 v1 w1... u4 v4 w4 3 7777777775 ⌘ 2 6666666664 q1 q2 q3... q10 q11 q12 3 7777777775 8< : 3i2 xdisp . 3i1 ydisp . 3iz disp . ,i =1 ,2,3 Tetrahedral element faces:Face 1: 1, 2, 3Face 2: 1, 4, 2Face 3: 2, 4, 3Face 4: 3, 4, 14Linear Tetrahedral Elementi= 1, 2, 3, 4 Natural Element System u= Nu (e) The interpolation functions can be defined in terms of a local system r-s-t. Isoparametricrepresentation: N1= r, N 2= s, N 3= t, N4=1 rst rt134(0,1,0)(0,0,0)(0,0,1)s2(1,0,0) x= Nx (e) (Nis a 3x12 matrix) N = 2 6666666666666666664 N1 00 0 N1 0 00 N1 N2 00 0 N2 0 00 N2 N3 00 0 N3 0 00 N3 N4 00 0 N4 0 00 N4 3 7777777777777777775 T 5Linear Tetrahedral ElementShape functions To calculate the strains we need to derive the displacement field with respect to the Cartesian Coordinates. Let us consider the displacement along the x-direction, u.6Four node Quadrilateral Element IsoparametricRepresentation. Derivatives of the displacement field!"!#=!"!%!%!#+!"!'!'!#+!"!(!(!#!"!)=!"!%!%!)+!"!'!'!)+!"!(!(!)!"!*=!"!%!%!*+!"!'!'!*+!"!(!(!*!"!#!"!)!"!*=!%!#!'!#!(!#!%!)!'!)!(!)!%!*!'!*!(!*!"!%!"!'!"!(J xij= xixj yij= yiyj zij= zizj J=!"!#!$!#!%!#!"!&!$!&!%!&!"!'!$!'!%!'="()#()$()"*)#*)$*)"+)#+)$+)!"!%!"!'!"!(=+!"!"!#!"!)!"!*J-1 = A =(,-./#*)$+)−#+)$*)#+)$()−#()$+)#()$*)−#*)$()$*)"+)−$+)"*)$+)"()−$()$+)$()"*)−$*)"()"*)#+)−"+)#*)"+)#()−"()#+)"()#*)−"*)#()The Jacobian is a constant matrix as well as the strain within the element. The element stressis found as = D✏= DBu (e) Following the same steps as for the triangle element, the strainfield can be expressed as: ✏= Bu (e) B = 2 66666664 A11 00 A12 00 A13 00 ˆA1 00 0 A21 00 A22 00 A23 00 ˆA2 0 00 A31 00 A32 00 A33 00 ˆA3 0 A31 A21 0 A32 A22 0 A33 A23 0 ˆA3 ˆA2 A31 0 A11 A32 0 A12 A33 0 A13 ˆA3 0 ˆA1 A21 A11 0 A22 A12 0 A23 A13 0 ˆA2 ˆA1 0 3 77777775 ˆA1= A11 +A12 +A13 ˆA2= A21 +A22 +A23 ˆA3= A31 +A32 +A33 A = J1 7Linear Tetrahedral ElementIsoparametricRepresentation. Element strain and stress8Linear Tetrahedral ElementIsoparametricRepresentation. Element stress and volumeNote that the element volumeis found to be: Ve= Z1 0 Z1r 0 Z1rs 0 det Jdtdsdr = 1 6det J = D✏= DBu (e) In case of isotropic linear elastic material=,1+ν1−2ν1−ννν000ν1−νν000νν1−ν000000212−ν000000212−ν000000212−ν3=24+555000524+550005524+5000000400000040000004=G=#21+ν(=#ν1+ν1−2ν K(e)is a symmetric 12x12 matrix K(e)= VeBTDB Force termBody force: Z e⌘TbdV b(e)= Ve 4[bx by bz ··· bx by bz]T Traction force: T(e)= Ae 3[Tx Ty Tz Tx Ty Tz Tx Ty Tz 000 ]T Z Ae ⌘TTdA Face 1Face 4 T(e)= Ae 3[Tx Ty Tz 000 Tx Ty Tz Tx Ty Tz]T IW e = Z e(⌘(e))TBTDBu (e)dV 9Linear Tetrahedral ElementElement stiffnessTrilinear Hexahedron local numberingHexahedral element faces:Face 1:1, 2, 3, 4Face 2:5, 8, 7, 6Face 3:1, 5, 6, 2Face 4:2, 6, 7, 3Face 5:3, 7, 8, 4Face 6:4, 8, 5, 14xy123z5678Cartesian System(-1,-1,-1) 1(1,-1,-1) 23 (1,1,-1)4 (-1,1,-1)(1,-1,1) 67 (1,1,1)8 (-1,1,1)5 (-1,-1,1)strNatural Element System10Trilinear Hexahedral Element (-1,-1,-1) 1(1,-1,-1) 23 (1,1,-1)4 (-1,1,-1)(1,-1,1) 67 (1,1,1)8 (-1,1,1)5 (-1,-1,1)str N1 = 1 8(1 r)(1 s)(1 t) N2 = 1 8(1 + r)(1 s)(1 t) N3 = 1 8(1 + r)(1 + s)(1 t) N4 = 1 8(1 r)(1 + s)(1 t) N5 = 1 8(1 r)(1 s)(1 + t) N6 = 1 8(1 + r)(1 s)(1 + t) N7 = 1 8(1 + r)(1 + s)(1 + t) N8 = 1 8(1 r)(1 + s)(1 + t) Ni = 1 8(1 + rir)(1 + sis)(1 + tit) ri, si, tiare the coordinates of node iin the (r,s,t)system 11Trilinear Hexahedral Element IsoparametricRepresentation. Shape functions u= Nu (e) x= Nx (e) Nis a 3x24 matrix = D✏= DBu (e) ✏= Bu (e) Strain and stressBis a non constant 6x24 matrixFollowing the steps for the quadrilateral element, the strain and stress within the element are expressed asStiffness matrixFollowing the steps for the quadrilateral element, the strain and stress within the element are expressed asK(e)is a symmetric 24x24 matrix K(e)= Z eBTDB dV 12Trilinear Hexahedral Element IsoparametricRepresentation 133D elements in AbaqusNaming convention143D elements in AbaqusLinearQuadraticTe t r a h e d r o nWedgeHexahedronC3D4C3D6C3D8C3D10C3D15C3D20 153D elements in AbaqusIntegration pointsWedge elementsTe t r a h e d r a l e l e m e n t sNote: only the layer of integration points close to 1-2-3 face is shown163D elements in AbaqusIntegration pointsHexahedral elements2x2x2 integration point scheme1x1x1 integration point scheme123745812375846§C3D8§C3D8R(full integration linear element)(reduced integration linear element) 173D elements in AbaqusIntegration pointsHexahedral elements3x3x3 integration point scheme2x2x2 integration point scheme§C3D20§C3D20R(full integration quadratic element)(reduced integration quadratic element)αMMα18Shear lockingIn an ideal situation, a solid under pure bending moment experiences the curved shape shown below and the angle αremains 90 degrees. MMααTo correctly model the ideal shape change, an element should be able to assume the curved shape. However, the edges of a fully integrated linear elementare not able to bend to curves.The top surface experiences tensile stress, the lower compressive stress. All dotted lines remain straight. Angle αdoes not remain 90 degrees.An incorrect artificial shear stress has been introduced. The element becomes locked or overly stiff, resulting in wrong displacements and false stresses. 19Hourglass effectTo address the shear locking: §quadratic elements§linear element with reduced integration (increase of computational efficiency)Reduced-integration linear elements suffers from hourglass effectbecause they tend to be excessively flexible αMMαVertical and horizontal dotted lines, and the angle α remain unchanged. Consequently, normal stresses and shear stresses are zero at the integration point and no strain energy is generated by the deformation. This zero-energy mode is a non-physical response, which may propagate when a coarse mesh is used. The propagation of such a mode may produce wrong results (indicating that the structure is too flexible)203D elements in AbaqusTr i l i n e a r h e x a h e d r a l e l e m e n t s w i t h i n c o m p a t i b l e m o d e s ( C 3 D 8 I )This element is a fully integrated element with added internal degrees of freedom (incompatible deformation modes) eliminating parasitic shear stress, thereby eliminating the shear locking phenomenon. It also removes artificial stiffening due to Poisson’s effect in bending.Pro:§No hourglass effect (because it is fully integrated)§It can capture bending as accurately as 2ndorder elements, with reduced runtimeCons:§The shape of the element has to be nearly perfect. Parallelogram or trapezoidal shapes result in less accurate results Numericalintegration: GaussianquadraturemethodDarioGastaldiComputationalBiomechanicsLaboratory E. è / " E ' ? ?.dk/i/iBI9iHEEG-ii)JH;Ytedgdy...=IIfH;)JG;y)d9o4 - I-1 MetododiGauss :µ, rifletti /Jtgi ; 4) Wi Ii! ! -→ ask.ie il •• (1,11Polinomi ordine ✗× (21-1) > g gg××@ kit Numericalintegration: GaussianquadraturemethodDarioGastaldiComputationalBiomechanicsLaboratory E. è / " E ' ? ?.dk/i/iBI9iHEEG-ii)JH;Ytedgdy...=IIfH;)JG;y)d9o4 - I-1 MetododiGauss :µ, rifletti /Jtgi ; 4) Wi Ii! ! -→ ask.ie il •• (1,11Polinomi ordine ✗× (21-1) > g gg××@ kit •As we have seen, computing the stiffness matrix and load force in isoparametricelements require the use of numerical integration.•The numerical integration used in finite element applications is the Gaussian Quadrature.Gaussian QuadratureThe idea behind the Gaussian Quadrature is the following I= Z1 1f(s)ds = X i !if(⇠i) We want to find wi(weight) and xi(Gauss point) in order to integrate f(s)exactly in the case f(s)is a polynomial. Depending on the order to the polynomial, different quadraturesare defined, i.e., one point, three point quadrature, etc.The basic assumption is that the quadrature so define can be applied for a general form of f(x), not necessarily a polynomial.2Numerical IntegrationConsider the one point quadrature I= Z1 1f(s)ds = !1f(⇠1) We need to determine two parameters w1and x1. The most general polynomial than can be characterized using two parameters is the linear polynomial f(s)= a·s+b Substituting the expression for f(s)in the expression for the quadrature we obtain: 2b= a!1⇠1+b!1 0= a!1⇠1+b(!12) Since this must hold for any first order polynomial !1=2 ,⇠1=0 3Numerical Integration % !•I -yy ] The two point quadrature formula is found as:We need to determine four parameters w1, w2, x1 and x2. The most general polynomial than can be characterized using four parameters is a cubic polynomialSubstituting the expression for f(s)in the expression and collecting all terms multiplying the coefficients a, b and c, the following system of nonlinear equations is defined: I= Z1 1f(s)ds = !1f(⇠1)+ !2f(⇠2) f(s)= as3+bs2+cs +d !1+!2 =2 !1⇠1+!2⇠2 =0 !1⇠21+!2⇠22 = 2 3 !1⇠31+!2⇠32 =0 This system has a unique solution: !1=1 , !2=1 ⇠1= 1p3, ⇠1= 1p3 4Numerical Integration ÷ : - Il's Note that an n-point Gauss Quadrature integrates exactly a 2n-1degree polynomial Number of Points, n Location, ⇠i Weights, !i 1 0.0 2.0 2 ±0.5773502691896257 1.0 30 .0 0.8888888888888888 ±0.7745966692414834 0.5555555555555556 4 ±0.3399810435848563 0.6521451548625461 ±0.8611363115940526 0.3478548451374538 50 .0 0.5688888888888889 ±0.5384693101056831 0.4786286704993665 ±0.9061798459386640 0.2369268850561891 6 ±0.2386191860831969 0.4679139345726910 ±0.6612093864662645 0.3607615730481386 ±0.9324695142031521 0.1713244923791704 70 .0 0.4179591836734694 ±0.4058451513773972 0.3818300505051189 ±0.7415311855993945 0.2797053914892766 ±0.9491079123427585 0.1294849661688697 I= Z1 1f(s)ds = X i !if(⇠i) 5Numerical Integration Multidimensional Integrals: I= Z1 1 Z1 1f(r, s )drds Z1 1 Z1 1f(r, s )drds = Z1 1 X i !if(⇠i,s)ds = X j X i !j!if(⇠i,⇠j) For a two point quadrature: Z1 1 Z1 1f(r, s )drds =1 ·1·f(1/p3,1/p3) + 1 ·1·f(1/p3,1/p3) +1 ·1·f(1/p3,1/p3) + 1 ·1·f(1/p3,1/p3) Extension to three dimensions follows readily6Numerical IntegrationTr i a n g u l a r e l e m e n t s Number of Points, n Location, ri Location, si Weights, !i 1 0.3333333333333333 0.3333333333333333 0.500000000000000 3 0.6666666666666666 0.1666666666666666 0.1666666666666666 0.1666666666666666 0.6666666666666666 0.1666666666666666 0.1666666666666666 0.1666666666666666 0.1666666666666666 4 0.3333333333333333 0.3333333333333333 -0.281250000000000 0.6000000000000000 0.2000000000000000 0.260416666666666 0.2000000000000000 0.6000000000000000 0.260416666666666 0.2000000000000000 0.2000000000000000 0.260416666666666 6 0.8168475729804590 0.0915762135097710 0.054975871827661 0.0915762135097710 0.8168475729804590 0.054975871827661 0.0915762135097710 0.0915762135097710 0.054975871827661 0.1081030181680700 0.4495484909159650 0.111690794839006 0.4495484909159650 0.1081030181680700 0.111690794839006 0.4495484909159650 0.4495484909159650 0.111690794839006 Z ⌦e f(r, s )drds = X i !if(ri,si) 7Numerical Integration Linear Tetrahedron n Location, ri Location, si Location, ti Weights, !i 1 0.2500000000000000 0.2500000000000000 0.2500000000000000 0.1666666666666666 4 0.5854101966249690 0.1381966011250110 0.1381966011250110 0.0416666666666666 0.1381966011250110 0.5854101966249690 0.1381966011250110 0.0416666666666666 0.1381966011250110 0.1381966011250110 0.5854101966249690 0.0416666666666666 0.1381966011250110 0.1381966011250110 0.1381966011250110 0.0416666666666666 5 0.2500000000000000 0.2500000000000000 0.2500000000000000 -0.1333333333333333 0.5000000000000000 0.1666666666666666 0.1666666666666666 0.0750000000000000 0.1666666666666666 0.5000000000000000 0.1666666666666666 0.0750000000000000 0.1666666666666666 0.1666666666666666 0.5000000000000000 0.0750000000000000 0.1666666666666666 0.1666666666666666 0.1666666666666666 0.0750000000000000 Z ⌦e f(r, s, t )drds = X i !if(ri,si,ti) 8Numerical Integration Non linear problems: iterativetechniquesDarioGastaldiComputationalBiomechanicsLaboratory2The actual relationship between loads and displacements (shown with dotted line) is not known beforehand.Consequently, a series of linear approximations with corrections is performed.Newton-Raphson method (shown as solid red lines)In the N-R method the total load Fais applied in iteration 1. The results is x1.From the displacements, the internal forces F1can be calculated.If Fa ≠ F1 , the system is not in equilibrium. Hence, a new stiffness matrix (slope of dotted line) is calculated based on the current conditions.The difference of (Fa -F1) is the out-of-balance or the residual forces. The residual forces must be small enough for the solution to converge.The process is repeated until Fa = Fiwhereiis the ithiteration of the process for which the solution is said to be converged (i=4 in thiexample).Non linear problems: iterative techniquesNewton-Raphson method 3Non linear problems: iterative techniquesNewton-Raphson method4Non linear problems: iterative techniquesNewton-Raphson methodQuadratic asymptotic rate of convergence (!) but some drawbacks... µ(un +. ) = [ +, - Plum)