$$\varepsilon^{\rm vK}\left(\mathbf{v},\mathbf{R}\right) = \int\_{\Omega} \left[\frac{1}{24}\mathcal{Q}\_2\left(\nabla^2\mathbf{v}\right) - \mathbf{R}\_{33}\mathbf{f}\_3 \cdot \mathbf{v}\right] d\mathbf{x}.\tag{34}$$

*Moreover, the non-linear strain satisfies*

() () () () 1 1 = ,= , *h hh h h h*

**Theorem 2**. *(iv) (Linearized isometry constraint). Suppose* 2 < < 3 *and set*  = 2 2*,*

( ) = 0. *T*

2 33 3

é ù Ñ- × ê ú

*(v) (von Kármán theory). Suppose that*  = 3 *and set*  = 4*,*  = 2*,*  = 1*. Then and for a subsequence of a -minimizing sequence, (1) and (2) hold and the limit triple* ,, *minimizes*

2 2

*Q Q d*

é ù æ ö é ù Ñ + Ñ +Ñ ÄÑ + Ñ - ê ú ç ÷ ë û ë û è ø

( ) ( ) 2

1 E ,= , 24

*lin Q d* W

( ) ( ) ( ) 2

**uvR u u vv v**

1 1 1 E ,, = 2 2 24

.

*x*

*d*

 d

**u Uv V** (27)

2*, suppose in addition that* Ω *is simply connected. Then,*

( ) , **R R** *h* ® (28)

Ñ + Ñ +Ñ ÄÑ **u u vv** (31)

ë û ò **v R v Rf v** *x* (32)

*x*

(29)

(30)

(33)

(ℎ) *is a -minimizing sequence then there exist constant* (ℎ) ∈ SO(3) *and c(h)*

g

respectively.

146 Perusal of the Finite Element Method

Then, we have

ℎ ≤ inf

 *such that*

*subject to* det∇2 = 0.

*von Kármán functional*

∈ ℝ*3*

= 2( 2)*,* 2 = *. If* 2 < < 5

ℎ ≤ 0*. If*

*Moreover, the pair* , *minimizes the functional*

*vK*

33 3

ò

W

**R fv**

*vK T*

W

ò

$$\frac{1}{h^{1-\alpha}} \left( \left[ \left( \nabla\_h \mathbf{r}^{(h)} \right)^T \nabla\_h \mathbf{r}^{(h)} \right]^{\frac{1}{2}} - \text{Id} \right) \to \text{sym} \, \mathbf{G}. \tag{35}$$

*In all cases, we have convergence of the rescaled energy to the minimum of the limit functional* . *Moreover, for* 3 ≡0 *we have* 33 = 1 *or* 1.

**Remark 1. (***i) Practically, Theorems 1 and 2 establish -convergence and give the -limits of threedimensional non-linear energy functional in the range*  0,5 3 ∪ 2,∞ *. The range*  5 3,2 *remains unexplored*.

*(ii) Even in the case* = 5 3, *it is still possible to find a sequence* (ℎ) *for which the functional is bounded when* ℎ 0. *The limiting deformations are composed of finitely many affine isometries separated by sharp folds and are called origami deformations [46]. It is proved ibid that the class of origami deformations is the closure of the class of short deformations with respect to the uniform convergence*.

*(iii) In particular, when*  = = 2*, for isotropic homogeneous bodies we obtain*

$$\mathcal{Q}\_z(\mathbf{II}) = 2\mu \left| \mathbf{II} \right|^2 + \frac{\lambda \cdot 2\mu}{\lambda + 2\mu} \text{tr}^2 \mathbf{II},\tag{36}$$

*in which and are Lamé constants, coinciding with plate energy derived by Kirchhoff much earlier in [49]*.

## **4. Relaxation of epithelial elastic energy by Pipkin's procedure and by adding the bending contribution**

As was mentioned in Section 1, there are two ways to relax the epithelial elastic energy. The first way is to relax the stretching energy for deformation regimes corresponding to compressive stresses. For that reason, Pipkin's procedure [41] seems to be the most common tool.

Suppose that the stretching energy density function = ∇ is frame indifferent, isotropic and positive except the cases

$$
\nabla \mathbf{r} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ 0 & 0 \end{pmatrix} \text{and } \nabla \mathbf{r} = \begin{pmatrix} \mathcal{S} & 0 \\ 0 & \mathcal{S} \\ 0 & 0 \end{pmatrix}, \quad 1 \le \mathcal{S} \in \mathbb{R}, \tag{37}
$$

where it turns to zero, that is, attains its minimum. On elastic energies with two minima, see [50]. This assumption is motivated by bistability of cell shapes during apical constriction (see **Figure 4**).

We restrict attention to deformations with

$$
\nabla \mathbf{r} = \begin{pmatrix} \sigma\_1 & 0 \\ 0 & \sigma\_2 \\ 0 & 0 \end{pmatrix}, \tag{38}
$$

and denote ∇ : = (1,2). Here = , 1; 2 , are the principal stretches, the squares of which are the eigenvalues of the first fundamental form. Therefore

$$\mathcal{A}\_{\mathbf{r},2}\left(\nabla\mathbf{r}\right) = \sqrt{\frac{1}{2}\left[\operatorname{tr}\left(\mathbf{I}\left(\nabla\mathbf{r}\right)\right) \pm \sqrt{\left\|\mathbf{I}\left(\nabla\mathbf{r}\right)\right\|^2 - 2\det\left(\mathbf{I}\left(\nabla\mathbf{r}\right)\right)}\right]}\,\tag{39}$$

In view of frame indifference of , and thence, , we conclude that ( ± 1, ± 2) = (1,2).

From quasi-convexity of , in particular, follows its rank-one convexity, therefore since the matrix

Improved Reduced Order Mechanical Model and Finite Element Analysis of Three-Dimensional Deformations... http://dx.doi.org/10.5772/65027 149

$$\mathbf{A} = \begin{pmatrix} \delta\_{11} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} \end{pmatrix}, \quad \delta\_{11} \neq \mathbf{0}, \tag{40}$$

is rank one, we have

**4. Relaxation of epithelial elastic energy by Pipkin's procedure and by**

As was mentioned in Section 1, there are two ways to relax the epithelial elastic energy. The first way is to relax the stretching energy for deformation regimes corresponding to compressive stresses. For that reason, Pipkin's procedure [41] seems to be the most common tool.

Suppose that the stretching energy density function = ∇ is frame indifferent, isotropic

= 0 1 and = 0 , 1 < R,

d

where it turns to zero, that is, attains its minimum. On elastic energies with two minima, see [50]. This assumption is motivated by bistability of cell shapes during apical constriction (see

1

( ) ( ( )) ( ) ( ( )) 2

é ù Ñ Ñ± Ñ - Ñ ê ú

= tr 2det .

In view of frame indifference of , and thence, , we conclude that

From quasi-convexity of , in particular, follows its rank-one convexity, therefore since the

ë û **r Ir Ir Ir** (39)

æ ö ç ÷ Ñ ç ÷ ç ÷ è ø

s

squares of which are the eigenvalues of the first fundamental form. Therefore

2 0 =0 , 0 0

s

d

 d

**r r** (37)

**r** (38)

, 1; 2 , are the principal stretches, the

1 0 0

æö æ ö ç÷ ç ÷ ÑÑ Î

00 00

èø è ø

**adding the bending contribution**

We restrict attention to deformations with

and denote ∇ : = (1,2). Here =

1

2

1,2

l

( ± 1, ± 2) = (1,2).

matrix

and positive except the cases

148 Perusal of the Finite Element Method

**Figure 4**).

therefore is convex with respect to 1. Repeating the argument when

$$\mathbf{A} = \begin{pmatrix} 0 & 0 \\ 0 & \delta\_{22} \\ 0 & 0 \end{pmatrix}, \quad \delta\_{22} \neq 0,\tag{41}$$

we see that (1,2) is separately convex.

Under assumptions made, we have (1,2)=0 if 1 = 2 = 1 and 1 = 2 = . Since > 1 and 0 ≤ are convex along lines parallel to 1 = 0 and 2 = 0, therefore = 0 on edges connecting the vertexes ±, ± in the 12 plane. Repeating the argument, we conclude that = 0 inside the square Ω0 = 1,2 ; 1 ≤ , 2 ≤ . So, for -short maps ≡ 0.

Since (1,2) is even and convex in 2, it will attain its minimum with respect to 2 only when 2 = 0. Denote by 0(1) = min 2 (1,2) when 1 > . Suppose that 0(1) is the largest value of 2 for which (1,0) = 0(1), 1 > . The convexity of in 2 implies that (1,2) = 0(1) for 2 ≤ 0(1), 1 > . In the same way, using the symmetry of in its arguments, we will arrive at (1,2) = 0(2) for 1 ≤ 0(2), 2 > . Evidently, for the rest of deformation regimes, is convex.

If we denote by Ω1 = 1,2 ; 1 > , 2 ≤ 0(1) and Ω2 = 1,2 ; 1 ≤ 0(2), 2 > , the ranges of principal stretches, corresponding to uniaxial tensions in directions of those stretches, respectively [51], we finally arrive at

$$
\tilde{\boldsymbol{W}}\_{q}(\boldsymbol{\lambda}\_{1}, \boldsymbol{\lambda}\_{2}) = \begin{cases} \boldsymbol{0}, & \text{in } \Omega\_{0}, \\ \tilde{\boldsymbol{W}}^{0}(\boldsymbol{\lambda}\_{1}), & \text{in } \Omega\_{1}, \\ \tilde{\boldsymbol{W}}^{0}(\boldsymbol{\lambda}\_{2}), & \text{in } \Omega\_{2}, \\ \tilde{\boldsymbol{W}}(\boldsymbol{\lambda}\_{1}, \boldsymbol{\lambda}\_{2}), & \text{else}. \end{cases} \tag{42}
$$

Besides opportunity of explicit relaxation, in [41] the derivation of criteria for convexity and quasi-convexity is described. In view of rank-one convexity, satisfies Legendre-Hadamard inequality from which it follows

$$
\tilde{W}\_q^1 \ge 0, \quad \tilde{W}\_q^2 \ge 0, \quad \tilde{W}\_q^{11} \ge 0, \quad \tilde{W}\_q^{22} \ge 0,\tag{43}
$$

$$G \ge 0, \quad \sqrt{\tilde{W}\_q^{11} \tilde{W}\_q^{22}} - \tilde{W}\_q^{12} \ge H - G, \quad \sqrt{\tilde{W}\_q^{11} \tilde{W}\_q^{22}} + \tilde{W}\_q^{12} \ge -H - G. \tag{44}$$

For convexity of in corresponding domain instead of the last two equalities, the stronger conditions

$$\left(\tilde{W}\_q^{11}\tilde{W}\_q^{22} - \left(\tilde{W}\_q^{12}\right)^2\right) \ge 0, G \succeq \left| H \right|, \tag{45}$$

must be satisfied. Above

$$
\tilde{W}\_q^{'} = \frac{\partial \tilde{W}\_q}{\partial \lambda\_j}, \quad \tilde{W}\_q^{'k} = \frac{\partial^2 \tilde{W}\_q}{\partial \lambda\_j \partial \lambda\_k}, \quad G = \frac{\lambda\_l \tilde{W}\_q^{1} - \lambda\_z \tilde{W}\_q^{2}}{\lambda\_1^2 - \lambda\_2^2}, \quad H = \frac{\lambda\_z \tilde{W}\_q^{1} - \lambda\_l \tilde{W}\_q^{2}}{\lambda\_1^2 - \lambda\_2^2}. \tag{46}
$$

The other way is the adding of pure bending contribution to the elastic energy of the tissue. Since actin-myosin network contraction leads to compressive or non-stretching stresses, we have to incorporate the elastic energy mainly for such deformations. According to Section 3 for non-stretching stresses, we have

**i.** (short deformations and compressive stresses)

$$\Gamma - \lim\_{h \to 0} \frac{1}{h^{\beta}} E^{(h)} = 0 \text{ in the range } \beta \in \left[0, \frac{\mathfrak{F}}{\mathfrak{Z}}\right), \tag{47}$$

**ii.** (isometric deformations, neither stretch nor compressive stresses)

$$
\Gamma - \lim\_{h \to 0} \frac{1}{h^2} E^{(h)} = E\_b. \tag{48}
$$

The deformation regime Id corresponding to tensile stresses leads to , which is already lower semi-continuous. In the case of uniaxial tension, there is no bending and the energy expression will look like the middle rows of (42). Thence, the total energy might be given as follows:

$$\mathbf{E}\left[\mathbf{r}\right] = -\int\_{\Omega} \mathbf{f} \cdot \mathbf{r} \, d\mathbf{x} + \frac{1}{2} \int\_{\Omega} \tilde{W}^{0} \left(\lambda\_{i}(\nabla \mathbf{r})\right) d\mathbf{x}, \text{ for unit-axial tension in the direction of } \lambda\_{i},$$

$$\mathbf{E}\left[\mathbf{r}\right] = -\int\_{\Omega} \mathbf{f} \cdot \mathbf{r} \, d\mathbf{x} + \frac{1}{2} \int\_{\Omega} \tilde{W}^{0} \left(\lambda\_{z}(\nabla \mathbf{r})\right) d\mathbf{x}, \text{ for unit-axial tension in the direction of } \lambda\_{z}, \tag{49}$$

$$\begin{aligned} \frac{1}{24} \int\_{\Omega} Q\_{z}(\mathbf{H}) \, d\mathbf{x}, \text{ I = Id,} \\ \frac{1}{2} \int\_{\Omega} W\_{\neq}(\nabla \mathbf{r}) \, d\mathbf{x}, \text{ otherwise.} \end{aligned} \tag{40}$$

We assume that the tissue is homogeneous and isotropic, so according to Remark 1

Besides opportunity of explicit relaxation, in [41] the derivation of criteria for convexity and quasi-convexity is described. In view of rank-one convexity, satisfies Legendre-Hadamard

*qq q q* ³³ ³ ³ (43)

*qq q* - ³³ (45)

 ll

> l l

*qq q* + ³- - (44)

1 2 11 22 0, 0, 0, 0, *WWW W* %% % %

11 22 12 11 22 12 0, , . *G WW W H G WW W H G* ³ %% % *qq q* - ³- %% %

> ( ) 2 11 22 12 0, | |, *WW W G H* %% %

= , = ,= , = . *j q jk q q q qq*

¶¶ - - ¶ ¶¶ - -

*W W WW WW*

 l l

The other way is the adding of pure bending contribution to the elastic energy of the tissue. Since actin-myosin network contraction leads to compressive or non-stretching stresses, we have to incorporate the elastic energy mainly for such deformations. According to Section 3

> 1 5 lim = 0 in the range 0, , 3

> > ( )

The deformation regime Id corresponding to tensile stresses leads to , which is already lower semi-continuous. In the case of uniaxial tension, there is no bending and the energy expression will look like the middle rows of (42). Thence, the total energy might be given

