Modeling the Local Viscoelastic Behavior of Living Cells Under Nanoindentation Tests

Abstract

Aiming at multiscale numerical investigations of fibrous connective tissues, the present work proposes a finite strain viscoelastic model suitable to represent the local mechanical response of living cells in conjunction with finite element simulations of nanoindentation tests. The material model is formulated in a thermodynamically consistent framework based on a variational constitutive approach. As a case of study, a numerical investigation of the local compressible response of fibroblast cells is addressed. Moreover, a set of cyclic, stress relaxation and creep experiments are simulated in order to investigate the capability of the model to predict rate dependence of cells, showing sound agreement with experimental data. The proposed model may be used as a suitable tool for a better understanding of stiffness and energetic dissipation of cells.

Keywords:
Living cells; Viscoelasticity; Constitutive model; AFM nanoindentation; Multiscale

1 INTRODUCTION

Withstanding mechanical loadings is the primary function of many biological connective tissues such as tendons, ligaments and cartilage. The mechanical response to external loadings are highly dependent on their hierarchical microstructure, cellular organization and interactions between one another (Lavagnino et al., 2015).

On its hierarchy inside the tissue, cells are subjected to a complex mechanical environment (Qi et al., 2006). These mechanical stimuli are transduced by cells into biochemical signals leading to tissue adaptation. The particular mechano-chemical pathways that induce changes in tissue morphology are usually known as mechanotransduction mechanisms (Lavagnino et al., 2015; Wang, 2006).

Regardless their physiological functions, living cells are structural units holding particular nonlinear behaviors, including viscous and plastic effects (Bonakdar et al., 2016; Haase and Pelling, 2015; Nawaz et al., 2012). Accordingly, cells may present important mechanical contributions in the macroscopic response of tissues, where multiscale analysis can be employed to investigate these micromechanical behaviors.

In multiscale numerical analyses based on representative volume elements (RVE), the homogenized behavior relies strongly on how accurate the behavior of each material phase is represented by corresponding specific material models. In addition, these models must show sound agreement with experimental data at that scale, specifically in cases where these phases present nonlinear behavior and are subjected to finite strains.

Several measurement techniques have been employed to assess the mechanical response of living cells, such as micropipette aspiration, magnetic and optical tweezers and atomic force microscopy (AFM) indentation, also known as nanoindentation test (Haase and Pelling, 2015; Janmey and McCulloch, 2007; Kollmannsberger and Fabry, 2011). Among these experimental methods, the later provides a convenient structure to finite element simulations. For instance, whereas the cell sample can be considered locally as a homogenous medium, the tip of the indenter is modeled as a rigid body.

With these motivations in mind and aiming future numerical multiscale investigations of fibrous soft tissues, the present work proposes a finite strain viscoelastic model suitable to represent the local mechanical response of living cells in conjunction with finite element simulations of nanoindentation tests. The proposed model is formulated in a thermodynamically consistent framework based on the variational formalism addressed in Ortiz and Stainier (1999) and Radovitzky and Ortiz (1999), particularized in Fancello et al. (2006) and Vassoler et al. (2012) to viscoelastic materials. In addition, an alternative incremental solution strategy is proposed.

This manuscript is organized as follows. Section 2 presents the theoretical background employed to model the viscoelastic behavior of cells. Aiming finite element simulations, the continuum constitutive equations are discretized in time, leading to a variational constitutive update algorithm, shown in section 3. In order to verify if the material model is able to represent the available experimental data, section 4 addresses a set of numerical simulations of nanoindentation experiments of fibroblast cells. In addition, the compressible and viscoelastic responses of the model are studied. Particularities of the finite element simulations and further discussions about the mechanical behavior of cells are highlighted in section 5. Finally, appendices provide operational details and the analytical expression for the consistent material tangent modulus.

2 CONTINUUM FORMULATION

2.1 Constitutive Modeling of Dissipative Materials in a Variational Framework

