An Introduction to the Finite Element Method

AE1007 Finite Element Method Unit - I An Introduction to the Finite Element Method The Finite Element Method Defined The Finite Element Method (FEM) is a weighted residual method that uses compactly-supported basis functions. Brief Comparison with Other Methods Finite Difference (FD) Method: FD approximates an operator (e.g., the derivative) and solves a problem on a set of points (the grid) Finite Element (FE) Method: FE uses exact operators but approximates the solution basis functions. Also, FE solves a problem on the interiors of grid cells (and optionally on the gridpoints as well). Brief Comparison with Other Methods Spectral Methods: Finite Element (FE) Method: FE methods use compact basis functions to approximate a solution on individual elements. Spectral methods use global basis functions to approximate a solution across the entire domain. Overview of the Finite Element Method S   W   G  M  Strong form Weak form Galerkin approx. Matrix form Sample Problem Axial deformation of a bar subjected to a uniform load (1-D Poisson equation) L d 2u EA 2 = p0 dx u 0  = 0 du EA dx =0 xL px  = p0 x = 0, L u = axial displacement E=Young’s modulus = 1 A=Cross-sectional area = 1 Strong Form The set of governing PDE’s, with boundary conditions, is called the “strong form” of the problem. Hence, our strong form is (Poisson equation in 1-D): d 2u = p0 2 dx u 0  = 0 du =0 dx x  L Weak Form We now reformulate the problem into the weak form. The weak form is a variational statement of the problem in which we integrate against a test function. The choice of test function is up to us. This has the effect of relaxing the problem; instead of finding an exact solution everywhere, we are finding a solution that satisfies the strong form on average over the domain. Weak Form d 2u = p0 2 dx d 2u  p0 = 0 2 dx L  d 2u    dx 2  p0 vdx = 0    0 Strong Form Residual R=0 Weak Form v is our test function We will choose the test function later. Weak Form Why is it “weak”? It is a weaker statement of the problem. A solution of the strong form will also satisfy the weak form, but not vice versa. Analogous to “weak” and “strong” convergence: strong: lim xn  x n  weak : lim f xn  f x f n  Weak Form Choosing the test function: We can choose any v we want, so let's choose v such that it satisfies homogeneous boundary conditions wherever the actual solution satisfies Dirichlet boundary conditions. We’ll see why this helps us, and later will do it with more mathematical rigor. So in our example, u(0)=0 so let v(0)=0. Weak Form Returning to the weak form:  d 2u    dx 2  p0 vdx = 0    0 L L d 2u  dx 2 vdx = 0 p0vdx 0 L Integrate Integrate LHS by parts: du dv du     dx  v( x)  dx dx dx  x 0  0 L xL du dv L du   dx v dx dx dx 0 L 0du v xL dx x 0 Weak Form Recall the boundary conditions on u and v: u 0 = 0 du =0 dx x  L v(0)  0 Hence, H du dv du  dx vL  dx dx dx 0 L du dv  dx dx dx  0 p0vdx 0 L L du  v0  dx xL x 0 The weak form satisfies Neumann conditions automatically! Weak Form Why is it “variational”? L du dv  dx dx dx = 0 p0vdx 0 L Variational statement: 1 Find u  H 1 such thatBu, v   F (v) v  H 0 B a bilinear functional F a linear functional , u and v are functions from an infinite-dimensional function space H Galerkin’s Method We still haven’t done the “finite element method” yet, we have just restated the problem in the weak formulation. So what makes it “finite elements”? Solving the problem locally on elements Finite-dimensional approximation to an infinitespace → Galerkin’s Method dimensional Galerkin’s Method Choose finite basis   Then, u  x    c j j  x , j 1 N N N i i i c j unkowns to solve for b j arbitraril y chosen v x    b j j  x , j 1 Insert these into our weak form : L du dv 0 dx dx dx  0 p0vdx N N L N L d j d i 0  c j dx x  bi dx x dx  0 p0  b j j x dx j 1 i 1 j 1 L Galerkin’s Method  L N 0 N L d i  c j dx x  bi dx x dx  0 p0  b j j x dx j 1 i 1 j 1 N d j Rearranging : b c  i 1 i j 1 j N N L 0 N L d j d i dx  bi  p0i dx 0 dx dx i 1 Cancelling : c  j 1 j N L 0 L d j d i dx   p0 i dx 0 dx dx Galerkin’s Method c  j 1 j N L 0 L d j d i dx   p0 i dx 0 dx dx We now havea matrix problemKc  F, where c j is a vectorof unknowns, K ij   L 0 d j d i dx, dx dx L 0 and Fi   p0 i dx We can already see K ij will be symmetric since we can interchang i, j without effect. e Galerkin’s Method So what have we done so far? 1) Reformulated the problem in the weak form. 2) Chosen a finite-dimensional approximation to the solution. Recall weak form written in terms of residual:  L 0 L L  d 2u   2  p0 vdx   Rvdx   bi  Ri dx  0  dx  0 0 i   This is an L2 inner-product. Therefore, the residual is orthogonal to our space of basis functions. “Orthogonality Condition” Orthogonality Condition  L 0 L L  d 2u   2  p0 vdx   Rvdx   bi  Ri dx  0  dx  0 0 i   The residual is orthogonal to our space of basis functions: u H φi Hh uh Therefore, given some space of approximate functions Hh, we are finding uh that is closest (as measured by the L2 inner product) to the actual solution u. Discretization and Basis Functions Let’s continue with our sample problem. Now we discretize our domain. For this example, we will discretize x=[0, L] into 2 “elements”. Ω1 0 h Ω2 2h=L In 1-D, elements are segments. In 2-D, they are triangles, tetrads, etc. In 3-D, they are solids, such as tetrahedra. We will solve the Galerkin problem on each element. Discretization and Basis Functions For a set of basis functions, we can choose anything. For simplicity here, we choose piecewise linear “hat functions”. Our solution will be a linear combination of these functions. φ1 φ2 φ3 x1=0 x2=L/2 x3=L Basis functionssatisfy: i x j    ij  Our solution w be ill interpolat ory. Also,theysatisfy th partition of unity. e Discretization and Basis Functions To save time, we can throw out φ1 a priori because, since in this example u(0)=0, we know that the coefficent c1 must be 0. φ2 φ3 x1=0 x2=L/2 x3=L Basis Functions  2x L   2x  2  2  L  0    2x  1 3   L 0  if x  0, L  2 if x   L , L  2 otherwise x1=0 x2=L/2 x3=L φ2 φ3 if x   L , L  2 otherwise Matrix Formulation Given our matrix problemKc  F : L d  j d i p  dx  c j 0 dx dx dx  0  Kc  F     0 i j 1   N L c K F We can insert the i chosenon thepreviousslide and arrive at a linear algebra problem. Differentiating the basis functions,then evaluatingthe integrals,we have: p0  1  1  4  2 2 K  , F   L  2 2  L 1  4 In a computercode,differentiating the basis functionscan be donein advance,since the basis functionsare known,and integration would be performed numerically by quadrature . It is standardin FEM to use Gaussian quadrature sinceit is exact , for polynomial s. Notice K is symmetric as expected. Remarks on Variational Problem L du dv  dx dx dx = 0 p0vdx 0 L u ' s first derivativeis integrable u  H 1 v also in H 1 , with addedrequirement that it vanishon boundaries 1 Find u  H 1 such thatBu , v   F (v) v  H 0 1  v  H0 H 1 is a Sobolevspace,which is a subspace a Hilbert space. of H 1 is a spaceof functionsthat can include discontinu ous functionsand functionswith singularities.  singularities and discontinu ities in our solutionOK! Remarks on Variational Problem 1 Find u  H 1 such thatBu, v   F (v) v  H 0 Analyticalmathematics that apply tobilinear functional also apply s to the finite element method. For example : Lax - MilgramTheorem A solutionexists to the variational problemaboveif : 1) B is continuous 2) B satisfiesthe inf  sup condition:   0, inf sup Bx, y    x y 3) sup B( x, y )  0 y, y  0 Also,Babusk a- Brezzi Condition stability for Bar Element example u1 Deformed shape u2 f1 x Element Node (a hinge) f2 Bar Element example (i) Conjecture a displacement function u(x) x  a1  u x   a1  a2 x  1 x    N a (1)    a2  Bar Element example (ii) Express u(x) in terms of nodal displacements by using boundary conditions. Deformed shape u(0) = u1 u(L) = u2  u1  1 0   a1  u   1 L  a    2   2  u  Aa (2) Bar Element example Sub (2) into (1) 1 0  1  u u x   N A  u  1, x    1 L  1  x u  x   1  ,  L x  u  C u L (3) Displacement polynomial that satisfies boundary conditions Bar Element example (iii) Derive strain-displacement relationship by using mechanics theory du d  1  x     Cu   B u   dx dx  L 1 u L ( 4) Axial Strain Bar Element example (iv) Derive stress-displacement relationship by using elasticity theory  x  E x  EBu Elastic Modulus Axial Stress (5) Bar Element example (v) Use principle of Virtual Work Work = Stress x Strain x Volume Internal work Bar cross-sectional area A WI   x . x dxdydz   x . x dx dydz  A  x . x dx  EA B u.B u dx  EAu T BT Bdxu  External work WE  u1  f1     uT f u2     f2  Bar Element example Equate internal and external work WE  WI u f  EAu T T  B Bdxu T f  k  u , Stiffness matrix k   EA BT Bdx (6) Bar Element example Resultant stiffness matrix k   EA  EA L  1   1L  0   L   1 L2  1 0  2  L L 1 1  L L dx   1 L2 dx  1 L2   EA  1  1   1 1  (7) L   Assembling issue (1); Element coords. Element axes are not all the same. So there is a need for a coordinate transformation Z Global axes Y X Assembling issue (1); Element coords. for forces FZ q  FX  cos q       FZ   sin q FX Fz2 fv2  Fx1  cosq        Fz1   sin q    Fx 2   0        Fz 2   0    F  R  f  sin q   f u      cos q   f v   sin q cosq 0 0 0 0 cosq sin q   f u1      0   f v1     sin q   f u 2      cosq   f v 2    0 Fz1 fv1 fu1 fu2 Fx2 q Fx1 (8) Assembling issue (1); Element coords. Similarly for displacement Z2 v2 Z1 v1 u2 X2  sin q cosq 0 0 0 0 cosq sin q   u1      0   v1     sin q  u 2      cosq   v2    0 u1 q X1  X 1  cosq        Z1   sin q   X2  0        Z2   0    U  R u (9) Assembling issue (1); Element coords. Element force-displacement in global coordinate Local coordinates Element nodal forces and displacements in local coordinates  f u1  1 f   EA  0  v1    f u 2  L  1    fv2  0  0  1 0  u1  0 0 0  v1    0 1 0 u2    0 0 0   v2   f  k  u Element stiffness matrix in local coordinates (extended) Global coordinates F  KU , K  Rk R , R  R1 Element stiffness matrix in global coordinates Element nodal forces and displacements in global coordinates Assembling issue (1); Element coords. Element stiffness matrix in global coordinates  cos2 q  EA  cosq sinq 1 K  R k R    L  cos2 q   cosq sinq  cosq sinq sin 2 q  cosq sinq  sin 2 q  cos2 q  cosq sinq cos2 q cosq sinq  cosq sin q    sin 2 q  cosq sinq   sin 2 q   (10) Assembling issue (2); Structure matrix  1 0 1 0 Element and nodal numbering Node number 2 45 0 0 K  1  EA  L  1 0  0 0 0 1 0  2 4 2 4 2 4 2 4 0  0  0 3 90 2 1 P [kN] 1 3 Lm Element number  42  42  2 2  4 K  2  EA  2 42 L  4  24  42  4  0 0 0 1 K  3  EA  L 0 0  0  1  0 0 0 0    42   42   2 4   0  1  0  1 2 4 Assembling issue (2); Structure matrix Create Structure stiffness matrix from element stiffness matrices  Fx1  F   z1    Fx3     Fz 3   1 0  1 0  X 1   0 1,1 0 0 1,3 0  Z  EA   1 L  1 0 1 0  X 3  3,3     3,1  0 0 0 0  Z 3   Fx1  F   z1   Fx 2     Fz 2   Fx3     Fz 3    1 0   EA L    1  0  0 0 0 0  1 0  X 1  0 0  Z1     X 2     Z2  1 0  X 3    0 0  Z 3     Fx 2  F   z2    Fx 3     Fz 3  Assembling issue (2); Structure matrix  42  42  22,2 2  4 4 EA  L 2  2  24 3,2 4 2  4  4    2 4 2 4 2 4 2 4  X2 2,3     42   Z 2   42   X 3   3,3 2   Z3  4   2 4  Fx1  F   z1   Fx 2     Fz 2   Fx 3     Fz 3    1 0   EA L    1  0  0 0 2 4 1 0   2 4 2 4 2 4 2 4 2 4 2 4 2 4  0  0 1   2 4 2 4 2 4 2 4 0  X1  0   Z1    2  X 2  4    42   Z 2   42   X 3    2 4   Z3     Fx1  F   z1    Fx 2    Fz 2   Assembling issue (2); Structure matrix 0 0 0 1,1 1 EA  L 0 0  2,1 0  1  Fx1  F   z1   Fx 2     Fz 2   Fx 3     Fz 3    0   X1  2,1   0  1  Z1   0 0  X2 2,2    0 1   Z2  0 0 0 1 1 0 0 1 0 1 0  2 0 0  42  42 4 EA L  2 0  1  42 1  42 4  2  1 0  42 1  42 4  2  42  42 0 0 4  0  X1  0   Z1    2  X 2  4    42   Z 2   42   X 3    2 4   Z3    Assembling issue (2); Structure b matrix ith Element i a ith Element in Structure Stiffness Matrix with Node numbers a and b a a k aa b kba b k ab kbb a 0   0  0 b    0 0 a  k aa  k ab  ith Element’s Stiffness Matrix with Node numbers a and b b  kba  kbb  Assembling issue (2); Structure matrix  Fx1  F   z1   Fx 2     Fz 2   Fx 3     Fz 3    0 0 1 1 0 0 1 0 1 0  2 0 0  42  42 4 EA  L 2 0  1  42 1  42 4  2  1 0  42 1  42 4  2  42  42 0 0 4  0  X1  0   Z1    2  X 2  4    42   Z 2   42   X 3    2 4   Z3    vector of nodal forces F s  K s U s vector of nodal displacements Structure stiffness matrix Assembling issue (3); Supports • How to deal with the problem of supports (restraints) ? • These are nodes where the displacements are known, zero in the perfectly rigid support case. unknown support  Fx1  Reactions  Fz1  R   0 0 1 1 0 0 1 0 1 0  2  Fx 2  EA  0 0  42  42 4   L  2 0  1  42 1  42 Fz 2  4   2  1 0  42  0  1  42 4    2  42  42  P  0 0   4  known applied nodal loads P 0  0  0  0    2  0  4    42   0   42   X 3    2 4   Z3    known support Displacements Ds unknown nodal Displacements D Assembling issue (3); Supports • This system of equations needs to be partitioned to determine the unknown nodal displacements and support reactions.  R K11   P  K     21 K12  D s  K22   D    1 R  K11 D s  K12 D P  K 21 D s  K 22 D Solve system of algebraic equations key numerical process Determine D  K22  nodal displacements  X 3  P  K21Ds   1  42 1 L  Z   D  K 22  P  EA  2  3  4  Fx1  F   z1   K D  12  Fx 2     Fz 2   1  0 EA  L  42  2  4 2 4 2 4   0       P  1 PL EA  1  12 2    Determine support reactions 0   P  0  PL  1   0     EA     P  2  1 2 2  4     42   P  • Calc. nodal displacements in local coordinates, eqn (9) For Element 1 0   1   0 0 PL   EA   1  0     1  2 2  0  0 0 0  u1  1 0 0  v1    0 1 0 u2    0 0 1   v2  Finishing off; calculate element internal stresses, strains and actions.  u1   0  u    PL   2   EA  • Calc. element strains, equation (4) For Element 1  x   B u   1 L 1 L  0  P  PL    EA  EA  Negative sign indicates Compression Finishing off; calculate element stresses, strains and actions. • Calc. element stresses, equation (5) For Element 1 P  x   E x    A Negative sign indicates Compression • Calc. element axial forces For Element 1 F x   xA  P