2 0 1 lim = . *h b h E E*

® *h*

ll

For convexity of in corresponding domain instead of the last two equalities, the stronger

2 1 2 12

% % % % %% % % (46)

b

é ö G - Î ê ÷ ë ø (47)

G - (48)

12 21 2 2 2 2 1 2 1 2

inequality from which it follows

150 Perusal of the Finite Element Method

conditions

as follows:

must be satisfied. Above

*q q*

l

for non-stretching stresses, we have

*j j k*

**i.** (short deformations and compressive stresses)

0

®

*h*

 ll

*WW G H*

( )

*E h*b

*h*

**ii.** (isometric deformations, neither stretch nor compressive stresses)

$$\mathcal{Q}\_z(\mathbf{II}) = 2\mu \left| \mathbf{II} \right|^2 + \frac{\lambda \cdot 2\mu}{\lambda + 2\mu} \text{tr}^2 \mathbf{II}. \tag{50}$$

Furthermore, the membrane energy density function we take in the following form [52, 53]:

$$W\left(\nabla \mathbf{r}\right) = \gamma\_1 \left(I\_1 - 3\right) + \gamma\_2 \left(I\_1 - 3\right)^2,\tag{51}$$

which is valid for large deformations of incompressible hyper-elastic membranes. Above, 1 = tr is the first strain invariant, 1 > 0 and 2 > 0 are material parameters, evaluated experimentally in [53]. For a neo–Hookean material 21 = and 22 = , in which is the shear modulus and is the bulk modulus. and 0 are determined according to the procedure described above. For that reason, 1 must be represented in terms of the principal stretches 1 and 2 as follows:

$$I\_1(\boldsymbol{\lambda}\_1, \boldsymbol{\lambda}\_2) = \boldsymbol{\lambda}\_1^2 + \boldsymbol{\lambda}\_2^2 + \boldsymbol{\lambda}\_3^2 = \boldsymbol{\lambda}\_1^2 + \boldsymbol{\lambda}\_2^2 + \frac{1}{\boldsymbol{\lambda}\_1^2 \boldsymbol{\lambda}\_2^2},\tag{52}$$

in which the incompressibility condition 123 = 1 is taken into account.

For any fixed 1 > 1 (uniaxial tension in the first principal direction), the corresponding 1,2 attains its minimum with respect to 2 at the maximal positive root of

$$\frac{\partial \tilde{W}}{\partial \lambda\_2} = 2\lambda\_2 \left( 1 - \frac{1}{\left(\lambda\_1^2 + \lambda\_2^2\right)^2} \right) \left[ \gamma\_1 + 2\gamma\_2 \left( I\_1 - 3\right) \right] = 0. \tag{53}$$

Since in the uniaxial tension regime 1 2 + 2 2 > 1, we finally get

$$
\mathcal{N}^0(\mathcal{A}\_1) = \frac{1}{\sqrt{\mathcal{A}\_1}}.\tag{54}
$$

Then,

$$W^0(\mathcal{J}\_\mathfrak{l}) = \frac{1}{\mathcal{J}\_\mathfrak{l}^2} (\mathcal{J}\_\mathfrak{l}^3 - 3\mathcal{J}\_\mathfrak{l} + 2) \left[ \mathcal{Y}\_\mathfrak{l}\mathcal{J}\_\mathfrak{l} + \mathcal{Y}\_2(\mathcal{J}\_\mathfrak{l}^3 - 3\mathcal{J}\_\mathfrak{l} + 2) \right]. \tag{55}$$

The correspondent tension is defined by

$$\pi\_1(\boldsymbol{\lambda}\_t) = \frac{\partial W^0(\boldsymbol{\lambda}\_t)}{\partial \boldsymbol{\lambda}\_t} = 2 \frac{\boldsymbol{\lambda}\_1^3 - 1}{\boldsymbol{\lambda}\_1^3} [\boldsymbol{\gamma}\_1 \boldsymbol{\lambda}\_t + 2 \boldsymbol{\gamma}\_2 (\boldsymbol{\lambda}\_t^3 - 3 \boldsymbol{\lambda}\_1 + 2)]. \tag{56}$$

Since is obviously symmetric in its arguments, in the case of uniaxial tension in the second principal direction, we would have

$$W^0(\mathcal{J}\_2) = \frac{1}{\lambda\_2^2} (\mathcal{\lambda}\_2^3 - 3\mathcal{\lambda}\_2 + 2)[\mathcal{\chi}\_1\mathcal{\lambda}\_2 + \mathcal{\chi}\_2(\mathcal{\lambda}\_2^3 - 3\mathcal{\lambda}\_2 + 2)],\tag{57}$$

and

$$\pi\_2(\boldsymbol{\lambda}\_2) = \frac{\partial W^0(\boldsymbol{\lambda}\_2)}{\partial \boldsymbol{\lambda}\_2} = 2 \frac{\boldsymbol{\lambda}\_2^3 - 1}{\boldsymbol{\lambda}\_2^3} [\boldsymbol{\gamma}\_1 \boldsymbol{\lambda}\_2 + 2 \boldsymbol{\gamma}\_2 (\boldsymbol{\lambda}\_2^3 - 3 \boldsymbol{\lambda}\_2 + 2)].\tag{58}$$

It is evident, that when we substitute 1 = 1 in the first case and 2 = 1 in the second case, we arrive at undeformed configuration of the membrane, so the relaxed energy is fully consistent.

Forces driving tissue deformation are strongly localized and in general are compressive. Force must possess those properties.

## **5. Finite element analysis of three-dimensional shape of some characteristic deformations of tissues during morphogenesis**

In this section, we summarize main results of finite element analysis of a single layer tissue model, the elastic energy of which is given by

$$E\left[\mathbf{r}\right] = hE\_{\text{w}}\left[\mathbf{r}\right] + h^{\text{3}}E\_{b}\left[\mathbf{r}\right] \tag{59}$$

in which and are defined in Eqs. (49)–(51). Since the discretization procedure is standard, we omit the details and bring only the main results.

In **Figure 6**, we bring the model of a single-cell (element) and cell-cell junction (in red). All structures (plate and shell) considered in this section entirely consist of such cell groups. In all tissues considered below, the height of a cell ℎ=4 and ℎ = 200, in which is its side and is the diameter of actin fibres. The ratio of Young's moduli of a cell and actin fibres is act = 3cell, and of actin fibres and cell-cell links act = 1.5link.

We consider

Since in the uniaxial tension regime 1

152 Perusal of the Finite Element Method

Then,

and

2 + 2

0 1

0 3 3 1 1 1 11 2 1 1 2

 l

0 3

¶ -

l

l

2

l

 l

l

1 1

 l

0 3 3 2 2 2 12 2 2 2 2

 l

0 3

¶ -

 l

> l

l

l

1 *W* ( ) ( 3 2)[ ( 3 2)],

 gl g l

2 2 3 2 2 3 12 2 2 2 2 2

gl

It is evident, that when we substitute 1 = 1 in the first case and 2 = 1 in the second case, we arrive at undeformed configuration of the membrane, so the relaxed energy is fully consistent.

Forces driving tissue deformation are strongly localized and in general are compressive. Force

 g l

() 1 ( ) 2 [ 2 ( 3 2)]. *W*

= = + -+

 l

1 *W* ( ) ( 3 2) [ ( 3 2)].

 gl g l

1 1 3 1 1 3 11 2 1 1

gl

Since is obviously symmetric in its arguments, in the case of uniaxial tension in the second

 g l

() 1 ( ) 2 [ 2 ( 3 2)]. *W*

= = + -+

1

l

 l

l

The correspondent tension is defined by

t l

principal direction, we would have

t l

must possess those properties.

l l() . l

2 > 1, we finally get

= (54)

 l

= -+ + -+ (55)

 l

 l

> l

¶ (58)

= -+ + -+ (57)

¶ (56)

1 1

**Figure 6.** Single-cell model (left) and adherens junction model (right): diagonals of the top hexagon imitate actin fibres, red areas between cells imitate junction bonds.

**Figure 7.** Apical constriction of initially flat tissue. It consists from 5050 elements or cells and has 529,245 DOFs.

**Figure 8.** von Mises stress distribution in deformed configurations.

**Figure 9.** Apical constriction of initially cylindrical tissue. It consists from 5975 elements or cells and has 621,375 DOFs.

**Figure 10.** von Mises stress distribution in deformed configurations.

Improved Reduced Order Mechanical Model and Finite Element Analysis of Three-Dimensional Deformations... http://dx.doi.org/10.5772/65027 155

**Figure 11.** Apical constriction of initially spherical tissue. It consists from 3126 elements or cells and has 325,245 DOFs.

Elements of the middle part of the rectangular tissue are compressed in apical sides to imitate apical constriction in cells. Increasing the compressing stresses, the tissue is bent and a blastershaped pattern is formed as shown in **Figure 7** (right). The quantitative picture of the stresses arising in the tissue is drawn in **Figure 8**.

**Figure 12.** von Mises stress distribution in deformed configurations.

Next, we consider a cylindrical shell to imitate ventral furrow (generally all tubular patterns) formation. Elements of the top part of the cylinder are constrained in the apical sides by compressing the links standing for apical fibres (see **Figure 9**). Increasing the compressing stresses, the ventral furrow formation is simulated similar to stages presented in **Figure 1**. The quantitative picture of the stresses arising in the tissue is drawn in **Figure 10**.

Finally, a spherical shell is simulated. Cells at the top of a semi-sphere are constrained in apical sides and depending on values of compressing stresses various stages of blastopore formation in archenteron can be described (see **Figure 11** (right)). The quantitative picture of the stresses arising in the tissue is drawn in **Figure 12**.

## **6. Conclusions**

**Figure 7.** Apical constriction of initially flat tissue. It consists from 5050 elements or cells and has 529,245 DOFs.

**Figure 9.** Apical constriction of initially cylindrical tissue. It consists from 5975 elements or cells and has 621,375 DOFs.

**Figure 8.** von Mises stress distribution in deformed configurations.

154 Perusal of the Finite Element Method

**Figure 10.** von Mises stress distribution in deformed configurations.

Analysis based on lower semi-continuous energy functionals reveals real three-dimensional deformations of soft epithelial tissues. Having different forms for different deformation regimes, such as compressive stresses (short deformations), uniaxial tensions collinear to principal directions of the first fundamental form, no stretching or pure bending stresses (isometric deformations) and stretching stresses (large deformations), the resulting total energy is lower semi-continuous, so the existence of its minimizers, that is, real deformations, is ensured. Particular energy density functions are chosen and the explicit form of the total energy functional is obtained, thereby the discretization is made easy.

On the basis of obtained energy functional, a three-dimensional discretized model of epithelial tissues undergoing combined stretching and bending deformations is constructed. Discretization elements correspond to single cells forming the tissue. Actin fibres and cell-cell adhesion links, mainly contributing on the tissue energy, are explicitly embedded in elements. Deformations characteristic to specific embryonic tissues (ventral furrow, neutral tube, neurosphere) observed earlier are described quantitatively increasing contractile stresses in fibres.

## **Dedication**

We thankfully dedicate the chapter to the blessed memory of our good fellow and colleague, a candidate of physical and mathematical sciences, Hamlet V. Hovhannisyan (1956–2016), who unexpectedly died before he could realize his best scientific ideas.

## **Acknowledgements**

The theoretical part of the chapter is investigated under the guidance of Doctor, Professor Benedikt Wirth, Institute for Computational and Applied Mathematics, University of Münster, whom we are heartily thankful. The work of As. Kh. and S. O. was made possible in part by a research grant from the Armenian National Science and Education Fund (ANSEF) based in New York, NY, USA.

## **Nomenclature**

## **Author details**

principal directions of the first fundamental form, no stretching or pure bending stresses (isometric deformations) and stretching stresses (large deformations), the resulting total energy is lower semi-continuous, so the existence of its minimizers, that is, real deformations, is ensured. Particular energy density functions are chosen and the explicit form of the total

On the basis of obtained energy functional, a three-dimensional discretized model of epithelial tissues undergoing combined stretching and bending deformations is constructed. Discretization elements correspond to single cells forming the tissue. Actin fibres and cell-cell adhesion links, mainly contributing on the tissue energy, are explicitly embedded in elements. Deformations characteristic to specific embryonic tissues (ventral furrow, neutral tube, neurosphere)

We thankfully dedicate the chapter to the blessed memory of our good fellow and colleague, a candidate of physical and mathematical sciences, Hamlet V. Hovhannisyan (1956–2016), who

The theoretical part of the chapter is investigated under the guidance of Doctor, Professor Benedikt Wirth, Institute for Computational and Applied Mathematics, University of Münster, whom we are heartily thankful. The work of As. Kh. and S. O. was made possible in part by a research grant from the Armenian National Science and Education Fund (ANSEF) based in

observed earlier are described quantitatively increasing contractile stresses in fibres.

energy functional is obtained, thereby the discretization is made easy.

unexpectedly died before he could realize his best scientific ideas.

**Dedication**

156 Perusal of the Finite Element Method

**Acknowledgements**

New York, NY, USA.

**Nomenclature**

⋅ Inner (dot) product , Lamé coefficients

1, 2 Principal stretches

ℎ Rescaled three-dimensional total energy

(Ω) Measure of bounded open set Ω

∇ Gradient (nabla) operator

Ara S. Avetisyan, Asatur Zh. Khurshudyan\* and Sergey K. Ohanyan

\*Address all correspondence to<EMAIL_ADDRESS>
Department of Dynamics of Deformable Systems and Connected Fields, Institute of Mechanics, National Academy of Sciences of Armenia, Yerevan, Armenia

## **References**

[14] Du. X., Osterfield M., Shvartsman S. Y., Computational analysis of three-dimensional epithelial morphogenesis using vertex models. *Physical Biology*, 2014, vol. 11, issue 6, 15 p.

**References**

158 Perusal of the Finite Element Method

pp. 139–156.

446–462.

247–261.

3400.

2012, vol. 338, pp. 201–203.

*Biology*, 2010, vol. 341, pp. 5–19.

*Development*, 2010, vol. 137, pp. 1407–1420.

esis. *Developmental Biology*, 2010, vol. 341, pp. 114–25.

*Part C: Embryo Today*, 2007, vol. 81, issue 4, pp. 253–269.

[1] Rhumbler L., Zur Mechanik des Gastrulationsvorganges insbesondere der Invagination Eine entwickelungsmechanische Studie. *Archiv für Entwicklungsmechanik der*

[2] Lewis W. H., Mechanics of invagination. *The Anatomical Record*, 1947, vol. 97, issue 2,

[3] Odell G. M., Oster G., Alberch P., Burnside B., The mechanical basis of morphogenesis: I. Epithelial folding and invagination. *Developmental Biology*, 1981, vol. 85, issue 2, pp.

[4] Polyakov O., Bing H., Michael S., Joshua W. S., Matthias K., Wieschaus E., Passive mechanical forces control cell-shape change during *Drosophila* ventral furrow forma-

