The Finite Element Method

1 / 64

# The Finite Element Method - PowerPoint PPT Presentation

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

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

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

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

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
• 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)

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

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

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

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
• 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 Equation

[K]U = F

[K] is thestiffness matrix and F is theload (RHS) vector

[K] is symmetric and positive definite

• 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

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)

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’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’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

U4=9

8

6

u

U3 =?

4

2

U2=?

U1=0

1

2

3

4

x

1. Formulate the weighted residual (weak) form

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

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

1

f1

f2

0.5

element basis functions

0

0

0.5

1

x

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)

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

[k]u = f

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

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
• 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
• 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

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
• 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.

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.

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.

x=0

x=1

x=0

x=0

s3

s2

s1

Scaling Factors

arc length

Global to local mapping:

Scaling Factors  arc lengths

Two-DimensionalTensor-Product Elements

Bilinear interpolation can be constructed

where

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

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

Tri-Cubic Basis Functions

In each node we define:

x3

7

5

x2

6

3

1

2

x1

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
• 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

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

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

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

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

knife

BASE

plunger

APEX

Rabbit Ventricular Anatomy





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

endo>0

epi<0

RIGHT VENTRICLE

LEFT VENTRICLE

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)

X

, longitudinal

2

X

, crossfiber

c

X

, circumferential

1

X

r

X

3

X

, fiber

f

Strain Analysis

A/P View

Lateral View

Reconstructed

3D Coordinates

Transform

0.04

0.00

-0.04

-0.07

End-Systolic Circumferential Strain

Baseline

2 minutes ischemia

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

-0.05 0.00 0.05

0.0 1.5 3.0

mL/min/g

Myocardial Blood Flow

Fiber Strain

Cross-fiber Strain

Control