Finite Element Analysis - 4

Lecture 4 Special Topics of FEA Somenath Mukherjee Gangan Prathap Scientist E2, Structural Technologies Division, National Aerospace Laboratories (NAL), Bangalore. Director, National Institute of Science and Information Resources (NISCAIR), New Delhi. 1 Lecture 4 Special Topics of FEA Chapters 1. Shear Locking in Timoshenko beam elements. 2. Error Analysis in Computational Elastodynamics. 3. Rank Deficiency in elements. 2 Lecture 4 Special Topics of FEA Chapter 1 Shear Locking in Timoshenko beam elements (A Pathological Problem) 3 1.1 The Pathological problem of locking Locking is a pathological problem encountered in formulating a certain class of elements for structural analysis, although these elements satisfy completeness and continuity requirements. • Locking causes slow convergence even for very fine mesh. • Locking is manifested as Spurious Stiffening and Stress Oscillations. Explanations: (1) Locking is caused by ill conditioning of the stiffness matrix due to the very large magnitude of the shear stiffness terms as compared to the those of bending stiffness (Tessler and Hughes). (2) Locking occurs due to coupling between the shear deformation and bending deformation, and that it can be eliminated by appropriate decoupling (Carpenter et al). (3) Elements lock because they inadvertently enforce spurious constraints that arise from inconsistencies in the strains developed from the assumed displacement functions. (Prathap et al). 4 1.2 The Shear-flexible beam (Timoshenko) In the classical Euler beam (meant only for thin beams), it has been shown that despite the presence of shear stress in the beam sections, the shear strain is ignored. The Euler beam is of infinite shear rigidity (!) For thick beams (of wider webs), the Euler beam theory is not valid. Shear deformation of the web requires shear-flexible formulations. Euler beam Prof S P Timoshenko Timosheko beam x dw/dx dw/dx x z,w =dw/dx NA zx= - dw/dx 5 Elementary beam theory as constrained media problem Euler beam model L x 1 dθ δΠ = δ EI 2 dx 0 x z,w =dw/dx NA The Euler beam has infinite shear rigidity But the practice of using a large shear rigidity for thin beams creates a problem called Shear Locking in shear-flexible beam elements. L L 1 dw dx − qwdx + κ θ− 2 dx 0 0 2 =0 Equilibrium Equations EI dw/dx 2 d 2θ dx 2 −κ θ − dw =0 dx dθ d 2 w =q − κ dx dx 2 Combining i.e. ....(i ) d dw =q κ θ− dx dx ....(ii ) (i ) & (ii ) d 2θ d ...(iii ) EI 2 − q = 0 dx dx Boundary conditions at x = 0 & x=L dθ = 0 or δθ = 0 dx dw Either = 0 or δw = 0 κ θ− dx dw As κ → ∞ θ → dx Either EI Equation (iii ) reduces to EI d 4w dx 4 −q →0 6 Equilibrium equations of the Shear flexible (deep) beams L dθ 1 δΠ = δ EI dx 2 0 dw/dx 2 L 1 dw dx + kGA θ − dx 2 0 Equilibrium Equations d 2θ dw/dx V dw EI 2 − kGA θ − =0 dx dx dθ d 2 w kGA − =q dx dx 2 Combining V Shear strain = - dw/dx=V/(kGA) (1.1) Shear rigidity: kGA i.e. L − qwdx = 0 0 -dM/dx-V=0 ....(i ) (1.2) d dw kGA θ − =q dx dx (1.3) ...(iii ) at ....(ii ) dV/dx=q (i ) & (ii ) d d 2θ EI 2 − q = 0 dx dx Boundary conditions dθ Either EI =0 dx dw Either kGA θ − dx 2 x=0 & x=L or δθ = 0 =0 or δw = 0 7 Shear rigidity of deep beams dw/dx γ τ= dw/dx V VQ Ib U shear γ= τ G = VQ GIb L h/2 V A I2 h/2 Q2 dy b −h / 2 L U shear h L V 1 1 V = dx = V γ dx kGA 2 2 0 0 L U shear b 1 1 Q2 dy = kGA G −h / 2 I 2b 1 V V dx = 2 kGA 0 1 k= y h/2 h/2 L U shear L 1 V 2Q 2 τ2 = .τ .γ .dV = . .(bdy)dx = . .(bdy)dx 2 2 2G V 0 y =− h / 2 0 y = − h / 2 2G ( Ib) L γ= L V dw =θ − kGA dx 2 dw 1 1 1 = V γ dx = kGA(γ ) 2 dx = kGA θ − dx 2 2 2 dx 0 0 0 (1.4) k is called the shear correction factor k=5/6 for a rectangular section 8 Example 1. Find the tip deflection of a cantilever subjected to a concentrated tip load P. (Include shear deformation) L P Deflection at the free end : PL3 PL3 P PL3 3EI L= 1+ δ= +γL = + 3EI 3EI kGA 3EI kGAL2 For thin beams, kGAL2 → ∞, EI EI 2 kGAL →0 9 1.3 Formulation of the two-noded Timoshenko Beam Element (Using Linear Lagrangian C0 Shape Functions) Element displacement and geometry (iso-parametric): wh = N1w1 + N 2 w2 θ h = N1θ1 + N 2θ 2 x = N1 x1 + N 2 x2 with N1 = e x1 = 0, x2 = L 1−ξ 2 N2 = 1+ ξ 2 Le x = (ξ + 1), 2 Le dx = dξ 2 θ h = N1 0 N2 0 θ1 0 N1 0 N2 w2 θ2 = [ N ]{δ e } (ε ) = dθ h dx θ h − dwh dx {ε } = [B]{δ } h e = 2x e L −1 (1.5) w2,F2 Le 2,M2 1,M1 ξ = −1 Element Strain vector: h ξ= w1,F1 w1 wh −1 ≤ ξ ≤ 1 ξ =1 − 1/ L 0 0 1/ L {δ e } 1 / L (1 − ξ ) / 2 − 1 / L (1 + ξ ) / 2 (1.6) 10 Element stress resultants : M = V dθ / dx 0 EI = [ D][ B]{δ e } kGA θ − dw / dx 0 (1.7) Element potential energy: L 1 dθ Π= EI 2 dx 0 2 L 1 dw dx + kGA θ − 2 dx 0 ( 2 w1,F1 L w2,F2 − qwdx − {δ e }T R e } Le 0 ) 1 Π = {δ e }T [ K e ]{δ e } − {δ e }T {F e } + {R e } 2 1,M1 Equilibrium [ K e ]{δ e } = {F e } + {R e } δΠ = 0 Element Stiffness matrix 1 Le [ K ] = [ B ] [ D ][ B ] dξ 2 −1 e 2,M2 (1.8) Element Force vector 1 T Le {F } = [ N ] q dξ 2 −1 e (1.9) T (1.10) 11 Using a 2 point Gauss integration the stiffness matrix is 1 Le [ K ] = [ B ] [ D ][ B ] dξ = [ K eb ] + [ K e s ] 2 −1 e T 0 EI 0 [ K e ] = [ K eb ] + [ K e s ] = e L 0 0 1 0 0 0 −1 1 Le / 2 −1 Le / 2 kGA Le / 2 ( Le ) 2 / 3 − Le / 2 ( Le ) 2 / 6 + e 0 0 0 1 −1 − Le / 2 − Le / 2 L 0 −1 0 −1 Le / 2 ( Le ) 2 / 6 − Le / 2 ( Le ) 2 / 3 (1.11) 12 FE results of analysis of deep beam cantilever beam under tip load L E=1000 G=375 b=1, h=1 L=4 Locked results y Observations h b P No of elements Normalized tip Displacement (Locked) 1 0.2 (10)-5 2 0.8 (10)-5 •A pattern in the error. 4 0.32 (10)-4 •Slow convergence 8 0.128 (10)-3 16 0.512 (10)-3 •Large errors 13 Antidote for shear locking. Use a 1 point (instead of 2 point) Gauss integration scheme for the stiffness matrix is 1 Le [ K ] = [ B ] [ D ][ B] dξ = [ K eb ] + [ K e s ] 2 −1 e T 0 EI 0 [ K e ] = [ K eb ] + [ K e s ] = e L 0 0 1 0 0 0 −1 1 Le / 2 −1 Le / 2 kGA Le / 2 ( Le ) 2 / 4 − Le / 2 ( Le ) 2 / 4 + e 0 0 0 1 −1 − Le / 2 − Le / 2 L 0 −1 0 −1 Le / 2 ( Le ) 2 / 4 − Le / 2 ( Le ) 2 / 4 (1.12) 1 Normalized displacement No. of elements 1 2 3 4 5 Magic! An error in the integration eliminates locking! WHY ? 14 FE results of analysis of deep beam cantilever beam under tip load L E=1000 G=375 b=1, h=1 L=4 Normalized tip Displacement (Locked) Normalized tip Displacement (Lock-free, Reduced Int.) 1 0.2 (10)-5 0.75 2 0.8 (10)-5 0.75+0.75/4 =0.9375 Number of elements y h b P 15 Example problems solved using a single Timoshenko beam element Observations: Spurious shear oscillations and bending stiffening for the locked case. 16 1.4 Explanations for the origin of locking (The field-consistency paradigm) q L Linear displacements: θ h = ax wh = bx h dw Shear strain γ = θ − = ax − b dx Bending strain h Rayleigh-Ritz procedure EI a= Element locks L 0 0 0 +κ L3 / 3 − L2 / 2 a − L2 / 2 L b − 3qL2 12 EI + κL2 As κ → ∞, a → 0, dθ h =a dx = 0 qL2 / 2 1.5qL3 qL − b=− 2κ 12 EI + κL2 b→0 dθ h →0 dx The parameter a effects both bending and shear strains. It is17 a spurious constraint that stiffens bending as well as shear strains 1.5 Explanation of shear locking in the element by FieldConsistency Theory The shear strain in the element is θ h - dwh dx = α + βξ where α = (θ + θ ) 2 − (w − w ) L and β = (θ − θ ) 2 2 1 2 1 2 1 • For thin beams, the shear strain energy term vanishes, leading to two α→0 β→0 constraints: (First constraint is physically meaningful in terms of the equivalent Euler beam model, but the second constraint is a spurious one. The spurious term β effectively enhances the element's bending stiffness to EI*=EI+kGA(Le)2 /12, where EI and kGA are the bending and shear rigidities respectively of the actual beam, leading to locking. wLF wL = I * I = 1 + kGAL2 (12EI ) = 1 + e e=kGA(Le)2/(12EI)=K/n2 , (l=total beam length, n=total number of equal elements, Le= element length=l/n). The parameter e becomes larger for thinner beams, leading to spuriously high bending stiffness, and spurious shear strain oscillations in the elements. 18 1.6 How shear locking is eliminated by reduced integration The integrand in the element stiffness matrix [Ke] is quadratic, so we need a 2 point Gauss rule for exact integration. This element suffers shear locking. 1 Le [ K ] = [ B ] [ D ][ B] dξ = [ K e b ] + [ K e s ] 2 −1 e T 0 0 0 EI 0 [K ] = [K b ] + [K s ] = e L 0 1 0 −1 0 0 e e e 0 0 0 −1 0 −1 1 + Le / 2 −1 Le / 2 kGA Le / 2 ( Le ) 2 / 3 − Le / 2 ( Le ) 2 / 6 Le −1 − Le / 2 1 − Le / 2 Le / 2 ( Le ) 2 / 6 − Le / 2 ( Le ) 2 / 3 A reduced integration actually eliminates (ignores) the spurious term of the shear strain (associated with linear variation in ) so that only constant terms are needed to be integrated. This elimination of the spurious constraint is done by a 1 point Gaussian rule for integration. 0 0 EI 0 1 [ K e ] = [ K eb ] + [ K e s ] = e L 0 0 0 −1 −1 1 Le / 2 Le / 2 kGA Le / 2 ( Le ) 2 / 4 − Le / 2 ( Le ) 2 / 4 + e 0 0 −1 − Le / 2 − Le / 2 L 1 0 −1 Le / 2 ( Le ) 2 / 4 − Le / 2 ( Le ) 2 19 /4 0 0 0 −1 If one uses a Reduced Integration scheme with a one-point rule of Gauss Quadrature, instead of the two-point rule necessary for accurate integration in the shear strain energy, it leads to · Elimination of shear locking by releasing the stiffening constraint β. · Elimination of spurious shear stress oscillations. Reduced integration effectively drops the Second Legendre Polynomial from the shear strain, α + βξ → α 20 The Function Space Approach to Locking Problems 21 1.7 Definition of the Inner product The inner product for the Timoshenko beam element is defined through the symmetric bilinear forms: h {ε } h T e a(u , u ) = EI 0 0 kGA e 1 = {ε } EI h T 0 −1 e = 1 −1 {ε } h T Le {ε } dξ =< ε h , ε > kGA 2 0 {ε } h T a (u h , u h )e = EI 0 {ε }dx EI 0 0 {ε h }dx kGA e 0 h L {ε } dξ =< ε h , ε h > kGA 2 EI=Flexural Rigidity, (1.13) kGA=Shear Rigidity 22 1.8 The B Subspace The B subspace is the space in which the column vectors of the strain-displacement matrix [B] lie. [B ] = − 1/ L 0 0 1/ L (1.14) 1 / L (1 − ξ ) / 2 − 1 / L (1 + ξ ) / 2 The Gram-Schmidt Algorithm for getting the orthogonal basis vectors spanning the B Space: {v1} = {b1} k bk +1, v j j =1 vj ,vj {vk +1} = {bk +1} − {v j } (1.15)) After scaling, only TWO NON-ZERO orthogonal basis vectors are obtained that span the B Space (of 2 dimensions, m=N-R=4-2=2) {v1} = B ⊂ Pnr==22 ; P22 0 1 , {v2 } = 2 = {p} : {p} = 2/ L (1.16) ξ {ai }ξi −1 ,−1 ≤ ξ ≤ 1 , { ai } ∈ R 2 i =1 dim( B) = 2 < dim Pnr==22 = 2 × 2 = 4 (1.17) 23 1.9 Strain projections on the B Subspace; Shear Locking Orthogonal Projection of the Analytical Strain onto the B Subspace yields the FEA computed element strains (best-fits). {ε } = {ε } = h < ε ,vj > {v j }, j =1 < v , v > j j m=2 < v1 , v2 >= 0 (1.18) However, we have problems for thin beams: 1. 2. 3. The bending strain is a lot smaller than the analytical one, showing that spurious bending stiffness has been introduced through FEA . There is spurious shear strain oscillation in FEA results. Slow Convergence even with many elements. These are the symptoms of locking 24 Locked FEA solutions agree with the best-fit strain vector at the element level. Thus locked solutions are variationally correct ε h = ε = Best − fit < ε ,vj > {ε } = {ε } = j =1 {v j }, < vj ,vj > m=2 h < v1 , v2 >= 0 A best fit satisfies the Projection Theorem (Pythagoras) ε −ε Thus i.e. ε −ε 2 h 2 = ε 2 = ε 2 −ε −ε 2 h 2 The Energy of the Error= Error of the Energies 25 26 TABLE 1 Analytical strains and their locked projections as finite element strains e=kGAL2/(12EI). Cantilever with tip moment M0 Analytical strain vector Locked strain vector {ε } = {ε } = h {ε } = {ε } = ( M 0 / EI ) /( 1 + e ) {ε } = PL( 1 + ξ ) /( 2 EI ) P / kGA 0 − < ε ,vj > {v j }, j =1 < v , v > j j m=2 M 0 / EI Cantilever with tip load P ( PL / 2 EI ) /( 1 + e ) − 6 e M 0ξ ( 1 + e ) LkGA < v1 , v2 >= 0 {ε } = {v1} = 3eξ P ) (1+ 1+ e kGA 0 1 , {v2 } = 2/ L ξ 27 TABLE 2 Error norm square for locked strain projections with the linear two noded Timoshenko beam element. e=kGAL 2/(12EI) 1 q 2 L = {q}T [D] {q} dξ 2 −1 Case Cantilever with tip moment, Mo Cantilever with tip transverse load P {q} = {ε } − {ε } Locked Solution 2 L 2M 0 e . 2 EI 1 + e L ( PL ) . 2 2EI 2 e 1 + 1+ e 3 28 1.10 The Function Space explanation of locking and its elimination The original field-inconsistent [B] matrix is [B ] = 0 − 1/ L 0 1/ L 1 / L (1 − ξ ) / 2 − 1 / L (1 + ξ ) / 2 Locking occurs because the 2-dimensional B subspace is field-inconsistent, which cannot be spanned by the standard basis vectors of its 4-dimensional parent space P 22 (linear in ξ), {L1} = [0, 1]T , {L2} = [1, 0]T, {L3} = [0, ξ]T, {L4} = [ξ, 0]T., (1.19) Actually, the field-inconsistent B space is spanned by non-standard basis vectors, {v } = 0 , {v } = 2 / L 1 1 2 ξ 29 1.11 Elimination of shear locking Reduced Integration effectively sets the highest order Legendre Polynomial ξ in the [B] matrix to zero. It replaces [B] by a (modified) [B*]. Lock-free strain vector is expressed as, (1.20) {ε h *} = 0 − 1 / L 0 1/ L {δ *e } = [B * ] {δ *e } 1/ L 1/ 2 1/ L 1/ 2 A new field-consistent space B* emerges from [B*]. This lockfree, field-consistent space B* is two-dimensional, and can be spanned by the standard orthogonal basis vectors, 1 0 {v1*} = , {v2 *} = 0 1 (1.21) B* ⊂ Pnr==12 ; P12 = {{p} : {p} = { ai } , { ai } ∈ R 2 } dim( B) = 2 = dim Pnr = 2 × 1 = 2 (1.22) 30 Lockfree stiffness matrix for the Timoshenko beam is obtained from the field-consistent (lockfree) strain-displacement matrix [B*] with exact integration 1 e L e e [ K ] = [ B*] [ D][ B*] dξ = [ K b ] + [ K s ] 2 −1 e T 0 0 0 0 Le / 2 Le / 2 −1 1 kGA Le / 2 ( Le ) 2 / 4 − Le / 2 ( Le ) 2 / 4 EI 0 1 0 − 1 + e [K ] = [K b ] + [K s ] = e 0 0 0 0 L L 0 −1 0 −1 e e e −1 − Le / 2 1 − Le / 2 Le / 2 ( Le ) 2 / 4 − Le / 2 ( Le ) 2 / 4 (1.23) 31 1.12 Orthogonal Projection on B* space In general, Reduced Integrated FEA results are NOT variationally correct. (RI is a variational crime !) Reduced Integrated FEA strains will agree with the best-fit solution, provided the following rule holds good, {F e E }= − [[B] − [B *]]T [D]{ε }dx = 0 (1.24) e Then: {ε *} = {ε *} = h m i =1 < ε ,ν i * > {ν i *}, < ν i *,ν j * >= 0 < ν i *,ν i * > for i≠ j (1.25) When this extraneous force {F e E }does not vanish, then the best-fit solution (on the B* space) will suffer additional strain from this extraneous force vector, over the lockfree (reduced integrated) FEA solution. 32 TABLE 3 Analytical strains and their locked and lockfree projections as finite element strains e=kGAL2/(12EI). Cantilever with tip moment M0 Analytical strain vector M 0 / EI {ε } = 0 − M 0 / EI { ε *} = 0 < ε ,vj > {v j }, j =1 < v , v > j j m=2 < ε,v j > j =1 < v j , v j 6 e M 0ξ ( 1 + e ) LkGA − m=2 {ε *}= {ε *}= h {ε } = P / kGA > {v * j }, ( PL / 2 EI ) /( 1 + e ) − {ε } = Lockfree strain vector h PL( 1 + ξ ) /( 2 EI ) ( M 0 / EI ) /( 1 + e ) Locked strain vector {ε } = {ε }= Cantilever with tip load P < v1 , v2 >= 0 v1 v2 {ε } = 3eξ P ) (1+ 1+ e kGA PL /( 2 EI ) − { ε *} = {v1} = 0 1 {v1*} = P / kGA {v2 } = , 0 1 , 2/ L ξ {v2 *} = 1 0 B B* 33 TABLE 4 Error norm square for strain projections with the linear two noded Timoshenko beam element. e=kGAL 2/(12EI) 1 L q = {q}T [D] {q} dξ 2 −1 2 Case Cantilever with tip moment, Mo Cantilever with tip transverse load P {q} = {ε }− {ε } Locked Solution 2 L 2M 0 e . 2 EI 1 + e L ( PL ) . 2 2EI 2 e 1 + 1+ e 3 Lockfree Solution 0 L 2 (. PL )2 6EI 34 A case of variational incorrectness through reduced integration A cantilever beam with uniformly distributed loading ρ FI: Field inconsistent, Locked, but variationally correct FE results. FC: Field consistent, Lock free, Reduced Integrated FE results. Note that FC (by FEA) deviates from the field-consistent best-fit results. For this case : The extraneous force vector (a non-zero vector) from Reduced Integration consists of self-equilibrating moments, that shift the FC Best-fit from the FC-FEM results. {F }= − e [[B ] − [B *]]T [D]{ε }dx E e 0 {FEe } = ρL2 / 12 0 − ρL2 / 12 {ε *}= {ε *}+ δε * h 35 1.13 Lockfree an-isoparametric formulation h Displacements: 3 w = _ h θ = N i wi i =1 {ε } h Strain vector: {ε }= h 2 _ (1.26) N iθ i i =1 dθ dx = [B ] δ e = θ − dw / dx { } 0 −1/ L 0 1/ L 0 δe − (2ξ − 1) / L (1 − ξ ) / 2 − (2ξ + 1) / L (1 + ξ ) / 2 4ξ / L { } (1.27) Standard basis vectors spanning 3-dimensional B space: {v } = [0 , ξ ] , {v } = [1 , 0 ] and {v } = [0 , 1] (1.28) T 1 T 2 T 3 36 Summary Shear locking in Timoshenko’s Shear Flexible beam element occurs from spurious constraints that arise from reducing the discretized domain into an Euler beam (of infinite shear rigidity). Shear locking is displayed through slow convergence, Spurious bending stiffening and shear oscillations. The field consistency paradigm identifies the spurious constraints related to locking, and suggests methods to eliminate field inconsistency by eliminating the spurious constraints (thereby enforcing field consistency). Reduced integration (RI) eliminates shear locking by eliminating the spurious constraint in the strain. The function space approach shows that locked strain vector in an element (through FEA) is actually the orthogonal projection of the analytical strain vector onto a fieldinconsistent subspace B, arising from a field-inconsistent [B] matrix (straindisplacement matrix). B cannot be spanned by standard orthogonal basis vectors. FEA through reduced integration (RI) effectively projects the analytical strain vector onto a field-consistent subspace B*. However, RI is variationally incorrect in general, and the FE strain vector agrees with the orthogonal projection on B* only when the spurious extraneous force vector vanishes. 37 Lecture 4 Special Topics of FEA Chapter 2 Error Analysis in Computational Elastodynamics A comedy of errors… 38 2.1 Finite Element Elastodynamic Equations using the Principle of Least Action 2 Action Lagrangian I = L(q,q, t )dt L = T −V 1 Hamilton’s Principle δI = 0 for q(t), q(t1)= q(t2)=0 Lagrange’s Equation for motion d ∂L ∂L − = Qi dt ∂q i ∂qi Qi = Non − conservative generalised force 39 In elastodynamics, the equations of motion are generally derived in a global sense (with element assembly) N L = T −V = e N T − e =1 e N U − e=1 W e =1 e N 1 e T = .{δ } [ M e ]{δ e } − e =1 2 Element Stiffness Matrix: N e T e e N {δ } [ K ]{δ } − e =1 {δ e }T {F e } e =1 [ K e ] = [ B ]T [ D][ B]dV e Element Consistent Mass Matrix: [ M e ] = [ N ]T [ ρ ][ N ]dV e Element Generalized Force Vector (time dependent) : {F e } = [ N ]T { f (t )}dV e With element assembly, we get the global form 1 1 L = T − V = {δ G }T [ M G ]{δ G } − {δ G }T [ K G ]{δ G } − {δ G }T {F G } 2 2 Equation of motion d ∂L ∂L − =0 dt ∂q i ∂qi [ M G ]{δ G } + [ K G ]{δ G } = {F G } (2.1) 40 2.2 Free Vibration Analysis [ M G ]{δ G } + [ K G ]{δ G } = {0} Let {[K G {δ } = {φ}.sin(ω t ) G n } ] − ωn 2 [ M G ] {δ G } = 0 { (2.2) (2.3) (2.4) } det [ K G ] − ωn 2 [ M G ] = 0 Eigenvalue ωi is ωi 2 , Eigenmode {φi }, natural circular frequency (rad / sec) 41 Orthogonality of the Eigen-modes (Normal modes) {φi }T [ K G ]{φ j } = 0 {φi }T [M G ]{φi } = 0 i ≠ j, i ≠ j, {φi }T [ K G ]{φi } = kii {φi }T [M G ]{φi } = mii (2.5) Natural Frequencies (rad/sec) kii ωi = mii kii : generalized modal stiffness for mode i mii:: generalized modal mass for mode i kii /\/\/\/\/\ mii qi(t) Qi(t) 42 Example 1. Free vibration analysis of a simple cantilever beam using 10 Euler beam elements. L=1m, b=0.1m, t=0.001m I=2.5×10-7m4, A=3×10-4m2 Density ρ=2722.77 kg/m3, Mass per unit length of the beam is ρA=0.816 kg/m E=7.1 ×1010 N/m2, m Table 3.1Comparison of the natural frequencies in bending of the uniform cantilever beam obtained by different methods Different methods Classical solution Natural Frequenc y fn (Hz) 82.4915 Natural Circular Frequency Different methods n (rad/sec) 518.31 Classical solution Natural Frequency fn (Hz) Natural Circular Frequency n (rad/sec) 516.935 3.248×103 Differen t methods Natural Frequency fn (Hz) Natural Circular Frequency n (rad/sec) Classical solution 1447.51 9.095×103 9.097×103 FE result 82.4836 518.26 FE result 516.935 3.248×103 FE result 1447.83 43 Example 2. Free vibration analysis of an aircraft wing using Euler beam elements. 7.105 m Leading edge 2.1469m Trailing edge Fig (a) A typical subsonic aircraft wing 0.849 m Tip Root Fig (b) Length of the wing divided in 22 elements for FE formulation 44 Example 3. Dynamic Characterization of an aircraft using a Stick Model Stick Stick Model Model ofofSARAS SARASfor forDynamic Dynamic Characterization Characterization •Beam Model with provision for Bending-Torsion Coupling (Shear Center offset). •Results for components from in-house code benchmarked with those from NASTRAN. Mode Results from Beam Model Results from GVT Results from 3D Model % Error (Betwee n Beam & GVT) Rudder Rotation 6.93 7.09 7.08 -2.25 Wing Symm. Bending 6.01 5.52 6.76 8.8 Wing 1st Anti-symm. Bending 7.69 6.85 5.36 6.27 Elevator Rotation 11.82 11.02 11.03 7.2 VT long. Bending 17.83 16.0 14.75 11.43 HT 1st Anti-symm. Bending 10.87 11.53 10.4 -5.72 VT lateral Bending 19.14 17.37 19.09 10.18 HT 1st Symm. Bending Aileron Symm. Rotation 21.39 21.68 20.25 1.43 22.78 24.0 22.18 -1.36 23.75 24.0 25.61 -1.04 1st Fuselage Lateral Bending Fuselage Longitudinal Bending Wing First Symmetric Mode 6.71 Hz (Stick) and 6.72 Hz Wing Second Symmetric Mode 18.89 Hz (Stick) and 19.51 Hz 45 30.29 30.22 31.69 7.00 HT Anti-Symmetric Mode 10.4 Hz (Stick) and 9.1 Hz 2.3 Definitions of Inner Products in Elastodynamics [D]= element elastic rigidity matrix [ ρ ]= element inertia density matrix. Stiffness-inner product < a, b >= Ne ele=1 ele {a}T [ D]{b}dx (2.6) Stiffness-norm squared value of the vector {a} is given as 2 a =< a, a > Inertia-inner product ( c, d ) = Ne ele =1 ele {c}T [ ρ ]{d }dx (2.7) (2.8) Inertia-norm squared value of the vector {c} is given as 2 c = (c, c ) (2.9) 46 2.4 The Rayleigh Quotient Free vibration of a system in a given mode can be expressed as {U ( x, t )} = {u ( x)}eiωt (2.10) Rayleigh Quotient from exact solution 2 ω = 2 ε u (2.11) 2 Let u and ε be the approximate modal vector and the strain vector. Rayleigh Quotient from FEA solution 2 ω = ε u 2 2 (2.12) But interestingly (!) for a variationally correct solution, < ε, ε > ω = (u , u ) 2 (2.13) 47 2.5 The Error Statements of Elastodynamics Combining equations (5.6) and (5.7), we get one rule 2 2 2 2 ε − ε = ω2 u −ω u (2.14) Error of global strain energy = Error of global kinetic energy Combining equations (5.7) and (5.8), we get another rule, valid for variationally correct solutions only, 2 < ε , ε − ε >= (u , ω 2u − ω u ) (2.15) Observation: The Errors in Elastodynamics are decided by both displacements and strains. 48 2.6 The Frequency-Error Hyperboloid It can be shown that for variationally correct formulations, u −u u 2 2 2 + ε −ε 2 ω − =1 ω2 ω2 u 2 (2.16) Z = ε −ε H F 1 A E (ω(ω == ωω) ) Y= ω ω X = u −u Fig 1. Geometric interpretation of eigenvalue analysis of the variationally correct formulation using Frequency-Error Hyperboloid. Approximate eigenvalues obtained form a variationally correct formulation lie in the shaded portion of the Hyperboloid. 49 2.7 Why is the approximate Rayleigh Quotient higher than the Exact one ? (Valid only for variationally correct formulations) ε u θu θε u ε Fig 2. Included angles (u , u ) cos(θ u ) = , uu cos(θε ) = cos θ u (u , u ) ε ε 1 ε ε = = cos θ ε u u < ε ,ε > ω 2 u u < ε ,ε > ε ε Hence for variationally correct formulations 2 2 2 cos θ u 1 ε ε 1 2 2 ω = 4 2 2 = 4ω ω = 2 2 cos θε ω u u ω ω 2 (2.17) Geometrically, the modal displacement vector suffers less deviation 2 than that of the modal strain vector. Hence ω >1 ω2 • 50 Example 4: Free Vibration of a Simply Supported Beam L x Analytical Modal Disp. Modal Strain Approximate w = a sin (πx / L ) ε = (− d w / dx ) = a 2 Eigenvalue 2 w=b 2 π d2w 2b ε= − 2 = 2 dx L sin (πx / L ) L ω 2 = π 4 EI /( ρAL4 ) 2 ω = 120 * EI /( ρAL4 ) EI ε − ε = ω u −ω u = 3 L 2 2 2 2 2 x x 1− L L 2 2 < ε , ε − ε >= (u , ω 2u − ω u ) = 4 π4 2 a 2 − 4b 2 EI 2 ( π ab − b ) L3 51 2.8 Replacement of Consistent Mass by Lumped Mass; A variational crime L eigen value Example 5. Free Vibration Analysis of a Fixed-Fixed Bar using 2 elements (First Mode Only shown) 36 33 30 27 24 21 18 15 12 9 6 3 0 αL (1-α)L consistent mass A B exact solution lumped mass 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 element length ratio Any variationally incorrect formulation (with Lumped Mass, Reduced Integration etc.) that does not conform exactly to the Weak form is variationally incorrect. Variationally incorrect formualtions - Do not satisfy the Hyperboloid Rule - Cannot guarantee and upper bound of the frequency. 52 Lecture 4 Special Topics of FEA Chapter 3 Rank Deficiency in elements 53 3.1 What is rank deficiency ? The rank of the stiffness matrix is the dimension of the B subspace that emerges from the strain-displacement matrix [B], i.e. Rank [Ke]= dim (B) (3.1) In the dimension of the B subspace is given by dim (B)=N-R (3.2) N= Number of degrees of freedom of the element R= Number of rigid body motions To eliminate locking, a reduced order integration effectively coverts the Field-inconsistent [B] matrix into a Field-consistent [B*] matrix, by simply removing the highest Legendre Polynomial in the fieldinconsietnt spurious term of [B]. Using Gram Schmidt algorithm for orthogonal basis vector spanning B* it can be shown that for some elements dim (B*) < dim (B) or (N-R*) < (N-R), i.e. R* > R (3.3) Rank [Ke*] < Rank [Ke] because of introduction of spurious rigid body motions Reduced integration may introduce rank deficiency 54 Rank deficiency of the plane stress Quad 4 element 1,2,3 are rigid body modes 4,5,6 are constant strain modes 7,8 are bending modes, but cannot be sensed by a 1x1 reduced integration 1x1 reduced integration (with sampling point at element center of zero strain) effectively considers modes 7 and 8 as zero energy hour-glass modes (spurious rigid body motions) Rank deficiency of this plane stress Quad 4 element is thus 2. 55 56 57 The Mindlin Plate element 58 Appendix The basis vectors for spanning the B subspaces of the elements 59 60 61 62 63 Some Thoughts A burning question: Does Mesh Optimisation Maximize Numerical Entropy? Total Strain Energy Analytical Strain Energy with the analytical solution u remains Invariant (Maximum entropy) Optimized mesh corresponds to maximum FEA Strain Energy (with highest possible entropy with the approximations made). A FEA Strain Energy with the approximate solution uh depends on meshing (the position of the middle node). Lower entropy than at A. 2 ε ≥ εh 2 Position of Middle Node. 64 Cui bono ? (For whose good ?) How the best-fit paradigm helps (a) Gives the exact, but hidden, mechanism of the way the Finite Element Method works. It shows that computations in FEM are actually determined in a best-fit manner of the strains (and stresses), instead of the existing myth that they are based on displacements. (b) Helps one to make a priori error estimates for bench mark problems easily. ( c) Helps one to evaluate the quality of the element that he/she develops. The origins of the pathological problems of elements can now be understood, diagnosed and eliminated by appropriate methods. 65 When Arts and Science met at the crossroads… An extract from “ Sanchaita ” by Rabindranath Tagore. A 900 A−A v u A P Subspace 66 Bibliography Variational Principles of Classical & Computational Solid Mechanics and Linear Algebra 1. J Taylor. Classical Mechanics. University Science Books. 2. L Meirovitch. Analytical Mechanics. Dover. 3. C Lanczos. The Variational Principles of Mechanics. Dover. 4. C L Dyms, I H Shames. Solid Mechanics- A Variational Approach. McGraw Hill. 5 T H Richards. Energy Methods in Stress Analysis. Ellis Horwood. 6 S Rajashekaran, G Sankarasubramanian. Computational Structural Mechanics. Prentice Hall of India. 7. R Solecki, J R Conant. Advanced Mechanics of Materials. Oxford University Press. 8. L H Edwards, D E Penny. Elementary Linear Algebra. Prentice Hall. 9. G Strang. Linear Algebra. Thomson Brooks/Cole. 67 Bibliography (Finite Element Method) Elementary FEM 1. T R Chandrapatla and A D Belegundu. Finite Elements in Engineering. Eastern Economy Edition. 2. R D Cook, D S Malkus, M E Plesha. Concepts and Applications of Finite Element Analysis. John Wiley & Sons. 3. J N Reddy. An introduction to the Finite Element Method. McGraw Hill. 4. K H Huebner, D L Dewherst, D E Smith, T G Byrom. The Finite Element Method for engineers. John Wiley & sons. 5. O C Zienkiewicz. The Finite Element Method. McGraw Hill Advanced FEM 6. G Strang, G J Fix. An analysis of the Finite element Method. Prentice Hall, NJ. . 7. G Prathap. The Finite Element Method in Structural Mechanics. Kluwer Academic Press, Dordrecht. 8. T J R Hughes. The Finite Element Method. Dover. 68 Recent Publications from NAL 1. S. Mukherjee and G. Prathap 2001 17 (6), pp 385-393. Communications in Numerical Methods in Engineering. Analysis of Shear Locking in Timoshenko beam elements using a function space approach. 2. S. Mukherjee and G. Prathap 2002 Sadhana.27(5) 507-526. Analysis of delayed convergence in the three noded isoparametric Timoshenko beam element using the function space approach. 3. G. Prathap and S. Mukherjee 2003 Current Science, 85(17), pp 989-994. The engineer grapples with Theorem 1.1 and Lemma 6.3 of Strang and Fix. 4. H. Mishra and S. Mukherjee 2004 Sadhana 29(6), pp 573-588. Examining the best-fit paradigm in FEM at element level. 5. S. Mukherjee, P. Jafarali and G. Prathap 2005 Journal of Sound and Vibration. Vol 285(3), pp 615-635. A variational basis for error analysis in finite element elastodynamics. 6. K. Sangeeta , Somenath Mukherjee, and Gangan Prathap 2005 Structural Engineering and Mechanics Vol 21(5), pp 539-551. A function space approach to study rank deficiency and spurious modes in finite elements. 7. K. Sangeeta, Somenath Mukherjee and Gangan Prathap 2006 International Journal for Computational Methods in Engineering Science & Mechanics, Vol 7, pages 1-12. Conservation of the Best-Fit Paradigm at Element Level. 8. P. Jafarali, M Ammen, S. Mukherjee, G. Prathap, 2007 Journal of Sound and Vibration, Vol 299(2), pp 196-211. Variational Correctness in Timoshenko beam finite element elastodynamics. 9. S. Mukherjee, P. Jafarali, (Accepted in October 2008 by Communications in Numerical Methods in Engineering; available on-line). Prathap’s best-fit paradigm and optimal strain recovery points in indeterminate tapered bar analysis using linear element. 69 Thank You The Blind Men and the Elephant And so these men of Indostan, Disputed loud and long, Each in his own opinion, Exceeding stiff and strong, Though each was partly in the right, And all were in the wrong! - John Godfrey Saxe (1816-1887) 70