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

a02v7n4

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

Carlos E. L. Pereira^*; Marco L. Bittencourt

Department of Mechanical Design, Faculty of Mechanical Engineering, University of Campinas, P.O. Box 6122, Campinas, SP, 13083-970 – Brazil

ABSTRACT

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.

Keywords: Topological Optimization, Finite Elements, Nonlinear Elasticity, Hyperelasticity, Topological Derivative.

1 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–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–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.

2 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–3].

Consider a new domain , that has a boundary denoted by is a ball of radius e and centered at the material point ∈ , 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,t) with ∈ and τ ∈ ℜ. In this way, the sequence of domains and respective perturbed reference boundaries can be defined as

and where where and . In this way, the domain , perturbed by a smooth expansion δε de of the ball , and its respective boundary can be written with relation to τ as .

The mapping between the non-perturbed reference domain and the perturbed domain 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 ∈ and ∀ X_τ ∈ )

Considering that the domain used in the shape variation sensitivity is the reference domain Ω⁰, then the topological derivative can be rewritten as

where f(ε) is the regularizing function chosen in such a way that 0 < |D_T ()| < ∞.

3 Nonlinear Hyperelasticity

In a large deformation problem, the final configuration of the body Ω ∈ ℜⁿ can differ greatly from the initial or reference configuration Ω⁰∈ ℜⁿ. A body is deformed through a one to one mapping f that relates each material point X ∈ Ω⁰ 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]

where p ∈ ℜ 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 () and the volumetric strain energy densities as

where = 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 density [10]

In this work, the distortional and volumetric densities are written in terms of the Mooney-Rivlin form with two-parameters, respectively, as [9–12,21]

where A₁₀ and A₀₁ are material constants and is 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]

with

where with dim ≤ 3 and Q = { p ∈ L₂( Ω⁰) } . It is assumed that the reference domain Ω⁰∈ ℜⁿ 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 Γ⁰, except in a finite set of zero measure.

4 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 , when the hole centered in X ∈ 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 non-perturbed domain , assuming that the domain is limited and open and sufficiently regular, is given by: find u_ε ∈ and p_ε ∈ such that

where , and 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 may be defined by the family of perturbed reference domains , remembering that , as: find u_τ ∈ such that

where are the spaces of kinematically admissible functions, their variations and pressures defined in , respectively.

Considering the strain energy functional per unit of non-deformed volume given in (10), the total potential energy functional Ψ_τ( u_τ) is written in as

To obtain the sensitivity of the cost function (19), the lagrangian function for the hyperelastic problem, in the non-perturbed configuration , is written as

where β_τ = m₁δu_τ and θ_τ = m₂δp_τ (m₁ and m₂∈ ℜ) 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 t, the derivative of the lagrangian function (20) will be the same as the derivative of the total potential energy functional

where

The directional derivatives in (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_τ ∈ such that

which represents the state equation for the considered problem. The adjoint problem is: find β_t∈ V_t and θ_t∈ Q_t such that

where the symmetry of the bilinear forms has been taken into account.

The directional derivative of the total potential energy functional in u_τ in the direction _τ, 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 , 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 is given by

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 , given in (39), equation (38), after a few manipulations, becomes

Equation (40) may be expressed in the following form

where is the Eshelby's energy moment tensor in the non-perturbed reference configuration 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

Replacing equation (46) in (45) and considering the velocity field (37), the following equation is obtained

From the definition of the Eshelby's tensor given in (43), the following relation is valid

Considering a homogeneous Neumann boundary condition in the hole, i.e., P_εn₀ = 0 on , and no body forces, b₀ = 0, the partial derivative of the lagrangian function in τ = 0 becomes

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,

5 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 e 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₀ along the edges of Ω, as illustrated in Figure 2. It is assumed plane strain and Mooney-Rivlin model with A₁₀ = 0.55 N/mm², A₀₁ = 0.138 N/mm² and = 666.66 N/mm². 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_ε 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_ε = 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₀ = { 6.25;12.50;25.00;37.50} ×10^-3 N/mm² and compression loads t₀ = -{ 6,25;12,50;25,00}×10^-3 N/mm² 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 e. A function f(e) which satisfies the condition 0 < | D_T()| < ∞ is f(ε) = -| B_ε| = -πε². Consequently, equations (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 e. 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(ε) = -πε², the following quotient is obtained

The plots in Figure 6 show the behavior of this quotient versus 1/ε for each value of t₀ 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 ( = 100). This behavior may be related to the distortion of the elements around the hole which is amplified as the values of t₀ increase. In the compressive case, this point was much more critical and made impossible to run the simulations for the case with t₀ ≤ -25.00×10^-3 N/mm². 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 = 1.0 N/mm² 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.

Tables 2 to ⁵ compare the values of C and W. It may be observed that the variation of the 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].

6 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.

6.1 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₁₀ = 0.55 N/mm², A₀₁ = 0.138 N/mm² and = 666.66 N/mm². The beam is clamped on the regions indicated by a = 5 mm and subjected to the distributed load t₀ = -0.4444 N/mm² 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₀, with A₀ 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.

6.2 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 is = 0.36V₀. 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.

6.3 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 is = 0.20V₀ 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.

7 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].

