Muscle-Tendon Complex Model

44 Chapter 4: Constitutive Model

The advantage of choosing this formulation over the other ones is the fact that one obtains for the single component a C¹-continuous function. The qualitative stress contributions are depicted for uniaxial constraints in Figure 4.3.

0.6 0.8 1 1.2 1.4 1.6

−0.1 0 0.1 0.2 0.3 0.4 0.5 0.6

fibre stretch λ [−]

stress [N/mm2 ]

iso passive active total

Figure 4.3:A qualitative representation of the Cauchy stress contributions for the uniaxial case.

If a stress free initial configuration is assumed, the initial pressure can be determined by

S_muscle(χ, t= 0) = 2∂Ψ_muscle(I₁, I₂, I₄)

∂C −pC⁻¹ =0

= 2 (c₁+I₁c₂)I−2c₂I+0−pI.

(4.31) Using I₁(t = 0) = 3, (4.31) reduces to the initial, scalar-valued pressure

p(χ, t= 0) = 2c₁+ 4c₂. (4.32)

4.4 Muscle-Tendon Complex Model 45

For γ_ST = 1, the anisotropic part of the material behaviour vanishes. The soft tissue is assumed to behave as the isotropic part of the muscle tissue. Hereby, a relatively soft material is defined, where compressive forces will dominate due to the incompressibility behaviour of the tissue. For 0 < γ_ST < 1, the material exhibits anisotropic behaviour.

For γ_ST = 0 and γ_M = 1, the material behaves to 100% like skeletal muscle tissue and coincides with the material behaviour introduced in Section 4.3. For γ_ST = 0 and γ_M = 0, the material reduces to tendon tissue. For any other values between zero and one, a homogenised material behaviour of the mixture is obtained by linear interpolation.

Further, the different tissue material parameters are enforced by a linear interpolation:

c₁ = γ_Mc_1M + (1−γ_M)c_1T , c₂ = γ_Mc_2M + (1−γ_M)c_2T , c₃ = γ_Mc_3M + (1−γ_M)c_3T , c4 = γMc4M + (1−γM)c4T .

(4.35)

Herein, c_iM and c_iT are the material parameters defined in Section 4.3 for the skeletal muscle tissue and tendon tissue, respectively.

5 Finite Element Method in Space

Chapter 3 and 4 provide the fundamentals for the initial boundary value problem. It consists of the

• kinematic relations, providing the suitable strain measures (Section 3.1),

• definitions of stress (Section 3.1.3),

• local balance relations (Section 3.2), and

• constitutive equations for the skeletal muscle tissue (Chapter 4).

As quasi-static conditions are assumed, time contributions can be neglected. However, it is common to successively change the boundary conditions. Hereby, successively new BVP are declared to pretend a pseudo time-transient problem. Hence, all derivations within this chapter are carried out for one instant in time and the variabletindicating the current deformation is omitted. While body forces are assumed to be small in comparison to the internal forces, the strong formulation of the governing system reduces to:

divT=0. (5.1)

As it will be convenient for the contact formulation later on, Equation (5.1) is pulled back to the reference configuration, see Holzapfel et al. (2005). For the pull-back opera-tion, (5.1) is reformulated to its global form and the area and volume elements are pulled back using Equation (3.7). As the obtained formulation is independent of the choice of the volume element, the moment balance reads:

DivP(X) =0, where Div (·) := ∂(·)

∂X ·I. (5.2)

By including the constitutive behaviour and a complex geometry, the resulting system leads to a second-order partial differential equation for the current position, x = χ(X).

As this system cannot be analytically solved for most BVPs, an adequate numerical framework needs to be introduced. Here, the finite element method is used to discretise the governing equations and the geometry. More information can be found in standard text books of, for example, Bonet (1997) or Zhu et al. (2005).

5.1 Fundamentals

When the finite element method is applied to a boundary value problem, the procedure can be classified into several distinct steps. Each of them will be briefly outlined. The numerical treatment is realised within CMISS and outlined in Section 5.2.

(i) Form the weak formulation based on the governing equations, here given in Equa-tion (5.2),

47

48 Chapter 5: Finite Element Method in Space

(ii) chose a suitable finite element ansatz (basis functions, Galerkin weights),

(iii) evaluate element integrals (calculate the element stiffness matrices and right-hand side vectors),

(vi) assemble the global matrix, (v) apply boundary conditions, and (vi) solve the global system.

The first step is to formulate the weak formulation of the current problem. This step is important, as it is often impossible to find the correct solution of the strong formulation for the primary variables of the PDE using efficient and suitable discretisation schemes. The primary variables are the unknowns, i.e. in our case the current position,x=χ(X) and the hydrostatic pressure, p, arising from the incompressibility condition of Equation (4.20).

The generalised assembled strong or local form of the PDEs can be represented by:

DivP(X) = 0,

pI₃(X)−1 = 0. (5.3)

To derive the weak formulation, the PDEs are multiplied by test functionsδχandδpand integrated over the domain Ω₀ resulting in the following equivalent global system,