[5] Keller R., Developmental biology. Physical biology returns to morphogenesis. *Science*,

[6] Mammoto T., Ingber D. E., Mechanical control of tissue and organ development.

[7] Martin A. C., Goldstein B., Apical constriction: themes and variations on a cellular mechanism driving morphogenesis. *Development*, 2014, vol. 141, pp. 1987–1998.

[8] Sawyer J. M., Harrell J. R., Shemer G., Sullivan-Brown J., Roh-Johnson M., Goldstein B., Apical constriction: a cell shape change that can drive morphogenesis. *Developmental*

[9] Martin A. C., Pulsation and stabilization: contractile forces that underlie morphogen-

[10] Martin A. C., Kaschube B., Wieschaus E. F., Pulse contractions of an actin-myosin network drive apical constriction. *Nature*, 2009, vol. 457, issue 22, pp. 495–501.

[11] Jones G. W., Chapman S. J., Modelling apical constriction in epithelia using elastic shell theory. *Biomechanics and Modeling in Mechanobiology*, 2010, vol. 9, issue 3, pp.

[12] Ciarletta P., Ben Amar M., Labouesse M., Continuum model of epithelial morphogenesis during *Caenorhabditis elegans* embryonic elongation. *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 2009, vol. 367, pp. 3379–

[13] von Dassow M., Davidson L., Variation and robustness of the mechanics of gastrulation: the role of tissue mechanical properties during morphogenesis. *Birth Defects Research*

tion. *Biophysical Journal*, 2014, vol. 107, issue 4, pp. 998–1010.

*Organismen*, 1902, vol. 14, issue 3, pp. 401–476.

[45] Le Dret H., Raoult A., The nonlinear membrane model as variational limit of nonlinear three-dimensional elasticity. *Journal de Mathématiques Pures et Appliquées*, 1995, vol. 74, pp. 549–578.

[28] Dacorogna B., *Direct Methods in the Calculus of Variations*. Springer, 2009.

[29] Dal Maso G., *An Introduction to -Convergence*. Birkhäuser, Boston, 1993.

*Mechanics and Applied Mathematics*, 1992, vol. 45, issue 1, pp. 57–75.

*Applied Mathematics*, 1992, vol. L, issue 2, pp. 389–400.

*of Applied Mathematics*, 1986, vol. 36, pp. 85–99.

180, pp. 183–236.

*Journal of Applied Mathematics*, 1994, vol. 52, pp. 297–308.

[43] Ciarlet P. G., *Mathematical Elasticity*, in 3 vol. Elsevier, 1993.

*of Applied Mathematics*, 1996, vol. LIV, issue 2, pp. 307–316.

bending stiffness. *Journal of Elasticity*, 1993, vol. 31, pp. 125–139.

order. *IMA Journal of Applied Mathematics*, 1992, vol. 48, pp. 265–269.

integrals of any order. *Transactions of AMS*, 1965, vol. 119, pp. 125–149.

353-377 (III).

160 Perusal of the Finite Element Method

*Journal of Mechanics and Applied Mathematics*, 1989, vol. 42, pp. 427–440.

[30] Steigman D. J., Pipkin A. C., Finite deformations of wrinkled membranes. *Quarterly*

[31] Kohn R. V., Strang G., Optimal design and relaxation of variational problems I, II, III. *Communications on Pure and Applied Mathematics*, 1986, vol. 39, pp. 1–25 (I), 139–182 (II),

[32] Hilgers M. G., Pipkin A. C., Elastic sheets with bending stiffness. *Quarterly Journal of*

[33] Hilgers M. G., Pipkin A. C., Bending energy of highly elastic membranes. *Quarterly of*

[34] Hilgers M. G., Pipkin A. C., Bending energy of highly elastic membranes II. *Quarterly*

[35] Hilgers M. G., Pipkin A. C., Energy–minimizing deformations of elastic sheets with

[36] Hilgers M. G., Pipkin A. C., The Graves condition for variational problems of arbitrary

[37] Schmidt B., Fraternali F., Universal formulae for the limiting elastic energy of membrane networks. *Journal of the Mechanics and Physics of Solids*, 2012, vol. 60, pp. 172–180.

[38] Meyers N. G., Quasi-convexity and lower semi-continuity of multiple variational

[39] Dal Maso G., Fonseca I., Leoni G., Morini M., Higher-order quasi-convexity reduces to quasi-convexity. *Archive for Rational Mechanics and Analysis*, 2004, vol. 171, pp. 55–81.

[40] Cagnetti F., k-Quasi-convexity reduces to quasi-convexity. *Proceedings of the Royal Society of Edinburgh: Section A Mathematics*, 2011, vol. 141, issue 4, pp. 673–708.

[41] Pipkin A. C., The relaxed energy density for isotropic elastic membranes. *IMA Journal*

[42] Pipkin A. C., Relaxed energy densities for large deformations of membranes. *IMA*

[44] Friesecke G., James R. D., Muller S., A hierarchy of plate models derived from nonlinear elasticity by Gamma-convergence. *Archive for Rational Mechanics and Analysis*, 2006, vol.

#### **Simulating Contact Instability in Soft Thin Films through Finite Element Techniques Simulating Contact Instability in Soft Thin Films through Finite Element Techniques**

Jayati Sarkar, Hemalatha Annepu and Jayati Sarkar, Hemalatha Annepu and Satish Kumar Mishra

Satish Kumar Mishra

Additional information is available at the end of the chapter Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/65357

## **Abstract**

When a thin film of soft elastic material comes in contact with an external surface, contact instability triggered by interaction forces, such as van der Waals, engenders topologi‐ cally functionalized surfaces. Innumerable technological applications such as adhesives; microelecromechanical systems (MEMS), and nanoelectromechanical systems (NEMS) demand understanding of the physics behind the mechanical contact, relationship between the morphologies, and detachment forces in such films. Indentation tests are important experimental approach toward this; there also exist many simulation procedures to model the mechanical contact. Both atomistic level and analytical continuum simulations are computationally expensive and are restricted by the domain geometries that can be handled by them. Polymeric films also particularly demonstrate a rich variety of nonlinear behavior that cannot be adequately captured by the aforementioned methods. In this chapter we show how finite element techniques can be utilized in crack opening and in contact‐instability problems.

**Keywords:** ABAQUS, thin films, patterned substrate, VDLOAD, contact instability

## **1. Introduction**

In nature, there are several examples such as lotus leaves, peacock's feathers, butterfly wings, gecko feet, and moths' eyes where meso‐nano‐scaled physical patterns on the surface impart characteristic properties such as hydrophobicity, structural color, adhesiveness, and optical properties, which are absent otherwise, to the surfaces. To engineer such surfaces artificially and to fabricate such features, there are mainly two most commonly used approaches: one a

© 2016 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2016 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

top‐down approach and the other a bottom‐up approach. The top‐down approach consists of the various lithographic techniques such as electron beam lithography, dip‐pen lithography, embossing lithography, laser lithography [1], laser photoablation [2], ion implantation [3], micromachining [4], etc. The bottom‐up approach, on the other hand, consists of self‐assembly of molecules to produce bigger building blocks to a growing structure as in chemical vapor deposition [5]. While the self‐assembly procedure involves numerous steps, lithographic methods such as dip‐pen lithography, which are used for creating patterns on hard silicon surfaces, are time‐consuming, involve great cost and may not at all be suitable to create patterns on soft surfaces such as polymers. Self‐organization methods, on the other hand, are for bulk production of desired structures involving lower cost and therefore are more potentially applicablefabricationtechniquesonsoft surfaces.Inself‐organizationtechnique,materials such as polymeric films, which have highly tunable mechanical properties, can be used to obtain a desired mesoscopic structure due to the influence of internal or external forces. The mechanical and interaction properties of these films can be controlled by the amount of cross‐linker added during the synthesis step itself. The instabilities during self‐organization can be harnessed to engineer desired periodic arrangements of sub‐micron features in technological applications such as in flexible electronics [6], optoelectronics [7], functional coatings [8], sensors [9], microfluidic devices, microelecromechanical systems (MEMS) and nanoelectromechanical systems (NEMS) [10], pressure‐sensitive adhesives [11], high‐efficiency light‐emitting diodes, thin‐filmtransistors,etc.Thus,thesethinfilmsaremodel‐mesoscopicsystemsforunderstanding physicalphenomenasuchas tunableadhesion,adhesiondebonding,patternformation,wetting dewetting and friction at soft interfaces, etc.

When a thin elastic film (shear modulus of <1 MPa) is brought in contact proximity of an external contactor, it forms miniaturized surface patterns because of a competition between the destabilizing interaction forces between the film and the contactor and the restoring elastic forces present in the bulk of the elastic film. The wavelengths of instabilities scale as 4 h in the case of peeling geometry [12–14] (whereas glass slide rests on a rigidly bonded film in cantilever configuration) and application of normal force at the end of the glass slide leads to well‐defined one‐dimensional (1D) finger‐like patterns at the crack opening and 3 h in the case of adhesive geometry [15–26]. In the adhesive‐contact geometry, the free surface of the film in contact proximity with an external contactor spontaneously roughens to form a two‐dimen‐ sional (2D) labyrinthic pattern when the stiffness of the attractive forces due to van der Waals (VDW) interactions (arising because of the contactor) exceeds the elastic stiffness of the film. This pattern remained unchanged in experiments either when the upper glass slide was silanized with a chemisorbed layer of fluorinated monochlorodecyl silane [15] or when the interactions were switched to electrostatic ones [21, 24, 25]. In recent past, much effort has been devoted to understand how further miniaturization of length scales and hence enhanced functional properties can be made technologically feasible in these already short wavelength exhibiting adhesive soft elastic thin films in contact geometry. From a fabrication point of view, the smallest patterns realized till date at such soft interfaces are ~ 0.1 h and are obtained with elastic bilayers when the surface energies of the films are very low and the top film is considered to be much thinner and of very small shear modulus compared to the bottom film [27, 28]. But pattern formation in elastic bilayers through squeezing modes of instabilities has the disad‐ vantage of delaminating at the film-film interface. To alleviate this problem, one has found that films cast on patterned substrates can give rise to surface patterns, which are at length scales that are about an order of magnitude less than that found in films cast on flat substrates [29–31].

Apart from direct experiments, these systems are well analyzed both theoretically and numerically. In both cases, atomistic scale models can be considered but are computationally expensive and are restricted by the domain geometries that can be handled by them. For elastic films, linear stability analysis is generally carried out to understand the critical length scales at the incipience of instability by perturbing Navier's equations with bifurcation modes and obtaining the nontrivial solutions [16, 17]. In nonlinear analysis [18–20], on the other hand, it is possible to find the surface morphology by trying to find the equilibrium displacements that minimize the total energy, which is composed of the interaction, and the elastic energies. Finding the optimized Fourier coefficients, which define the top surface, with the help of conjugate-gradient method, does this. However, for elastic films cast on generally patterned substrates the obtained equations show highly nonlinear behavior, which cannot be solved by the simplistic analytical methods or even by the energy-minimizing schemes developed for the base case. On top of these, polymeric films particularly demonstrate a rich variety of nonlinear behavior ranging from viscoelasticity to plasticity that cannot be adequately captured by the aforementioned methods. Finite element analysis (FEA) methods can be used to alleviate these problems. FEA methods have been used extensively in recent years to study crack opening, crack propagation [32–34], delamination [35], film instability over soft substrates, and bifurcation analysis of elastic contact instability [36–39].

In this chapter we show how to implement cohesive-zone modeling (CZM) for simple shear and normal modes of crack opening in the finite element software and how to simulate contact instability and adhesion debonding in soft thin films cast on complex geometries.

## **2. Cohesive-zone modeling**

top‐down approach and the other a bottom‐up approach. The top‐down approach consists of the various lithographic techniques such as electron beam lithography, dip‐pen lithography, embossing lithography, laser lithography [1], laser photoablation [2], ion implantation [3], micromachining [4], etc. The bottom‐up approach, on the other hand, consists of self‐assembly of molecules to produce bigger building blocks to a growing structure as in chemical vapor deposition [5]. While the self‐assembly procedure involves numerous steps, lithographic methods such as dip‐pen lithography, which are used for creating patterns on hard silicon surfaces, are time‐consuming, involve great cost and may not at all be suitable to create patterns on soft surfaces such as polymers. Self‐organization methods, on the other hand, are for bulk production of desired structures involving lower cost and therefore are more potentially applicablefabricationtechniquesonsoft surfaces.Inself‐organizationtechnique,materials such as polymeric films, which have highly tunable mechanical properties, can be used to obtain a desired mesoscopic structure due to the influence of internal or external forces. The mechanical and interaction properties of these films can be controlled by the amount of cross‐linker added during the synthesis step itself. The instabilities during self‐organization can be harnessed to engineer desired periodic arrangements of sub‐micron features in technological applications such as in flexible electronics [6], optoelectronics [7], functional coatings [8], sensors [9], microfluidic devices, microelecromechanical systems (MEMS) and nanoelectromechanical systems (NEMS) [10], pressure‐sensitive adhesives [11], high‐efficiency light‐emitting diodes, thin‐filmtransistors,etc.Thus,thesethinfilmsaremodel‐mesoscopicsystemsforunderstanding physicalphenomenasuchas tunableadhesion,adhesiondebonding,patternformation,wetting

When a thin elastic film (shear modulus of <1 MPa) is brought in contact proximity of an external contactor, it forms miniaturized surface patterns because of a competition between the destabilizing interaction forces between the film and the contactor and the restoring elastic forces present in the bulk of the elastic film. The wavelengths of instabilities scale as 4 h in the case of peeling geometry [12–14] (whereas glass slide rests on a rigidly bonded film in cantilever configuration) and application of normal force at the end of the glass slide leads to well‐defined one‐dimensional (1D) finger‐like patterns at the crack opening and 3 h in the case of adhesive geometry [15–26]. In the adhesive‐contact geometry, the free surface of the film in contact proximity with an external contactor spontaneously roughens to form a two‐dimen‐ sional (2D) labyrinthic pattern when the stiffness of the attractive forces due to van der Waals (VDW) interactions (arising because of the contactor) exceeds the elastic stiffness of the film. This pattern remained unchanged in experiments either when the upper glass slide was silanized with a chemisorbed layer of fluorinated monochlorodecyl silane [15] or when the interactions were switched to electrostatic ones [21, 24, 25]. In recent past, much effort has been devoted to understand how further miniaturization of length scales and hence enhanced functional properties can be made technologically feasible in these already short wavelength exhibiting adhesive soft elastic thin films in contact geometry. From a fabrication point of view, the smallest patterns realized till date at such soft interfaces are ~ 0.1 h and are obtained with elastic bilayers when the surface energies of the films are very low and the top film is considered to be much thinner and of very small shear modulus compared to the bottom film [27, 28]. But pattern formation in elastic bilayers through squeezing modes of instabilities has the disad‐

