20.1 The Finite Element Methodology

Let us consider for demonstration purposes the mechanical response of a femur under moderate loads, which after idealization can be considered to behave linearly elastic having isotropic but inhomogeneous material properties. One would be interested in the displacements, strains, and stresses at any point in the femur to determine, for example, if there is a risk of fracture. For this purpose, an FE analysis may be performed that requires three ingredients: (1) the geometric description, (2) a “constitutive model” and material properties that describe the physical phenomenon and the behavior of the material, and (3) boundary conditions (forces and constrains) that are applied.

20.2 The Geometry and Finite Element Mesh

For a patient-specific analysis, the geometry description of the femur may be obtained from a computed tomography (CT) scan, so that the outer and inner contours of the femur can be segmented (Fig. 20.2). After segmentation, the data is manipulated using a computer-aided design software and a solid model of the bone is generated, call it Ω. This solid model of a complex topology can be divided into a collection of smaller, simple geometric shapes such as hexahedral, pentahedral, or tetraherdral elements—each called a finite element. The collection of these elements is called the FE “mesh”. Each element is denoted by Ω_⍳, so that ⋃ Ω_⍳ = Ω.

20.3 The Mathematical Basis of the Finite Element Method

To compute the elastic response of the femur, a set of partial differential equations (the Navier-Lame system,1 known also as the elasticity system, without body forces) has to be solved in the domain of the femur denoted by Ω:

Here u_Simp = [u_x(x,y,z), u_y(x,y,z), u_z(x,y,z)]^T is the sought displacement vector, E(x,y,z) (Young′s modulus), and v(x; y; z) (Poisson ratio) are the material properties at each point. The elasticity system (20.2) is complemented by traction T or displacement boundary conditions. The analytical solution that solves this problem would be the “classical” or “strong” solution u_Simp. However, due to the complex domain of the femur, it is impossible to obtain analytically u_Simp. To obtain an approximated solution by FE methods, we first transform the strong form into a “weak formulation” by multiplying the set of equations by a “test function” v = (v_x(x,y,z),v_y(x,y,z),v_z(x,y,z))^T, then integrate over the entire femur and apply the Green′s lemma (integration by parts) so we finally obtain the formulation:

The weak formulation of equation (20.3) can also be cast into an equivalent formulation, the “minimum potential energy”:

Seek u_Simp that minimizes the “potential energy”:

Once u_Simp is found, one can easily compute the stress vector by σ = [E] [D] u_Simp.

To solve equation (20.4), a major difficulty still exists, namely, one is required to find one displacement field among an infinite number of displacements that when inserted into ∏(u) will result in a minimum value. Let us denote by ℰ(Ω) the space of all possible displacements, so we seek for u_Simp ∊ ℰ(Ω). From the practical point of view, this is of course impossible, and instead of an infinite number of possible functions, we seek a displacement field within a finite number of functions (we denote this finite set of functions S (Ω) ⊂ ℰ(Ω) that satisfies equation (20.4). This function is the FE solution, u_FE. Of course, that u_FE may not provide the minimum value to ∏(u), among all possible functions in ∏(Ω), but it is the “closest” to u_Simp among all functions in S (Ω). The discretization error is introduced because we seek for a solution in S(Ω) instead of in ℰ(Ω). If one uses hierarchical spaces, that is a family of increasingly larger spaces, each containing the smaller spaces, S₁(Ω) ⊂ S₂(Ω) ⊂ S₃(Ω) … ⊂ ℰ(Ω), then the FE solutions u_FE1, u_FE2, u_FE3 … obtained are progressively closer to u_Simp. Of course, as the space S (Ω) is enriched by more and more functions, then the approximated solutions become closer to the exact solution. The systematic enrichments of the subspaces are called extensions, and these are mandatory for estimating the discretization error to quantify how close are u_FE1, u_FE2, u_FE3 … to u_Simp. We will get back to the estimation of the discretization later on.

Having a mesh, the integral over the entire domain in equation (20.4) is the sum of integrals over the elements, so that equation (20.4) can now be stated as:

One may observe that we partitioned the overall problem into “elemental problems,” and have now to compute the integrals in equation (20.5) over each element. Integrating over different-shaped hexahedra elements, for example, is a complicated procedure. To overcome this difficulty, a standard element is introduced Ω_ST = {(ζ, ŋ, ς) −1 ≤ ζ ≤ 1, −1 ≤ ŋ ≤ 1, −1 ≤ ς ≤ 1,} and an associated “mapping” is possible so to “map” the standard hexahedral element into any “real” hexahedral element with curved faces (see Fig. 20.3) by using the blending function method.2 In this manner, each integral in equation (20.5) is evaluated over the standard element introducing the Jacobian of the mapping. For example:

Because the minimum potential energy formulation has been split into a sum over all elements, and furthermore, all integrations are performed over the standard element, it is only natural to define a basis function on the standard element. In each element, the displacements are represented by a set of basis functions—these are also called shape functions. Let u_FE in each element be expressed in terms of elemental M basis functions, which are constructed by monomials, N_i(ζ, ŋ, ς) (spanning a finite-dimensional subspace):

where a_i^(l) are the amplitudes of the basis functions in element l, and for example N1=18(1−ζ)(1−η)(1−ς)N_1 = {1 \over 8}(1 – \zeta )(1 – \eta )(1 – \varsigma ). Seeking for u_FE reduces to the determination of the unknowns a_i^(l). These are collected from all elements in the FE model and stored in a vector denoted a. The number of entries in the vector a represents the dimension of the space in which u_FE is sought and is called the number of degrees of freedom (DOFs).

Substituting equation (20.7) in equation (20.6) and thereafter in equation (20.5), then collecting the contribution of all elements (a process named the assembly process), one obtains the set of algebraic equations, which may be stated as:

Find the vector a that minimizes

where [K] is denoted stiffness matrix and r is the load vector. The FE solution is obtained by inverting the stiffness matrix and multiplying it by the load vector. In most FE methods, the stiffness matrix is symmetric and sparse, therefore efficient numerical methods exist that invert it (e.g., a problem with a million DOFs can be solved nowadays on a personal computer within less than half an hour).