1. Introduction

The mechanical behavior for solids is significantly affected by irreversible evolution of multiple distributed defects in the immediate neighborhood of a material point. In the past decades, a variety of homogenization methods have been presented for developing equivalent homogeneous material models that effectively represent the mean mechanical constitutive response of the defected solids. The appropriateness of a homogenization method indeed relies on the purpose of simulation. In order to derive a sufficiently accurate model, it is in general necessary that the local defect interactions be properly captured. Also, in order for easy implementation, one of the major factors in concern is the computational effort of the respective homogenization procedure.

The concept of using micromechanics constitutive models has been extensively employed as an efficient homogenization method in linear elasticity. A large number of analytical and numerical schemes of micromechanical modeling have been presented for problems containing either periodically distributed cracks (e.g., Nemat-Nasser et al. ^[1]), or irregularly distributed cracks (e.g., Kachanov ^[2], Petrova et al. ^[3], Shen and Li ^[4], etc.), and inclusions (e.g., Jasiuk ^[5], Nozaki and Taya ^[6], Tsukrov and Novak ^[7], etc.). In many of these models, the complex geometry of heterogeneous materials is simplified by considering the multi-defected solid as one inhomogeneity embedded in a homogeneous infinite matrix medium. By imposing a set of uniform loads on the remote boundaries of the matrix, the corresponding effective constitutive model is then derived by integrating and averaging the governing kinematic, stress, and energetic quantities over the respective material neighborhood. In this sense, the macroscopic effective properties lead to equivalent mechanical behavior as that of the material with the defected microstructure. Nevertheless, solution to mutual interaction of the defects in most of these models typically costs significant computational effort, especially when the number of defects becomes relatively large. A scalar energy parameter was presented (Chang and Liu ^[8]) to efficiently characterize the interactive behavior. With this energy parameter, the effective linear elastic moduli can thus be more easily obtained by direct use of numerical schemes such as finite element method.

For problems involving more complicated geometry and mechanics, the multi-scale modeling methods are more appropriate for alternatively estimating the effective mechanical properties. A variety of multi-scale schemes have been presented for this purpose. One of the most commonly-used approaches is the computational homogenization technique (e.g., Miehe et al. ^[9], Kouznetsova et al. ^[10], etc.). This technique has also been extensively exploited in many mechanical and engineering aspects including, e.g., masonry works (Mistler et al. ^[11]), thermomechanical applications (Shabana and Noda ^[12]), material layer models (Matous et al. ^[13], Hirschberger et al. ^[14]), and dynamic analysis (Pham et al. ^[15]), etc. The applications for the multi-scale methods also include the following remarkable works, e.g., Belytschko and Xiao ^[16], Liu et al. ^[17], Budarapu et al. ^[18], Budarapu et al. ^[19], Talebi et al. ^[20], Yang et al. ^[21], etc. Through the concept of scale transitions, the numerical solutions to the coupled multi-scale boundary value problems can be accurately provided. Thereby, large deformations, nonlinearity, plasticity, damage, fracture, etc., can be properly taken into account in the modeling process. Nevertheless, the nested iterative algorithm in the calculation iscomputationally intensive and sometimes inconvenient to be implemented in practice.

In this paper, a micromechanics constitutive model is proposed to characterize the effective elastic behavior of multi-defected solids subjected to large deformation. The model is formulated in terms of anisotropic hyperelasticity and allows for a straightforward incorporation of material and geometrical nonlinearities. By employing an energy parameter that evaluates the released energy due to presence of the defects, a correspondence relation between the effective hyperelastic strain energy density function and this energy parameter is established in the context of large deformation. With combined use of the "direct difference approach" and the extended "modified compliance contribution tensor" method, the constitutive model and the associated material coefficients can be properly constructed by direct use of solutions from numerical schemes such as finite element method. The presented constitutive model can be easily implemented in macroscopic structural analysis in practice.

---

2. Linear Effective Elastic Moduli

In this section, two commonly-used homogenization procedures for determining linear effective elastic moduli, both based on the concept of micromechanics constitutive modeling, are briefly reviewed. Then, a modified method is proposed in the third subsection. All the three procedures are formulated in conjunction with a scalar energy parameter. Each of them thus represents a correspondence relation between the effective moduli of a multi-defected solid and the energy parameter. Finally in the fourth subsection, the numerical approaches for calculation of this energy parameter are illustrated.