dewetting and friction at soft interfaces, etc.

164 Perusal of the Finite Element Method

Most of the failures in engineering structures are the result of cracks that are preexisting (during production) or the cracks that are formed in service. Continuum mechanics calculations fail to handle such cracks and these are handled by fracture mechanics, which primarily studies resistance of a material to fracture. When a crack is formed, stress is usually not transmitted between the two "virtual" surfaces that form the following crack formation. We can consider a cohesive zone, of zero thickness, to bridge these two virtual surfaces and it can be suggested that fracture happens when material strength in this zone decays. The basic premise for cohesive-zone model is that infinite stresses (stress singularity) at crack tip are unrealistic. Dugdale [40] and Barenblatt [41] developed the first models to address this shortcoming and model brittle fracture. Their models divide the crack tip into stress-free region and cohesively stressed regions. More than a decade later, Hillerborg [42] was the first one to use finite element analysis for cohesive-zone modeling of brittle fracture. Tvergaard's [43–44] modification to Dugdale and Barenblatt's model ensured that continuum elements remain intact in the cohesive model. It is their first decohesion model where traction is dependent on both normal and tangential separation, whereas in Barenblatt's model, it is dependent on the distance from the crack tip. They do this by defining the cohesive elements between the continuum elements such that damage initiation is characterized by opening of the cohesive elements (and governed by some material law) and at failure these elements lose their stiffness, resulting in the continuum elements disjointing from each other. This model is also called the cohesive‐zone model or damage zone model (DZM). The initiation and evolution of the damage are defined according to material‐softening law/cohesive law/bridging law/traction‐ separation law (TSL). TSL can be defined based on the properties of the cohesive zone. Traction‐ separation laws are not universal, and these are not physical laws to characterize fracture at bonded interfaces. Rather, they are highly material dependent. We can define the TSLs in many ways but all separation laws have one thing in common: they all suggest that the traction initially increases with separation and beyond a critical separation distance it diminishes to zero (material softening) at which point these cohesive surfaces detach from each other and a real crack tip is defined.

For implementing TSL in CZM, two constitutive relations are considered: one for the contin‐ uum or the "bulk behavior" of the blocks and another for the gap between the surfaces that is defined by the traction‐displacement relationships. If the cohesive zone is a thin adhesive layer bonding the blocks, we can define the TSL as the traction at this interface with respect to the separation at the interface. In the modeling by commercial finite element software ABAQUS, the cohesive zone is represented by a single layer of elements that connect the two blocks/ surfaces and is constrained by one of the blocks as can be seen from **Figure 1**. Each element on the surface of one block is "tied" with its corresponding neighboring element in the other block through a series of "inelastic springs" which govern the interaction between these two surfaces. This tie constraint prevents the propagation of singular modes in deformation that can arise due to lack of membrane stiffness. The debonding/crack‐opening process is akin to the failure of these springs.

**Figure 1.** Meshing with tie constraints for cohesive elements in ABAQUS.

#### **2.1. Cohesive zone implementation in ABAQUS with an example problem**

To show the implementation of cohesive zone model with traction‐separation law in ABAQUS, we have provided an example problem [43] below. In this, we study the propagation of a crack in a material block. A debonding model is taken from Tvergaard's work [43] and is described in brief below. This model does not consist of a debonding potential, but is designed to address the effect of fiber debonding from a metal‐matrix composite on the tensile properties (as this requires one to consider both normal and tangential separations in the debonding). For this, first, a function *F*(*λ*) was chosen such that

$$F\left(\mathcal{X}\right) = \frac{27}{4}\sigma\_{\text{max}}\left(1 - 2\mathcal{X} + \mathcal{X}^2\right) \text{ and} \\ 0 \le \mathcal{X} \le 1 \tag{1}$$

The nondimensional parameter *λ* is defined as

both normal and tangential separation, whereas in Barenblatt's model, it is dependent on the distance from the crack tip. They do this by defining the cohesive elements between the continuum elements such that damage initiation is characterized by opening of the cohesive elements (and governed by some material law) and at failure these elements lose their stiffness, resulting in the continuum elements disjointing from each other. This model is also called the cohesive‐zone model or damage zone model (DZM). The initiation and evolution of the damage are defined according to material‐softening law/cohesive law/bridging law/traction‐ separation law (TSL). TSL can be defined based on the properties of the cohesive zone. Traction‐ separation laws are not universal, and these are not physical laws to characterize fracture at bonded interfaces. Rather, they are highly material dependent. We can define the TSLs in many ways but all separation laws have one thing in common: they all suggest that the traction initially increases with separation and beyond a critical separation distance it diminishes to zero (material softening) at which point these cohesive surfaces detach from each other and a

For implementing TSL in CZM, two constitutive relations are considered: one for the contin‐ uum or the "bulk behavior" of the blocks and another for the gap between the surfaces that is defined by the traction‐displacement relationships. If the cohesive zone is a thin adhesive layer bonding the blocks, we can define the TSL as the traction at this interface with respect to the separation at the interface. In the modeling by commercial finite element software ABAQUS, the cohesive zone is represented by a single layer of elements that connect the two blocks/ surfaces and is constrained by one of the blocks as can be seen from **Figure 1**. Each element on the surface of one block is "tied" with its corresponding neighboring element in the other block through a series of "inelastic springs" which govern the interaction between these two surfaces. This tie constraint prevents the propagation of singular modes in deformation that can arise due to lack of membrane stiffness. The debonding/crack‐opening process is akin to the failure

real crack tip is defined.

166 Perusal of the Finite Element Method

of these springs.

**Figure 1.** Meshing with tie constraints for cohesive elements in ABAQUS.

**2.1. Cohesive zone implementation in ABAQUS with an example problem**

To show the implementation of cohesive zone model with traction‐separation law in ABAQUS, we have provided an example problem [43] below. In this, we study the propagation of a crack

$$\mathcal{A} = \sqrt{\left(\frac{\mu\_n}{\mathcal{S}\_n}\right)^2 + \left(\frac{\mu\_t}{\mathcal{S}\_t}\right)^2} \tag{2}$$

For a monotonically increasing value of the tractions are = and = . Here, the normal and tangential tractions are and respectively, normal and tangential displace‐ ments are and respectively, δ and δ are the maximum separation distances in normal and traction modes, respectively. Here, we have considered that the two virtual surfaces that form after a crack growth are of unit dimension. For solving this, we have developed our own code using MATLAB and ABAQUS. MATLAB is excellent to call external commands to run the PYTHON scripts for generating the input files for ABAQUS.

ABAQUS has vast libraries of constitutive material models but this library is neither a complete one nor a completely flexible one. However, fortunately, user‐defined constitutive material models can be implemented in ABAQUS through either UMAT or VUMAT. UMAT is an implicit time integration method (used with ABAQUS/Standard) and it requires the material stiffness matrix as input for forming the consistent Jacobian (*C*), which is also called algorithmic tangential stiffness or ATS, of the nonlinear equations. For linearly elastic materials, it is defined as

$$C = \frac{\partial \Delta \sigma}{\partial \Delta \varepsilon} \tag{3}$$

where is the increment in Cauchy stress and is the increment in strain.

VUMAT, on the other hand, uses explicit time integration and hence does not require Jacobian matrix. It is implemented in ABAQUS/Explicit and is easier to write. In the case of cohesive zone modeling, we can include element failure and progressive damage failure in the VUMAT subroutine. The VUMAT subroutine consists of the strain decomposition of the total strain into elastic (computed from a linear elastic or hyperelastic constitutive model) and plastic compo‐ nents and flow rule to define the material's plastic behavior. A sample of the input file (filename.inp) and VUMAT subroutine (vumat.f) are annotated and provided at the end of this section. These are run in ABAQUS through the following command:

abaqus job=*filename.inp* user=*vumat.f*

There are a minimum of five parameters that have to be prescribed for the CZM. They are as follows:

Cohesive strength (related to yield stress) and cohesive energy signify the maximum resistance to fracture and dissipation following material separation, respectively. For the cohesive zone modeling with traction‐separation law, we need to specify the initial elastic stiffness (*K*).

**Figure 2.** Cohesive zone model (CZM) in ABAQUS. Loads and boundary conditions for the case of purely normal trac‐ tion (T*t* = 0).

**Figure 3.** Normal traction (*Tn*) versus separation (*un*) obtained after CZM in ABAQUS.

We cannot use the exact value of the material's Young's modulus (*E*) for this as this stiffness here refers to a "penalty stiffness," that is, we are free to choose a value that is necessary to prevent interpenetration of the surfaces and for the convergence of our model. One rule of thumb is to take a value of *K* that is equal to *E/h* (*h* is the cohesive layer thickness). Hence, *K* can be taken to be a very large value as the cohesive elements are considered to be a layer of zero thickness. Also, for simplification, identical values of *K* are taken for the crack‐opening mode (mode I) and crack‐shearing modes (mode II). The chosen damage initiation and damage evolution criteria affect how these two modes of crack opening interact with each other.

The variables we output (stress, strain, etc.) depend on the kind of problems being modeled by the cohesive elements. This is indicated in the .inp file itself in the definition of the section properties by specifying the response type required. The stress components are as follows:

11—Direct membrane stress

nents and flow rule to define the material's plastic behavior. A sample of the input file (filename.inp) and VUMAT subroutine (vumat.f) are annotated and provided at the end of

There are a minimum of five parameters that have to be prescribed for the CZM. They are as

Cohesive strength (related to yield stress) and cohesive energy signify the maximum resistance to fracture and dissipation following material separation, respectively. For the cohesive zone modeling with traction‐separation law, we need to specify the initial elastic stiffness (*K*).

**Figure 2.** Cohesive zone model (CZM) in ABAQUS. Loads and boundary conditions for the case of purely normal trac‐

We cannot use the exact value of the material's Young's modulus (*E*) for this as this stiffness here refers to a "penalty stiffness," that is, we are free to choose a value that is necessary to prevent interpenetration of the surfaces and for the convergence of our model. One rule of

, the peak values of the traction‐separation curve).

, usually the separation value corresponding to the cohesive

this section. These are run in ABAQUS through the following command:

**2.** Cohesive energy (area under the traction‐separation curve).

**Figure 3.** Normal traction (*Tn*) versus separation (*un*) obtained after CZM in ABAQUS.

abaqus job=*filename.inp* user=*vumat.f*

168 Perusal of the Finite Element Method

