1 / 64

The Finite Element Method

NBCR Summer Institute 2006: Multi-Scale Cardiac Modeling with Continuity 6.3 Wednesday: Finite Element Discretization and Anatomic Mesh Fitting Andrew McCulloch and Fred Lionetti. The Finite Element Method. Evolved first from the matrix methods of structural analysis in the early 1960’s

blaine-everett
Download Presentation

The Finite Element Method

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. NBCR Summer Institute 2006:Multi-Scale Cardiac Modeling with Continuity 6.3Wednesday:Finite Element Discretization and Anatomic Mesh FittingAndrew McCulloch and Fred Lionetti

  2. The Finite Element Method • Evolved first from the matrix methods of structural analysis in the early 1960’s • Uses the algorithms of linear algebra • Later found to have a more fundamental foundation • The essential features are in the formulation • There are two alternative formulations that are broadly equivalent in most circumstances • Variational formulations, e.g. the Rayleigh-Ritz method • Weak or weighted residual formulations, e.g.the Galerkin method • Both approaches lead to integral equations instead of differential equations (the strong form)

  3. 24 23 15 14 22 21 12 13 Finite differences: R = 0   Finite elements: R = 0 The Finite Element Method • Solution is discretized using a finite number of functions • Piecewise polynomials (elements) • Continuity across element boundaries ensured by defining element parameters at nodes with associated basis functions, • FE equations are derived from the weak form of the governing equations

  4. 8 7 24 23 15 14 3 4 22 21 12 13 5 6 Global equations 2 1 Element equations S The Finite Element Method • Integrate governing equations in each element • Assemble global system of equations by adding contributions from each element

  5. Integral Formulations Consider the strong form of a linear partial differential equation, e.g. 3-D Poisson’s equation with zero boundary conditions: On region R on boundary S Strong FormLu = f Variational Principle, e.g. minimum potential energy Weighted Residual (weak) form, e.g. virtual work

  6. Weak Form for 2-D Poisson’s Equation 0 0 Where, u and w vanish at the boundary On region S on boundary C Weak form Integrate by parts

  7. Galerkin’s Method for 2-D Poisson’s Equation • Choose a finite set of approximating (trial) functions, i(x,y), i = 1, 2, …, N • Allow approximations to u in the formU(x,y) = U11 + U22 + U33 + … + UNN (that can also satisfy the essential boundary conditions) • Solve N discrete equations for U1, U2, U3, …, UN

  8. Galerkin’s Method for 2-D Poisson’s Equation [K]U = F [K] is thestiffness matrix and F is theload (RHS) vector [K] is symmetric and positive definite

  9. Comments on Galerkin’s Method • Galerkin is more general than Rayleigh-Ritz. If we add u/x, symmetry & the variational principle are lost, but Galerkin still works • If w is chosen as Dirac delta functions at N points, weighted residuals reduces to the collocation method • If w is chosen as the residual functions Lu-f, weighted residuals reduces to the least squares method • By choosing w to be the approximating functions, Galerkin’s method requires the error (residual) in the solution to be orthogonal to the approximating space. • The integration by parts (Green-Gauss theorem) automatically introduces the Neumann (natural) boundary conditions • The Dirichlet (essential) boundary conditions must be satisifed explicitly when solving [K]U=F • Since discretized integrals are sums, contributions from many elements are assembled into the global stiffness matrix by addition. • The Ritz-Galerkin FEM finds the approximate solution that minimizes the error in the energy

  10. Steps in the Finite Element Method 1. Formulate the weighted residual (weak form) 2. Integrate by parts (or Green-Gauss Theorem) • reduces derivative order of differential operator • naturally introduces derivative (Neumann) boundary conditions, e.g. flux or traction. Hence called that natural boundary condition 3. Discretize the problem • discretize domain into subdomains (elements) • discretize dependent variables using finite expansions of piecewise polynomial interpolating functions (basis functions) weighted by parameters defined at nodes

  11. Steps in the Finite Element Method (…cont’d) 4. Derive Galerkin finite element equations • substitute dependent variable approximation in weighted residual integral • Choose weight functions to be interpolating functions — the Galerkin assumption (Galerkin, 1906) 5. Compute element stiffness matrices and RHS • integrate Galerkin equations over each element subdomain • integrate right-hand side to obtain element load vectors which also include any prescribed Neumann boundary conditions

  12. Steps in the Finite Element Method (…cont’d) 6. Assemble global stiffness matrix and load vector • Add element matrices and RHS vectors into global system of equations • Structure of global matrix depends on node ordering 7. Apply essential (i.e. Dirichlet) boundary conditions • at least one is required (essential) for a solution • prescribed values of dependent variables at specified boundary nodes, e.g. prescribed displacements • eliminate corresponding rows and columns from global stiffness matrix and transfer column effects of prescribed values to Right Hand Side •  the constraint reduced system

  13. Steps in the Finite Element Method (…cont’d) 8. Solve global equations • for unknown nodal dependent variables • using algorithms for Ax = b or Ax = x 9. Evaluate element solutions • interpolate dependent variables • evaluate derivatives, e.g. fluxes • derived quantities, e.g. stresses or strain energy • graphical visualization; post-processing 10. Test for convergence • refine finite element mesh and repeat solution

  14. Galerkin FEM: Simple 1-D Example U4=9 8 6 u U3 =? 4 2 U2=? U1=0 1 2 3 4 x

  15. 1. Formulate the weighted residual (weak) form 2. Integrate by parts (or Green-Gauss Theorem)

  16. 3. Discretize the problem • 4 global nodal parameters U1, U2, U3, U4 • 3 linear elements each with 2 element nodal parameters u1, u2. • Adjacent elements share global nodal parameters, e.g., global parameter U2 is element parameter u2 of element 1 and u1 of element 2. • Two (linear) element interpolation functions for each element, i(x), i = 1, 2 • Allow element approximations to u in the form • u(x) = u1 1 + u2 2 = ui i i=1,2

  17. Element Basis Functions 1 f1 f2 0.5 element basis functions 0 0 0.5 1 x

  18. 4. Derive Galerkin equations for each element In each element, let u(x)  u1 1 + u2 2 = ui i (x) and w(x)i (x)

  19. 4. Derive Galerkin equations for each element (… cont’d) e.g. for Element 1(no derivative boundary conditions): [k] = [(kij)] is the element stiffness matrix f = (fi) is the element load vector

  20. 5. Compute element stiffness matrices [k]u = f Element stiffness matrix, [k] and load (RHS) vector, f

  21. 5. Compute element RHS matrices In this problem, each element is the same size and thus: [k](ele 1) = [k](ele 2) = [k](ele 3) and: f(ele 1) = f(ele 2) = f(ele 3)

  22. 6. Assemble global stiffness matrix and load vector

  23. 7. Apply essential (i.e. Dirichlet) boundary conditions That leaves global equations 2 and 3

  24. 8. Solve global equations (constraint-reduced) Exact!

  25. Representing a One-Dimensional Field • Polynomials are convenient, differentiated and integrated readily • For low degree polynomials this is satisfactory • If the polynomial order is increased further to improve the accuracy, it oscillates unacceptably • Divide domain into subdomains and use low order piecewise polynomials over each subdomain – called elements

  26. Making Piecewise Polynomials Continuous • constrain the parameters to ensure continuity of u across the element boundaries • or better, replace the parameters a and b in the first element with parameters u1 and u2, which are the values of u at the two ends of that element: • where is a normalized measure of distance along the curve

  27. u = u(x) u u = c + dx u = e + fx u = a + bx + + + + + + + + + + + + + x u x x u = (1- ) u + u 1 2 u 1 + nodes u 0 + 3 + + + + + x + + + + + u 1 2 u + 4 element 1 element 2 element 3 x

  28. Linear Lagrange Interpolation

  29. Global-Element Mapping • Associate the nodal quantity un with element node n • Map the value U defined at global node  onto local node n of element e by using a connectivity matrix  (n, e), Thus, in the first element withu1=U1 and u2=U2.. In the second element u is interpolated by Withu1=U2andu2=U3.

  30. u u1 u u1 u2 x 1 u2 x x x1 x2 x2 x1 x 1 Isoparametric Interpolation We haveu () but to defineu (x) we needx (). Definexas an interpolation of nodal values, e.g.

  31. Quadratic Lagrange Basis Functions Use three nodal parametersu1, u2 and u3 2 1 1.0 1.0 x x 1.0 0 1.0 0.5 0 0.5 3 1.0 x 1.0 0 0.5 are the quadratic Lagrange basis functions.

  32. Cubic Hermite Basis Functions 1 x 0 1 x 1 0 1 x 0 1 x 1 0

  33. x=0 x=1 x=0 x=0 s3 s2 s1 Scaling Factors arc length Global to local mapping: Scaling Factors  arc lengths

  34. Two-DimensionalTensor-Product Elements Bilinear interpolation can be constructed where

  35. 2 1 u =  u n n u  1 0 1 y  1 1 x  2 x =  x n n y =  y n n 0 1  1 Bilinear Tensor-Product Basis Functions

  36. A Six-NodedQuadratic-Linear Element x2 1.0 x1 0 0 0.5 1.0

  37. Three-dimensional Linear Basis Functions e.g. trilinear element has eight nodes with basis functions: x3 7 8 5 x2 6 3 4 1 2 x1

  38. Tri-Cubic Basis Functions In each node we define: x3 7 5 x2 6 3 1 2 x1

  39. Tri-Cubic Basis Functions (Cont’d)

  40. x=0 x=1 x=0 x=0 s3 s2 s1 Scaling Factors arc length Global to local mapping: Scaling Factors  arc lengths

  41. Coordinate Systems • Rectangular Cartesian global reference coordinate system • Orthogonal curvilinear coordinate system to describe geometry and deformation • Curvilinear local finite element coordinates • Locally orthonormal body coordinates define material symmetry and structure, related to the finite element coordinates by a rotation about the -normal axis through the "fiber angle" , From Costa et al, J Biomech Eng 1996;118:452-463

  42. Curvilinear World Coordinates D) Prolate Spheroidal Coordinates (L,M,Q)

  43. Fiber/Sheet Coordinates

  44. Coordinate System Notations

  45. Fitting with Linear Lagrange 1-D Elements Two linear Lagrange elementsfit the data with a root-mean-squared-error (RMSE) of 0.614892. Result of twice refining the mesh (yielding 8 elements) andre-fitting: RMSE = 0.0930764

  46. Least Squares Fitting The least squares fit minimizes the objective function: where is measured coordinate or field variable; is the interpolated value at are weights applied to the data points are smoothing weights

  47. Fitting a Coronary Vascular Tree with Quadratic Lagrange 1-D Elements

  48. heart cast in rubber knife tube plunger Rabbit Ventricular Anatomy • anesthetized & ventilated New Zealand White rabbit • heart arrested in diastole, excised • pulmonary vessels removed, aorta cannulated • heart suspended in Ringers lactate, perfused in unloaded state with buffered formalin at 80 mm Hg for 4 minutes • heart cast in polyvinylsiloxane

  49. knife BASE plunger APEX Rabbit Ventricular Anatomy

More Related