---

2.1. Eshelby's Equivalent Inclusion Method

We consider a homogeneous infinite matrix medium, modeled by linear elasticity, containing a set of irregularly distributed defects and subjected to uniform remote loads σ^∞ = (σ₁^∞, σ₂^∞, τ^∞) (Figure 1).

In order to describe the damage level due to presence of the set of defects, we delimit a cutoff area Ω, as shown in the shaded region in Figure 1. The size of Ω is chosen to be large enough to enclose the whole set of defects. Also, the shape of Ω can be specified either according to the statistical distribution of defects or by using any other convenient selections, e.g., the rectangle in Figure 1.

In order to determine the effective elastic moduli of Ω, we take the whole cutoff area as a single inclusion and postulate that the effective mechanical properties of Ω be equivalent to those of this inclusion. An energy parameter Δ∏, which represents the released strain energy due to presence of the inclusion in the infinite matrix, is employed here to characterize the homogenized fracture state of Ω. Based on the concept of the Eshelby's equivalent inclusion theory (Eshelby ^[22]), the correspondence relation between the linear effective elastic moduli of Ω and Δ∏can then be expressed as

${{\sigma }^{\infty }}:{{[{{({{\mathbf{H}}_{eff}}-{{\mathbf{H}}_{o}})}^{-1}}:{{\mathbf{H}}_{o}}+\mathbf{S}]}^{-1}}:{{\mathbf{H}}_{o}}^{-1}:{{\sigma }^{\infty }}=-2\Delta \Pi /{{A}_{\Omega }}$

(1)

where H_eff is defined as the effective elastic moduli tensor of Ω(i.e., the inclusion), H_o is the elastic stiffness tensor of the matrix material, S is the Eshelby tensor, and A_Ω is the area of Ω. Details for derivation of the left hand side of Eqn. (1) were presented by, e.g., Nemat-Nasser et al. ^[1]. Also, expressions for S associated with various shapes of inclusion can be found in, e.g., Walpole ^[23].

To find each component of the anisotropic lineareffective elastic tensor H_eff various combinations for the remote loads σ^∞ need to be considered. By evaluating the released strain energy Δ∏ associated with these loading conditions, the effective material moduli can then be determined.

---

2.2. Compliance Contribution Tensor Method

We again consider a homogeneous infinite linearly elastic matrix medium subjected to uniform remote loads σ^∞ and take the whole cutoff area Ω (Figure 1) as a single inclusion. Based on the concept of thecompliance contribution tensor (CCT) C_inc (e.g., Kachanov et al. ^[24], Tsukrov and Novak ^[7], etc.), the correspondence relation between the linear effective elastic moduli of the inclusion Ω and the released strain energy Δ∏ can alternatively be written as

${{\sigma }^{\infty }}:{{\mathbf{C}}_{inc}}:{{\sigma }^{\infty }}\approx 2\Delta \Pi /{{A}_{\Omega }}$

(2)

where

${{\mathbf{C}}_{inc}}\equiv {{\mathbf{H}}_{eff}}^{-1}-{{\mathbf{H}}_{o}}^{-1}$

(3)

Still, by applying various conditions for the remote loads σ^∞ and evaluating the respective released strain energy Δ∏, the corresponding component of the tensor H_eff can be properly determined through combined use of Eqns. (2) and (3).

---

2.3. Modified Compliance Contribution Tensor Method

Analytically, Eqn. (1), the correspondence relation based on Eshelby's theory, is applicable for use in solving H_eff under all combinations of inclusion and matrix. On the other hand, it is observed that the CCT actually represents the increment of the compliance tensor. This means that the relation between C_inc and Δ∏ shown in Eqn. (2) is only an approximation. Thus, equivalence of both sides of this equation holds asymptotically as the value of Δ∏ becomes infinitesimally small, which indicates that the CCT method is feasible only for the conditions when the inclusion and the matrix possess very similar material behavior. Here, by combing Eqns. (2) and (3), and partly replacing the remote loads σ^∞ with σ^itf, a modified compliance contribution tensor (MCCT) can be written as