**1.** Cohesive strength (*Tn* and *Tt*

**3.** Characteristic length (*un* and *ut*

strength).

follows:

tion (T*t* = 0).

**Figure 4.** Cohesive zone model (CZM) in ABAQUS. Loads and boundary conditions shown for the case of purely shear traction (T*n*= 0).

**Figure 5.** Tangential traction (*Tt* ) versus shear separation (*ut* ) obtained from CZM in ABAQUS.

During a purely normal crack opening ( = 0), according to our TSL, the normal traction initially increases with the normal separation , reaches a maximum value max and reduces to 0 when = that is, when the final separation occurs (see **Figures 2** and **3**). Similarly, during a purely shear crack opening ( = 0), the shear traction initially increases with tangential separation , reaches a maximum value max and then at = the final separation occurs (see **Figures 4** and **5**). The negative values in **Figure 5** indicate compressive stresses.

### **2.2. Writing ABAQUS .inp file and VUMAT subroutine**

## *2.2.1. Snippet of input file*

Snippet of input file (filename.inp) containing the geometry, meshing, material parameters, and output variables

\*Heading

```
** Job name: cohesive‐zone Model name: Model‐1
```

```
** Generated by: Abaqus/CAE 6.9‐1
```

```
*Preprint, echo=NO, model=NO, history=NO, contact=NO
```
\*\* PARTS

```
*Part, name=Part‐1
```
\*End Part

\*\*

```
** ASSEMBLY
```
\*\*

```
*Assembly, name=Assembly
```
\*\*

```
*Instance, name=Part‐1‐1, part=Part‐1 ***(two unit blocks/surfaces, each with 4 nodes)
```
\*Node

⋮

```
*Element, type=CPE4R***(4 noded continuum plane-strain elements for the blocks)
```
⋮

```
*Element, type=COH2D4 ***(4 noded cohesive elements for the cohesive zone)
```
⋮

\*Nset, nset=\_PickedSet*i* , internal, generate

⋮

```
*Elset, elset=_PickedSeti
                       , internal
⋮
** Section: Section‐1
*Solid Section, elset=_PickedSeti
                               , material=steel
,
** Section: Section‐2
*Cohesive Section, elset=_PickedSeti
                                     , material=coh, response=TRACTION SEPARATION
***(specifying response as TRACTION SEPARATION to implement TSL)
*End Instance
**
⋮
*Surface, type=ELEMENT, name=_PickedSurf21, internal
__PickedSurf21_S3, S3
*End Assembly
*Amplitude, name=Amp‐1 ***(time dependent loading conditions)
0., 0., 0.5, 0.5, 1., 1.
t=0 amp_at_t=0 t=0.5 amp_at_t=0.5 t=1 amp_at_t=1
**
** MATERIALS
**
*Material, name=coh ***(cohesive zone parameters)
*User material, Constants = 5 ***(user defined subroutine with 5 parameters)
200.0,0.25,0.25,1.0,0.5 * * *(, ,
                                     ,  .   these are read through props( ) matrix in vumat.f)
*Density
5.4e‐7
*Depvar
7, ***(number of solution-dependent variables SDVs)
*Initial conditions, type=solution
Part‐1‐1._PickedSet4,0.0,0.0,0,0.0,0.0,0.0,0.0 ***(initial values of the 7 SDVs)
*Material, name=material‐1 ***(material properties of the continuum surfaces)
```
to 0 when = that is, when the final separation occurs (see **Figures 2** and **3**). Similarly,

separation occurs (see **Figures 4** and **5**). The negative values in **Figure 5** indicate compressive

Snippet of input file (filename.inp) containing the geometry, meshing, material parameters,

\*Instance, name=Part‐1‐1, part=Part‐1 *\*\*\*(two unit blocks/surfaces, each with 4 nodes)*

\*Element, type=CPE4R*\*\*\*(4 noded continuum plane-strain elements for the blocks)*

\*Element, type=COH2D4 *\*\*\*(4 noded cohesive elements for the cohesive zone)*

, internal, generate

, reaches a maximum value max and then at =

initially increases with

the final

during a purely shear crack opening ( = 0), the shear traction

**2.2. Writing ABAQUS .inp file and VUMAT subroutine**

\*\* Job name: cohesive‐zone Model name: Model‐1

\*Preprint, echo=NO, model=NO, history=NO, contact=NO

\*\* Generated by: Abaqus/CAE 6.9‐1

tangential separation

170 Perusal of the Finite Element Method

*2.2.1. Snippet of input file*

and output variables

stresses.

\*Heading

\*\* PARTS

\*End Part

\*\* ASSEMBLY

\*Assembly, name=Assembly

\*Nset, nset=\_PickedSet*i*

\*\*

\*\*

\*\*

⋮

⋮

⋮

⋮

\*Node

\*Part, name=Part‐1

```
*Density
7.8e‐06,
*Elastic
2.5e+8, 0.25***(Young's modulus E and Poisson's ratio of Material-1)
** ‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐*
** STEP: Step‐1
**
*Step, name=Step‐1
*Dynamic, Explicit
, 1.
*Bulk Viscosity
0.06, 1.2
** Mass Scaling: Semi‐Automatic
** Whole Model
*Fixed Mass Scaling, factor=50.
**
** BOUNDARY CONDITIONS
**
** Name: BC‐2 Type: Displacement/Rotation
*Boundary
_PickedSet17, 2, 2
** Name: BC‐3 Type: Displacement/Rotation
*Boundary
_PickedSet18, 1, 1
** Name: BC‐4 Type: Displacement/Rotation
*Boundary, amplitude=Amp‐1
_PickedSet19, 2, 2, 1.
**
** LOADS
**
```
\*\* Name: Load‐1 Type: Pressure

\*Dsload, amplitude=Amp‐3

\_PickedSurf21, P,50. *\*\*\*(Surface-2 of and permitted to have maximum displacement 0.5 in xdirection i.e. =0.5)*

\*\*

\*Density 7.8e‐06, \*Elastic

172 Perusal of the Finite Element Method

\*\* STEP: Step‐1

\*Bulk Viscosity

\*\* Whole Model

\*\* Mass Scaling: Semi‐Automatic

\*Fixed Mass Scaling, factor=50.

\*\* BOUNDARY CONDITIONS

\*\* Name: BC‐2 Type: Displacement/Rotation

\*\* Name: BC‐3 Type: Displacement/Rotation

\*\* Name: BC‐4 Type: Displacement/Rotation

\*Boundary, amplitude=Amp‐1

\*Step, name=Step‐1 \*Dynamic, Explicit

\*\*

, 1.

\*\*

\*\*

\*Boundary

\*Boundary

\*\*

\*\*

\*\* LOADS

\_PickedSet17, 2, 2

\_PickedSet18, 1, 1

\_PickedSet19, 2, 2, 1.

0.06, 1.2

2.5e+8, 0.25*\*\*\*(Young's modulus E and Poisson's ratio of Material-1)*

\*\* ‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐\*

\*\* OUTPUT REQUESTS

\*\*

\*Restart, write, number interval=1, time marks=NO

\*\*

\*\* FIELD OUTPUT: F‐Output‐2

\*\*

\*Output, field,op=new,number interval=50

\*Node Output

RF, U *\*\*\*(output variables at the nodes: reaction forces and displacements)*

\*Element Output, directions=YES

S, SDV *\*\*\*(output variables over the elements: stresses and SDVs which are defined in the .f file of VUMAT)*

\*\*

\*\* FIELD OUTPUT: F‐Output‐1

\*\*

\*Output, field, variable=PRESELECT

\*\*

\*\* HISTORY OUTPUT: H‐Output‐1

\*\*

\*Output, history, variable=PRESELECT

\*End Step

*2.2.2. Snippet of vumat.f subroutine*

The headers of this subroutine are common to every VUMAT and can be obtained from ABAQUS documentation. We have listed below only portions of the user‐defined calculations for the properties.

*! Reading material properties specified in abaqus input file through props() array* = props(1) *! max nominal stress at which damage is initiated* = props(2) *! max normal displacement at which the material fails and Tn=0* = props(3) *! max tangential displacement at which the material fails and Tt=0*

 $a = \text{props(4)}$   $\left[ \text{pparameter in } T\_t = a \frac{u\_t}{\delta\_t} F(\lambda) \right]$

mu = props(5)

do ii = 1,nblock

*! Reading the old state variables or SDVs*.

= stateOld(ii,1) *! normal displacement*

= stateOld(ii,2) *!tangential displacement*

counter = stateOld(ii,3)

 = stateOld(ii,4)*!tangential traction*

= stateOld(ii,5)*!normal traction*

*! Reading the strain Increments*.

*! In VUMAT shear components are stored as tensor components and not as engineering components*.

*! so* Δ  *will be twice of tensorial shear strain*.

```
Δ = strainInc(ii,1)
```
Δ = 2.0\*strainInc(ii,2)

*!Depending on pure normal traction or pure shear traction or mixed mode of debonding, calculate the following:*

$$\begin{array}{rcl} \mathcal{A} & = & \cdots & \cdots \\\\ T\_{\kappa, \text{new}} & = & \frac{\mathcal{U}\_{\kappa}}{\delta\_{\kappa}} F(\mathcal{A}) \end{array}$$

$$T\_{\iota, \text{new}} \quad = \quad \alpha \frac{\mathcal{U}\_{\iota}}{\delta\_{\iota}} F(\mathcal{A})$$

Simulating Contact Instability in Soft Thin Films through Finite Element Techniques http://dx.doi.org/10.5772/65357 175

$$\begin{split} \boldsymbol{\sigma}\_{11} &= \left(\bigvee\_{\boldsymbol{\beta}\boldsymbol{\lambda}\_{n}} \right) \Bigg( \boldsymbol{F}\left(\boldsymbol{\lambda}\right) + \boldsymbol{u}\_{n} \frac{\partial \boldsymbol{F}\left(\boldsymbol{\lambda}\right)}{\partial \boldsymbol{u}\_{n}} \Bigg); \frac{\partial \boldsymbol{F}\left(\boldsymbol{\lambda}\right)}{\partial \boldsymbol{u}\_{n}} = \frac{\partial \boldsymbol{F}\left(\boldsymbol{\lambda}\right)}{\partial \boldsymbol{\lambda}} \frac{\partial \boldsymbol{\lambda}}{\partial \boldsymbol{u}\_{n}} \\ \boldsymbol{\sigma}\_{12} &= \left(\bigvee\_{\boldsymbol{\beta}\boldsymbol{\lambda}\_{n}} \right) \frac{\partial \boldsymbol{F}\left(\boldsymbol{\lambda}\right)}{\partial \boldsymbol{u}\_{i}}; \frac{\partial \boldsymbol{F}\left(\boldsymbol{\lambda}\right)}{\partial \boldsymbol{u}\_{i}} = \frac{\partial \boldsymbol{F}\left(\boldsymbol{\lambda}\right)}{\partial \boldsymbol{\lambda}} \frac{\partial \boldsymbol{\lambda}}{\partial \boldsymbol{u}\_{i}} \\ \boldsymbol{\sigma}\_{13} &= \boldsymbol{\alpha} \Big(\boldsymbol{u}\_{\boldsymbol{\beta}} \Big/ \boldsymbol{F}\left(\boldsymbol{\lambda}\right) \Bigg) \\ \boldsymbol{\sigma}\_{14} &= \left(\boldsymbol{\alpha} \Big{\boldsymbol{\beta}}\_{1}\right) \Big(\boldsymbol{F}\left(\boldsymbol{\lambda}\right) + \boldsymbol{u}\_{i} \frac{\partial \boldsymbol{F}\left(\boldsymbol{\lambda}\right)}{\partial \boldsymbol{u}\_{i}} \Bigg) \end{split}$$

*!Increment in stress*

*! Reading material properties specified in abaqus input file through props() array*

= props(2) *! max normal displacement at which the material fails and Tn=0*

= props(3) *! max tangential displacement at which the material fails and Tt=0*

*! In VUMAT shear components are stored as tensor components and not as engineering*

*!Depending on pure normal traction or pure shear traction or mixed mode of debonding,*

,new

*u T F*

=

=

*u T F*

=

*n*

l

*t*

,new

( )

l

*n*

LLL

*n t*

d

a

*t*

d

( )

 l

= props(1) *! max nominal stress at which damage is initiated*

*components*.

Δ = strainInc(ii,1)

= 2.0\*strainInc(ii,2)

*calculate the following:*

*! so* Δ

Δ

mu = props(5)

do ii = 1,nblock

counter = stateOld(ii,3)

= props(4) *!parameter in*  =

174 Perusal of the Finite Element Method

*! Reading the old state variables or SDVs*.

= stateOld(ii,1) *! normal displacement*

= stateOld(ii,4)*!tangential traction*

 *will be twice of tensorial shear strain*.

= stateOld(ii,5)*!normal traction*

*! Reading the strain Increments*.

= stateOld(ii,2) *!tangential displacement*

$$\begin{aligned} \Delta T\_{\boldsymbol{u}} &= (\sigma\_{11} \, ^\ast \Delta \boldsymbol{u}\_{\boldsymbol{n}}) + (\sigma\_{12} \, ^\ast \Delta \boldsymbol{u}\_{\boldsymbol{\iota}}), \\ \Delta T\_{\boldsymbol{\iota}} &= (\sigma\_{21} \, ^\ast \Delta \boldsymbol{u}\_{\boldsymbol{n}}) + (\sigma\_{22} \, ^\ast \Delta \boldsymbol{u}\_{\boldsymbol{\iota}}). \end{aligned}$$

*! stress update*

```
stressNew(ii,1) = stressOld(ii,1)+
stressNew(ii,2) = stressOld(ii,2)+
! state variable update
stateNew(ii,1) = stateOld(ii,1) +
stateNew(ii,2) = stateOld(ii,2) +
stateNew(ii,3) = stateOld(ii,3) + 1
stateNew(ii,4) =
stateNew(ii,5) =
enddo
return
end
```
## **3. Elastic contact instability through FEM simulations**

**Figure 6A** and **B** shows the schematic setup of a soft incompressible initially stress‐free elastic film of shear modulus <1 MPa, mean thickness *h*, and lateral length *L*. It is cast over a rough substrate. To demonstrate the effect of roughness here, we have chosen a substrate with a uniform sawtooth profile characterized by its amplitude ℎ and wavelength = / where is the number of substrate troughs or crests. In the asymptotic case when ℎ 0 or / ∞, we get the classical case of a smooth substrate, though in real systems it is almost impossible to come across a substrate that is atomistically smooth and almost all surfaces have inherent roughness. Again in real situations mechanical/thermal noises and vibrations almost always disturb the film surface at the air interface. To incorporate them for numerical simulations, the film surface is perturbed by sinusoidal disturbances having different frequencies and random amplitudes of the order of 10−3 nm. The mean distance from the free surface of the film and the contactor is termed the gap distance and is represented by the symbol *d* (refer to **Figure 6A**). In simulations, the contactor, as in experiments, is brought near the film surface from a finite distance at an approach velocity ( 1 /) maintained below the critical distance where the instabilities first initiate (adhesion) for around ~15 s. This is necessary to allow the film to have sufficient time to arrange itself under the interplay of imposed force and elasticity of the film. Once the patterns get mature, the debonding phase is started. In the debonding phase, the contactor is pulled off at different velocities (debonding: ). The contactor positions at different times can be visualized from **Figure 6B**. The thickness of the film considered for the simulations is 500 nm and length *L* of the film considered is 8ℎ. The *β* is varied from 0.1 to 0.7 and *n* is varied from 1 to 12 (for *L = 8h*) in the simulations. Though in reality the film is three‐ dimensional (3D) as shown in **Figure 6A**, for ease of computation the simulations are carried out considering a plane‐strain model (as because of normal traction from the contactor the displacements can be considered to be confined in the plane of the paper) and a 2D geometry as shown in **Figure 6B** has been mostly considered. To understand the morphology in greater detail, 3D simulations are also carried out in certain cases with the 3D geometry.

**Figure 6.** Schematic diagram of an elastic thin film case on a (A) 3D sawtooth and (B) 2D sawtooth‐patterned substrate and in contact proximity to an external contactor. The position of the contactor at different time is also shown.

If one considers an infinitesimal control volume inside the elastic film, a force balance will give rise to the following governing thin‐film (Navier's) equation:

$$
\vec{\nabla} \cdot \overleftarrow{\sigma} + \vec{f} = 0 \tag{4}
$$

where, = + ∇ + ∇ is the stress developed inside the film, is the displacement vector, *p* is the pressure across the film, and **I** is the identity matrix. The second term, , represents the body force (per unit volume) that the film experiences because of the contactor. The equation as can be seen is a partial differential equation in the unknown and does not yield an analytical solution in the case when the substrate is patterned or when the film material displays geometric or material nonlinearity. A numerical solution can, however, be obtained when the differential equation is discretized by considering a trial/shape function for the displacement and converted into a set of algebraic equations. In the Galerkin finite element formulation, the residual obtained because of the trial function can now be minimized in a weighted fashion according to

$$\int\_{\Omega} \left(\vec{\nabla} \cdot \vec{\sigma} + \vec{f}\right) \phi d\Omega = 0 \tag{5}$$

where Ω is the discretized domain inside the film and Φ is the shape function. The desired discretization and subsequent simulations to realize the film deformation and morphology due to influence of the external contactor are carried out with the help of the commercial finite element software ABAQUS

#### **3.1. Geometry and mesh creation**

is the number of substrate troughs or crests. In the asymptotic case when ℎ 0 or / ∞, we get the classical case of a smooth substrate, though in real systems it is almost impossible to come across a substrate that is atomistically smooth and almost all surfaces have inherent roughness. Again in real situations mechanical/thermal noises and vibrations almost always disturb the film surface at the air interface. To incorporate them for numerical simulations, the film surface is perturbed by sinusoidal disturbances having different frequencies and random amplitudes of the order of 10−3 nm. The mean distance from the free surface of the film and the contactor is termed the gap distance and is represented by the symbol *d* (refer to **Figure 6A**). In simulations, the contactor, as in experiments, is brought near the film surface from a finite distance at an approach velocity ( 1 /) maintained below the critical distance where the instabilities first initiate (adhesion) for around ~15 s. This is necessary to allow the film to have sufficient time to arrange itself under the interplay of imposed force and elasticity of the film. Once the patterns get mature, the debonding phase is started. In the debonding phase, the contactor is pulled off at different velocities (debonding: ). The contactor positions at different times can be visualized from **Figure 6B**. The thickness of the film considered for the simulations is 500 nm and length *L* of the film considered is 8ℎ. The *β* is varied from 0.1 to 0.7 and *n* is varied from 1 to 12 (for *L = 8h*) in the simulations. Though in reality the film is three‐ dimensional (3D) as shown in **Figure 6A**, for ease of computation the simulations are carried out considering a plane‐strain model (as because of normal traction from the contactor the displacements can be considered to be confined in the plane of the paper) and a 2D geometry as shown in **Figure 6B** has been mostly considered. To understand the morphology in greater

176 Perusal of the Finite Element Method

detail, 3D simulations are also carried out in certain cases with the 3D geometry.

**Figure 6.** Schematic diagram of an elastic thin film case on a (A) 3D sawtooth and (B) 2D sawtooth‐patterned substrate and in contact proximity to an external contactor. The position of the contactor at different time is also shown.

If one considers an infinitesimal control volume inside the elastic film, a force balance will give

*f* 0 r r (4)

Ñ× + = s

rise to the following governing thin‐film (Navier's) equation:

To mimic the physical geometry and to simulate several cases of substrates (having varying frequency and amplitudes) and to incorporate different perturbation conditions at the filmcontactor interface, the graphical user interface in the software proves to be difficult. Python scripting was used instead to generate the multitude of geometries. For the film's lateral edges and the substrate boundary, straight lines are used, while for achieving the perturbed top surface cubic splines are fitted to get a continuous line out of the random sinusoidal fluctuations at the film-contactor interface. Meshing forms an integral part of the preprocessing steps, which helps in discretizing the geometry. The accuracy that can be obtained from any finite element analysis is directly related to the finite element mesh that is used. Generally, for the large deformations in 2D, eight node elements are used but those are computationally intensive. In 2D as discussed earlier for the present scenario one can consider a plain-strain condition, therefore, continuum plain-strain 4 node-reduced integration elements (CPE4R) are used in 2D. For 3D geometry, continuum three-dimensional eight-node-reduced integration (C3D8R) elements are used. To circumvent the problem of artificial stiffness due to shear locking that can arise due to four node elements, a large number of elements are used while maintaining the aspect ratio of the elements close to 1, which also helps in avoiding problems associated with skewness. Further, for the large deformation problem, reduced integration elements are shown to be successfully implemented without impacting the results and reduce the computational load. On the top of the film, the force is much higher and decreases drastically with depth in the film due to the highly nonlinear dependence of force on distance. Thus, it becomes necessary to have high density of mesh toward the top of the film. The geometry of the film in case of a smooth substrate is seeded with a bias ratio to have smaller elements toward the top as shown in **Figure 7A** and **B** (for 2D and 3D geometries, respectively). For sawtooth profile, the geometry is divided into two blocks. The top block has similar seeding as in case of smooth substrate and the bottom is meshed as shown in **Figure 7C** and **D** (for 2D and 3D geometries, respectively). Since the force applied depends on the amplitude of substrate roughness as well, the meshed geometry is tested with the grid independence test for every case and the finalized meshes are then chosen.

**Figure 7.** The meshing done inside the film geometry for (A)&(B) a planar substrate, (C)&(D) for a sawtooth substrate, respectively. Figures on the left and right represent meshing in 2D and 3D respectively.

## **3.2. Boundary conditions**

To imitate the rigid boundary conditions of the base in the simulations, the bottom edge of the film is pinned to the substrate and movement toward either direction is restricted with the help of an Encastre boundary condition. To incorporate periodic boundary conditions on the lateral sides of the film, x‐symmetric boundary conditions are applied.

#### **3.3. Force calculations**

As discussed earlier, so far the failure at the interfaces was captured via a cohesive zone model. In the model, the contactor (a rigid body) and the elastic film were considered to be in intimate contact separated only by a very thin layer of cohesive zone. In the cohesive zone, a traction‐ separation law (inbuilt in ABAQUS) is defined in such a way that the traction initially increases with separation and beyond a critical separation distance it diminishes to zero (material softening) at which point these cohesive surfaces detach from each other. The energy release rate at debonding is given by the area under the curve. The problem with the method is except for the energy all the other parameters are not well characterized and the traction‐separation law is not based on any physical law. Also, while the model is good for peeling/debonding experiments and almost always considers the initial stage as one where the rigid contactor and the film is in intimate contact, it fails to capture the essential morphologies at critical separation distance.

geometry of the film in case of a smooth substrate is seeded with a bias ratio to have smaller elements toward the top as shown in **Figure 7A** and **B** (for 2D and 3D geometries, respectively). For sawtooth profile, the geometry is divided into two blocks. The top block has similar seeding as in case of smooth substrate and the bottom is meshed as shown in **Figure 7C** and **D** (for 2D and 3D geometries, respectively). Since the force applied depends on the amplitude of substrate roughness as well, the meshed geometry is tested with the grid independence test

**Figure 7.** The meshing done inside the film geometry for (A)&(B) a planar substrate, (C)&(D) for a sawtooth substrate,

To imitate the rigid boundary conditions of the base in the simulations, the bottom edge of the film is pinned to the substrate and movement toward either direction is restricted with the help of an Encastre boundary condition. To incorporate periodic boundary conditions on the

As discussed earlier, so far the failure at the interfaces was captured via a cohesive zone model. In the model, the contactor (a rigid body) and the elastic film were considered to be in intimate contact separated only by a very thin layer of cohesive zone. In the cohesive zone, a traction‐ separation law (inbuilt in ABAQUS) is defined in such a way that the traction initially increases with separation and beyond a critical separation distance it diminishes to zero (material softening) at which point these cohesive surfaces detach from each other. The energy release rate at debonding is given by the area under the curve. The problem with the method is except

respectively. Figures on the left and right represent meshing in 2D and 3D respectively.

lateral sides of the film, x‐symmetric boundary conditions are applied.

**3.2. Boundary conditions**

178 Perusal of the Finite Element Method

**3.3. Force calculations**

for every case and the finalized meshes are then chosen.

Thus, to capture both the adhesion and debonding of an elastic film due to the interaction between the film and the contactor, a body‐force approach has been considered, where each node in the film experiences a varying magnitude body force depending on its relative position with the contactor. The most generic interaction between two different materials when it comes in contact proximity is given by van der Waals type of potential. The attractive potential of VDW prevalent between two atoms or molecules separated at a distance *r* can be described as [45]

$$\Psi(r) = -M \,/ \, r^6 \tag{6}$$

where *M* is the coefficient of attraction. However, a molecule in the film experiences attraction from all the molecules in the contactor. Thus, if one considers a ring‐like element in the contactor as shown in **Figure 8A**, the total number of atoms/molecules present in the ring will be 2 12 where denotes the number density of atoms/molecules of the contactor. Hence, the total effective potential that a molecule in the film will experience due to a contactor that is semi‐infinitely thick is laterally unbounded and is situated at a surface‐to‐surface separation distance of *D* which is given by

$$\text{cov}(D) = -2\pi \rho\_c M \int\_{\mathbf{x}\_2 - D}^{\mathbf{x}\_2 - \epsilon} d\mathbf{x}\_2 \int\_{\mathbf{x}\_1 - 0}^{\mathbf{x}\_1 - \epsilon} \frac{\mathbf{x}\_1 d\mathbf{x}\_1}{\left(\sqrt{\mathbf{x}\_1^2 + \mathbf{x}\_2^2}\right)^6} = -\frac{\pi M \rho\_c}{6D^3} \tag{7}$$

If the single molecule of the film is now replaced by a solid body of unit volume having number density of atoms/molecules , the body force that the elastic film will experience is given by

$$f\left(D\right) = -\rho\_f \frac{dw}{dD} = -\frac{\pi M \,\rho\_f \rho\_c}{2D^4} = -\frac{A}{2\pi D^4} \tag{8}$$

where *A* is known as the Hamaker constant and has the value of 10−19 J for the present simulations. In the simulations, the above body‐force term is modified and a born‐repulsion term is also added to avoid extremely high forces at total contact such that

$$f\left(D\right) = -\frac{A}{2\pi D^4} + \frac{9B}{D^{10}}\tag{9}$$

where = 6/18 Jm6 (*le* is an equilibrium distance and has a value of 0.138 nm for the present case). The nature of the body force can be seen in **Figure 8B**.

**Figure 8.** (A) Schematic diagram of the calculation of force between a contactor surface and a molecule inside the elas‐ tic film. (B) It shows the van der Waals body force present between the contacting surface and any nodal point in the mesh inside the film geometry.

The contactor in the simulations is virtual in nature, and its presence is felt only through the body‐force term *f*. Since the contactor position keeps on changing with time (refer to **Figure 6B**) so also does the body force. Hence, the gap distance between a particular node and the surface of the contactor keeps on changing with time as per the following:

During adhesion: = initial − × − 2 for 0 ≤ t ≤ 1 where 1 = ( + 1)/ During resting: = initial − + 1 − 2 for 1 ≤ t ≤ 2 where 2( + 1)/ + 15 During debonding phase:

> ( initial ( )) ( 2 2 ) ( ) 2 1 For t t *D d v tt ut c d* x= - + + ´- - £ (10)

In the above set of equations, initial denotes the initial gap distance, 2 the normal position of the particular node at time *t*, the critical distance at which the film begins to experience substantial force due to the contactor and starts deforming. Its value can be roughly estimated analytically from linear stability analysis on sinusoidal substrates according to [30]

$$\frac{hA}{2\pi\mu d\_c^4} = \frac{6.22}{\left(1-\beta\right)}\tag{11}$$

In VDLOAD subroutine of ABAQUS, the body force for a particular node is thus calculated according to

Simulating Contact Instability in Soft Thin Films through Finite Element Techniques http://dx.doi.org/10.5772/65357 181

$$f\left(\xi\right) = -\frac{A}{2\pi\xi^4} + \frac{9B}{\xi^{10}}\tag{12}$$

The morphological evolution of the film due to the presence of the contactor is thus dependant on the dynamic contactor position and thus is a time‐dependent problem. The solution is obtained by marching forward in time using an explicit scheme such that the body force calculated at a particular time *t* is given by

$$\int f\left(t\right) = f\left(t - \Delta t\right) + \frac{\partial \tilde{\xi}}{\partial t} \times \frac{\partial f}{\partial \tilde{\xi}}\bigg|\_{\left(t - \Delta t\right)}\tag{13}$$

#### **3.4. Structure of VDLOAD subroutine**

The headers of this subroutine are common to every VDLOAD and can be obtained from ABAQUS documentation. We have listed below only portions of the user‐defined load subroutine

**\*\*\* Declare new variables used (if any)**

**\*\*\* Declare the values of contactor approach velocity, retracting velocity, van der Waals potential/force parameters, initial separation distance, cutoff distance**

do k = 1, nblock

where =

180 Perusal of the Finite Element Method

mesh inside the film geometry.

During debonding phase:

according to

case). The nature of the body force can be seen in **Figure 8B**.

6/18 Jm6 (*le* is an equilibrium distance and has a value of 0.138 nm for the present

**Figure 8.** (A) Schematic diagram of the calculation of force between a contactor surface and a molecule inside the elas‐ tic film. (B) It shows the van der Waals body force present between the contacting surface and any nodal point in the

The contactor in the simulations is virtual in nature, and its presence is felt only through the body‐force term *f*. Since the contactor position keeps on changing with time (refer to **Figure 6B**) so also does the body force. Hence, the gap distance between a particular node

1 = ( + 1)/

2( + 1)/ + 15

1 ≤ t ≤ 2 where

*d* = - (11)

= - + + ´- - £ (10)

and the surface of the contactor keeps on changing with time as per the following:

( initial ( )) ( 2 2 ) ( ) 2 1 For t t *D d v tt ut c d*

analytically from linear stability analysis on sinusoidal substrates according to [30]

2 1 *c hA* pm

In the above set of equations, initial denotes the initial gap distance, 2 the normal position of the particular node at time *t*, the critical distance at which the film begins to experience substantial force due to the contactor and starts deforming. Its value can be roughly estimated

> ( ) 4 6.22

In VDLOAD subroutine of ABAQUS, the body force for a particular node is thus calculated

 b

During adhesion: = initial − × − 2 for 0 ≤ t ≤ 1 where

During resting: = initial − + 1 − 2 for

x

u2= curCoords (k,2) **\*\*\* curCoords(k,1) and curCoords(k,2) represent the updated coordinates in** *x***1 and** *x***2 direction respectively)**

**\*\* Calculate gap distance according to Eq. 10**

**\*\*\* Calculate forces according to Eq. 12**

**\*\* Calculate force update according to Eq. 13 and input the magnitude in** the array value (*k*), which is the prescribed **ABAQUS declared array to contain forces**

end do

The FE results displaying the effect of the substrate roughness on the morphologies formed, the work of adhesion exhibited, are discussed next in "Results" section.

#### **3.5. Results**

Linear stability analysis performed with the help of a single Fourier mode [16, 17] and energy‐ minimizing nonlinear analysis [18–20] has already been able to identify the nontrivial bifur‐ cation modes for which instability is possible. For critical conditions at which instability first initiates, it has been observed that the critical wave length λ 2.96 h is at par with the exper‐ imental observations. In FEM to obtain the critical wave mode, the geometry is created with the help of a series of cosine waves having random amplitudes and frequencies in the vicinity of a mean value *j*. When the particular geometries are simulated using VDLOAD, only in the adhesive branch, the displacements at the corner‐most node, which has the highest amplitude (because of the cosine nature of the waves), is tracked. The displacement of the node thus obtained at any generic time near critical separation distance is plotted against the average wave modes in **Figure 9A**. It can be seen that the curve has a bell shape indicating that near a particular mean surface wave mode the deformation of the surface is highest. This marks the dominant wave mode, whose value matches with those from Linear stability analysis (LSA) and experiments involving flat substrates.

**Figure 9.** (A) The displacement of the top‐most corner node of the film for different perturbed film surfaces at a partic‐ ular time near critical. The dominant wave mode is ~3 for flat substrate. (B) The dominant wave mode and the wave‐ length for different substrate amplitude averaged over the substrate wavelength. With increased roughness in substrate miniaturization of the surface patterns occurs.

**Figure 10.** The trajectory of the corner‐most point with time.

To study whether the substrate roughness has any influence on the emerging surface patterns, similar studies have been conducted for sawtooth‐profiled substrates having different substrate parameters of and *n*. The results reveal that though the substrate wavelengths have negligible effect, the substrate amplitudes have a profound influence on the emerging wavelengths. The surface wave mode/wave lengths are found to be an increasing/decreasing function of the substrate amplitude averaged over the different substrate wave modes considered, indicating that the patterned substrates do engender miniaturized length scales as can be seen from **Figure 9B**.

of a mean value *j*. When the particular geometries are simulated using VDLOAD, only in the adhesive branch, the displacements at the corner‐most node, which has the highest amplitude (because of the cosine nature of the waves), is tracked. The displacement of the node thus obtained at any generic time near critical separation distance is plotted against the average wave modes in **Figure 9A**. It can be seen that the curve has a bell shape indicating that near a particular mean surface wave mode the deformation of the surface is highest. This marks the dominant wave mode, whose value matches with those from Linear stability analysis (LSA)

**Figure 9.** (A) The displacement of the top‐most corner node of the film for different perturbed film surfaces at a partic‐ ular time near critical. The dominant wave mode is ~3 for flat substrate. (B) The dominant wave mode and the wave‐ length for different substrate amplitude averaged over the substrate wavelength. With increased roughness in

To study whether the substrate roughness has any influence on the emerging surface patterns, similar studies have been conducted for sawtooth‐profiled substrates having different

and experiments involving flat substrates.

182 Perusal of the Finite Element Method

substrate miniaturization of the surface patterns occurs.

**Figure 10.** The trajectory of the corner‐most point with time.

The displacement profile of the corner node at different time is shown in **Figure 10**. It can be seen that at initial time *t* = 0 since the initial gap is much greater than the critical value, the film surface remains planar with zero displacement. Near critical distance (near 1), the dis‐ placements start to have nonzero values and at just greater than 1 when the contactor is in resting position the film comes in total contact. The node continues to be in contact with the contactor even in the debonding phase when the contactor is pulled up and it gets de‐ tached only at snap‐off distance dsnap. The trajectory of the corner node arising due to dif‐ ferent debonding velocity is found to be different in the debonding phase. It can be seen that the snap‐off also occurs at different times. Though the snap‐off time t snap for the lower velocity is found to be highest, the snap‐off distance dsnap = vd × tsnap is in reverse order 30, 16, 5, 10, 6.6 nm showing that in case of lower debonding velocity the detachment oc‐ curs at the lowest snap‐off distance.

**Figure 11.** (A1‐E1) Morphologies at different contactor positions. (A2‐E2) Stresses in the film at different contactor po‐ sitions for a film cast on a sawtooth‐profiled substrate of amplitude = 0.1 and = 4.

In ABAQUS, not only a single node but the whole evolution of the film can be traced throughout adhesion‐debonding. In **Figure 11**, the labels A1‐E1 indicate the morphologies and A2‐E2 the stress profiles that develop across a film cast on a sawtooth‐profiled substrate (of amplitude = 0.1 and = 4) for the different contactor positions. In the inset, a contactor displacement graph is provided for easy understanding of the positions. At position A1 where the contactor is just near critical distance the surface is seen to begin destabilizing where the structures formed have very less asperities. At stage B1, the film surface comes in physical contact with the contactor and the top surface takes the shape of columns bridging the gap with intervening cavities (which form the dispersed phase). At C1, the end of rest‐ ing phase, the columnar structures mature to form continuous ridges, which upon with‐ drawal to D1 shrinks in size due to peeling from the edges. Continuous peeling ultimately gives rise to isolated columns (which now forming the dispersed phase) as shown in E1 leading to snap‐off. Thus, from adhesion to debonding, the morphologies at the surface go through a phase inversion. The corresponding stress profiles can be seen from **Figure 11** A2‐ E2. It can be seen that the stresses lie concentrated at the column edges. Thus, more number or greater width of columns (B2 and C2) are seen to have higher stresses at the edges which are responsible for subsequent peeling and stress release as seen in **Figure 11** D2 and E2.

**Figure 12.** Contour plots of the work of adhesion (J/m2 ) for different substrate parameters for a film cast on a rough (sawtooth‐profiled) substrate. *(A) For debonding velocity* = 0.125 /. *(B) For debonding velocity*, = 0.25 /.

The sum of the stresses across the total top surface area is a signature of the forces re‐ quired by the contactor to pull the film up to a certain pull‐off distance. The area under the force displacement curve from the beginning of debonding phase to snap‐off yields the work of adhesion, which is required to separate the two surfaces of the film and the con‐ tactor. **Figure 12** reveals the work of adhesion exhibited by films cast on substrates having various degrees of roughness. It can be seen vividly that while the effect of the wavelength of the substrate is negligible, the higher the substrate amplitude, more is the substrate roughness more is the work of adhesion and better is the performance of the elastic film as an adherent. The work of adhesion is also seen to be a function of the debonding velocity. For a higher debonding velocity (**Figure 12B**), it can be seen that the work of adhesion is also high (compared to **Figure 12A**). The fact is intuitive considering in band‐aids when we pull at higher velocity higher adhesiveness is felt compared to that felt at lower velocities [46].

## **4. Conclusions**

In ABAQUS, not only a single node but the whole evolution of the film can be traced throughout adhesion‐debonding. In **Figure 11**, the labels A1‐E1 indicate the morphologies and A2‐E2 the stress profiles that develop across a film cast on a sawtooth‐profiled substrate (of amplitude = 0.1 and = 4) for the different contactor positions. In the inset, a contactor displacement graph is provided for easy understanding of the positions. At position A1 where the contactor is just near critical distance the surface is seen to begin destabilizing where the structures formed have very less asperities. At stage B1, the film surface comes in physical contact with the contactor and the top surface takes the shape of columns bridging the gap with intervening cavities (which form the dispersed phase). At C1, the end of rest‐ ing phase, the columnar structures mature to form continuous ridges, which upon with‐ drawal to D1 shrinks in size due to peeling from the edges. Continuous peeling ultimately gives rise to isolated columns (which now forming the dispersed phase) as shown in E1 leading to snap‐off. Thus, from adhesion to debonding, the morphologies at the surface go through a phase inversion. The corresponding stress profiles can be seen from **Figure 11** A2‐ E2. It can be seen that the stresses lie concentrated at the column edges. Thus, more number or greater width of columns (B2 and C2) are seen to have higher stresses at the edges which are responsible for subsequent peeling and stress release as seen in **Figure 11** D2 and E2.

(sawtooth‐profiled) substrate. *(A) For debonding velocity* = 0.125 /. *(B) For debonding velocity*,

The sum of the stresses across the total top surface area is a signature of the forces re‐ quired by the contactor to pull the film up to a certain pull‐off distance. The area under the force displacement curve from the beginning of debonding phase to snap‐off yields the work of adhesion, which is required to separate the two surfaces of the film and the con‐ tactor. **Figure 12** reveals the work of adhesion exhibited by films cast on substrates having various degrees of roughness. It can be seen vividly that while the effect of the wavelength of the substrate is negligible, the higher the substrate amplitude, more is the substrate roughness more is the work of adhesion and better is the performance of the elastic film as an adherent. The work of adhesion is also seen to be a function of the debonding velocity.

) for different substrate parameters for a film cast on a rough

**Figure 12.** Contour plots of the work of adhesion (J/m2

= 0.25 /.

184 Perusal of the Finite Element Method

Opening of crack surfaces through traction‐separation laws in cohesive‐zone modeling is demonstrated for cohesive blocks. In the case of adhesion‐debonding of elastic films, simula‐ tions are performed for topologically rough substrates to demonstrate the ability of the method to resolve complex geometries, which these analytical methods are incapable of. At incipience of instability when the film surface is perturbed by random fluctuations at the vicinity of prescribed wavelengths, the displacements made by those contacting surface nodes are found highest near the dominant wavelength. These dominant perturbations grow into columnar structures at contact. A thorough FE study reveals that film cast on rough substrate is ideal for fabrication of features/columns at miniaturized length scales and exhibits high work of adhesion. Thus, these films can act as better adhesives and can produce enhanced functional surfaces.

#### **Nomenclature**

*j* Mean wave mode of surface

## **Acknowledgements**

*L* Lateral length

*λs* Substrate wavelength *n* Substrate wave mode *d* Film‐contactor gap distance

186 Perusal of the Finite Element Method

*va* Approach velocity *vd* Debonding velocity Stress tensor in the film

 Displacement vector *P* Pressure across the film **I** The identity matrix

Φ The shape function

Body force (per unit volume)

*x1*‐*x2* Cartesian coordinate system *Ψ* Attractive van der Waals potential

*M* The coefficient of attraction

*w* Modified interaction potential

*A* The Hamaker constant *B* Born‐repulsion coefficient *le* An equilibrium distance

Dintial The initial gap distance

*t2* Adhesive phase end time

*j* Mean wave mode of surface

*dc* Critical gap distance *u2* Normal displacement λ Surface wavelength

*t2* Resting phase end time of contactor

*t* Time

Ω The discretized domain inside the film

*r* Separation distance between two atoms or molecules

*ρc* The number density of atoms/molecules of the contactor

*ρf* The number density of atoms/molecules of the film

*ξ* Gap distance between a particular node and the surface of the contactor

*D* Surface‐to‐surface separation distance of D

*β* Nondimensional substrate amplitude

Funding from the Department of Science & Technology, New Delhi, India, is highly acknowl‐ edged.

## **Author details**

Jayati Sarkar1\*, Hemalatha Annepu2 and Satish Kumar Mishra3

\*Address all correspondence to<EMAIL_ADDRESS>
1 Department of Chemical Engineering, Indian Institute of Technology Delhi, New Delhi, India

2 Institute of Advanced Simulations‐2, Forschungszentrum Gmbh, Jülich, Germany

3 Department of Chemical Engineering, Mississippi State University, MS, USA

## **References**

[23] J. Y. Chung, K. Kim, M. K. Chaudhury, J. Sarkar, and A. Sharma. Confinement-induced instability and adhesive failure between dissimilar thin elastic films. Eur. Phys. J. E. 2006;20:47–53.

[7] F. Zhang, T. Nyberg, and O. Inganäs. Conducting polymer nanowires and nanodots

[8] D. Yoo, A. Wu, J. Lee, M. F. Rubner. New electro‐active self‐assembled multilayer thin films based on alternately adsorbed layers of polyelectrolytes and functional dye

[9] Z. Li, P. Dharap, S. Nagarajaiah, E. V. Barrera, and J. D. Kim. Carbon nanotube film

[10] H. G. Craighead. Nanoelectromechanical systems. Science. 2000;290:1532–1535.

[11] R. A. Chivers. Easy removal of pressure sensitive adhesives for skin applications. Int.

[12] A. Ghatak, L. Mahadevan, J. Y. Chung, M. K. Chaudhury, and V. Shenoy. Peeling from a biomimetically patterned thin elastic film. Proc. R. Soc. Lond.

[13] A. Ghatak, M. K. Chaudhury, V. Shenoy, and A. Sharma. Meniscus instability in a thin

[14] A. Ghatak and M. K. Chaudhury. Adhesion‐induced instability patterns in thin

[15] W. Mönch and S. Herminghaus. Elastic instability of rubber films between solid bodies.

[16] V. Shenoy and A. Sharma. Pattern formation in a thin solid film with interactions. Phys.

[17] J. Sarkar, V. Shenoy, and A. Sharma. Spontaneous surface roughening induced by surface interactions between two compressible elastic films. Phys. Rev. E.

[18] J. Sarkar, V. Shenoy, and A. Sharma. Patterns, forces, and metastable pathways in

[19] J. Sarkar, A. Sharma, and V. Shenoy. Adhesion and debonding of soft elastic films: crack

[20] J. Sarkar, A. Sharma, and V. Shenoy. Adhesion and debonding of soft elastic films on

[21] N. Arun, A. Sharma, V. B. Shenoy, and K. S. Narayan. Electric‐field‐controlled surface

[22] M. Gonguntla, A. Sharma, J. Sarkar, S. A. Subramanian, M. Ghosh, and V. Shenoy. Contact instability in adhesion and debonding of thin elastic films. Phys. Rev. Lett.

patterns, metastable pathways, and forces. Langmuir. 2005;21:1457–1469.

debonding of elastic films. Phys. Rev. Lett. 2004;93:018302 (4).

rough and patterned surfaces. J. Adhes. 2005;81:271–295.

instabilities in soft elastic films. Adv. Mater. 2006;18:660–663.

made with soft lithography. Nano Lett. 2002;2(12):1373–1377.

molecules. Synth. Metals. 1997;85:1425–1426.

elastic film. Phys. Rev. Lett. 2000;85: 4329–4332.

confined elastic film. Langmuir. 2003;19:2621– 2631.

sensors. Adv. Mater. 2004;16:640–643.

J. Adhes. Adhes. 2001;21:381–388.

Europhys. Lett. 2001;53:525–531.

Rev. Lett. 2001;86:119–122.

2003;67:031607 (11).

2006;97:018303(4).

A. 2004;460:2725–2735.

188 Perusal of the Finite Element Method

## **The Role of Finite Element Analysis in Studying Potential Failure of Mandibular Reconstruction Methods The Role of Finite Element Analysis in Studying Potential Failure of Mandibular Reconstruction Methods**

Raymond C.W. Wong, John S.P. Loh and I. Islam Raymond C.W. Wong, John S.P. Loh and I. Islam Additional information is available at the end of the chapter

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/64890

#### **Abstract**

[37] Q. Wang and X. Zhao. Phase diagrams of instabilities in compressed film-substrate

[38] F. Xu, M. Potier-Ferry, S. Belouettar, and Y. Cong. 3D finite element modeling for instabilities in thin films on soft substrates. Int. J. Solids Struct. 2014;51:3619–3632. [39] F. Xu, Y. Koutsawa, M. Potier-Ferry, and S. Belouettar. Instabilities in thin films on hyperelastic substrates by 3D finite elements. Int. J. Solids Struct. 2015;69:7–1-85. [40] D. S. Dugdale, Yielding of steel sheets containing slits. J. Mech. Phys. Solids. 1960;8:100–

[41] G. I. Barenblatt, The mathematical theory of equilibrium cracks in brittle fracture. Adv.

[42] A. Hillerborg, M. Modeer, and P. E. Petersson. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cem. Concr.

[43] V. Tvergaard. Effect of fibre debonding in a whisker-reinforced metal. Mater. Sci. Eng.

[44] A. Needleman and V. Tvergaard. An analysis of ductile rupture in notched bars. J. Mech.

[45] J. N. Israelachvili. Intermolecular and Surface Forces, New York, NY: Academic, 2011. [46] S. Kundu and E. P. Chan. Adhesion Behavior of Soft Materials in Engineering Bioma-

systems. J. Appl. Mech. 2014;81:510041–5100410.

104.

190 Perusal of the Finite Element Method

Appl. Mech. 1962;7:55–129.

Res. 1976;6:773–781.

A. 1990;125:203–213.

Phys. Solids. 1984;32:461–490.

terials for Regenerative Medicine, Springer, 2012.

Defects of the mandible occur after trauma or resection after infection or tumours. There have been many methods espoused, but many methods can fail especially if the biomechanics of the mandible is not considered fully. As the only moveable, loadbearing bone of the skull, the mandible is subject to loads and stresses unique to it due to its shape, location and function. This chapter reviews the basic knowledge of the mandible necessary to perform finite element analysis, the challenges and then reviews several studies that have been done. The authors' personal research is detailed to illustrate how finite element analysis can be used to look at potential failure of a new method for mandibular reconstruction and implant evaluation.

**Keywords:** mandible reconstruction, finite element analysis, potential failure

## **1. Introduction**

The lower jaw or mandible is the only load bearing, moveable bone of the skull. Defects in the lower jaw or mandible can happen as a result of trauma, infection or after a resection for pathology, which can be benign or malignant. The field of medicine or surgery has never been able to satisfactorily reconstruct a mandible after the original tissue has been gone. An unrepaired defect in the mandible, depending on the location of the defect would lead to (1) constriction of the remaining tissue around the defect leading to a malocclusion due to the pull of the muscles, (2) collapse of the arches leading to inability to eat and function and (3) difficulty swallowing or even breathing especially when lying flat due to loss of attachment of the tongue, sometimes leading to aspiration or asphyxiation.

© 2016 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2016 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## **2. Methods for mandibular reconstruction**

There have been many methods advocated for reconstruction of the mandible, most of which have never fully considered the biomechanical forces acting on the mandible both in the short and long term. These methods have been "tried", most times with catastrophic results. Among the methods for mandibular reconstruction advocated are:

The current gold standard is still the vascularized free flap, which needs a long operation time and harvest of tissue from another part of the body [1, 2].

This brings the question of "what is the ideal method of mandibular reconstruction?" to the fore.

## **2.1. Ideal method**

The ideal method of mandibular reconstruction would:

No such ideal method exists.

#### **2.2. Challenges to mandibular reconstruction**

The challenges to replacing a section of the mandible lie in its form and function as well as its unique location in the body. Any hardware (plates and screws) used to fix the mandible would undergo unique stresses not seen in other parts of the body. There is a non-axial load in that the long axis of the teeth is about perpendicular to the long axis of the mandible. The overlying tissue in the mouth is thin; any break in the soft tissue could lead to exposure of the hardware to the bacteria and saliva in the mouth, leading to the formation of biofilms. Ingestion of hot and cold food and liquids can cause expansion and contraction of the hardware, which is dissimilar to that of the underlying bone.

## **3. Mandibular biomechanics**

**2. Methods for mandibular reconstruction**

the methods for mandibular reconstruction advocated are:

and harvest of tissue from another part of the body [1, 2].

The ideal method of mandibular reconstruction would:

**5)** Not need to take tissue from another part of the body

**8)** Allow the patient to eat and function early or immediately