References

[1] S. Amstutz. The topological asymptotic for the Helmholtz equation. SIAM Journal on Control and Optimization, 42(5):1523–1544, 2003.

[2] S. Amstutz. The topological asymptotic for the Navier-Stokes equations. ESAIM: Control, Optimisation and Calculus of Variations, 11(3):401–425, 2005.

[3] S. Amstutz. Topological sensitivity analysis for some nonlinear PDE systems. Journal de Mathématiques Pures et Appliquées, 85(4):540–557, 2006.

[4] S. Amstutz and A.A. Novotny. Topological optimization of structures subject to von mises stress constraints. Struc- tural and Multidisciplinary Optimization Journal, 41(3):407–420, 2010.

[5] S. Amstutz, T. Takahashi, and B. Vexler. Topological sensitivity analysis for time-dependent problems. ESAIM: Control, Optimisation and Calculus of Variations, 14(3):427–455, 2008.

[6] M. P. Bendsϕe. Optimization of Structural Topology, Shape and Material. Springer, Heidelberg, 1995.

[7] M. Burger, B. Hackl, and W. Ring. Incorporating topological derivatives into level set methods. Journal of Compu- tational Physics, 194(1):344–362, 2004.

[8] J. C´ea, Ph. Guillaume, and M. Masmoudi. The shape and topological optimizations connection. Technical Report, UFR MIG, Universit´e Paul Sabatier, Toulouse, French, 1998.

[9] J.S. Chen, W. Han, C.T. Wu, and W. Duan. On the perturbed lagrangian formulation for nearly incompressible and incompressible hyperelasticity. Computer Methods in Applied Mechanics and Engineering, 142:335–351, 1997.

[10] J.S. Chen and C. Pan. A pressure projection method for nearly incompressible rubber hyperelasticity, Part I: Theory. Journal of Applied Mechanics, 63:862–868, 1996.

[11] J.S. Chen, C.T. Wu, and C. Pan. A pressure projection method for nearly incompressible rubber hyperelasticity, part II: Applications. Journal of Applied Mechanics, 63:869–876, 1996.

[12] J.S. Chen, S. Yoon, H.P. Wang, and W. K. Liu. An improved reproducing kernel particle method for nearly incompressible finite elasticity. Computer Methods in Applied Mechanics and Engineering, 181:117–145, 2000.

[13] J. Rocha de Faria, A.A. Novotny, R.A. Feijo, and E. Taroco. First and second order topological sensitivity analysis for inclusions. Inverse Problems in Science and Engineering Journal, 17(5):665–679, 2009.

[14] H. Eschenauer, N. Olhoff, and (eds). Optimization methods in structural design. In EUROMECH-Colloquium, volume 164, Wien, 1982. B.I.-Wissenschaftsverlag.