${{\sigma }^{itf}}:{{\mathbf{H}}_{eff}}^{-1}:{{\sigma }^{itf}}-{{\sigma }^{\infty }}:{{\mathbf{H}}_{0}}^{-1}:{{\sigma }^{\infty }}\approx 2\Delta \Pi /{{A}_{\Omega }}$

(4)

where σ^itf is the average of the stress tensors along the interface betweden the inclusion and the matrix. In practice, σ^itf can be properly extracted from the numerical fields resulting from, say, FE analysis.

In the first term on the LHS of Eqn. (4), the remote load σ^∞ that originally contributes to H_eff in Eqn. (2) has been replaced by σ^itf. With such modification, the effect resulting from the internal strain energy due to presence of the inclusion is more appropriately taken into account. The feasibility of MCCT will be illustrated in the following numerical examples.

---

2.4. Evaluation of Δ∏

For accuracy, it is required that the released strain energy Δ∏ be computed with a high degree of precision. A path-independent contour integral termed the "M-integral" was proposed by Chang and Liu ^[8] for this purpose. For linear elasticity, the M-integral evaluates twice the value of Δ∏. Due to path-independence, the integration contour can be arbitrarily chosen as long as they contain the whole set of defects.

Alternatively, Δ∏can be evaluated by using a "direct difference approach" (DDA). In this approach, we consider the defected specimen in Figure 1, yet taking the state containing no defect as an original configuration. The values of the strain energy associated with both the defected and original configurations under the same loading condition are calculated by using finite elements. The value of Δ∏ can then be obtained by directly evaluating the corresponding energy difference between these two configurations.

It was numerically investigated by Chang and Liu ^[8] that, for linear elasticity, Δ∏ can be properly evaluated by either the M-integral or the DDA as long as the cutoff area Ω being delimited in an extended region that is large enough to be regarded as an infinite matrix medium. Also, it was observed that the calculations are rather insensitive to the local near-defect FE discretization so that a complicated FE model around the defects is not required. Further, by comparing the computational aspects of these two approaches, DDA appears to be more straightforward in practice in that Δ∏can be easily computed by direct use of any well-developed FE code with almost no need of extra coding cost in the post-processing stage. Such an advantage indicates the superiority of DDA particularly in its extended application to problems subjected to large deformation. Therefore, M-integral is not used in this work. Instead, Δ∏ is evaluated by using DDA, which appears to be more computationally efficient.

---

3. Effective Constitutive Model (Large Deformation)

We again consider a homogeneous infinite elastic matrix medium containing a cutoff area Ω at its undeformed state, as shown in Figure 1. By applying a system of uniform remote loads that lead to large deformation, the body then reaches its current deformed state (this state is not shown in Figure 1). In the context of large elastic deformation, the highly nonlinear mechanical behavior can be properly characterized by using hyperelastic strain energy density function, with all state variables reinterpreted with respect to its undeformed state. In this sense, it is more appropriate to employ the concept of "effective strain energy density", rather than the effective elastic moduli, for use as a micromechanics model for describing the constitutive behavior of Ω.

Further, by examining the Eshelby's and the MCCT methods, it is observed that MCCT be more suitable for extension to nonlinear analysis. To construct the correspondence relation between the effective strain energy density and Δ∏, we have Eqn. (4) extended as

${{W}_{eff}}|{{\sigma }^{itf}}-\text{ }{{W}_{0}}|{{\sigma }^{\infty }}=\Delta \Pi /{{A}_{W}}$

(5)

where W_eff and W₀ are the effective strain energy densities of the inclusion Ω and the matrix, respectively. Note that, while W₀ is evaluated at the stress level σ^∞, W_eff is taken at σ^itf so that the contribution from the internal strain energy of the inclusion is taken into account.

By applying various combinations of remote boundary conditions (either external loads or stretches) and evaluating the associated released strain energy Δ∏ with DDA, the effective strain energy density function W_eff can then be properly determined.

---

4. Anisotropic Hyperelastic Material Model

