6. Numerical Results 141
4.7. Chebyshev semi-iteration
Realizing a fixed point iteration we need that the spectrum of A is contained in [a−b, a+b]. In particular for efficient application of the Chebyshev semi-iteration good estimates for spectral bounds λmin and λmax are required. For the case of a mass matrix that is preconditioned by a one-step Jacobi preconditioner, sufficiently good estimates for the spectral bounds have been obtained for many discretiza-tions [273]. In this case, the Chebyshev semi-iteration provides a cheap and efficient
110
4.5 Approximation of operators
iterative solver and, if employed with a fixed number of iterations, it realizes a linear preconditioner.
Linearity of preconditioners is an important property. First because this is what we employ in theoretical considerations. Second, if nonlinear preconditioners such as a CG method – even with a fixed number of steps – are used as preconditioners for Krylov solvers this can significantly reduce their performance. This issue is nicely illustrated in [274].
4.5.2. Approximation of the stiffness matrix
Regarding the discretizations of the differential operators, occurring in the con-straint and the adjoint equation, hierarchical multigrid solvers and preconditioners are among the most promising methods for adaptive algorithms. The two main reasons are the h-independent convergence rate and linear complexity with respect to the number of unknowns. Both properties are not only provable, but also can be verified for reasonable implementations of multigrid algorithms [260].
These mainly exploit the observation that the classification of high- and low-frequency error components in finite element (FE)-computations is dependent on the resolu-tion of the spatial domain. High-frequency errors on a fine grid are not captured on coarser grids. In contrast low-frequency errors on fine grids may appear high-frequency on coarser ones, whereas smooth error components essentially can be represented on significantly coarser grids.
In order to exploit this insight, multigrid solvers use cheap smoothers, such as damped Jacobi- or Gauss-Seidel-iterations [79], to eliminate high-frequency error components. Repeated application on grid with different spatial resolution then admits to significantly reduce the oscillatory components of the algebraic error.
Eventually on the coarsest grid the remaining error can be eliminated by a direct solver.
To express this formally we consider the operator equation Au = b in X∗ and a sequence of hierarchical grids, resp. nested FE-spaces
S0 ⊂ · · · ⊂Sj ⊂X
with corresponding projected differential operators Ak = A|Sk. We denote the pro-jection operators by Ik−1k : Sk−1 → Sk and the corresponding restriction operators byIkk−1 : Sk →Sk−1. The essential ingredients of the multigrid method are captured in Alg. 4.8.
For adequate choices of the algorithmic parameters, ν1 applications of the smoother eliminate the high-frequency error components in Sk. Then, the remaining error is projected to a coarser grid and eliminated there. Assuming that the remaining error is “smooth” with respect to Sk we expect relatively small loss due to restriction and interpolation operators. After the computation of the coarse grid correction one
Chapter 4 Computation of steps for optimal control problems
Require: given smoothing parametersν1, ν2, grid levelk, operatorAkand right hand side bk, initial value uk.
1: mgCycle:
2: if k = 0 do
3: use direct solver to solveA0δu0 =b0
4: u0 ←u0+δu0
5: else do
6: apply ν1 steps of a smoother to the system Akuk=bk
7: compute the restriction of the residual to the coarse spacerk−1 ←Ikk−1rk= Ikk−1(bk−Akuk)
8: compute Ak−1 =Ik−1k AkIkk−1, set δuk−1 = 0
9: apply mgCycle(ν1, ν2, k−1, Ak−1, rk−1, δuk−1)
10: correct the fine grid solutionuk ←uk+Ik−1k δuk−1
11: apply ν2 steps of a smoother to the system Akuk=bk Algorithm 4.8: Two-grid correction scheme.
further applies a smoother ν2 times to eliminate possibly remaining high-frequency error contributions. On the coarsest grid a direct solver can be employed to compute a highly accurate coarse grid correction. Alg. 4.9 is essentially the V-cycle multigrid method. For the realization of a linear solver, that can satisfy prescribed accu-racy requirements, we repeatedly apply the above algorithm on the defect equation.
See Briggs et al. [45], Trottenberg et al. [260] for more details and the W-cycle as well as the full multigrid (FMG) scheme. Same as the Chebyshev semi-iteration the multigrid method provides a linear preconditioner if employed with a fixed number of iterations.
Require: given smoothing parametersν1, ν2, operator A and right hand side b.
1: while convergence test failed do
2: setr←b−Au and δu= 0
3: apply mgCycle(ν1, ν2,0, A, r, δu)
4: set u←u+δu
Algorithm 4.9: V-cycle multigrid method.
On optimized academic examples multigrid solvers are extremely fast. This is il-lustrated in Tab. 4.2 where the iteration numbers are given for a simple example of linear heat transfer on the unit cube, with right hand side f = 1. On real-world problems, where already small and coarse geometries may contain larger numbers of degrees of freedoms in a linear finite element space, their performance is less out-standing, but still impressive. This is illustrated in Tab. 4.3 for equations of linear elasticity with the same right hand side. On the left of this table the iterations
112
4.5 Approximation of operators
numbers and computational time on the unit cube are given. On the right we give the corresponding numbers on the geometry used in Sec. 6.2.2.1. Thus, we will have to apply several V-cycles to get a reasonable preconditioner for our KKT-systems.
If used to approximate the differential operators in the constraint preconditioner (4.2.4) we even need to actually solve the constraint and adjoint equation. For these reason the repeated application of multigrid solvers and preconditioners for problems in constraint preconditioners for PPCG methods still is quite expensive.
Laplace (2D, unit cube) Laplace (3D, unit cube)
dof #iter. time dof #iter. time
545 7 2.2 ms 369 8 2.6 ms
8 321 7 32 ms 2 465 9 21 ms
33 025 7 0.13 s 17 985 10 0.18 s
131 585 7 0.59 s 137 345 10 1.7 s
2 099 201 7 11 s 1 073 409 11 16 s
Table 4.2.: Computation times for simple test problems (rel. acc.: 10−9, smoothing steps: 10).
Linearized elasticity (unit cube) Linearized elasticity (real world geometry)
dof #iter. time dof #iter. time
53 955 43 3.6 s 213 687 25 18 s
412 035 49 35 s 1 523 604 48 4.4 min
3 220 227 53 5 min 11 461 518 72 49 min
Table 4.3.: Computation times for simple test problems from linearized elasticity (both 3D, rel. acc.: 10−9, smoothing steps: 10, resp. 20 for the second problem)
We note that in the context of pressure-type boundary conditions the lower left and upper right block in the KKT-system not only contain the differential operator Wϕϕ, but also non-symmetric contributions from the Piola-transformed pressure-type boundary conditions gcof(∇ϕ)n. Moreover, due to the polyconvexity of the stored energy function, the corresponding differential operator may not be elliptic.
Consequently, in our setting, the application of multigrid solvers and preconditioners is not backed by a solid theoretical basis. However, in our computations they seemed to work well.
Chapter 4 Computation of steps for optimal control problems
4.6. Summary
Let us summarize the step computations within the different settings. For the computations of δn, px and δs we can always assume positive definiteness ofM on kerc0(x), and thus unique solvability of the corresponding system. For moderately sized problems, or a low dimensional control space, the solution can be found by direct elimination methods. Otherwise a PPCG method, i.e. a conjugate gradient method combined with a constraint preconditioner, can be used.
The situation is different for tangential steps δt. As Lxx is in general indefinite we have to use one of the modifications from Sec. 4.3, the HCG method, to compute descent directions for the cost functional. The restriction to kerc0(x) is again in-corporated with the help of a constraint preconditioner. For problems of moderate size or low dimensional control space we can reuse the direct factorization which was computed for the determination of the normal step as preconditioner. If this approach is not admissible we will use the same block triangular preconditioner
Qsc =
as in the computation of the normal step. To efficiently evaluate this preconditioner we have to efficiently solve the state, variational and adjoint equation. For this we replace the inversion of bothA and AT by a multigrid solver. Since the constraint preconditioner has to project onto kerc0(x0), and in the absence of further analysis, it is necessary to solve the arising systems A0(y0)δny =rp+Bδnu to high accuracy.
Relaxing this condition on Qsc is subject to ongoing work. In contrast Mu−1 can be replaced by a fixed number of Chebyshev semi-iterations [110, 124, 274], which needs not to be overly accurate. The required spectral bounds for the Chebyshev semi-iterations for the preconditioned matrixQjacMu, whereQjacrepresents one step of the Jacobi iteration, are taken from [273].
For the purpose of error estimation and adaptive mesh refinement in an affine co-variant setting, a hierarchical estimator for the error of the primal variables in the Lagrange-Newton step was proposed. Both the block triangular preconditioner and the error indicator require the solution of equations involving the mass matrix resp. the differential operators of the state and adjoint equation. For the first again the Chebyshev semi-iteration is employed, whereas the differential operators are treated with multigrid preconditioners (25 V-cycles and 20 pre- and post-smoothing steps).
114
5. Mechanical behavior of biological soft tissues
In the last chapters a model for the description of the implant shape design problem and an algorithm for its solution have been presented. It leaves to specify de-scriptions for the occurring soft tissue types. In literature often models of linearized elasticity are employed. This is on the one hand due to their simple formulation that are straightforward to implement efficiently. On the other hand realistic models are challenging for numerical solvers. Thus, in order to perform relevant computations, we need to understand state-of-the-art models for biological soft tissues.
Compared with industrially manufactured materials, biological soft tissues exhibit a far more complex mechanical behavior. This includes anisotropy [168], thermo-and viscoelastic behavior [220, 280] as well as complex metabolic interactions [99, 135, 136, 139, 147, 255] and self-regulating mechanisms such as tissue growth [212, 253, 288].
Due to the various different observed phenomena, general accurate mathematical descriptions of the mechanical properties of biological soft tissues are not avail-able. Nonetheless in the described setting of polyconvex hyperelasticity reasonable descriptions of the static mechanical behavior have been proposed for many soft tissue types. These descriptions are mainly based on phenomenological continuum models1.
In most biological soft tissues the mechanical properties are essentially determined through the properties of the extracellular matrix (ECM) between the cells. The ECM itself mainly consists of ground substance, connecting the cells and its chemical processes, and of three different fiber types (elastic, collagen and reticulin fibers).
• Elastic fibers mainly consist of elastin surrounded by fibrillin and are respon-sible for the tensile elastic properties of a tissue, i.e. the ability to return to its initial configuration when being stretched and then released.
1In the context of modeling it is often distinguished between three different approaches. Con-tinuum models are derived from general considerations on the macroscopic material structure.
Phenomenological modelsrely on the fitting of heuristically chosen expressions to experimental data. Structural models use probabilistic descriptions of insights in the microscopic material structure for the description of macroscopic properties. We will that state-of-the-art models for biological soft tissues are in general based on a combination of all three of the mentioned strategies. Insight in a tissue’s micro-structure and functional properties are extensively used to derive structural models and as guidance for the development of phenomenological models, both incorporating basic requirements from continuum mechanics.
Chapter 5 Mechanical behavior of biological soft tissues
• Collagen fibers, often occurring as parallely aligned bundles, can sustain large tensile forces thus giving strength to biological tissues. “Its importance to man may be compared to the importance of steel to our civilization [...]” Fung [99, p. 251].
• Reticulin fibersconsist of collagen surrounded by a glycoprotein and form nets from fine collagen bundles. Similar to collagen fibers they are responsible for the ability to sustain forces in many tissues such as blood vessels or smooth muscle tissue.
The presence of elastic and collagen fibers is illustrated in Fig. 5.1 for the case of subcutaneous tissue such as found in the papillary dermis, a deeper layer of skin tissue.
Eventually muscles themselves are assembled in a fibrous structure. Fibers largely determine the tensile mechanical behavior of biological soft tissues. Under compres-sion they buckle and do not contribute directly to the mechanical response [136].
But they contribute indirectly, due to their influence on the permeability of liquids.
This indirect relation is usually neglected in biomechanical models.
Figure 5.1.: Fiber structure of subcutaneous tissue.
We begin with the introduction of general strategies for the description of constitu-tive relations in the hyperelastic setting and the commonly used setting for fiber-reinforced materials in Sec. 5.1. Then, specifications for the behavior of different tissue types with respect to tensile (Sec. 5.2) and compressive forces (Sec. 5.3) are discussed. Eventually we shortly address the attainment of patient-specific material parameters in Sec. 5.4.
116
5.1 Modeling framework
5.1. Modeling framework
In order to reduce the complexity in the development of constitutive laws for biolog-ical soft tissues it is commonly assumed that the polyconvex stored energy function W can be split into independent contributions2. Most of these splittings are of the form
W =Wt+Wc, (5.1.1)
where Wt describes the behavior with respect to isochoric (volume-preserving) de-formations , and Wc the behavior with respect to purely volumetric deformations.
The latter is assumed to be isotropic, whereas
Wt =Wt,iso+Wt,aniso, (5.1.2)
can be split into isotropic and anisotropic contributions Wt,iso, resp. Wt,aniso. The stored energy function Wt,aniso describes the influence of fibers. Models of this form are called fiber-reinforced models. Since fibers buckle under compression they are considered to be relevant only in the presence of tensile forces. The contrary holds for liquids that are dominate the elastic behavior under compressive forces. In biomechanical models, the splitting (5.1.1) often corresponds to different models for the description of tension and compression experiments.
5.1.1. Isotropic materials
Isotropic materials do not contain any directional information on a macroscopic level.
Definition 5.1. A response function ˆT : ¯Ω×M3+ →S3 is isotropic at x∈ Ω if it¯ satisfies
Tˆ(x, F Q) = ˆT(x, F) for all F ∈M3+ and all Q∈O3+.
We call an elastic material isotropic if the response function of the corresponding stress tensor is isotropic at every x∈Ω.¯
Isotropy implies that we can express the stress tensor in terms of the strain tensor instead of the deformation gradient.
Theorem 5.2. A response function Tˆ: ¯Ω×M3+ → S3 is isotropic at x ∈ Ω¯ if and only if there exists a mapping T˜(x,·) :S3+ →S3 such that for all F ∈M3+ holds
Tˆ(x, F) = ˜T(x, F FT).
2In fact all common splittings do not yield independent contributions. In particular for complex models it is rather unclear how the interaction between the different summands affects the overall model.
Chapter 5 Mechanical behavior of biological soft tissues Proof. See Ciarlet [56, Thm. 3.4-1].
Incorporating also the property of frame-indifference admits to specify the form of stress tensors for isotropic materials. First recall the invariants of a matrix.
Definition 5.3. Let A ∈ M3 with eigenvalues λi, i = 1,2,3. Then its principal invariantsιA = (ι1, ι2, ι3) are
ι1 = tr(A) = λ1+λ2+λ3, ι2 = 1
2
tr(A)2−tr(A2)
= tr(cof(A))
= λ1λ2+λ2λ3+λ3λ1, ι3 = det(A) =λ1λ2λ3.
Remark 5.4. If A = FTF ∈ S3+, for F ∈ M3+, is a left Cauchy-Green strain tensor, then its eigenvalues are also called principal stretches [201]. The corresponding invariants ι1, ι2 and ι3 are called principal strain invariants. Note that in this case we have
ι3 = det(A) = det(FTF) = det(F)2.
Thus for some deformation ϕ and F = ∇ϕ, the third principal strain invariant describes volume changes.
Next to the famous existence result of Ball [15], which is not restricted to isotropic materials, the Rivlin-Ericksen theorem is the most important theorem in the theory of isotropic elasticity.
Theorem 5.5 (Rivlin-Ericksen). A mapping Tˆ:M3+ →S3 satisfies
Tˆ(QF) = QTˆ(F)QT and Tˆ(F Q) = ˆT(F) (5.1.3) for all F ∈M3+ and all Q∈O3+, if and only if for all F ∈M3+ holds
Tˆ(F) = ˜T(F FT), where T˜: S3+→S3 is of the form
T˜(C) =f0(ιC)I +f1(ιC)C+f2(ιC)C2,
and fi, i= 1,2,3 are real-valued functions of the principal strain invariants ιC = (ι1, ι2, ι3) of the strain tensor C =FTF ∈S3+.
Proof. See [56, Thm. 3.6-1]. Note that no regularity assumptions on ˜T are required.
Instead symmetry of the strain and stress tensors, the axiom of frame-indifference as well as the assumption of isotropy (5.1.3) yield simultaneous diagonalizability of C and ˜T(C) as essential tool in the proof of the Rivlin-Ericksen theorem.
118
5.1 Modeling framework
The Rivlin-Ericksen theorem tells us that for hyperelastic materials the stress tensor, and consequently the stored energy function W, can be written in terms of the principal strain invariants. Since the invariant ι3 describes volume changes the splitting(5.1.1) is mostly assumed to be of the form
W =Wt(ι1, ι2) +Wc(ι3).
However, the first and second principal strain invariant are not isochoric. For this reason first expressions for the influence of the shear parts of the first and second principal invariant on the volumetric part had already been derived in [209]. How-ever, for sake of simplicity, it is commonly assumed that Wc only depends on the third strain invariant.
For a proper splitting into isochoric and volumetric contributions Penn [210] intro-duced a modified set of invariants ¯ιC.
Definition 5.6. Let A ∈ M3 with eigenvalues λi, i = 1,2,3. Then its modified principal invariants ¯ιA = (¯ι1,¯ι2,¯ι3) are
¯
ι1 = ι1ι−
1 3
3 ,
¯
ι2 = ι2ι−
2 3
3 ,
¯
ι3 = ι3.
Remark 5.7. Since we can recover the principal invariants from the modified principal invariants we may also use the latter in the Rivlin-Ericksen theorem.
Penn then investigated models of the form
W =Winc(¯ι1,¯ι2) +Wvol(¯ι3).
Assuming constant compressibility, Penn concluded that for large strains the relation between volume changes and stresses is no more physically reasonable and rejected this approach. Similar results in numerical experiments have been presented in [90].
Nonetheless this strategy is appealing as it facilitates the interpretation of the contri-butions of tensile and compressive responses and its fitting to experimental results, and is widely used [128, 135, 140, 137]. Moreover, premature rejection of the splitting based on modified invariants is challenged by Hartmann and Neff [128]. Stressing the fact that actually the interplay of complex nonlinear models for applied tensile and compressive forces is not well understood, an expression for the volumetric part is proposed that does not admit the previously observed deficiencies, see Sec. 5.3.
The isotropic parts of the models presented in Sec. 5.2 and Sec. 5.3 will all be of the form
W =Wt(ι1, ι2) +Wc(ι3) resp. W =Wt(¯ι1,¯ι2) +Wc(¯ι3),
Chapter 5 Mechanical behavior of biological soft tissues depending of the employed set of invariants. The fact that the model based on the principal invariants does not properly separate isochoric and volumetric contribu-tions is often neglected in the modeling process. In many cases this is not a problem since the subsequent fitting to measurement data may at least partially correct these inaccuracies in the modeling process.
Remark 5.8.
• Another idea was followed in the development of Ogden’s models, cf. Ogden [199, 200, 201]. These are expressed in terms of the principal stretches, the eigenvalues of the strain tensor. The most popular Ogden-type models, the neo-Hookean and the Mooney-Rivlin material law, can also be expressed in terms of the principal strain invariants. Models which do not allow a formu-lation in terms of its invariants require the additional solution of eigenvalue problems and thus are not widely used.
• An interesting alternative approach was followed in Criscione et al. [68, 69].
There, use of the “K-invariants”, based on the Hencky, “true”, strain tensor ln(∇ΦT∇Φ) and orthogonality requirements, was proposed. The use of orthog-onal, physically meaningful, invariants admits the construction of constitutive relations that are more robust with respect to measurement errors than for-mulations based on (modified) principal invariants. Yet further examinations have to reveal merits and shortcomings of this approach.
5.1.2. Fiber-reinforced materials
The presence of fibers leads to anisotropic mechanical behavior of biological soft tissues. In view of the splitting
Wt=Wt,iso+Wt,aniso,
these properties are incorporated into hyperelastic models by augmenting an isotropic model with an anisotropic one.
The standard strategy for the derivation of anisotropic model is isotropization. In order to incorporate directional information we need a projection of the strain tensor onto the fiber direction v ∈ R3, |v| = 1. For d ∈ R3 the orthogonal projection on span(v) is given through
v(v, d) =vvTd= (v⊗v)d,
where (·,·) denotes the Euclidean scalar product inR3and⊗the Kronecker product.
Thus, the projection is characterized by the structural tensor
M :=v⊗v. (5.1.4)
This tensor must be incorporated in the definition of strain invariants that describe anisotropic materials. This leads to a new set of invariants.
120
5.1 Modeling framework
Definition 5.9. Let A ∈ M3 and M ∈S3. The corresponding mixed invariants
Definition 5.9. Let A ∈ M3 and M ∈S3. The corresponding mixed invariants