**7)** Not need prolonged hospitalization and recovery

**2.2. Challenges to mandibular reconstruction**

**1)** Reconstruct the missing soft tissue and bone

**6)** Newer methods like the endoprosthesis or alloplastic replacement

**3)** Mandibular reconstruction/bridging plate

**4)** Cancellous bone in titanium mesh

**7)** Tissue engineered bone scaffold

**2)** Allow replacement of teeth **3)** Not need a long surgery

No such ideal method exists.

**1)** Soft tissue flap

192 Perusal of the Finite Element Method

**2)** Autogenous bone blocks

**5)** Vascularized free flap

fore.

**2.1. Ideal method**

**6)** Cost effective

There have been many methods advocated for reconstruction of the mandible, most of which have never fully considered the biomechanical forces acting on the mandible both in the short and long term. These methods have been "tried", most times with catastrophic results. Among

The current gold standard is still the vascularized free flap, which needs a long operation time

This brings the question of "what is the ideal method of mandibular reconstruction?" to the

**4)** Easy to learn without needing extensive training and skills that are hard to learn

**9)** Be able to withstand the forces of biting and chewing permanently or for a long time

The challenges to replacing a section of the mandible lie in its form and function as well as its unique location in the body. Any hardware (plates and screws) used to fix the mandible would The biomechanics of the mandible has not been very well studied. Mandibular biomechanics can be studied with finite element analysis, strain gauges or photoelastic models. Finite element analysis allows study of the forces throughout an entire structure but is limited to static forces. An understanding of the anatomy and biomechanics of the intact mandible is necessary prior to looking at how to conduct a finite element analysis of the mandible.

