Topological sensitivity analysis for a two-parameter Mooney-Rivlin hyperelastic constitutive model

The Topological Sensitivity Analysis (TSA) is represented by a scalar function, called Topological Derivative (TD), that gives for each point of the domain the sensitivity of a given cost function when an infinitesimal hole is created. Applications to the Laplace, Poisson, Helmoltz, Navier, Stokes and Navier-Stokes equations can be found in the literature. In the present work, an approximated TD expression applied to nonlinear hyperelasticity using the two parameter Mooney-Rivlin constitutive model is obtained by a numerical asymptotic analysis. The cost function is the total potential energy functional. The weak form of state equation is the constraint and the total Lagrangian formulation used. Numerical results of the presented approach are considered for hyperelastic plane problems.


INTRODUCTION
Topological optimization may be classified as continuous and discrete. The topological optimization of continuous structures are given in terms of material (micro) and geometrical (macro) approaches. The main technique of the material approach is the SIMP (Solid Isotropic Microstructures with Penalization) [6]. The methods ESO (Evolutionary Structural Optimization) and TSA (Topological Sensitivity Analysis) are techniques related to the geometrical approach [25,31,33,34]. For discrete structures, the optimal topology is determined by the optimal number, position and connectivity of the structural elements. A complete review of discrete topological optimization may be found in [6,14,17,26,[28][29][30].
The TSA technique calculates the sensitivity of the cost function when an infinitesimal hole is created in the domain of the problem [8,16,31,33]. The sensitivity is given by the topological derivative or gradient. In [25], the topological derivative was presented using shape sensitivity analysis and applied to the steady-state heat conduction and linear elasticity problems.
The topological derivative is obtained for the non-linear torsion problem of a primastic bar under stationary fluency described by the p-Poisson equation in [24]. Exact expressions for the topological derivative in the nonlinear Helmholtz and Navier-Stokes equations have been obtained in [1][2][3]. Topological derivative and level set methods were presented in [7]. Applications to time dependent problems and image processing are found in [5,22]. In [27], the TSA concept was applied to large deformation problems and an approximated expression to the topological derivative was obtained.
The purpose of this work is to extend the results obtained for large deformation with linear elastic material in [27] to nonlinear hyperelasticity using the 2 parameter Mooney-Rivlin constitutive model. In Sections 2 and 3, the topological derivative for the Lagrangian formulation and the nonlinear hyperelasticity model are reviewed. In Section 4, the expression of the topological derivative is obtained for nonlinear hyperelasticity and the total potential energy functional. In Section 5, an asymptotic numerical analysis is developed to determine the behavior of the analytical expression of the topological derivative. Finally, a heuristic topological optimization algorithm is applied to three hyperelastic two-dimensional structures and conclusions are addressed in Sections 6 and 7, respectively.

TSA FOR THE TOTAL LAGRANGIAN FORMULATION
The topological derivative [25] can be extended to large deformation problems using the total Lagrangian formulation as in [27]. A review of that application is presented below. It should be emphasized that this paper does not intend to present any formal proof of the topological derivative to non-linear elastic problems. The topological derivative in the cases of semilinear problem and linear elasticity is demonstrated in [19,23], based on the asymptotic expansions which results in the topological derivatives of general shape functionals. Asymptotic analysis is also presented for nonlinear Helmholtz and Navies-Stokes equations in [1][2][3].
Consider a new domain Ω 0 ε ∈ R n , Ω 0 ε = Ω 0 −B 0 ε , that has a boundary denoted by Γ 0 ε = Γ 0 ∪∂B 0 ε andB 0 ε = B 0 ε ∪∂B 0 ε is a ball of radius ε and centered at the material pointX ∈ Ω 0 ε , as illustrated in Figure 1. A total Lagrangian description for the topological derivative can be written as The action of increasing the hole can be interpreted as a sequence of perturbed configurations characterized by the parameter τ and described by a smooth and invertible function T (X, τ ) with X ∈ Ω 0 ε ⊂ R n and τ ∈ R. In this way, the sequence of domains Ω 0 τ and respective perturbed reference boundaries Γ 0 τ can be defined as In this way, the domain Ω 0 ε+δε , perturbed by a smooth expansion δε of the ball B 0 ε , and its respective boundary Γ 0 ε+δε can be written with relation to τ as Ω 0 The mapping between the non-perturbed reference domain Ω 0 ε and the perturbed domain Ω 0 τ can be written as where V n is the normal component to the hole of the velocity field V in the reference configuration. The perturbation δε is associated to the parameter τ in the following way (∀ X ∈ ∂B 0 Considering that the domain used in the shape variation sensitivity is the reference domain Ω 0 , then the topological derivative can be rewritten as where f (ε) is the regularizing function chosen in such a way that 0 < |D T (X)| < ∞.

NONLINEAR HYPERELASTICITY
In a large deformation problem, the final configuration of the body Ω ∈ R n can differ greatly from the initial or reference configuration Ω 0 ∈ R n . A body is deformed through a one to one mapping f that relates each material point X ∈ Ω 0 to a point x ∈ Ω in such a way that x = f (X) ∈ Ω and det ∇f > 0. The vector u (X) = f (X) − X represents the displacement of the material point X.
For a homogeneous incompressible hyperelastic material, the constitutive equation is given by [18] W where p ∈ R is the Lagrange multiplier related to the hydrostatic pressure and J = det F. Following [10], the strain energy density functional W for an isotropic hyperelastic material may be decomposed as the sum of the distortional (W ) and the volumetricW strain energy densities as and I 3 = det C is the third invariant of the right Cauchy-Green deformation tensor C. In addition, the following equation relates the hydrostatic pressure p and the volumetric strain energy densityW [10] In this work, the distortional and volumetric densities are written in terms of the Mooney-Rivlin form with two-parameters, respectively, as [9][10][11][12]21] where A 10 and A 01 are material constants andkis the bulk modulus, which is a numerical penalization term for nearly incompressible materials. After replacing equations (8) and (9) in (6), the strain energy density functional for an isotropic nearly incompressible hyperelastic material is given by The mixed variational form of the non-linear nearly incompressible hyperelastic problem is expressed by: find u ∈ U and p ∈ Q such that [32] { a (u, δu) with where It is assumed that the reference domain Ω 0 ∈ R n is open and bounded. Its boundary is sufficiently regular and admits the existence of a normal unitary vector n in almost all of the points of Γ 0 , except in a finite set of zero measure.

TSA IN NONLINEAR HYPERELASTIC PROBLEMS
In this section, the expression of the topological derivative for the isotropic nearly incompressible hyperelastic problem is obtained from equation (1). This expression requires the evaluation of a given cost function Ψ, defined on the non-deformed configuration Ω 0 ε , when the hole B 0 ε centered in X ∈ Ω 0 ε and ε → 0 increases in size according to the velocity field defined. The topological derivative is related to the shape sensitivity analysis as indicated in equation (4).
The mixed variational statement for the elastic problem with large deformation and homogeneous, isotropic and nearly incompressible hyperelastic material in the reference nonperturbed domain Ω 0 ε , assuming that the domain Ω 0 ε is limited and open and Γ 0 ε sufficiently regular, is given by: where U ε , V ε and Q ε are, respectively, the spaces of the admissible kinematically functions, their variations and pressures defined on the reference domain with the non-perturbed hole B ε . An equivalent mixed variational problem to the system given in (17) in the non-perturbed reference domain Ω 0 ε may be defined by the family of perturbed reference domains Ω 0 τ , remembering that Ω 0 τ | τ =0 = Ω 0 ε , as: find u τ ∈ U τ and p τ ∈ Q τ such that where U τ , V τ and Q τ are the spaces of kinematically admissible functions, their variations and pressures defined in Ω 0 τ , respectively. Considering the strain energy functional per unit of non-deformed volume given in (10), the total potential energy functional Ψ τ (u τ ) is written in Ω 0 τ as To obtain the sensitivity of the cost function (19), the lagrangian function for the hyperelastic problem, in the non-perturbed configuration Ω 0 τ , is written as where β τ = m 1 δu τ and θ τ = m 2 δp τ (m 1 and m 2 ∈ R) are, respectively, the Lagrange multipliers related to the first and second equations of the mixed variational system (18). However, considering that equation (18) is satisfied for all τ , the derivative of the lagrangian function (20) will be the same as the derivative of the total potential energy functional (21) are given, respectively, by The partial derivative of the lagrangian function (20) in τ is the same as the total derivative if the directional derivatives presented in (21) are zero. Therefore, from (22) to (30), the following variational problem is obtained: find u τ ∈ U τ and p τ ∈ Q τ such that which represents the state equation for the considered problem. The adjoint problem is: find β τ ∈ V τ and θ τ ∈ Q τ such that where the symmetry of the bilinear forms δa τ (s τ ;u τ , β τ ) and δg τ (p τ ;ṗ τ , θ τ ) has been taken into account. The directional derivative of the total potential energy functional in u τ in the directionu τ , according to (19), may be written as where it has been assumed that u τ and p τ satisfy the state equation (18). Consequently, the solution of the adjoint equation (32) is (β τ , θ τ ) = (0, 0) and the partial derivative of the lagrangian (20), written in the perturbed configuration Ω 0 τ , is Replacing the definition of the total potential energy for the non-linear nearly incompressible hyperelastic problem, given in (19), in (34) and using the Reynolds Transport Theorem, the partial derivative of the lagrangian function in the non-perturbed configuration Ω 0 where for u τ and β τ fixed and the velocity field V defined by Therefore, equation (35) is rewritten as where Using the definition of the partial derivative of the Green-Lagrange strain tensor E τ in the non-perturbed configuration Ω 0 ε , given in (39), equation (38), after a few manipulations, becomes Equation (40) may be expressed in the following form where Σ 0 ε is the Eshelby's energy moment tensor in the non-perturbed reference configuration Ω 0 ε for the total lagrangian formulation [15]. For the considered problem, it is defined by or in terms of the first Piola-Kirchhoff tensor From the divergence theorem and the tensorial expression equation (41) is rewritten as It is possible to show that DivΣ 0 ε = 0.
Replacing equation (46) in (45) and considering the velocity field (37), the following equation is obtained From the definition of the Eshelby's tensor Σ 0 ε given in (43), the following relation is valid Considering a homogeneous Neumann boundary condition in the hole, i.e., P ε n 0 = 0 on ∂B 0 ε , and no body forces, b 0 = 0, the partial derivative of the lagrangian function in The topological derivative, except for the limit ε → 0, is obtained after the substitution of (48) in the expression of the topological derivative for the total lagrangian formulation given in (4). Therefore,

NUMERICAL ASYMPTOTIC ANALYSIS
Equation (49) represents the topological derivative, except for the limit with ε → 0, for the case of large deformations and nonlinear nearly incompressible hyperelasticity. In order to obtain the topological derivative, it is necessary that the limit for ε → 0 in (49) be calculated either analytical or approximately. For the present non-linear problem, an analytic asymptotic analysis of equation (49) becomes impracticable. A procedure based on numerical experiments for the calculation of the limit with ε → 0 in (49) was used in [27] for large deformation and linear elastic material. The same asymptotic analysis is applied below to large deformation and Mooney-Rivlin material.
The topological sensitivity analysis aims to provide an asymptotic expansion of a shape functional with respect to the size of a small hole created inside the domain [3]. The asymptotic analysis makes possible to obtain the behavior of the integrand of the topological derivative expression (49) in the limit when the hole radius becomes small. Another point is to compare the quotient between the topological derivative in the domain with a small hole and the strain energy functional in the domain without the hole.
Consider the function d T (u ε ) defined by in such a way that A numerical study of the asymptotic behavior of the function d T (u ε ) with relation to the radius ε is developed. Consider a plane square domain, denoted by Ω, with size L = 2 mm and a hole of radius ε at the center of the domain, subjected to the distributed load cases t 0 along the edges of Ω, as illustrated in Figure 2. It is assumed plane strain and Mooney-Rivlin model with A 10 = 0.55 N /mm 2 , A 01 = 0.138 N /mm 2 andk = 666.66 N /mm 2 . Figure 3 shows the meshes for the original domain and the domain with the hole. The purpose of the first two model cases is to verify the influence of each of the displacement components u and v separately on the results of the asymptotic analysis of the topological derivative. The third and fourth load cases aim to study the influence of the two displacements components acting together. For all cases, traction and compression loads were applied. Other load conditions were considered and the conclusions of the asymptotic analysis were the same. These load cases constitute the minimum set to obtain the main results from the asymptotic analysis.   The finite element meshes of quadratic triangles are constructed in such a way that they have the same number of elements at the boundary of the hole, independently of the value of the radius ε. Consequently, the approximated size of the elements is calculated as being n e the number of elements required over the boundary of the hole. Table 1 shows the total number of elements N E of the meshes generated in the domain Ω for different values of the radius ε and n e = 60.   The plots of the asymptotic behavior of d T (u ε ) f ′ (ε) × ε are shown in Figures 4 and 5, respectively, for the traction loads t 0 = {6.25; 12.50; 25.00; 37.50} × 10 −3 N /mm 2 and compression loads t 0 = − {6, 25; 12, 50; 25, 00} × 10 −3 N /mm 2 loads. Based on these results, it is reasonable to assume that the integrand term in equation (50) behaves as a straight line trough the origin in relation to ε. A function f (ϵ) which satisfies the condition 0 (50) and (51) give  The analysis of the constant C is based on the asymptotic behavior of the function d T (u ε ) in relation to the radius ε. For that purpose, consider the quotient between the function d T (u ε ), given in (50), and the integrand calculated on the node of the mesh without hole, whose coordinates coincide with the coordinates of the center of the hole in the mesh with holes (see Figure 3). Therefore, for f (ε) = −πε 2 , the following quotient is obtained The plots in Figure 6 show the behavior of this quotient versus 1/ε for each value of t 0 and the four load cases. The plots for the first, second and fourth load cases miss a certain smoothness and have an increase in the quotient for ε = 0.01 ( 1 ε = 100). This behavior may be related to the distortion of the elements around the hole which is amplified as the values of t 0 increase. In the compressive case, this point was much more critical and made impossible to run the simulations for the case with t 0 ≤ −25.00 × 10 −3 N /mm 2 . However, it seems reasonable to assume an asymptotic behavior of equation (54), in function of ε, for all traction and compression load cases.  For the third load case, the values for C are much larger than for the other load cases. The C values depend on the combination of loads and material properties. For the third load case, C is larger because the strain energy density calculated on the central node of the mesh without hole (denominator of equation (54)) is very small. The distortional strain energy is almost zero. Due to the near-incompressibility behavior of the Mooney-Rivlin material, the volumetric part is also very small. Figure 7 shows the results for the third load case with k = 1.0 N /mm 2 which makes the material more compressible. The values for C are now much smaller and close to the values obtained for the linear elastic material given in [27]. However, the main point of the asymptotic analysis is that the variation range of the C values is more limited than the strain energy densities W for a given load and material properties.    strain energy density is much larger when compared to the variation of C for all the considered cases.
Based on the previous results, the expression of the topological derivative for the non-linear nearly incompressible hyperelastic problem may be written approximately as where C * is taken from the previous values obtained for C. The effective value of the constant C * becomes of little practical interest, because the topological derivative will be calculated for all the nodes of the mesh and the holes will be created where the topological derivative assumes the least values [27].

RESULTS
In this section, the heuristic algorithm used in [27] is applied to obtain the topology of three two-dimensional non-linear hyperelastic large deformation problems.

Clamped short beam
In this example, the initial domain is a square modeled as plane strain problem with sides L = 50 mm, Mooney-Rivlin constants A 10 = 0.55 N /mm 2 , A 01 = 0.138 N /mm 2 andk = 666.66 N /mm 2 . The beam is clamped on the regions indicated by a = 5 mm and subjected to the distributed load t 0 = −0.4444 N /mm 2 on the region indicated by b = 4.5 mm, as illustrated in Figure 8(a). A mesh of 3656 triangular quadratic finite elements showed in Figure 8(b) was used. The stopping criterion was based on the final areaĀ = 0.40A 0 , with A 0 the initial domain area, and 1% of the area was removed in each iteration. The final result was reached after 89 iterations of the algorithm presented in the last section and is showed in Figure 8(c). The maximum deflexion was −7.28 mm in the upper right side of the beam. Figure 8(d) shows the decreasing of the total potential energy functional versus the number of iterations. Figure 9 illustrates the final distribution of the strain energy density W .

Two bars
In this example, the initial domain is subjected to a concentrated load F = −1.0 N as illustrated in Figure 10(a). The other parameters are the same used in the previous example. The required final area isV = 0.36V 0 . The final topology, illustrated in Figure 10(c), is achieved after 99 iterations and the maximum deflexion is −3.75 mm on the point of load application. Figure  10(d) illustrates the behavior of the total potential energy functional along the iterations of the algorithm. Figure 11 shows the distribution of the strain energy density W on the domain after 99 iterations.  Iterations Ψ(u) Figure 11 Total potential energy versus iteration number.

Clamped-clamped short beam
In this third example, the clamped-clamped beam given in Figure 12 is considered for L = 20 mm and a concentrated force F = −0.04 N applied on the beam middle span. The same previous material constants are used. The final considered area isV = 0.20V 0 and the percentage of removed area is 3% for each iteration. The final topology is obtained after 52 iterations and illustrated in Figure 13(b). A similar result was obtained in [20,27] for large deformation and linear elastic material.

CONCLUSIONS
This paper presented the application of the TSA to large deformation and nearly incompressible hyperlastic Mooney-Rivlin material using the total Lagrangian formulation.
As in [27], an approximated expression for the topological derivative was obtained using a numerical asymptotic analysis for four different load conditions. The same behaviour was observed from the asymptotic analysis of the topological gradient for linear elastic and hyperelstic

materials.
A heuristic topological optimization algorithm was applied to three planar problems with deformation of about 20%. In the case of the third example, similar topologies were obtained with those ones calculated using large deformation and linear material [27].
Most of the computational cost is related to the solution of the state equation due to the nearly incompressible material model and the possibility of the locking phenomenon. Once the solution has converged, the topological derivative is given by the strain energy functional as indicated in expression (55).
The final solution in Figure 13(b) shows the formation of mechanisms. This is related to the way the optimization algorithm is implemented based on the removal of a fixed amount of material. In this sense, the implemented TSA procedure has the same limitations of most hard-kill methods.
Originally, the TSA approach did not define a classical minimization problem. In [4], the topological derivative is used as a descent direction in the optimization problem. Also a penalty method is used for the constraint imposition in terms of the von Mises stress.
The procedure presented in this paper allows only the removal of elements. The topology optimization of non-linear problems may also require the inclusion of material which is presented in [13].