For those materials that remain nonlinearly elastic under very large deformation, their mechanical behavior can be characterized by using hyperelastic models. The associated strain energy density function Ω can usually be expressed as a scalar function of the three principal invariants of the Cauchy-Green strain tensor V² (i.e., I₁, I₂ and I₃) as

$W\text{ }=\text{ }W\left( {{I}_{1}},\text{ }{{I}_{2}},\text{ }{{I}_{3}} \right)$

(6)

In the present study, the hyperelasitc behavior of the infinite matrix medium is modeled with a compressible isotropic NeoHookean strain energy density function as

$W\text{ =}=\Lambda (\frac{{{I}_{1}}}{I_{3}^{1/3}}-3)\text{ }+\frac{1}{2}K{{({{I}_{3}}^{1/2}-1)}^{2}}$

(7)

where ∧ and K are material coefficients related to the initial shear and bulk moduli respectively.

The material of the cutoff area Ω, on the other hand, exhibits highly anisotropic behavior due to presence of the defects. The microstructure may also be rearranged due to reorientation of the defects under large deformation. To simulate such an anisotropic feature, a compressible anisotropic strain energy density function, "Anisotropic NeoHookean" (ANH), is constructed by superposing Eqn. (7) with an anisotropic potential Φ as

$W=\Lambda (\frac{{{I}_{1}}}{I_{3}^{1/3}}-3)\text{ }+\frac{1}{2}K{{({{I}_{3}}^{1/2}-1)}^{2}}+\Phi ({{I}_{1}},\text{ }{{I}_{3}};{{I}_{p,}}_{\alpha })$

(8)

Various models of anisotropic potentials have been presented in the literature. A micromechanically-based form presented by Gasser, et al. ^[25] is employed here as

$\begin{align} & \Phi ({{I}_{1}},\text{ }{{I}_{3}};{{I}_{p,}}_{\alpha })\text{ }=\frac{{{k}_{1}}}{{{k}_{2}}}\sum\limits_{\alpha =1}^{n}{\{exp}[{{k}_{2}}(\zeta (\frac{{{I}_{1}}}{I_{3}^{1/3}}-3)+ \\ & (1-3\varsigma )(\frac{{{I}_{p,}}_{\alpha }}{I_{3}^{1/3}}-1){{)}^{2}}]-1\}\text{ } \\ \end{align}$

(9)

where k₁ and k₂ are material coefficients that characterize the anisotropic behavior, N is the number of families of preferred directions in the set of defects. The parameter $\zeta$ (0 ≦ $\zeta$ ≦ 1/3) describes the level of dispersion in the defect directions; when V = 0, the defects are perfectly aligned; when $\varsigma$ = 1/3, the defects are randomly distributed and the material becomes isotropic. The variables I_{p, α} (α=1, …, N) are pseudo-invariants of V² and defined as

${{I}_{p,}}_{\alpha }=\text{ }{{\mathbf{d}}_{\alpha }}{{\mathbf{V}}^{2}}{{\mathbf{d}}_{\alpha ~~~~~~}}(\alpha =\text{ }1,\text{ }\ldots ,\text{ }N)$

(10)

where d_α is a unit vector along the a-th preferred direction of the defects at its undeformed state. Also, in ANH, L + k₁ and Kare related to the initial shear and bulk moduli respectively.

---

5. Numerical Examples (Effective Material Models)

Three set of numerical examples are presented in the following three subsections. The problems are analyzed using finite elements. Quadratic elements are used for displacement interpolation in the calculation for both small and large deformation problems. No particular singular element is used throughout the study.

---

5.1. A Single Inclusion

In Problems 1.1 and 1.2, we consider a plane stressspecimen containing one square inclusion of size d (Figure 2). Both the specimen matrix and the inclusion are modeled by isotropic elastic materials. The inclusion is relatively small compared with the specimen so that the finite width effect is negligible. The single inclusion is taken as the whole cutoff area Ω. In such a case, the effective mechanical properties of Ω are analytically equivalent to those of the inclusion so that the accuracy of our numerical results can be properly examined.

Problem 1.1 A single inclusion (linear elasticity)