[15] J. D. Eshelby. The elastic energy-momentum tensor. Journal of Elasticity, 5(3-4):321–335, April 1975.

[16] S. Garreau, Ph. Guillaume, and M. Masmoudi. The topological gradient. Technical Report, Universit´e Paul Sabatier, Toulouse, French, 1998.

[17] W. Gutkowski and Z. Mroz. Second world congress of structural and multidisciplinary optimization. In WCSMO-2, Intitute of Fundamental Technological Research, volume 182, Warsaw, Poland, 1997.

[18] G. A. Holzapfel. Nonlinear Solid Mechanics: A Continuum Approach for Engineering. John Wiley & Sons, 2000.

[19] M. Iguernane, S.A. Nazarov, J.-R. Roche, J. Sokolowski, and K. Szulc. Topological derivatives for semilinear elliptic equations. Appl. Math. Comput. Sci., 19(2):191–205, 2009.

[20] D. Jung and C. Gea. Topology optimization of nonlinear structures. Finite Elements in Analysis and Design, 40:1417 – 1427, 2004.

[21] N. H. Kim. Shape Design Sensitivity Analysis and Optimization of Nonlinear Static/Dynamic Structures with Contact/Impact. PhD thesis, The University of Iowa, May 1999.

[22] I. Larrabide, R.A. Feij´oo, A.A. Novotny, E. Taroco, and M. Masmoudi. Topological derivative: A tool for image processing. Computers & Structures, 86:1386–1403, 2005.

[23] S.A. Nazarov and J. Soko lowski. Asymptotic analysis of shape functionals. Journal de Mathématiques Pures et Appliquées, 82(2):125–196, 2003.

[24] A. A. Novotny, R. A. Feij´oo, E. Taroco, M. Masmoudi, and C. Padra. Topological sensitivity analysis for a nonlinear case: The p-Poisson problem. In 6th World Congress on Structural and Multidisciplinary Optimization, Rio de Janeiro, Brazil, 2005.

[25] A. A. Novotny, R. A. Feij´oo, E. Taroco, and C. Padra. Topological sensitivity analysis. Computer Methods in Applied Mechanics and Engineering, 192:803–529, 2003.

[26] N. Olhoff and G. I. N. Rozvany. Structural and multidisciplinary optimization. Proc of WCSMO-1, Oxford/UK, 1995.

[27] C. E. L. Pereira and M. L. Bittencourt. Topological sensitivity analysis in large deformation problems. Stuctural Multidisciplinary Optimization, 37(2):149–163, 2008.

[28] G. I. N. Rozvany, editor. Topology Optimization in Structural Mechanics, Vienna/Austria, 1997. Springer. Vol 374 of CISM Course and Lectures.

[29] G. I. N. Rozvany. A critical review of established methods of structural topology optimization. Structural and Multidisciplinary Optimization, 37(3):217–237, 2009.

[30] G. I. N. Rozvany, M. P. Bendsϕe, and U. Kirsch. Layout optimization of structures. Applied Mechanics Reviews, 48:41–119, 1995.

[31] A. Schumacher. Topologieoptimierung Von Bauteilstrukturen unter Verwendung Von Lopchpositionierungkrieterien. PhD thesis, Universit¨at-Gesamthochschule-Siegen, 1995.

[32] C.A.C. Silva and M.L. Bittencourt. Structural shape optimization of 3d nearly-incompressible hyperelasticity problems. Latin American Journal of Solids and Structures, 5(2):129–156, 2008.

[33] J. Sokolowski and A. Z˙ ochowski. On topological derivative in shape optimization. Technical Report, INRIA-Lorraine, French, 1997.

[34] C. Zhao, G. P. Steven, and Y. M. Xie. Evolutionary optimization of maximizing the difference between two natural frequencies. Structural Optimization, 13(2-3):148–154, 1997.

Received 24 May 2010;

In revised form 19 Aug 2010

Publication Dates

• Publication in this collection
  06 Mar 2012
• Date of issue
  2010

History

• Received
  24 May 2010
• Accepted
  19 Aug 2010