The constitutive modeling of dissipative materials can be cast in a thermodynamically consistent framework by defining a potential function, in present case the Helmholtz free energy ψ(F, α), and a so-called dissipation pseudo-potential function ϒ(de Souza Neto et al., 2009; Jirásek and Bazant, 2002). The argument of the free energy, in the case of an isothermal process, is the thermodynamic state defined by the deformation gradient F and a set of internal variables α. The dissipation potential ϒ, on the other hand, depends also on the corresponding rates , and , where the semicolon in the argument list indicates the state variables play a role of parameters. For the sake of clarity, only the dependence on the rates is kept henceforward in the notation of this potential., F, α) (

It is hypothesized that = (F, α), which allows the following decoupling:

(2.1)

Based on these assumptions, it is shown in Ortiz and Stainier (1999) and Radovitzky and Ortiz (1999) that a potential 𝒫 of the form

(2.2)

may be defined such that the optimality condition of the minimization problem

(2.3)

provides a thermodynamically consistent equation for the evolution of the internal variables:

(2.4)

The first Piola-Kirchhoff stress tensor is computed from the first partial derivative of the rate potential with respect to also results in the first Piola-Kirchhoff stress tensor:. Moreover, the rate potential evaluated at the minimizer argument of (2.2) defines a reduced potential 𝒫red, whose total derivative with respect to

(2.5)

Potential 𝒫^red plays an important role in present framework, since its incremental counterpart presents a compelling mathematical structure for the calculation of the incremental updates of the internal variables.

2.2 Thermodynamic Potentials for a Class of Viscoelastic Models

The classic multiplicative decomposition of the deformation gradient F is used herein to separate elastic F ^e and viscous F ^v contributions:

(2.6)

Based on this, total, elastic and viscous right Cauchy-Green strain tensors are given, respectively, by

(2.7)

Similarly, total and viscous rate of deformations are formally defined as

(2.8)

where the assumption of null inelastic spin is used (Anand and Gurtin, 2003; Gurtin and Anand, 2005).

The proposed model follows the rheological assembly illustrated in Figure 1. According to this, the set of internal variables is reduced to α: = {F ^v}, and the Helmholtz free energy and dissipation potential take the form:

(2.9)

Figure 1:
Schematic representation of the rheological model.

The superscript ∞ represents the time-independent response. As can be noted, the dissipation potential φ of equation (2.1) is null. Moreover, the arguments of the potentials (2.9) are objective entities (Gurtin et al., 2010).

In view of (2.9), the rate potential (2.2) reduces to the expression:

(2.10)

with d ^v depending on by means of (2.8). Equation (2.10) provides a continuum background for the variational approach introduced in section 2.1. The incremental counterpart of the constitutive equations and numerical issues are addressed in the next section.

3 INCREMENTAL FORMULATION

3.1 Incremental Potential and Updating Rule for the Internal Variable

An incremental counterpart of the rate potential (2.10) is used to directly compute the incremental internal variable updates needed in finite element schemes. Considering a time increment ∆t = _{^tn +1} - _^tn , a possible general expression for this incremental potential is

(3.1)

where the variables d ^v and are discrete approximations for the rates . One can note that ϕv is evaluated at an intermediate time tn +ϑ inside the time increment, where the parameter ϑ ∈[0, 1] is closely related to the discretization rule employed for . The consistency of the incremental form (3.1) can be verified if ∆t → 0 ⇒𝒫inc = 𝒫. and

The most common option to obtain an update expression for the internal variable Fvn+1, in function of the discrete tensor dv, is the well-known exponential mapping (Weber and Anand, 1990):

(3.2)

This choice has been extensively used in combination with the spectral parameterization of Fancello et al., 2006; Fancello et al., 2008; Vassoler et al., 2012). Moreover, under an incompressible viscous flow assumption, the incremental relation tr(F ^v _{n +1}) =1 is automatically fulfilled.) = 0 ⇔ det ( to handle arbitrary isotropic hyperelastic potentials or in combination with Hencky’s isotropic model to obtain some numerical advantages (

Another alternative expression to obtain the updated value F ^v _{n +1} is based on the classical Euler integration rule:

(3.3)

(3.4)

When ϑ = 1, a fully implicit integration scheme is recovered. In this case, substituting equation (3.4) into (3.3) leads to the following updating rule:

(3.5)

where I is the second order identity tensor.

In view of the parameterization (3.5), the resulting variational problem is formulated in a fully tensorial format. While this implies the simultaneous evaluation of a set of independent components, the formulation remains able to handle anisotropic behaviors. Moreover, neither the spectral decomposition of , nor the first and second derivatives of these spectral quantities are needed. In addition, if the viscous flow is considered incompressible, it is enforced into the variational problem by means of a Lagrangian functional (further details in next section).

To end this section it is worth mentioning that, even both approaches are quite different, is not the goal of the present work to compare and quantify its numerical advantages or drawbacks.

3.2 Constitutive Algorithm

The local constitutive algorithm is composed by two main stages: the solution of the variational principle and the stress evaluation. Accordingly, the general concepts introduced in section 2 are particularized to the present class of viscoelastic models.

3.2.1 Solution of the Variational Principle

In view of the parameterization (3.5), the viscous rate of deformation becomes the primal variable of the variational problem. Moreover, in the present modeling, the viscous flow is considered incompressible. Therefore, the minimum principle (2.3) can be presented as

(3.6)

where ℐso: = {ym |det[F ^v _{n +1}(ym the space of symmetric second order tensors. The kinematic constraint J ^v _{n +1}:= det(F ^v _{n +1})=1 is taken into account by means of the Lagrangian functional)] =1} represents the isochoric viscous space and 𝒮 ∈ 𝒮

(3.7)

where ϒ_{n +1} is a Lagrange multiplier. Consequently, the minimum principle (3.6) can be rewritten as an unconstrained problem, such that,

(3.8)

The solution of (3.8) defines the internal variable updates algorithm. Once the optimal solution F ^v _{n +1} is updated by equation (3.5). The proposed solution strategy is based on the full Newton-Raphson procedure, where further operational details are addressed in ^{appendix A}. It is important to note that the symmetric nature of Itskov (2002) and Jog (2007) for further details about directional derivatives in symmetric spaces). In this particular case, the solutions of the principle (3.8) always return symmetric tensors is obtained from (3.8), the internal variable is not included as a constraint into the problem (3.8). However, its symmetry is considered throughout the directional derivatives required by the Newton’s procedure, resulting in suitable Fréchet operators (see .

3.2.2 Stress Evaluation and Consistent Material Tangent Modulus

Once the solution of (3.8) is obtained, the stress can be updated as follows. Taking into account the incremental potential (3.1), the incremental first Piola-Kirchhoff stress tensor is given by

(3.9)

The partial derivatives of the Helmholtz strain energies in relation to F _{n +1} result in

(3.10)

where S _{n +1} is the second Piola-Kirchhoff stress tensor. In addition, the time-independent and elastic second Piola stresses are defined as

(3.11)

It should be noted that, due to the directional derivative of equation (3.9), the elastic second Piola-Kirchhoff stress tensor in (3.10) is consistently pulled back to the referential configuration.

Within the framework of a conventional nonlinear finite element code, the consistent tangent modulus must be provided (Simo and Taylor, 1985). Therefore, a closed form for the consistent material tangent modulus is provided in ^{appendix B}.

3.3 Choice of Specific Potentials

A common assumption in the biomechanical modeling of cells is considering them as an incompressible material. However, under physiological conditions, cells can experience large volume changes (Hamann et al., 2002; Zlotek-Zlotkiewicz et al., 2015). Based on these observations, the following Neo-Hookean and quadratic potentials are chosen to represent the viscoelastic behavior of living cells:

(3.12)

where {µ ^∞ , κ ^∞ , µ ^e , η ^v} is the set of constitutive parameters to be identified according to experimental data. In the expression ψ ^∞, the parameter κ ^∞ controls changes in volume, while µ ^∞ accounts for distortional effects. To ensure an incompressible behavior, κ ^∞ typically takes values of ((10³ - 10⁴) µ ^∞ (Bonet and Wood, 2008).

4 NUMERICAL SIMULATIONS AND RESULTS

4.1 Compressibility of Fibroblast Cells

AFM indentation experiments were carried out by Nawaz et al. (2012) in fibroblast cells. Authors reported cyclic experiments under indentation depth of (0.38 ± 0.02 μm at 0.8 μm/s indentation rate, retrieving 72 ± 5 pN indentation forces. The relative energy dissipated by the hysteresis loop reached a value of (30±5% (Figure 3-a).

Figure 2:
Finite element model employed for the simulations of nanoindentation tests. (a) Perspective view of the model showing the spherical indenter and the cell sample in a quarter symmetry. (b) Top view of the numerical sample emphasizing the mesh refinement.

Figure 3:
(a) Comparison between experimental and numerical results obtained from the nanoindentation tests. The experimental points were plotted based on data provided by Nawaz et al. (2012). (b) Sensitivity of the numerical results in relation to the incompressible response of the model. The related constitutive parameter sets are shown in Table 1.

Numerical simulations of these experiments were run with a rigid 1.98 μm diameter spherical indenter and a numerical sample with diameter of 10 μm and height of 4 μm. As can be seen from Figure 2, the cell sample is considered locally homogeneous and it was discretized by 8, 130 quadratic hexahedral elements with reduced integration. Due to symmetry, only a quarter of the sample was modeled. Moreover, frictionless contact was considered in the interface between sample and indenter.

The finite element simulations presented herein were run in the software Abaqus, where the proposed viscoelastic model was implemented into the user-subroutine UMAT.

Using this finite element setup, a set of numerical simulations evaluated the capability of present material model to represent the experimental data, considering different compressibility levels. To this aim, for increasing values of κ ^∞, the parameter µ ^∞ was kept fixed and the viscoelastic parameters µ ^e and η ^v were adjusted in order to fit the numerical curves to the range of the experimental data (see Table 1).

Table 1:
Constitutive parameters related to the numerical curves shown in Figure 3.

Potential	Parameter	Sim. 1	Sim. 2	Sim. 3	Sim. 4
ψ ∞	μ∞ [Pa]	20	20	20	20
	κ∞ [Pa]	80	800	8,000	80,000
ψ e	μe[Pa]	60	48	40	40
ϕ v	ηv[Pa · s]	80	58	50	50

Figure 3-a compares the force-displacement curves of the indenter obtained from the experiments and Simulation 1. Figure 3-b show the corresponding results obtained with all simulations. Finally, Figure 4 depicts the displacement and the von-Mises (Cauchy) stress fields at maximum penetration of the indenter for Simulation 4.

Figure 4:
Displacement (a) and von-Mises stress (b) fields at the instant of the maximum penetration of the indenter for Simulation 4.

In order to gain a better understanding of how the parameters of Table 1 affect the incompressible response of the model, simulations of homogeneous uniaxial tensile tests were performed. Figure 5-a display the uniaxial first Piola-stretch curves, while the corresponding time history of the Jacobians are shown in Figure 5-b. All simulations were run under an engineering strain rate of 0.05s^-1.

Figure 5:
Stress-stretch curves (a) and time history of the volumetric jacobians (b) obtained from simulations of homogeneous uniaxial tensile tests. The related constitutive parameter sets are shown in Table 1.

4.2 Cyclic, Stress Relaxation and Creep Tests

Identification procedures often recover material parameters that successfully reproduce the used set of experimental data. However, it is not unusual these parameters fail to predict the behavior of the same material under a mechanical condition different from (but close to) that identified. In order to gain confidence on the behavior of the proposed model using the identified parameters, a brief sensitivity analysis considering different mechanical tests was carried out.

The same finite element model of the previous section with constitutive parameters of Simulation 1 (Table 1) was used.

A first set of numerical tests considering cyclic indentations with 0.5 μm depth at indentation rates of 0.1, 0.5, 1.0 and 10 μm/s, was run. The corresponding force-displacement curves are displayed in Figure 6-a. Different hystereses are seen, depending on the indentation rate. The relative percentage of the energy dissipated at each hysteresis loop is shown in Figure 6-b. Moreover, a quadratic trend curve is fitted in order to better illustrate how the dissipated energy changes within the range of indentation rates.

Figure 6:
(a) Cyclic nanoindentation tests at different rates. (b) Relative loss of energy by the hysteresis loops computed from cyclic curves.

Besides cyclic tests, stress relaxation and creep AFM indentation tests are also found in the literature as efficient experiments to get information on the viscoelastic properties of cells (Haase and Pelling, 2015; Moreno-Flores et al., 2010). Based on this motivation, a set of relaxation and creep simulations under different load conditions were tested, as shown in Figure 7. Concerning the relaxation tests, three indentation depths of 0.1, 0.3 and 0.5 μm were applied at a 10 μm/s rate, and then the indenter was kept fixed for 10 seconds. The corresponding relaxation curves are displayed in Figure 7-a. Concerning the creep tests, three indentation forces of 10, 25 and 50 pN were progressive and linearly applied in 0.05 seconds, and then kept constant for 16 seconds. The corresponding creep curves are shown in Figure 7-b.

Figure 7:
Mechanical responses predicted by the viscoelastic model when submitted to nanoindentation tests of stress relaxation (a) and creep (b) under different loading conditions.

5 DISCUSSIONS AND FINAL REMARKS

As can be seen from Figure 3, all numerical simulations were able to predict the local viscoelastic response of fibroblast cells under nanoindentation (in view of the standard deviations) using spherical indenter. In addition, the displacement and stress fields depicted in Figure 4 emphasize the numerical results seem to present minimum influence of the boundaries of the sample.

Concerning to the numerical curves plotted in Figure 3-b, it can be seen that when incompressibility is enforced, the force-indentation response changes its curvature within the firsts 0.1 μm displacement, diverging from the usual expected behavior. Based on these numerical results, one can hypothesize that cells seem to respond locally, under indentation tests, as a compressible material (see time histories of the volumetric jacobians plotted in Figure 5-b). However, further numerical and experimental investigations are required to support this hypothesis.

Taking into account the cyclic numerical experiments displayed in Figure 6-a one can observe that, for the identified material parameters, the dissipation decreases when the indentation rate increases. The energetic loss is quantified by the hysteresis loop as shown in Figure 6-b. Moreover, one can see the percentage of hysteresis presents a nearly quadratic decay over rates.

According to the relaxation and creep simulations shown in Figure 7, the material reaches equilibrium states beyond 5 seconds in stress relaxation and 10 seconds in creep. These time scales predicted by the model show sound agreement with experimental data reported in Moreno-Flores et al. (2010).

It is worth mentioning the use of analytical expressions based on Hertz contact model and indentation tests to estimate the cell stiffness is common in the literature (technical details may be found elsewhere, Vinckier and Semenza (1998)). This approach, however, is quite sensitive to the indentation depth, since the Hertz expressions were tailored for linear elastic materials and infinitesimal strains. In addition, it is unable to consider rate dependent effects. In contrast to these analytical approaches, the finite element simulation procedure presented in this work seems to provide a more accurate tool for a better understanding of stiffness and energetic dissipation of cells.

Within a multiscale numerical framework of biological tissues, the present constitutive model may be employed to represent the microscopic material phases composed of cells and/or other isotropic viscoelastic components as well, a necessary step for future developments in the mechanobiology field, where the knowledge of the complex mechanical responses of cells are issues of major concern.

Acknowledgments

The authors would like to thank the financial support provided by the Brazilian funding agencies CAPES - (Coordination for the Improvement of Higher Education Personnel) and CNPq - (National Council for Scientific and Technological Development).

References

• Anand, L., Gurtin, M.E., (2003). A theory of amorphous solids undergoing large deformations with applications to polymers and metallic glasses, International Journal of Solids and Structures 40, 1-29.
• Belytschko, T., Liu, W., Moran, B., (2000). Nonlinear Finite Elements for Continua and Structures, Chichester, New York, John Wiley. John Wiley & Sons, Ltd.
• Bonakdar, N., Gerum, R., Kuhn, M., Spörrer, M., Lippert, A., Schneider, W., Aifantis, K.E., Fabry, B., (2016). Mechanical plasticity of cells. Nature Materials 1-6.
• Bonet, J., Wood, R.D., (2008). Nonlinear continuum mechanics for finite element analysis, 2nd Editio. ed, Zhurnal Eksperimental’noi i Teoreticheskoi Fiziki. Cambridge University Press.
• de Souza Neto, E.A., Peric, D., Owen, D.R.J., (2009). Computational Methods for Plasticity: Theory and Applications, Engineering.
• Fancello, E., Ponthot, J.-P., Stainier, L., (2006). A variational formulation of constitutive models and updates in non-linear finite viscoelasticity. International Journal for Numerical Methods in Engineering 65, 1831-1864.
• Fancello, E., Vassoler, J.M., Stainier, L., (2008). A variational constitutive update algorithm for a set of isotropic hyperelastic-viscoplastic material models. Computer Methods in Applied Mechanics and Engineering 197, 4132-4148.
• Gurtin, M., Fried, E., Anand, L., (2010). The mechanics and thermodynamics of continua. Cambridge University Press.
• Gurtin, M.E., Anand, L., (2005). The decomposition F = FeFp, material symmetry, and plastic irrotationality for solids that are isotropic-viscoplastic or amorphous. International Journal of Plasticity 21, 1686-1719.
• Haase, K., Pelling, A.E., (2015). Investigating cell mechanics with atomic force microscopy. Journal of The Royal Society Interface 12.
• Hamann, S., Kiilgaard, J.F., Litman, T., Alvarez-Leefmans, F.J., Winther, B.R., Zeuthen, T., (2002). Measurement of Cell Volume Changes by Fluorescence Self-Quenching. Journal of Fluorescence 12, 139-145.
• Itskov, M., (2002). The derivative with respect to a tensor: some theoretical aspects and applications. ZAMM - Journal of Applied Mathematics and Mechanics 82, 535-544.
• Janmey, P.A., McCulloch, C.A., (2007). Cell mechanics: integrating cell responses to mechanical stimuli. Annual review of biomedical engineering 9, 1-34.
• Jirásek, M., Bazant, Z., (2002). Inelastic analysis of structures. John Wiley & Sons, Ltd.
• Jog, C.S., (2007). Foundations and Applications of Mechnics Volume I: Continuum Mechanics, 2 edition. ed. Alpha Science Intl Ltd.
• Kollmannsberger, P., Fabry, B., (2011). Linear and Nonlinear Rheology of Living Cells. Annual Review of Materials Research 41, 75-97.
• Lavagnino, M., Wall, M.E., Little, D., Banes, A.J., Guilak, F., Arnoczky, S.P., (2015). Tendon mechanobiology: Current knowledge and future research opportunities. Journal of Orthopaedic Research 33, 813-822.
• Moreno-Flores, S., Benitez, R., Vivanco, M. dM, Toca-Herrera, J.L., (2010). Stress relaxation and creep on living cells with the atomic force microscope: a means to calculate elastic moduli and viscosities of cell components. Nanotechnology 21.
• Nawaz, S., Sánchez, P., Bodensiek, K., Li, S., Simons, M., Schaap, I.A.T., (2012). Cell Visco-Elasticity Measured with AFM and Optical Trapping at Sub-Micrometer Deformations. PLoS ONE 7.
• Ortiz, M., Stainier, L., (1999). The variational formulation of viscoplastic constitutive updates. Computer Methods in Applied Mechanics and Engineering 7825, 419-444.
• Qi, J., Fox, A.M., Alexopoulos, L.G., Chi, L., Bynum, D., Guilak, F., Banes, A.J., (2006). IL-1beta decreases the elastic modulus of human tenocytes. Journal of applied physiology 101, 189-195.
• Radovitzky, R., Ortiz, M., (1999). Error estimation and adaptive meshing in strongly nonlinear dynamic problems. Computer Methods in Applied Mechanics and Engineering 172, 203-240.
• Simo, J.C., Taylor, R.L., (1985). Consistent tangent operators for rate-independent elastoplasticity. Computer Methods in Applied Mechanics and Engineering 48, 101-118.
• Vassoler, J.M., Reips, L., Fancello, E.A., (2012). A variational framework for fiber-reinforced viscoelastic soft tissues. International Journal for Numerical Methods in Engineering 89, 1691-1706.
• Vinckier, A., Semenza, G., (1998). Measuring elasticity of biological materials by atomic force microscopy. FEBS Letters 430, 12-16.
• Wang, J.H.C., (2006). Mechanobiology of tendon. Journal of Biomechanics 39, 1563-1582.
• Weber, G., Anand, L., (1990). Finite deformation constitutive equations and a time integration procedure for isotropic, hyperelastic-viscoplastic solids. Computer Methods in Applied Mechanics and Engineering 79, 173-202.
• Zlotek-Zlotkiewicz, E., Monnier, S., Cappello, G., Le Berre, M., Piel, M., (2015). Optical volume and mass measurements show that mammalian cells swell during mitosis. The Journal of Cell Biology 211, 765-774.

APPENDIX A - SOLUTION OF THE LOCAL CONSTITUTIVE PROBLEM

Herein are present the operational details related to the solution strategy employed to solve the variational constitutive principle (3.8). Furthermore, the nonlinear equations resulting from the stationarity condition are solved by the full Newton-Raphson procedure.

A.1 Stationarity Condition

In view of the variational principle (3.8), the stationarity condition is defined through the variation dℒ = 0, which results in the following equations:

(A.1)

The partial derivatives of (A.1) result in,

(A.2)

and,

(A.3)

where M _{^en+ 1} := C_{^{e n+} 1}S_{^{e n+} 1} are the elastic Mandel stress tensor.

A.2 Newton-Raphson Procedure

Concerning to the full Newton-Raphson procedure, one can define from the nonlinear equations (A.1) the quantities,

(A.4)

where r _{n+ 1} is a residual and x_{n+ 1} represents the set of unknown variables. Therefore, an increment δ related to a current iteration of the Newton’s algorithm, is computed from the linear system of equation,

(A.5)

which can be represented by the compact arrangement J _{n+ 1} ⋆ δx_{n+ 1}= -r _{n+ 1}. The symbol ⋆ represents the appropriate product between the Newton’s increments δx_{n+ 1} and the Jacobian matrix

(A.6)

In order to achieve the quadratic convergence rate provided by the Newton’s method, the derivatives of equation (A.6) must be computed precisely. Accordingly, the second partial derivatives of (A.6) result in

(A.7)

where it is introduced the fourth order operators,

(A.8)

In equations (A.8), 𝕀_sym is the fourth order symmetric identity tensor. In addition, are introduced the tensor products (A ⊗ B)_ijkl: =(A)_ij (B)_kl and (AB)_ijkl: =(A)_ik (B)_jl , for any second orders tensors A and B.

APPENDIX B - CONSISTENT MATERIAL TANGENT MODULUS

Taking into account a total Lagrangian formulation (Belytschko et al., 2000), the linearization of the equilibrium equations results in the material tangent modulus

(B.1)

where one can note that if 𝒫^red _inc is a convex potential of C _{n +1}, the tensor ℂ_{n +1} will present major symmetry. The total derivative of the second Piola-Kirchhoff stress tensor is given by

(B.2)

The partial derivatives of (B.2) result in

(B.3)

where,

(B.4)

in which the fourth order operators ∂S^e _n+1/∂^e _{n +1} were introduced in equations (A.8). The total derivative dC _n+1 of equation (B.2) is one of the solutions of the following linear system of equation:/d and ℂ

(B.5)

resulting from the total derivative of the residual (A.4). The linear system (B.5) can be presented in a convenient compact form:

(B.6)

One can note that from equation (B.6), the derivative ∂r _n+1/∂ x_n+1 is the Jacobian matrix previously defined in equation (A.6). Finally, the partial derivatives of the set ∂r _n+1/∂C _n+1are computed as

(B.7)

where,

(B.8)

Publication Dates

• Publication in this collection
  June 2017

History

• Received
  07 Feb 2017
• Reviewed
  07 Mar 2017
• Accepted
  13 Mar 2017