## **3.1. Anatomy of the mandible**

The mandible is a U-shaped bone, connected at the temporomandibular joints at both ends to the skull. It consists of a corpus (body), the symphysis, which connects both right and left at the midline, the alveolar bone, which supports the teeth, the ramus with the condylar as well as the coronoid processes (**Figure 1**).

**Figure 1.** Schema showing the gross anatomy of the mandible.

The teeth are connected to the alveolar bone by periodontal ligaments, which act as a gomphosis (allowing minute movement). These ligaments are inserted into the bone on one end and into the cementum of the teeth.

The muscular attachments can be divided into three groups:

The outer surface of the mandible consists of dense cortical bone, the thickness of which varies. In certain areas, there is only cortical bone throughout. The alveolar processes and in the middle of the mandibular body and part of the ramus consist of cancellous bone (bone marrow). There is a nerve coursing through the mandible in a canal, which usually does not play a role in terms of biomechanics.

The attachment of the mandible to the skull consists of the attachments of the muscles of mastication and the temporomandibular joint. The temporomandibular joint consists of two joint spaces, the superior and inferior joint spaces, surrounded by a capsule consisting of elastic collagen fibres and divided by the articular cartilage, which is a fibrocartilage. The movement of the joint consists of two distinct movements, which are Phase I (rotation) about a hinge for the first 20 mm followed by Phase 2 (translation), which is affected mainly by the action of the lateral pterygoid muscle pulling the entire condyle to the front and out of the glenoid fossa onto the part of the zygomatic process of the temporal bone. The forward movement of the condyle is limited somewhat by a protrusion called the articular eminence.