G^χ = Z

Ω0

DivP(X)·δχdV = 0, G^p =

Z

Ω0

(J(X)−1)δp dV = 0.

(5.4)

Since Equation (5.3) has to be true, Equation (5.4) is true for any choice of δχ and δp.

By integrating Equation (5.4)₁ by parts, the order of the operator of the primary variable is lowered by one and the order of the operator of the test function is increased by one:

G^χ = Z

Ω0

DivP(X)·δχdV

= Z

Ω0

P(X)·Grad δχdV − Z

Γ0

t₀·δχdA= 0.

(5.5)

As a result of the integration by parts, G^χ now contains a volumetric part and a part associated with the boundary, at which the initial traction vectort₀is applied asNeumann boundary condition. Equation (5.4)₂ of G^p is left untouched.

As the following procedure is in principle the same for both equations in Equation (5.4), the following procedure will be exemplarily explained using Equation (5.4)₁.

In order to obtain bounded results for the variational formulation, the first derivatives of the trial functions,S_χand test functionsS_p need to be square-integrable on the domain Ω₀. This requirement is satisfied, if the functions are assumed to be H¹(Ω).

S := { χ ∈ H¹(Ω)^d : χ(X) = χ(X)¯ on Γ},

T := {δχ ∈ H¹₀(Ω)^d : δχ(X) = 0 on Γ}, (5.6)

5.1 Fundamentals 49

Herein, d ∈ {1,2,3} is an integer denoting the dimensions in space and ¯χ denote the value of the Dirichlet boundary condition as the Dirichlet boundary conditions need to be satisfied exactly.

The defined trial spaces, S, and test spaces, T, have to be approximated by an N-dimensional subspace S^h and T^h, respectively. The discrete trial and test function are defined by:

χ(X) ≈ χ^h(X) = χ¯^h(X) +

N

X

i=1

φ^h_i(X)C_i ∈ S, δχ(X) ≈ δχ^h(X) =

N

X

i=1

ψ^h_i(X)δC_i ∈ T ,

(5.7)

where ¯χ^h denotes the Dirichlet boundary condition, φ^h_i and ψ^h_i are the global basis functions of the respective ansatz and test functions for the node i, which only depend on X, N indicates the total number of nodal points in the mesh, and C_i and δC_i are the vector-valued coefficients of dimension R^N×1 of the discretised geometrical primary variables.

The basis functions are linearly independent and chosen in such a way, that the nodal value is equal to one at the node itself, while the nodal value for all surrounding nodes are equal to zero.

φ^h_i(X_j) = δ_ij with i, j = 1, . . . , N . (5.8)

With this property, the coefficients of the discretised primary variables coincide with the nodal values of the primary variables, χ(X_i) and the nodal values of the test function, δχ_i. Hence,

C_i ≡χ(X_i), and δC_i ≡δχ(X_i). (5.9)

With this definition at hand, Equation (5.5) is discretised. As the number of unknowns equals the number of equations, the discretised PDE is, in principle, ready to be solved.

Yet, the integration of Equation (5.5) still poses a problem. The finite element method provides a widely used solution to this issue. By partitioning the solution space of the problem into a sum of small problems, Equation (5.5) is formulated with respect to a subdomain Ω_e of finite size. The domain Ω₀ is approximated by a spatially discretised domain Ω^h, which consists of E discretised, non-overlapping elements of finite size, Ω^h_e,

Ω₀ ≈Ω^h =

E

