slide1 n.
Download
Skip this Video
Loading SlideShow in 5 Seconds..
The Finite Element Method PowerPoint Presentation
Download Presentation
The Finite Element Method

Loading in 2 Seconds...

play fullscreen
1 / 64

The Finite Element Method - PowerPoint PPT Presentation


  • 226 Views
  • Uploaded on

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

loader
I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.
capcha
Download Presentation

PowerPoint Slideshow about 'The Finite Element Method' - blaine-everett


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.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.


- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -
Presentation Transcript
slide1

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

the finite element method
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)
the finite element method1

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
the finite element method2

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
integral formulations
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

weak form for 2 d poisson s equation
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

galerkin s method for 2 d poisson s equation
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
galerkin s method for 2 d poisson s equation1
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

comments on galerkin s method
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
steps in the finite element method
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
steps in the finite element method cont d
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
steps in the finite element method cont d1
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
steps in the finite element method cont d2
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
galerkin fem simple 1 d example
Galerkin FEM: Simple 1-D Example

U4=9

8

6

u

U3 =?

4

2

U2=?

U1=0

1

2

3

4

x

1 formulate the weighted residual weak form
1. Formulate the weighted residual (weak) form

2. Integrate by parts (or Green-Gauss Theorem)

3 discretize the problem
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
element basis functions
Element Basis Functions

1

f1

f2

0.5

element basis functions

0

0

0.5

1

x

4 derive galerkin equations for each element
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)

4 derive galerkin equations for each element cont d
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

5 compute element stiffness matrices
5. Compute element stiffness matrices

[k]u = f

Element stiffness matrix, [k] and load (RHS) vector, f

5 compute element rhs matrices
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)

representing a one dimensional field
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
making piecewise polynomials continuous
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
slide27

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

global element mapping
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.

slide30

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.

quadratic lagrange basis functions
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.

scaling factors

x=0

x=1

x=0

x=0

s3

s2

s1

Scaling Factors

arc length

Global to local mapping:

Scaling Factors  arc lengths

two dimensional tensor product elements
Two-DimensionalTensor-Product Elements

Bilinear interpolation can be constructed

where

slide35

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

slide37

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

slide38

Tri-Cubic Basis Functions

In each node we define:

x3

7

5

x2

6

3

1

2

x1

scaling factors1

x=0

x=1

x=0

x=0

s3

s2

s1

Scaling Factors

arc length

Global to local mapping:

Scaling Factors  arc lengths

coordinate systems
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

curvilinear world coordinates
Curvilinear World Coordinates

D) Prolate Spheroidal Coordinates (L,M,Q)

slide45

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

slide46

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

slide48

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
slide49

knife

BASE

plunger

APEX

Rabbit Ventricular Anatomy

slide51





4

data point projects onto surface at

ddd 

{

i

1

(

,

)

(

,

)

1

2

i

1

2

i

1

i



2

(

,

)

i

1

2



1

i



3

(

,

)

i

1

2



2

2

i

}

4

(

,

)

i

1

2





1

2

x = d cosh  cos 

y = d sinh  sin  cos 

z = d sinh  sin  sin 

Bicubic Hermite

isoparametric interpolation

slide55

endo>0

epi<0

RIGHT VENTRICLE

LEFT VENTRICLE

slide56

Anatomic Model

  • 8,351 geometric points
  • 14,368 fiber angles
  • 36 elements
  • 552 geometric DOF RMSE = ±0.55 mm
  • 184 Fiber angle DOF RMSE = ±19°

Vetter & McCulloch Prog Biophys & Mol Biol 69(2/3):157 (1998)

strain analysis

X

, longitudinal

2

X

, crossfiber

c

X

, circumferential

1

X

, radial

r

X

, radial

3

X

, fiber

f

Strain Analysis
slide58

A/P View

Lateral View

Reconstructed

3D Coordinates

Transform

slide59

0.04

0.00

-0.04

-0.07

End-Systolic Circumferential Strain

Baseline

2 minutes ischemia

slide60

0

.

7

0

.

6

0

.

5

RMS Fitting Error (mm)

0

.

4

0

.

3

0

.

2

0

.

1

0

-7

-6

-5

-4

-3

-2

-1

0

10

10

10

10

10

10

10

10

a

Smoothing Weight

slide62

-0.05 0.00 0.05

0.0 1.5 3.0

mL/min/g

Myocardial Blood Flow

Fiber Strain

Cross-fiber Strain

Control

LAD Occlusion

slide64

3months post-surgery

Pre-surgery

LATERAL

SEPTAL