In this problem, the validity of our proposed approach in linear elasticity is illustrated. Here, we have the matrix and the inclusion modeled by Young's moduli E₀ and E_inc respectively. Poisson's ratios for both materials are the same and denoted as n₀. By applying two far-field stress conditions along its exterior boundaries, including uniaxial tension σ^∞and pure shear τ^∞, the two effective moduli E_eff and G_eff can then be independently calculated.

The study in this problem is organized as follows. First, the effect of the local finite element modeling in the near-defect area is investigated. Next, the numerical results associated with varying size of the specimen are observed. Finally, the effect due to the stiffness ratio E₀/E_inc is examined.

In order to accurately evaluate the effective moduli, it is required that the released strain energy Δ∏be computed with a high degree of precision. To this end, three FE models with very different near-defect local meshes are used to conduct a mesh convergence study. Details of the near-defect local portions for the three FE meshes are shown in Figure 3, where the multi-point constraint is applied to tie the fine and coarse elements in the second mesh so that displacement continuity is ensured. The normalized results of Δ∏, for both loading conditions, from the three meshes are tabulated in Table 1. Here, Δ∏is evaluated by using both the M-integral and the DDA. Although the mesh in the first model is quite coarse, the results from the three models show very good consistency, with deviations less than 1%. As observed, this energy parameter Δ∏ appears to be rather insensitive to the local near-defect FE models.

Next, three specimens (with aspect ratio w/d equal to 25, 60, and 240 respectively) are used to examine the finite width effect. Here, Δ∏ is evaluated by using the DDA. The normalized results of Δ∏ for both loading conditions, under two material combinations (with E₀/E_inc equal to 2 and 5), are shown in Table 2. The results obtained from the three aspect ratios show very good agreement. This indicates that, in our calculation, the extent of the specimen is chosen sufficiently large so that it can be regarded as an infinite matrix medium.

After Δ∏ being properly evaluated, the effective elastic moduli of the inclusion, E_eff and G_eff, are then determined by using the three procedures described in Section 2. The numerical results associated with the two loading conditions, under various material combinations (with E₀/E_inc ranging from 1 to ∞), are shown in Table 3, where E_eff and G_eff are normalized with respect to the analytical values E_inc and G_inc respectively. We can see that, while both the Eshelby's and the MCCT methods are analytically applicable for all material combinations, they may yield inaccurate results when E_o/E_inc becomes too large due to inevitable numerical truncation error.This can be illustrated in the limiting case when the inclusion corresponds to a void (i.e., E₀/E_inc → ∞), where the numerical results of E_eff and G_eff are relatively small yet not vanishing. As to the CCT method, it is observed that this approach is valid only when E₀/E_inc is very small, say, 1.5, as anticipated.

Problem 1.2 A single inclusion (hyperelasticity)

The matrix and inclusion are modeled by the NeoHookean strain energy density function, with the material coefficients denoted as ($\Lambda$₀, K₀) and ($\Lambda$_inc, K_inc) respectively. Here ($\Lambda$_inc, K_inc) are taken as ($\Lambda$₀/m, K₀/m). Two sets of far-field loads, i.e., uniaxial and biaxial tensile loads, are considered. The magnitude of the applied tension σ^∞is varying, with the resulting maximum principal stretch λ_max along the boundaries of the specimen ranging from 1 to 10. After Δ∏ being properly evaluated by using DDA, the effective coefficients ($\Lambda$_eff, K_eff) are obtained by using Eqn. (5) and fitting the data within the specified range of deformation.More details about the numerical procedure will be illustrated in the next example problem.

By considering various material combinations (with m ranging from 1 to ∞), the effective coefficients associated with the two loading conditions are shown in Table 4, where $\Lambda$_eff and K_eff are normalized with respect to $\Lambda$_inc and K_inc respectively. It is observed that the extended MCCT method yields very accurate solutions when m is smaller than, say, 5. Still, the results appear to deviate from the analytical values as the value of m becomes too large.

---

5.2. Parallel Cracks

In Problems 2.1 and 2.2, we consider a plane stress specimen containing a family of N parallel one-sized cracks (Figure 4).