[

e=1

Ω^h_e. (5.10)

The discretised small domains, Ω^h_e, are called finite elements. The union of Ω^h_e forms the FE mesh. The elements are spanned by the nodes N_e, belonging to the element Ω^h_e. As the choice of the basis functions to represent the different primary variables does not have to be the same, the number of nodes, N_e, where the primary variable are defined, can vary. The number of nodes for the spatial variable and for the hydrostatic pressure is N_e^χ and N_e^p, respectively.

50 Chapter 5: Finite Element Method in Space

To discretise the nodal values within an element, the isogeometric mapping is employed.

The basic idea of the isogeometric mapping is to map a discretised, arbitrary element, Ω^h_e, into a reference element ˆΩ. If the element geometry is transformed by a function of the same polynomial order as the ansatz function, the mapping describes the isopara-metric concept. If the order of the geoisopara-metric transformation is higher, one refers to it as super-parametric. In general, for each independent variable a different parametric trans-formation can be chosen. Within the scope of this work, the geometry and the current configuration use quadratic ansatz functions and an isoparametric concept is chosen for the placement function. The pressure is assumed to change linearly within an element resulting in a super-parametric concept for the pressure. The reference element is defined by

N_e : ˆΩ→[0,1]^d, (5.11)

where N_e is the set of nodal basis functions for an element e, N_e :=

N₁(ξ), . . . , N_N^DoF

e (ξ)T

, ∈R^N

eDoF×1

, (5.12)

whereξ are the local coordinates of the reference element. Furthermore, the nodal values of the nodes belonging to an element e are defined by

X_e :=

X_e1, . . . ,X_eN^DoF

e

, ∈R^3×N

eDoF . (5.13)

to utilise the isogeometric mapping,

Φe: ˆΩ→Ω^h_e , (5.14)

ξ 7→X_e N_e(ξ), with e= 1, . . . , E . (5.15)

The mapping satisfies,

Ω^h_e = Φ_e( ˆΩ), with e= 1, . . . , E . (5.16)

and has the Jacobian J_e=∇Φ_e(ξ) =X_e







∇N₁(ξ) ...

∇N_N^DoF

e (ξ)





=X_e ∇N_e, ∈R^3×3. (5.17)

Note, the basis functions, φ^h_i, consist of contributions from N_e defined on different elements.

With the definitions (5.11) to (5.15), the ansatz functions of an element Ω^h can now be specified by choosing the appropriate nodal contributions belonging to this particular element.

Separating the variational formulation into a sum over small domains, the problem is defined by

f : Ω₀ →R with f see G^χ,G^p from (5.4), (5.18)

Z

Ω0

f dV ≈ Z

Ω^h

f dV ^(5.10)=

E

X

e=1

Z

Ω^h_e

dV_e

=

E

X

e=1

Z

Ω^h_e^(5.16)= Φe( ˆΩ)

f dV_e =

E

X

e=1

Z

Ωˆ

f(Φ_e)|det(Φ_e)|dV_e. (5.19)

5.1 Fundamentals 51

Herein, the first part is the exact statement and f is a placeholder for the variational formulations. With the first step, the solution is approximated by the discretised solution space. The integration of the discretised domain is split into a sum of many small domains using Equation (5.10). Each of theE elements is mapped by the element routine from Ω^h_e to the reference element ˆΩ by the isogeometric mapping, which is visualised in Figure 5.1.

ξ₂

ξ₁ x₂

x₁

Φ_e Φ⁻_e¹

Ω^h_e Ωˆ

global

Figure 5.1: Geometric transformation of an arbitrary rectangular 2-d element with quadratic basis functionsΩ^h_e to the reference elementΩ, defined in local coordinatesˆ ξ.

Within the reference element, a numerical integration scheme, such as the Gaussian quadrature scheme, is employed to compute the respective integrals. The numerical inte-gration, within one element, is carried out usingK Gauss points. Each pointkis weighted by theGaussian quadrature weights,w_k, and the sum of the weighting factors is one. The approximation of the integration for one element reads:

(5.19)≈

E

X

e=1 K

X

k=1

w_k f(Φ_e)|det(Φ_e)| (5.20)

Due to the fact that this integration is carried out on the element level with respect to the spatial reference element, the integration points and the weighting factors are always the same.

Dirichlet and Neumann boundary conditions are considered in the global system. For Dirichlet boundary conditions, the nodes are known. Hence, its contribution can be transferred to the right-hand side and the matrix can be reduced. For the Neumann boundary conditions, the respective force contribution is directly added to the right-hand side.

To solve the global problem, the element independent variable vectors are defined in 3 dimensions as:

u_e :=

χ^h(X₁), . . .χ^h(X_Ne^χ)T

∈ R^N

χ e×3, v_e :=

p^h(X₁), . . . p^h(X_Ne^p)T

∈ R^N

p e , w_e :=

u^T_e,v^T_eT

,

(5.21)

Even though the dimension of thew_eis not consistent, when a mixed-element formulation is used, the primary unknowns are stringed together to define the vectors, w_e on the

52 Chapter 5: Finite Element Method in Space

element level, and w on the global level. The global solution variable vectors are defined to:

u := [u₁, . . . ,u_Ne^χ]^T ∈ R^Nχ×3, v := [v₁, . . . ,v_Ne^p]^T ∈ R^N

p,

w := h

u^T₁, . . . ,u^T_N^χ

e, v^T₁, . . . ,v^T_N^p

e

iT

.

(5.22)

After the PDE is spatially discretised, the emerging system is of size 3N^χ+N^p and the emerging global system of equations can be written as

G_χ = K(u)−f = 0,

G_p = g(v) = 0, (5.23)

where (5.23)₂is the algebraic side condition emerging from the incompressibility condition.

Within the framework of finite elasticity, the emerging equations are still nonlinear.

Hence, linearisation is necessary. The Newton-Raphson method has proven to efficiently solve the nonlinear set of equations. The Newton-Raphson method requires the first derivative or tangent of the nonlinear functions. As the functions depend on more than just one variable, the directional derivative has to be applied, i.e. the observed function is derived with respect to each variable keeping the other variable constant at a distinct lo-cation. The tangent can be computed numerically or analytically. Within this framework, the constitutive equations are linearised numerically while the contact formulation uses an analytic tangent, cf. Section 6.1.3. As the numerical tangent is, within this framework, always determined in the reference configuration, the balance of momentum is used in the reference configuration.

Im Dokument A 3D continuum-mechanical model for forward-dynamics simulations of the upper limb (Seite 65-73)