## **3.2. Biomechanics of the intact mandible**

The mandible functions as a Class III lever. During function, there is a zone of tension on the alveolar part of the mandible and a zone of compression on the lower border. Studies by Meyer [3, 4] showed bone deformation in the mandibular condyle region, with tensile stress along the anterior ramus as well as along the sigmoid notch area and compressive stress along the posterior ramus border. This suggested that there is a tendency of the mandible to straighten during function. This somewhat simplistic model holds true when there is bilateral and equal function and bite forces (**Figure 2**).

Upon contraction of the muscles of mastication, the mandible is bent in a sagittal plane; this is produced by the vertical component of the muscle forces, the joint reaction forces and the reaction forces from chewing motions. During asymmetrical loading (biting on one side), the largest shear forces happen between the bite force and the muscle force on the working side (the biting side) and between the muscle force and joint force on the balancing side (non-biting side). This produces then a converse load distribution with a zone of tension on the lower border and a zone of compression on the alveolar portion in the working side and vice versa on the balancing side. This means that during incisal biting (biting on the front teeth), there is an equal amount of sagittal bending on both sides but a different deformation on the working and balancing sides during molar biting [5–7].

**Figure 2.** Zones of tension and compression in the mandible.

The teeth are connected to the alveolar bone by periodontal ligaments, which act as a gomphosis (allowing minute movement). These ligaments are inserted into the bone on one end

**1)** Muscles of mastication—medial pterygoid, lateral pterygoid, masseter and temporalis; only the lateral pterygoid muscle assist in opening the mandible, the rest closes the

**2)** Suprahyoid muscles that assist in some opening of the mouth and swallowing—mylo-

**3)** Muscles of facial expression that insert into the mandible—buccinators, depressor anguli

The outer surface of the mandible consists of dense cortical bone, the thickness of which varies. In certain areas, there is only cortical bone throughout. The alveolar processes and in the middle of the mandibular body and part of the ramus consist of cancellous bone (bone marrow). There is a nerve coursing through the mandible in a canal, which usually does not

The attachment of the mandible to the skull consists of the attachments of the muscles of mastication and the temporomandibular joint. The temporomandibular joint consists of two joint spaces, the superior and inferior joint spaces, surrounded by a capsule consisting of elastic collagen fibres and divided by the articular cartilage, which is a fibrocartilage. The movement of the joint consists of two distinct movements, which are Phase I (rotation) about a hinge for the first 20 mm followed by Phase 2 (translation), which is affected mainly by the action of the lateral pterygoid muscle pulling the entire condyle to the front and out of the glenoid fossa onto the part of the zygomatic process of the temporal bone. The forward movement of the

The mandible functions as a Class III lever. During function, there is a zone of tension on the alveolar part of the mandible and a zone of compression on the lower border. Studies by Meyer [3, 4] showed bone deformation in the mandibular condyle region, with tensile stress along the anterior ramus as well as along the sigmoid notch area and compressive stress along the posterior ramus border. This suggested that there is a tendency of the mandible to straighten during function. This somewhat simplistic model holds true when there is bilateral and equal

Upon contraction of the muscles of mastication, the mandible is bent in a sagittal plane; this is produced by the vertical component of the muscle forces, the joint reaction forces and the reaction forces from chewing motions. During asymmetrical loading (biting on one side), the largest shear forces happen between the bite force and the muscle force on the working side (the biting side) and between the muscle force and joint force on the balancing side (non-biting side). This produces then a converse load distribution with a zone of tension on the lower

condyle is limited somewhat by a protrusion called the articular eminence.

and into the cementum of the teeth.

play a role in terms of biomechanics.

**3.2. Biomechanics of the intact mandible**

function and bite forces (**Figure 2**).

mandible.

194 Perusal of the Finite Element Method

oris, mentalis.

The muscular attachments can be divided into three groups:

hyoid, hyoglossus, genioglossus, digastric muscles.

There is also a tendency for narrowing of the mandibular arch from parasagittal and transverse deformation upon clenching and incisal biting. This is caused by bilateral torsion of both mandibular bodies and bending at the symphyseal region leading to compression at the superior margin of the symphysis and tension at the inferior margin.

Hylander showed that the mandibular symphysis undergoes three distinct patterns of stress and deformation, that is, corporal rotation (relative outward rotation of both halves of the mandible), medial convergence (change in mandibular width during function) and dorsoventral shear (movement of both mandibular halves relative to one another in the vertical plane) [8] (**Figure 3**).

There is also some difference in the deformation between the outer (buccal) surface of the mandible and the inner (lingual) surface. Lateral transverse bending occurs and the bending moment increases from back to front during the late power stroke of biting/clenching. The maximum magnitude of the bending occurs near the symphysis. This bending produces compressive stresses at the buccal cortex and tensile stress at the lingual cortex. The deformation has been calculated to be as large as 0.6 mm in a simulated molar bite of 526N using finite element analysis [6]. The mandible deformed in a helical pattern upwards and towards the working side, with regions of high tensile stress (15–25 MPa) from the coronoid process and ramus towards the lingual side of the symphysis. The highest value of compressive stress (15– 25 MPa) was at the bite point and bilateral sigmoid notches, at the working side angle and in an area from the posterior surface of the balancing side ramus running to the lower border of the body till the symphysis. This then runs up to the buccal side from the inferior until the bite point. Overall, the shear stresses were larger on the working side with the exception of the balancing side condyle (peak shear stress of 25 MPa). In a nutshell, this means that the mandible changes in dimensions during function as a result of its shape and muscle pull with the greatest change in dimension and deformation in the midline.

What does all this information mean and what is the practical application? It means that application of any hardware to the mandible must take into account these forces and change in dimension. This has led to the creation of Champy's ideal lines of internal fixation for fixation of a mandibular fracture [9] (**Figure 4**). Placement of bone plates and screws in the area between the zones of tension and zones of compression will tend neutralize the forces and stabilize the bony fragments enough for healing. This however uses the principle of cross bracing and load sharing, that is, using inter-fragmentary bone friction to help stabilize the bony segments.

**Figure 3.** (A) Tendency for mandible to straighten as well as undergo torsion. (B) Forces acting about the midline of the mandible. CR – corporal rotation; MC – medial convergence; DVS – dorso-ventral shear.

In comminuted fractures or in a segmental mandibular defect, the principle of load sharing cannot be applied. We are then dependent on using load bearing bone plates, which means the material strength of the plate is the only thing keeping the mandibular segments together. The method of fixation of the bone plate to the bone also plays a factor.

**Figure 4.** Champy's lines for ideal internal fixation – placement of bone plates along these lines will stabilize the bone fragments enough for healing of fractures.

## **4. Finite element model of the mandible**

A very accurate finite element model is practically impossible to create due to the complex anatomy of the mandible. A true to accurate model that takes everything into account would need several supercomputers. Assumptions will have to be made to simplify the model and reduce the need for computing power.

Several authors have created finite element models of the mandible, each with different levels of complexity [10, 11]. Following are the steps taken to create an accurate finite element model [12].

## **4.1. Creation of a 3D mesh**

point. Overall, the shear stresses were larger on the working side with the exception of the balancing side condyle (peak shear stress of 25 MPa). In a nutshell, this means that the mandible changes in dimensions during function as a result of its shape and muscle pull with the greatest

What does all this information mean and what is the practical application? It means that application of any hardware to the mandible must take into account these forces and change in dimension. This has led to the creation of Champy's ideal lines of internal fixation for fixation of a mandibular fracture [9] (**Figure 4**). Placement of bone plates and screws in the area between the zones of tension and zones of compression will tend neutralize the forces and stabilize the bony fragments enough for healing. This however uses the principle of cross bracing and load sharing, that is, using inter-fragmentary bone friction to help stabilize the bony segments.

**Figure 3.** (A) Tendency for mandible to straighten as well as undergo torsion. (B) Forces acting about the midline of the

mandible. CR – corporal rotation; MC – medial convergence; DVS – dorso-ventral shear.

change in dimension and deformation in the midline.

196 Perusal of the Finite Element Method

The information about the external shaped of the mandible is needed to be able to create a mesh to input into the FEA software. Several possible methods that have been used to get the geometry of the mandible:

Once the 3D geometry has been obtained, it is subdivided into a finite, large number of geometrically simplified elements, connected together at the nodes. The mesh is a contiguous collection of these simple-shaped elements. Most FE software have automated mesh generation features, creating relatively dense meshes, which can be refined in different regions.

## **4.2. Input of material properties and boundary conditions**

The material properties of the elements, namely the elastic modulus and Poisson's ratio must be defined. Since the mandible consists of cortical and cancellous bone together with the teeth, several assumptions need to be made to simplify the model further. It is well known that bone is anisotropic in different dimensions. For purposes of simplifying the calculations, most authors have tended to assume that bone is isotropic and that the mandible is purely cortical bone. Any teeth present, which in real life would have periodontal ligaments and allows minute movements would tend to be assumed to be ankylosed to the bone, that is, fused to the bone. The teeth are composed of enamel on the outer surface of the crown followed by dentine and the pulp (a hollow cavity which contains the nerve fibres and blood supply). The outer surface of the root is covered by cementum. All these tissues have different material properties. The values for the material properties are have already been determined by studies. Some errors abound when it comes to the values of the cancellous bone. Most studies have tended to remove a block of cancellous bone and then subject the block to mechanical testing to ascertain its material properties. There is evidence to suggest that this is not entirely accurate. Misch et al. (1999) showed that the presence of cortical bone increases the elastic modulus of cancellous bone [14]. When the cortical bone was present, the elastic modulus ranged from 24.9 to 240 MPa (mean 96.2 MPa). When cancellous bone only was tested, the elastic modulus reduced dramatically (3.5–125.6 MPa). This means that the values from the literature for cancellous bone are not accurate.

Boundary conditions are important to prevent movement of the individual units so that the model can be loaded and deformed as a rigid structure, allowing computations to be performed. It can be divided into essential boundary conditions (displacement constraints to anchor the model and the non-essential boundary conditions or loading conditions, which are the forces to be applied to the model). Decisions will need to be made also about the insertion of the muscle forces and the force of each muscle. Some muscles like the masseter have three distinct types of fibres with different vectors. Most of these values are already in the literature.

## **4.3. Solution**

**1)** Digital creation of 3D model—not accurate representation of the anatomical detail. The question then arises: "At what level of intricacy would there be detriment to the level of

**2)** Conversion of digitized slices or sections of a human mandible into a whole 3D structure