These cracks, each of length l, are distributed in a doubly periodic manner, where s and a are the spacings of neighboring cracks in the parallel (x₁-) and perpendicular (x₂-) directions respectively. The crack distribution results in orthotropic mechanical behavior in the loading plane. The crack system is enclosed in a square cutoff area Ω of size d. The value of d is relatively small compared with the size w of the specimen. For such a crack system, a commonly-used measure of crack density parameter f in micromechanics is defined as

$f=~~\frac{N{{l}^{2}}}{4{{A}_{\Omega }}}$

(11)

In the present study, the crack system is taken to consist of 11 × 36 cracks in the cutoff area. Also, in the following calculations, the area of Ω and the spacing d are fixed, while crack length l is varying.

Problem 2.1 Parallel cracks (linear elasticity)

The original uncracked material in Ω is modeled by Young's modulus E₀ and Poisson's ratio 0.25. Three sets of far-field loads (i.e., uniaxial tension, pure shear, and biaxial tension) are considered. After Δ∏ being properly evaluated by using DDA with finite elements, the effective moduli are then determined by using MCCT.Note that accurate solutions can be achieved by appropriately adjusting the material properties of the matrix. The results of the normalized effective Young's modulus E_{22, eff}/E₀, shear modulus G_{12, eff/}G₀, and bulk modulus B_eff/B₀ versus the crack density parameter f are shown in Figure 5(a) and 5(b) respectively. Also included in the figures are the solutions of E_{22, eff} and G_{12, eff} from two other approaches (Nemat-Nasser, et al. ^[1], Shen and Li ^[4]). As can be seen, our results are in good agreement with those from these methods. As an aside, although not depicted in detail here, our calculations show that a very fine FE mesh is not required in the near-tip area to have accurate results. The feasibility for combined use of DDA and MCCT in linear elasticity is thus illustrated.

Problem 2.2 Parallel cracks (hyperelasticity)

The original uncracked material in Ω is modeled by the NeoHookean strain energy density function with the material coefficients denoted as ($\Lambda$₀, K₀). Three types of far-field loads (including uniaxial (x₁- and x₂-directions), and biaxial tensile loads) are considered, with the resulting maximum stretch λ_max ranging from 1 to 16. The orthotropic effective strain energy density W_eff of Ω is characterized by the ANH model as

$\begin{align} & W=\Lambda (\frac{{{I}_{1}}}{I_{3}^{1/3}}-3)\text{ }+\frac{1}{2}{{K}_{eff}}{{({{I}_{3}}^{1/2}-1)}^{2}}+ \\ & \frac{{{k}_{1,eff}}}{{{k}_{2,eff}}}\{exp[{{k}_{2}}eff{{(\frac{{{I}_{p,1}}}{I_{3}^{1/3}}-1)}^{2}}]-1\} \\ \end{align}$

(12)

where I_{p, 1} is determined from Eqn. (10) by specifying d₁ in x₁-direction, and k_2.eff is taken as 0.0001.

After Δ∏ being evaluated by using DDA, W_eff is determined by using extended MCCT (Eqn. (5)). The results of W_eff for, say, f = 0.367 under the three loading conditions are scaled with respect to $\Lambda$_o and depicted versus λ_max in Figure 6. By substituting these data into Eqn. (12), the material coefficients ($\Lambda$_eff, K_eff, k_{1, eff}) can then be solved with a least squares solution scheme. Also shown in Figure 6 are the resulting fitted curves for the three loading conditions. Note that, typically in hyperelasticity, the material data can be well fit within only a rather small range of deformation, usually for Δλ_max < 3 ^[26]. In this sense, we can see that our given model yields very good results, which are will fit within almost the whole range of deformation for the uniaxial-x₂ data, and in the range of very large deformation (12< λ_max < 17) for the uniaxial-x₁ and biaxial data.

The normalized results of the effective coefficients $\Lambda$_eff/$\Lambda$₀, K_eff/K₀, and k_{1, eff}/$\Lambda$₀ versus the crack density parameter f are shown in Figure 7. It can be seen that $\Lambda$_eff decreases monotonically with respect to f. On the other hand, as f increases, the anisotropic coefficient k_{1, eff} increases and thus makes more significant contribution to the initial shear modulus. Further, by comparing the results of $\Lambda$_eff + k_{1, eff} and K_eff with those of G_{12, eff} and B_eff for the linear case (Figure 5(b)), it is observed that they both coincide closely with each other. The physical meaning of L + k₁ and Kas the initial shear and bulk moduli is thus evident.

---

5.3. An Elliptical Void

In Problems 3.1 and 3.2, we consider a plane stress specimen containing a central elliptical void (Figure 8).

This void is enclosed in a square cutoff area Ω of size d and results in orthotropic mechanical behavior in the plane. The size d is relatively small compared with w. Three geometric instances are considered. The first two (Cases Ⅰ and Ⅱ) are circular holes with area of 5% and 25% of A_Ωrespectively. The third (Case Ⅲ) is an elliptical void with aspect ratio d₁/d₂ = 1.8 and area of 35% of A_Ω.

Problem 3.1 An elliptical void (linear elasticity)

The original unvoided material in Ω is modeled by Young's modulus E₀ and Poisson's ratio 0.3. Four sets of far-field loads are considered, including uniaxial tension in x₁- and x₂-directions, pure shear, and biaxial tension. The effective Young's moduli (E_{11, eff} and E_{22, eff}), shear modulus G_{12, eff}, and bulk modulus B_eff are evaluated by combined use of these data. The resulting normalized values for the three geometric cases are shown in Table 5. It is observed that, for Cases Ⅰ and Ⅱ, the effective moduli due to the circular voids turn out to slightly deviate from isotropy. Such feature is actually anticipated because the calculations are based on a square cutoff area Ω. As to the orthotropic case (Case Ⅲ), E_{11, eff} appears to be higher than E_{22, eff} and G_{12, eff}, as anticipated.

Problem 3.2 An elliptical void (hyperelasticity)

The original unvoided material in Ω is modeled by the NeoHookean strain energy density function with the material coefficients ($\Lambda$₀, K₀). Three types of far-field loading conditions are considered, including uniaxial (x₁- and x₂-directions) and biaxial tensile loads, with the resulting λ_max ranging from 1 to 16. The deformed finite element mesh in Figures 9(a) and 9(b) show the response of the specimen for Case Ⅲ under the action of uniaxial (x₁-direction) and biaxial tension, respectively, for λ_max = 3, where the undeformed configuration is depicted in heavy lines. Figures 10(a) and 10(b) show the distribution of maximum principal Cauchy stress forthe two scenarios in the near-void local area under the deformed state.

By taking the ANH strain energy density function (Eqn. (12)) with d₁ specified in x₁-direction, the results of the normalized effective coefficients $\Lambda$_eff/$\Lambda$₀, K_eff/K₀, and k_{1, eff}/$\Lambda$₀ for the three geometric cases are shown in Table 6. It is observed that, for Cases Ⅰ and Ⅱ, the "nearly"-isotropic behavior can be characterized by using only the first two effective coefficients $\Lambda$_eff and K_eff, which appear to be very close to the values of G_{12, eff} and B_eff in the linear case (Table 5). Also, by comparing the results of $\Lambda$_eff + k_{1, eff} and K_eff for Case Ⅲ with those of G_{12, eff} and B_eff in the linear case (Table 5), the physical meaning of L + k₁ and Kis again evident. Further, the results of W_eff/$\Lambda$₀ versus λ_max, along with the fitted curves for the three loading conditions, are shown in Figure 11. Still, we can see that our given model yields very good results, which are well fit within almost the whole range of deformation for both the uniaxial-x₂ and biaxial data, and in the range of very large deformation (10 < λ_max< 15) for the uniaxial-x₁ data.

---

6. Numerical Example

In this section, we consider a plane strain body of height h₀ and with infinite extent in x₁-direction. The body consists of the bulk and a thin material layer with height d (Figure 12(a)). Both the materials are modeled with hyperelasticity. The initial Young's moduli for the bulk and the layer are E₀ and E_m respectively, and Poisson's ratios are both 0.3. Perfect bonding between the materials is assumed. In many engineering applications (e.g., adhesive bonding layer, masonry, biomedical tissues, etc.), the material layer typically possesses substantially weaker mechanical properties than the surrounding bulk. More possibly, its material integrity can further be degraded if the layer contains heterogeneities such as microvoids or microcracks. In the present example, the material layer with a microstructure composed of a series of periodically-distributed identical elliptical voids is considered. A square "representative area element" (RAE) of width d (Figure 12(b)), containing a central elliptical void, is thus delimited for modeling the effective properties on the macro level. The following study is conducted by considering three instances of microstructures, including no void, a circular hole of 25% void ratio, and an elliptical void of 35% area ratio. The body is subjected to mixed-mode displacements (u₁, u₂) along its top and bottom boundaries. Both the local stress intensity in this layer and the global reaction at the boundaries are thus significantly governed by the geometric and material properties.

The FE analysis is performed by choosing a "representative column" of width w₀, as shown in Figure 12(a), where the infinite extent of the specimen is represented by imposing periodic boundary conditions along the two lateral sides. The FE mesh for the geometric instance of d= h₀/20 is shown in Figure 13(a). The effective ANH material coefficients listed in Table 6 are used for modeling the three microstructure instances of the material layer, with E_m = E₀/10. A shear-dominated pair of displacements (u₁, u₂) = (u₁, 0.1u₁) are applied incrementally until a final state of deformation due to u₁ = 0.1h₀ is reached. Mesh studies have been performed to ensure convergence of the solutions. The deformed FE meshes in Figure 13(b)-13(d) show the responses of the three microstructures at their final deformed states, where the local deformation of the material layer depends substantially on its stiffness. Figure 14(a) and 14(b) illustrate the results of the reactions R₁ and R₂ (in x₁- and x₂-directions, respectively) at the upper boundary versus the applied displacements. The values are normalized with respect to R_{1, ref}, the reaction in x₁-direction at u₁ = 0.1h₀ for the case of no void. Also included in the figures are the results from Hirschberger et al. ^[14], where a multiscale-based computational homogenization scheme is used for the solution. It appears that they are very well consistent for most instances in their macroscopic responses and the feasibility of our material constitutive model is therefore illustrated.

In order to further investigate the effect of the displacement mode, we define a dimensionless parameter "mode mixity" $\eta$ as

$\eta =\text{ }{{u}_{1}}/\left( {{u}_{1}}+\text{ }{{u}_{2}} \right)~$

(13)

The variations of R₁/R_{1, ref} and R₂/R_{1, ref} for the two microstructures (i.e., the circular and the elliptic voids) are depicted with respect to hin Figure 15. It is observed that, while R₁ varies almost linearly with respect to h, R₂ bears a more nonlinear trend as it decreases more significantly as hdecreases. As an aside, the reactions under varying material layer height d are also evaluated. The results for the microstructure with circular voids, deformed at (u₁, u₂) = (0.1, 0.01) h₀, are shown in Figure 16. We can see that the reactions decrease as d increases, as anticipated. It is further observed that such weakening trend of the global stiffness becomes less significant as the layer becomes thinner.

Although the presented numerical examples are all two-dimensional cases, it is straightforward to make the model extended to general three-dimensional cases. To this end, by applying combinations of remote boundary conditions in the third dimension and evaluating the associated released strain energy Δ∏, the three-dimensional effective strain energy density function W_eff can then be properly constructed.

---

7. Conclusions

An anisotropic hyperelastic constitutive model ANH is presented to quantitatively characterize the effective mechanical response of multi-defected solids under the action of large deformation. The corresponding effective strain energy density function can be properly constructed by combined use of the DDA and the extended MCCT methods. With our proposed approach, the associated material coefficients can then be easily determined by using a scalar energy parameter Δ∏, which can be directly extracted from regular FE solutions.

The feasibility of the proposed approach is illustrated by considering a series of numerical examples. By comparing with the results from other computational schemes, it is observed that our presented approach provides quite good approximations with sufficient accuracy in engineering applications and can be easily implemented in the FE model with no extra computational cost. When compared with other multi-scale computational homogenization techniques, the approach thus appears to be a very efficient computational device and is straightforward for use in practice.

---

Acknowledgments

This work has been partially supported by Ministry of Science and Technology Grant No. MOST 104-2221-E-008-068 to National Central University.

---

Conflict of Interest

The authors declare that there is no conflict of interest regarding the publication of this manuscript.

---