Finite element solution of crack-tip fields for an elastic porous solid with density-dependent material moduli and preferential stiffness

Abstract

In this paper, the finite element solutions of crack-tip fields for an elastic porous solid with density-dependent material moduli are presented. Unlike the classical linearized case in which material parameters are globally constant under a small strain regime, the stiffness of the model presented in this paper can depend upon the density with a modeling parameter. The proposed constitutive relationship appears linear in the Cauchy stress and linearized strain independently. From a subclass of the implicit constitutive relation, the governing equation is bestowed via the balance of linear momentum, resulting in a quasi-linear partial differential equation (PDE) system. Using the classical damped Newton’s method, the sequence of linear problems is then obtained, and the linear PDEs are discretized through a bilinear continuous Galerkin-type finite element method. We perform a series of numerical simulations for material bodies with a single edge-crack subject to a variety of loading types (i.e. the pure mode-I, II, and mixed-mode). Numerical solutions demonstrate that the modeling parameter in our proposed model can control preferential mechanical stiffness with its sign and magnitude along with the change of volumetric strain. This study can provide a mathematical and computational foundation to further model the quasi-static and dynamic evolution of cracks, utilizing the density-dependent moduli model and its modeling framework.

Introduction

Mechanical stiffness of a porous solid is dependent not only on the mechanical property of solid grain but also on the bulk skeleton composed of connected or disconnected pores (or microstructures). For example, the shear stiffness of sandy soil can be lost and act like liquid when soil liquefaction occurs. ¹ It is also well known that the mechanical properties of metal matrix composite or alloy are largely dependent on the defects, that is, porosity. ² Furthermore, the pore space can be partially or fully saturated with fluid (e.g. gas or liquid) inside, which can induce the poromechanical effects.^3–5 In the field of petroleum engineering and rock mechanics, ⁶ hydraulic fracturing^7,8 (or commonly known as fracking) is a well-known stimulation technology for the development of shale gas with high pressure gradient, which aims to generate tensile fractures in the rock. The success of hydraulic fracturing is highly dependent on the mechanical properties of rock,^7,9 as it is much more efficient for sand which is stiffer than soft clay or mud for tensile failure. Considering the heterogeneity in the composites, simple averaging of rock properties might mislead, but rigorous mathematical formulation (e.g. the multiple porosity model^4,10) is required.

Material properties and mechanical strength can also be reliant on certain modes of loading preferentially. For instance, brittle materials are generally weak in tension, and their tensile strength is known to be only one-tenth of compressive strength. ¹¹ Meanwhile, ductile materials are known to possess approximately equal strength in tension and compression, but weak in shear. Experimental results regarding hydraulic fracturing also substantiate that shear failure may occur in advance before the tensile failure around the crack-tip, ¹² where mixed-mode of loading exists. Regarding this stiffness anisotropy, we particularly pay attention to the density change of a porous material. Microstructures of material including pore space are required to bond with one another to transmit the tensile loading. ¹¹ Meanwhile, the fracture propagation is demanding under a compressive loading, as transverse cracks can emerge and tend to congregate, resulting in the increased density of material. Representatively, concrete or ceramics is weak under tension but strong against compression, thus composite (or reinforced) material can be synthesized with them to have higher tensile strength. Rock is another exemplary material such that it is stiffer against compressive loading but weaker against tensile or shear loading. On the other end, two-dimensional (2D) material ¹³ (e.g. the graphene) is considered to be one of the strongest materials able to withstand tensional loading, due to the 2D nature of the structure. In essence, the preferential strength can be related to material moduli that are dependent upon density of the material or change in volumetric strain: the dilation versus compaction.

Furthermore, this density-dependent stiffness concept is physically reasonable in that elastic properties of a porous material are known to have nonlinear relations along with the porosity,^14,15 even though it is not straightforward to model these in the continuum scale. Conventional practice of modeling porous solid is based on the elastic and infinitesimal strain regime, wherein the Cauchy stress is expressed as a function of deformation gradient and density at material points. This classical relationship, where material moduli are generally constants, has been successfully used in various applications. However, this classical linearized theory of elasticity cannot accommodate the phenomena found in many metallic alloys (see Saito et al. ¹⁶ and Hou et al. ¹⁷ ) and concrete (see Grasley et al. ¹⁸ ) that exhibit nonlinear mechanical responses well within the range of “small strains,” nor aforementioned realistic responses.

The objective of this study to investigate elastic porous solids of preferential stiffness with an explicit type of the density-dependent moduli model. Recently, Rajagopal^19–23 has established a new class of constitutive relations for the description of non-dissipative elastic bodies (for more details about “novel” elastic constitute relationships). These implicit constitutive relations have been explored in many studies to describe the state of stress-strain in the neighborhood of crack-tip,^24–33 for the behavior of thermoelastic bodies,^34,35 viscoelastic materials,^36–42 and quasi-static crack evolution.^43,44 Utilizing the implicit constitutive relation in Rajagopal^45,46 and Rajagopal and Saccomandi ⁴⁷ and the density-dependent moduli derived with a simplified parameter of the model, we address the preferential or anisotropic stiffness of a material with in the sense of tension/compaction that is, change in void: dilation and compaction. Another goal is to weigh up the applicability of the density-dependent material moduli model for the quasi-static or dynamic propagation of a network of cracks under mechanical/thermal loading.

To this end, we assume a homogeneous porous solid with a crack inside, where the mechanical regime is under small strain and pure elasticity without any failure. Particularly for the stiffness variations, we focus on the crack (or damaged pores) and its tip under different modes of loadings: mode-I, II, and mixed-mode. Instead of an explicit porosity, that is, one of the state variables based on the mixture theory^3,48,49 that can also be saturated with fluid inside, a pure solid matrix with implicit porosity is considered through the change of volumetric strain. We employ a computational model that is physically and mathematically consistent with the framework developed in Yoon et al. ⁴³ and Lee et al. ⁴⁴ a universal approach based on the Newton’s method and standard finite element method (FEM). This approach yields a numerically stable and efficient algorithm against the serious nonlinear problem of our interest. In the numerical experiments, our iterative algorithm displays an optimal order convergence for a BVP with a manufactured solution. Comparing this nonlinear model with the conventional linearized elasticity model, distinct variations of stress and strain distributions are demonstrated, especially near the crack-tip area under different modes of loading, that is, parallel and perpendicular to the crack. We also compare the stress intensity factor, drained bulk modulus, and volumetric strain change to elucidate the preferential stiffness of the elastic porous solid. Furthermore, we find that the intensity of preferential stiffnessis controllable through the parameter of the model, relating the void to density variation of a porous solid. Upon the efficacy of the density-dependent material moduli model and its modeling framework, current study can further be expanded to the topics such as the fracture evolution (hydraulic fracture), fluid transport, or heat transfer in deformable porous media.

Formulation of the density-dependent material moduli model

The current investigation is prompted by a recent study (see Rajagopal^45,46 and Rajagopal and Saccomandi ⁴⁷ for more details) on developing new nonlinear constitutive relationships between the linearized strain and Cauchy stress. The density-dependent material moduli model is worthwhile in characterizing the cracks or damaged pores in porous solids such as rocks and concrete.

Basic notations

In this section, we provide a brief introduction to the geometrical and mechanical description of the elastic porous solid body, highlighting its mass balance. Let $\mathrm{\Omega }\subset {\mathbb{R}}^2$ be a closed and bounded Lipschitz domain in the Euclidean space that represents the material body under consideration in the reference configuration. Let $\partial \mathrm{\Omega }$ be the Lipschitz continuous boundary and $\mathbf{\eta }=({\eta }_1,{\eta }_2)$ be the outward unit normal. Further, we consider a boundary part consists of two disjoint parts such that:

$$\partial \mathrm{\Omega }=\overline{{\mathrm{\Gamma }}_D}\cup \overline{{\mathrm{\Gamma }}_N}\mathrm{and}{\mathrm{\Gamma }}_D\cap {\mathrm{\Gamma }}_N=Ø,$$

where the partition consists of a Neumann boundary ${\mathrm{\Gamma }}_N$ and a nonempty Dirichlet boundary ${\mathrm{\Gamma }}_D$ . Let ${\mathrm{\Gamma }}_c\subseteq \mathrm{\Omega }$ be a geometrical boundary of an interface and we assume that the interface ${\mathrm{\Gamma }}_c$ is completely contained inside of $\mathrm{\Omega }$ , not on the boundary $\partial \mathrm{\Omega }$ . Let $\mathit{X}=(X_1,X_2)$ and $\mathit{x}=(x_1,x_2)$ denote typical points in the reference and deformed configurations of the body. Let $\mathit{u}:\mathrm{\Omega }\to {\mathbb{R}}^2$ denote the displacement field, and

$$\mathit{u}=\mathit{x}-\mathit{X}.$$

(1)

In the rest of this paper, we use the usual notations of Lebesgue and Sobolev spaces.^50,51 We denote $L^p(\mathrm{\Omega })$ as the space of all Lebesgue integrable functions with $p\in [1,\infty )$ . In particular, when $p=2$ , $L^2(\mathrm{\Omega })$ denote the space of all square-integrable (Lebesgue) functions on $\mathrm{\Omega }$ together with its inner product $(v,w):={\int }_{\mathrm{\Omega }}vwd\mathit{x}$ and the norm $\parallel v\parallel :=\parallel v{\parallel }_{L^2(\mathrm{\Omega })}={(v,v)}^{1/2}$ . Let $C^m(\mathrm{\Omega }),m\in {\mathbb{N}}_0$ denote the linear space of continuous functions on $\mathrm{\Omega }$ and let $H^1(\mathrm{\Omega })$ denote the standard Sobolev space:

$$H^1(\mathrm{\Omega }):=W^{1,2}=\{v\in L^2(\mathrm{\Omega });{\partial }_{j}v\in L^2(\mathrm{\Omega })j\in \{1,2,\dots n\}\},$$

(2)

with the inner (scalar) product and the norm defined, respectively, as:

$${(f,g)}_{H^1}:=(f,g)+\sum _{j=1}^n({\partial }_{j}f,{\partial }_{j}g),$$

(3a)

$$\parallel f{\parallel }_{H^1}:={(\parallel f^2{\parallel }^2+\sum _{j=1}^n\parallel {\partial }_j{f}^2{\parallel }^2)}^{1/2}.$$

(3b)

The closure of $C_0^{\infty }(\mathrm{\Omega })$ by the norm of $H^1(\mathrm{\Omega })$ is denoted as $H_0^1(\mathrm{\Omega })$ , that is,

$$H_0^1(\mathrm{\Omega })={\overline{C_0^{\infty }(\mathrm{\Omega })}}^{\parallel ·{\parallel }_{H^1}}.$$

(4)

The space of displacements satisfying the homogeneous and non-homogeneous Dirichlet boundary conditions are defined as:

$${\hat{V}}_0:=\{\mathit{u}\in {(H^1(\mathrm{\Omega }))}^2:\mathit{u}=0\mathrm{on}{\mathrm{\Gamma }}_D\},$$

(5)

$${\hat{V}}_{\mathit{g}}:=\{\mathit{u}\in {(H^1(\mathrm{\Omega }))}^2:\mathit{u}=\mathit{g}\mathrm{on}{\mathrm{\Gamma }}_D\}.$$

(6)

Let $\mathrm{Sym}({\mathbb{R}}^{2\times 2})$ is the space of $2\times 2$ symmetric tensors equipped with the inner product $\mathit{A}:\mathit{B}=\sum _{i,j=1}^d{\mathit{A}}_{\mathrm{ij}}{\mathit{B}}_{\mathrm{ij}}$ and for all $\mathit{A}=(\mathit{A}{)}_{\mathrm{ij}}$ and $\mathit{B}=(\mathit{B}{)}_{\mathrm{ij}}$ in $\mathrm{Sym}({\mathbb{R}}^{2\times 2})$ , the associated norm $\parallel \mathit{A}\parallel =\sqrt{\mathit{A}:\mathit{A}}$ . Further, let $\mathit{F}:\mathrm{\Omega }\to {\mathbb{R}}^{2\times 2}$ denote the deformation gradient, $\mathit{C}:\mathrm{\Omega }\to {\mathbb{R}}^{2\times 2}$ denote the right Cauchy-Green stretch tensor, $\mathit{B}:\mathrm{\Omega }\to {\mathbb{R}}^{2\times 2}$ denote the left Cauchy-Green stretch tensor, $\mathit{E}:\mathrm{\Omega }\to {\mathbb{R}}^{2\times 2}$ denote the Lagrange strain, $\mathbf{ϵ}:\mathrm{\Omega }\to \mathrm{Sym}({\mathbb{R}}^{2\times 2})$ denote the linearized strain tensor, respectively defined as:

$$\mathit{F}:={\nabla }_r\mathit{x}=\mathit{I}+\nabla \mathit{u},$$

(7a)

$$\mathit{B}:=\mathit{F}{\mathit{F}}^{\mathrm{T}},\mathit{C}:={\mathit{F}}^{\mathrm{T}}\mathit{F},\mathit{E}:=\frac{1}{2}(\mathit{C}-\mathit{I}),$$

(7b)

$$\mathbf{ϵ}:=\frac{1}{2}(\nabla \mathit{u}+\nabla {\mathit{u}}^{\mathrm{T}}),$$

(7c)

where ${(·)}^{\mathrm{T}}$ denotes the transpose operator for the second-order tensors, $\mathit{I}$ is the two-dimensional identity tensor, ${\nabla }_r$ , and ∇ are the gradient operators in the reference and current configuration, respectively. Under the standard assumption of linearized elasticity, we have the following

$$\underset{\mathit{X}\in \mathrm{\Omega }}{max}\parallel \nabla \mathit{u}\parallel <<\mathcal{O}(\delta ),\delta \le 1.$$

(8)

Henceforth, we shall not distinguish the dependence of the quantities on $\mathit{X}$ for the notational convenience. The premise of infinitesimal strains (8) implies

$$\mathit{B}\approx \mathit{I}+2\mathbf{ϵ},\mathit{C}\approx \mathit{I}+2\mathbf{ϵ},\mathit{E}\approx \mathbf{ϵ},$$

(9a)

$$det\mathit{F}=1+\mathrm{tr}(\mathbf{ϵ}).$$

(9b)

Let $\mathit{T}:\mathrm{\Omega }\to \mathrm{Sym}({\mathbb{R}}^{2\times 2})$ be the Cauchy stress tensor in the current configuration and it satisfies the linear momentum balance:

$$\rho \stackrel{··}{\mathit{u}}=\mathrm{div}\mathit{T}+\rho \mathit{b},$$

(10)

where $\rho$ is the density in the current configuration and $\mathit{b}:\mathrm{\Omega }\to {\mathbb{R}}^2$ is the body force in the current configuration and the notation $\stackrel{·}{\left(·\right)}$ denotes the time derivative. The first and second Piola-Kirchhoff stress tensor tensors, $\mathit{S}:\mathrm{\Omega }\to {\mathbb{R}}^{2\times 2}$ and $\overline{\mathit{S}}:\mathrm{\Omega }\to {\mathbb{R}}^{2\times 2}$ , are defined by

$$\mathit{S}=\mathit{T}{\mathit{F}}^{-\mathrm{T}}det\mathit{F},\overline{\mathit{S}}:={\mathit{F}}^{-1}\mathit{S}.$$

(11)

The principle of angular momentum balance implies that the Cauchy stress tensor is symmetric, that is,

$$\mathit{T}={\mathit{T}}^{\mathrm{T}}.$$

(12)

The balance of mass (or continuity equation) in the material description follows that

$${\rho }_0=\rho det\mathit{F},$$

(13)

where ${\rho }_0$ is the reference density. In the view of (9b), the balance of mass (13) reduces to

$${\rho }_0=\rho (1+\mathrm{tr}(\mathbf{ϵ})).$$

(14)

More details about the kinematics and kinetics can be found in Truesdell and Noll. ⁵²

Implicit constitutive relations

In this study, our focus is to study the behavior of elastic (nondissipative) porous solids whose material moduli depend on the density. The response of such materials to the mechanical loading can be straightforwardly described by implicit constitutive relations, introduced by Rajagopal^45,46 and the references therein.

A generalization for the elastic body defined through an implicit type constitutive relation^19,20 between the Cauchy stress and deformation gradient is of the form:

$$\mathit{f}(\rho ,\mathit{T},\mathit{F},\mathit{X})=0,$$

(15)

where $\mathit{f}$ is a tensor-valued function. In the case of isotropic bodies, the above implicit response relation reduces to

$$\mathit{f}(\rho ,\mathit{T},\mathit{B})=0.$$

(16)

Assuming that the function $\mathit{f}$ is isotropic, then one can write the most general form of the implicit constitutive relation:

$$\begin{array}{c}\multicolumn{1}{c}{0={\delta }_0\mathit{I}+{\delta }_1\mathit{T}+{\delta }_2\mathit{B}+{\delta }_3{\mathit{T}}^2+{\delta }_4{\mathit{B}}^2+{\delta }_5(\mathit{TB}+\mathit{BT})}\\ \multicolumn{1}{c}{+{\delta }_6({\mathit{T}}^2\mathit{B}+\mathit{B}{\mathit{T}}^2)+{\delta }_7({\mathit{B}}^2\mathit{T}+\mathit{T}{\mathit{B}}^2)+{\delta }_8({\mathit{T}}^2{\mathit{B}}^2+{\mathit{B}}^2{\mathit{T}}^2),}\end{array}$$

(17)

where the material moduli ${\delta }_i$ for $i=0,1,\dots ,8$ are scalar functions that depend on the density and the invariants for the pair $\mathit{T}$ and $\mathit{B}$ , that is,

$$\begin{array}{c}\multicolumn{1}{c}{\{\rho ,\mathrm{tr}(\mathit{T}),\mathrm{tr}(\mathit{B}),\mathrm{tr}({\mathit{T}}^2),\mathrm{tr}({\mathit{B}}^2),\mathrm{tr}({\mathit{T}}^3),\mathrm{tr}({\mathit{B}}^3),\mathrm{tr}(\mathit{T}\mathit{B}),}\\ \multicolumn{1}{c}{\mathrm{tr}({\mathit{T}}^2\mathit{B}),\mathrm{tr}(\mathit{T}{\mathit{B}}^2),\mathrm{tr}({\mathit{T}}^2{\mathit{B}}^2)\}.}\end{array}$$

(18)

Using the linearization assumption given in (8) and the subsequent asymptotic results for the classical terms (9), we obtain a special subclass of the implicit constitutive relations of the form:

$$\begin{array}{c}\multicolumn{1}{c}{{\hat{\delta }}_0\mathit{I}+{\hat{\delta }}_1\mathbf{ϵ}+{\hat{\delta }}_2\mathit{T}+{\hat{\delta }}_3{\mathit{T}}^2+{\hat{\delta }}_4(\mathbf{ϵ}\mathit{T}+\mathit{T}\mathbf{ϵ})}\\ \multicolumn{1}{c}{+{\hat{\delta }}_5(\mathbf{ϵ}{\mathit{T}}^2+{\mathit{T}}^2\mathbf{ϵ})=0,}\end{array}$$

(19)

where the terms ${\hat{\delta }}_i,i=0,\dots ,5$ are the functions of scalar-valued invariants of $\mathbf{ϵ}$ and $\mathit{T}$ . But the terms ${\hat{\delta }}_i,i=0,2,3$ are the functions linearly depend upon the invariants of $\mathbf{ϵ}$ and arbitrarily upon the invariants of the $\mathit{T}$ , while the terms ${\hat{\delta }}_i,i=1,4,5$ depends on the invariants of $\mathit{T}$ . An important subclass of the above general class of implicit constitutive relations (19) is

$$\mathbf{ϵ}={\hat{\alpha }}_0\mathit{I}+{\hat{\alpha }}_1\mathit{T}+{\hat{\alpha }}_2{\mathit{T}}^2,$$

(20)

where the material moduli ${\hat{\alpha }}_i,i=0,1,2$ depend on $\rho ,\mathrm{tr}\mathit{T},\mathrm{tr}{\mathit{T}}^2,\mathrm{tr}{\mathit{T}}^3$ , and by virtue of the mass balance (as in (14)). The model (20) has been studied in several investigations involving cracks in elastic bodies exhibiting strain-limiting behavior,^{24,28–30,53} quasi-static crack evolution,^43,44 thermo-elastic bodies,^34,35 quasi-linear viscoelastic bodies,^54–56 and nonlinear constitutive model for rock.^57,58

Another subclass of models of the above general class of relations (19) wherein the constitutive relation is linear in both $\mathbf{ϵ}$ and $\mathit{T}$ is given by (see also Itou et al. ⁵⁹ )

$$(1+{\kappa }_3\mathrm{tr}\mathit{T})\mathbf{ϵ}=C_1(1+{\kappa }_1\mathrm{tr}\mathbf{ϵ})\mathit{T}+C_2(1+{\kappa }_2\mathrm{tr}\mathbf{ϵ})\left(\mathrm{tr}\mathit{T}\right)\mathit{I}.$$

(21)

In the above relation (21), the moduli ${\kappa }_1,{\kappa }_2,{\kappa }_3,C_1,C_2$ are all constants. Note that the above model is both linear in stress and strain, and it can be used to describe the mechanical response of porous solids. The classical model can be recovered from the above model by choosing ${\kappa }_1,{\kappa }_2,{\kappa }_3$ as zero, and then we can identify

$$C_1=\frac{1+\nu }{E}=\frac{1}{2\mu }>0,C_2=-\frac{\nu }{E}<0,$$

(22)

where $E$ is the Young’s modulus and $\nu$ is the Poisson’s ratio, and these are related to the Lamè constants $\lambda$ and $\mu$ :

$$\lambda =\frac{E\nu }{(1+\nu )(1-2\nu )},\mu =\frac{E}{2(1+\nu )}.$$

(23)

Taking ${\kappa }_3=0$ in (21), we can obtain a special constitutive relation

$$\mathbf{ϵ}=C_1(1+{\kappa }_1\mathrm{tr}\mathbf{ϵ})\mathit{T}+C_2(1+{\kappa }_2\mathrm{tr}\mathbf{ϵ})\left(\mathrm{tr}\mathit{T}\right)\mathit{I}.$$

(24)

By the balance of mass in the reference configuration (as in (14)), one can express the above relation as

$$\mathbf{ϵ}={\mathrm{\Xi }}_1^{C_1}(\rho ,\mathrm{tr}\mathbf{ϵ})\mathit{T}+{\mathrm{\Xi }}_2^{C_2}(\rho ,\mathrm{tr}\mathbf{ϵ})\left(\mathrm{tr}\mathit{T}\right)\mathit{I},$$

(25)

with both functions ${\mathrm{\Xi }}_1^{C_1}(\rho ,\mathrm{tr}\mathbf{ϵ}),{\mathrm{\Xi }}_2^{C_2}(\rho ,\mathrm{tr}\mathbf{ϵ})$ being linear in the density variable $\rho$ , hence linear in $\mathrm{tr}\mathbf{ϵ}$ .

Remark 2.1. It is an interesting and challenging open topic about characterizing the crack-tip fields in the density-dependent model (25). The question is about whether the above model shares the same feature of crack-tip strain (as well as crack-tip stress) singularities as the classical linearized model. However, the general strain-limiting models^{24,26,27,29,31,32,60,61} lead to the bounded crack-tip strain behavior. The uniform upper bound on the strain in the entire body can be fixed a priori to a value, as small and bounded as modeled. In all these models the material parameters are globally constant. Our model (25) is certainly useful in describing the behavior of porous elastic solids, as the functions ${\mathrm{\Xi }}_1^{C_1}(\rho ,\mathrm{tr}\mathbf{ϵ})$ and ${\mathrm{\Xi }}_2^{C_2}(\rho ,\mathrm{tr}\mathbf{ϵ})$ can be recognized as the density-dependent material moduli.

Remark 2.2. We can also introduce a new class of implicit constitutive model to describe the response elastic material via splitting the stress in terms of its deviatoric and volumetric parts. Within the framework of models of type (21), such stress splitting leads to a three-field (displacement, volumetric, and deviatoric stresses) mixed formulation of the nonlinear elasticity model. Interesting issues of the static-crack field, such as the non-penetrating cracks in bodies, ⁵³ nucleation, and damage evolution, can then be addressed appropriately.

Boundary value problem and numerical method

In this section, we develop a boundary value problem and the numerical approach for computing the stress and strain fields in the elastic porous solid body which includes a crack-tip.

Mathematical model

The elastic material body under consideration is homogeneous, isotropic, initially unstrained, and unstressed. Let $\mathrm{\Omega }$ be an open, bounded, Lipschitz, and connected domain with the boundary $\partial \mathrm{\Omega }$ consisting of two smooth disjoint parts ${\mathrm{\Gamma }}_N$ and ${\mathrm{\Gamma }}_D$ such that $\partial \mathrm{\Omega }=\overline{{\mathrm{\Gamma }}_N}\cup \overline{{\mathrm{\Gamma }}_D}$ where ${\mathrm{\Gamma }}_D\cap {\mathrm{\Gamma }}_N=Ø$ . To describe the state stress of the material under investigation, we consider the balance of linear momentum for the static body in the absence of body forces, which reduces to

$$-\nabla ·\mathit{T}=0\mathrm{in}\mathrm{\Omega }.$$

(26)

To model the stress response of the material whose material moduli are dependent upon the density, for example, the elastic porous solids, we make use of implicit constitutive relations (see Rajagopal and Wineman ⁶² for details). A special subclass of such elastic bodies can be obtained if one starts with (24) and use ${\kappa }_1={\kappa }_2=\beta$ as:

$$\mathbf{ϵ}=\frac{(1+\nu )(1+\beta \mathrm{tr}(\mathbf{ϵ}))}{E}\mathit{T}-\frac{\nu (1+\beta \mathrm{tr}(\mathbf{ϵ}))}{E}\mathrm{tr}(\mathit{T})\mathit{I}.$$

(27)

We note that the above constitutive relationship is invertible, and the inverted relation is written for the Cauchy stress as:

$$\mathit{T}=\frac{\mathbb{E}[\mathbf{ϵ}]}{1+\beta \mathrm{tr}(\mathbf{ϵ})},$$

(28)

which is the definition for the stress of the proposed density-dependent model. Plugging (28) into (27), we note that the stress and strain relation reduces to the same form of the linearized elasticity:

$$\mathbf{ϵ}=\frac{(1+\nu )}{E}\mathbb{E}[\mathbf{ϵ}]-\frac{\nu }{E}\mathrm{tr}(\mathbb{E}[\mathbf{ϵ}])\mathit{I},$$

(29)

where the fourth-order tensor $\mathbb{E}[·]$ are defined as:

$$\mathbb{E}[\mathbf{ϵ}]:={\overline{c}}_1\mathbf{ϵ}+{\overline{c}}_2\mathrm{tr}(\mathbf{ϵ})\mathit{I},$$

(30)

and both the constants ${\overline{c}}_1$ and ${\overline{c}}_2$ depend on the linearized material parameters, that is, the conventional Young’s modulus $\left(E\right)$ and Poisson’s ratio $\left(\nu \right)$ , defined as:

$${\overline{c}}_1:=\frac{E}{1+\nu },{\overline{c}}_2:=\frac{\nu E}{(1+\nu )(1-2\nu )}.$$

(31)

Note that ${\overline{c}}_1$ and ${\overline{c}}_2$ are identical to linearized Lamé coefficients such that ${\overline{c}}_1=2\mu$ and ${\overline{c}}_2=\lambda$ . Unlike (22) and (23), material moduli are variants, thus starting from the linearized value of $E$ and $\nu$ , the nonlinear Lamé coefficients for the density-dependent model are:

$$\begin{array}{c}\multicolumn{1}{c}{\lambda :=\frac{\nu E}{(1+\nu )(1-2\nu )(1+\beta \mathrm{tr}(\mathbf{ϵ}))},}\\ \multicolumn{1}{c}{\mu :=\frac{E}{2(1+\nu )(1+\beta \mathrm{tr}(\mathbf{ϵ}))},}\end{array}$$

(32)

and the generalized (drained) bulk modulus $\left(K_{\mathrm{dr}}\right)$ is then derived from (31) and (32) as:

$$\begin{array}{c}\multicolumn{1}{c}{K_{\mathrm{dr}}:=\frac{1}{1+\beta \mathrm{tr}(\mathbf{ϵ})}({\overline{c}}_2+\frac{{\overline{c}}_1}{3})}\\ \multicolumn{1}{c}{=\frac{1}{1+\beta \mathrm{tr}(\mathbf{ϵ})}(\frac{\nu E}{(1+\nu )(1-2\nu )}+\frac{E}{3(1+\nu )}).}\end{array}$$

(33)

With $\mathrm{tr}(\mathbf{ϵ})$ , the density-dependent material moduli are set in (32) and (33). Note that we can also calculate the nonlinear Young’s modulus and Poisson’s ratio for this model from (32). For the preferential stiffness, we see the intensity with its direction can also be designed using the model parameter, $\beta$ .

Combining the balance of linear momentum (26) and the constitutive relationship (27), we obtain the following governing system of equations:

$$-\nabla ·\mathit{T}=0\mathrm{in}\mathrm{\Omega },$$

(34a)

$$\begin{array}{c}\multicolumn{1}{c}{\mathbf{ϵ}=\frac{(1+\nu )(1+\beta \mathrm{tr}(\mathbf{ϵ}))}{E}\mathit{T}}\\ \multicolumn{1}{c}{-\frac{\nu (1+\beta \mathrm{tr}(\mathbf{ϵ}))}{E}\mathrm{tr}(\mathit{T})\mathit{I}\mathrm{in}\mathrm{\Omega }.}\end{array}$$

(34b)

The above system of partial differential equations (PDEs) needs to be supplemented by appropriate boundary conditions with ${\mathrm{\Gamma }}_N$ and ${\mathrm{\Gamma }}_D$ :

$$\mathit{Tn}=\mathit{g},\mathrm{for}\mathrm{all}\mathit{x}\in {\mathrm{\Gamma }}_N,$$

(35a)

$$\mathit{u}={\mathit{u}}^0,\mathrm{for}\mathrm{all}\mathit{x}\in {\mathrm{\Gamma }}_D,$$

(35b)

where $\mathit{g}:\mathrm{\Omega }\to {\mathbb{R}}^2$ is the given traction and ${\mathit{u}}^0:\mathrm{\Omega }\to {\mathbb{R}}^2$ is the given boundary displacement data. Finally, we have the following formulation for the equilibrium problem for a two-dimensional elastic porous solid body:

Formulation 1. Given the material parameters, find $\mathit{u}=(u_1,u_2)$ , $\mathit{T}=\left\{{\mathit{T}}_{\mathrm{ij}}\right\}$ , and $\mathbf{ϵ}=\left\{{\mathbf{ϵ}}_{\mathrm{ij}}\right\},i,j=1,2$ such that

$$-{\mathit{T}}_{\mathrm{ij},j}(\mathit{u})=0,\mathrm{in}\mathrm{\Omega },\mathrm{and}i=1,2,$$

(36a)

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{ϵ}}_{\mathrm{ij}}(\mathit{u})=\frac{(1+\nu )(1+\beta \mathrm{tr}(\mathbf{ϵ}))}{E}{\mathit{T}}_{\mathrm{ij}}-\frac{\nu (1+\beta \mathrm{tr}(\mathbf{ϵ}))}{E}}\\ \multicolumn{1}{c}{\mathrm{tr}(\mathit{T}){\delta }_{\mathrm{ij}}\mathrm{in}\mathrm{\Omega },\mathrm{and}i=1,2,}\end{array}$$

(36b)

$${\mathit{u}}_i={\mathit{u}}_i^0,\mathrm{on}{\mathrm{\Gamma }}_D,\mathrm{and}i=1,2,$$

(36c)

$${\mathit{T}}_{\mathrm{ij}}(\mathit{u}){\mathit{n}}_j={\mathit{g}}_i,\mathrm{on}{\mathrm{\Gamma }}_N,\mathrm{and}i=1,2.$$

(36d)

The above boundary value problem can also be formulated in variational form as a problem of minimizing the total strain-energy density functional. This approach will be studied in future communication. As (36) is nonlinear, it is not tractable in the current form to any well-known analytical or numerical methods. In the following section, we propose Newton’s method for the linearization at the differential equation level, followed by the bilinear finite element method as a discretization technique.

Newton’s iterative method and variational formulation

To construct the numerical approximation to the boundary value problem (BVP) (36), we first linearize the differential equation to obtain a sequence of linear problems. Consider a mapping $\mathcal{L}(·):\mathrm{Sym}({\mathbb{R}}^{2\times 2})\mapsto \mathrm{Sym}({\mathbb{R}}^{2\times 2})$ , with $\mathcal{L}(0)=0$ , defined as:

$$\mathcal{L}(\mathbf{ϵ}(\mathit{u})):=\frac{\mathbb{E}[\mathbf{ϵ}(\mathit{u})]}{(1+\beta \mathrm{tr}(\mathbf{ϵ}))}.$$

(37)

The assumption of smoothness that we made for the Sobolev space (4) and the quasi-linearity of the operator in (37) bestow us to apply the Newton’s method to obtain the linearized version of the PDE model (34).

Given the initial guess ${\mathit{u}}^0\in {(C^2(\mathrm{\Omega }))}^2$ with $\mathbb{E}[\mathbf{ϵ}({\mathit{u}}^0)]\in \mathcal{C}\subseteq \mathrm{Sym}({\mathbb{R}}^{2\times 2})$ , for each $n\in {\mathbb{N}}_0$ denoting the iteration number, find ${\mathit{u}}^{n+1}\in {(C^2(\mathrm{\Omega }))}^2$ with $\mathbb{E}[\mathbf{ϵ}({\mathit{u}}^{n+1})]\in \mathcal{C}$ such that

$$D\mathcal{L}(\mathbf{ϵ}({\mathit{u}}^n))({\mathit{u}}^{n+1}-{\mathit{u}}^n)=-\mathcal{L}(\mathbf{ϵ}({\mathit{u}}^n)),$$

(38)

where $D\mathcal{L}(\mathbf{ϵ}(\mathit{u}))$ denotes the Fréchet derivative, defined as

$$D\mathcal{L}(\mathbf{ϵ}(\mathit{u})):=\underset{\xi \to 0}{lim}\frac{\mathcal{L}(\mathbf{ϵ}(\mathit{u})+\xi \mathbf{ϵ}(\mathit{v}))-\mathcal{L}(\mathbf{ϵ}(\mathit{u}))}{\xi },$$

(39)

for each $\mathit{v}\in {(C^2(\mathrm{\Omega }))}^2$ . Using (37) and (39), we obtain

$$D\mathcal{L}(\mathbf{ϵ}({\mathit{u}}^n))\delta {\mathit{u}}^n=[\frac{({\overline{c}}_1/2)(\nabla \delta {\mathit{u}}^n+{(\nabla \delta {\mathit{u}}^n)}^T)+{\overline{c}}_2(\nabla ·\delta {\mathit{u}}^n)\mathit{I}}{(1+\beta \mathrm{tr}(\mathbf{ϵ}({\mathit{u}}^n)))}-\frac{\beta \mathbb{E}[\mathbf{ϵ}({\mathit{u}}^n)](\nabla ·\delta {\mathit{u}}^n)}{{(1+\beta \mathrm{tr}(\mathbf{ϵ}({\mathit{u}}^n))}^2}].$$

(40)

Then, combining the equations (40), (34a), and (34b), we obtain the updated model equations:

$$\begin{array}{c}-\nabla ·[\frac{({\overline{c}}_1/2)(\nabla \delta {\mathit{u}}^n+{(\nabla \delta {\mathit{u}}^n)}^T)+{\overline{c}}_2(\nabla ·\delta {\mathit{u}}^n)\mathit{I}}{(1+\beta \mathrm{tr}(\mathbf{ϵ}({\mathit{u}}^n)))}-\frac{\beta \mathbb{E}[\mathbf{ϵ}({\mathit{u}}^n)](\nabla ·\delta {\mathit{u}}^n)}{{(1+\beta \mathrm{tr}(\mathbf{ϵ}({\mathit{u}}^n))}^2}]\hfill \\ =-\nabla ·\left[\frac{\mathbb{E}[\mathbf{ϵ}({\mathit{u}}^n)]}{(1+\beta \mathrm{tr}(\mathbf{ϵ}({\mathit{u}}^n)))}\right].\hfill \end{array}$$

(41a)

$$\mathbf{ϵ}(\mathit{u})=\frac{(1+\nu )(1+\beta \mathrm{tr}(\mathbf{ϵ}(\mathit{u})))}{E}\mathit{T}(\mathit{u})-\frac{\nu (1+\beta \mathrm{tr}(\mathbf{ϵ}(\mathit{u})))}{E}\mathrm{tr}(\mathit{T}(\mathit{u}))\mathit{I}.$$

(41b)

Then, the solution at the next iteration level is constructed by using

$${\mathit{u}}^{n+1}={\mathit{u}}^n+{\alpha }^n\delta {\mathit{u}}^n,$$

(42)

where ${\alpha }^n\in (0,1]$ is the “damping parameter” that controls the convergence rate. In the Formulation 2, we will address the procedure of computing the damping parameter with the line search algorithm. After the linearization, we have the formulation of the BVP at the continuous level.

Formulation 2. Given the linearized material parameters, Young’s modulus $E$ and Poisson’s ratio $\nu$ , initial guess ${\mathit{u}}^0\in {(C^2(\mathrm{\Omega }))}^2$ , find $\delta {\mathit{u}}^n$ for $n=0,1,\dots$ , $\mathit{T}=\left\{{\mathit{T}}_{\mathrm{ij}}\right\}$ , and $\mathbf{ϵ}=\left\{{\mathbf{ϵ}}_{\mathrm{ij}}\right\},i,j=1,2$ such that

$$\begin{array}{c}-\nabla ·[\frac{({\overline{c}}_1/2)(\nabla \delta {\mathit{u}}^n+{(\nabla \delta {\mathit{u}}^n)}^T)+{\overline{c}}_2(\nabla ·\delta {\mathit{u}}^n)\mathit{I}}{(1+\beta \mathrm{tr}(\mathbf{ϵ}({\mathit{u}}^n)))}-\frac{\beta \mathbb{E}[\mathbf{ϵ}({\mathit{u}}^n)](\nabla ·\delta {\mathit{u}}^n)}{{(1+\beta \mathrm{tr}(\mathbf{ϵ}({\mathit{u}}^n))}^2}]\hfill \\ =-\nabla [\mathcal{L}(\mathbf{ϵ}({\mathit{u}}^n))],\hfill \end{array}$$

(43a)

$${\mathit{u}}_i^n={\mathit{u}}_i^0\mathrm{on}{\mathrm{\Gamma }}_D,\mathrm{and}i=1,2,$$

(43b)

$${\mathit{T}}_{\mathrm{ij}}({\mathit{u}}^n){\mathit{n}}_j={\mathit{g}}_i\mathrm{on}{\mathrm{\Gamma }}_N,\mathrm{and}i=1,2,$$

(43c)

$${\mathit{u}}^{n+1}={\mathit{u}}^n+{\alpha }^n\delta {\mathit{u}}^n,$$

(43d)

$$\mathbf{ϵ}({\mathit{u}}^{n+1})=\frac{(1+\nu )(1+\beta \mathrm{tr}(\mathbf{ϵ}({\mathit{u}}^{n+1})))}{E}\mathit{T}({\mathit{u}}^{n+1})-\frac{\nu (1+\beta \mathrm{tr}(\mathbf{ϵ}({\mathit{u}}^{n+1})))}{E}\mathrm{tr}(\mathit{T}({\mathit{u}}^{n+1}))\mathit{I}.$$

(43e)

From equation (43d), note that it is clear that we need to impose the zero Dirichlet boundary conditions at each iteration level for Newton’s update $\delta {\mathit{u}}^n$ .

Continuous weak formulation

Here, we provide a variational formulation for the linearized version of the nonlinear BVP derived in the previous section. Also, the function spaces used to define the variational formulation have already been defined in the beginning, and the same setting is applied. To pose a weak formulation, we multiply the equations in the strong formulation (43a) with the test function from ${\hat{V}}_0$  as in (5), then via integrating by parts using Green’s formula together with the boundary conditions given in (43b) and (43c), we arrive at the following variational formulation.

Formulation 3. Given ${\mathit{u}}^0\in {\hat{V}}_{\mathit{g}}$ , for $n=0,1,2,\cdots$ , find ${\mathit{u}}^{n+1}:={\mathit{u}}^n+{\alpha }^n\delta {\mathit{u}}^n\in V$ , such that

$$A({\mathit{u}}^n;\delta {\mathit{u}}^n,\mathit{v})=L({\mathit{u}}^n;\mathit{v})\forall \mathit{v}\in {\hat{V}}_0,$$

(44)

where the bilinear term $A({\mathit{u}}^n;\delta {\mathit{u}}^n,\mathit{v})$ and the linear term $L({\mathit{u}}^n;\mathit{v})$ are given by

$$A({\mathit{u}}^n;\delta {\mathit{u}}^n,\mathit{v})={\int }_{\mathrm{\Omega }}\left[[\frac{({\overline{c}}_1/2)(\nabla \delta {\mathit{u}}^n+{(\nabla \delta {\mathit{u}}^n)}^T)+{\overline{c}}_2(\nabla ·\delta {\mathit{u}}^n)\mathit{I}}{(1+\beta \mathrm{tr}(\mathbf{ϵ}({\mathit{u}}^n)))}-\frac{\beta \mathbb{E}[\mathbf{ϵ}({\mathit{u}}^n)](\nabla ·\delta {\mathit{u}}^n)}{{(1+\beta \mathrm{tr}(\mathbf{ϵ}({\mathit{u}}^n))}^2}]\right]:\mathbf{ϵ}(\mathit{v})]d\mathit{x},$$

(45)

$$L({\mathit{u}}^n;\mathit{v})=-{\int }_{\mathrm{\Omega }}[\left[\frac{({\overline{c}}_1/2)(\nabla {\mathit{u}}^n+{(\nabla {\mathit{u}}^n)}^T)+{\overline{c}}_2(\nabla ·{\mathit{u}}^n)\mathit{I}}{(1+\beta \mathrm{tr}(\mathbf{ϵ}({\mathbf{u}}^n)))}\right]:\mathbf{ϵ}(\mathit{v})]d\mathit{x}+{\int }_{{\mathrm{\Gamma }}_N}\mathit{Tn}·\mathit{v}\mathrm{dx}.$$

Finite element discretization

Finally, we first recall some basic notions and structure of the classical finite element method (FEM) to discretize the weak formulation (44). The meshes used for all computations in the numerical examples presented in this paper are quadrilaterals. Let ${\left\{{\mathcal{T}}_h\right\}}_{h>0}$ be a conforming, shape-regular (in the sense of Ciarlet ⁵⁰ ) family of triangulation of the domain $\mathrm{\Omega }$ ; ${\mathcal{T}}_h$ is a finite family of sets such that $K\in {\mathcal{T}}_h$ which implies $K$ is an open simplex with the mesh size $h_K:=\mathrm{diam}(K)$ for each $K$ . Furthermore, we denote the largest diameter of the triangulation by

$$h:=\underset{K\in {\mathcal{T}}_h}{max}{h}_K.$$

For any $K_1,K_2\in {\mathcal{T}}_h$ , we have that ${\overline{K}}_1\cap {\overline{K}}_2$ is either a null set or a vertex or an edge or the whole of ${\overline{K}}_1$ and ${\overline{K}}_2$ , and $\bigcup _{K\in {\mathcal{T}}_h}\overline{K}=\overline{\mathrm{\Omega }}$ . We now define the classical piece-wise affine finite element space to approximate the displacement variable $\left(\mathit{u}\right)$ ,

$${\hat{V}}_h=\{{\mathbf{u}}_h\in {(C(\overline{\mathrm{\Omega }}))}^2:{{\mathbf{u}}_h|}_K\in {\mathbb{Q}}_k^2\forall K\in {\mathcal{T}}_h\},$$

(46)

where ${\mathbb{Q}}_k$ is a set containing the tensor-product of polynomials in two variables up to order $k$ on the reference cell $\hat{K}$ . Then, the discrete approximation space is:

$$V_h:={\hat{V}}_h\cap H^1(\mathrm{\Omega }).$$

(47)

The discrete counterpart of the continuous formulation (Formulation 3) is then as follows:

Formulation 4. Given ${\mathit{u}}_h^0\in V_h$ , and the $n^{\mathrm{th}}$ Newton’s iterative solution, that is, ${\mathit{u}}_h^n\in V_h$ , for $n=0,1,2,\cdots$ , find ${\mathit{u}}_h^{n+1}\in V_h$ , such that

$$A({\mathit{u}}^n;\delta {\mathit{u}}_h^n,{\mathit{v}}_h)=L({\mathit{u}}^n;{\mathit{v}}_h),\forall {\mathit{v}}_h\in {\hat{V}}_h,$$

(48)

where the linear and bilinear terms are given by:

$$A({\mathit{u}}_h^n;\delta {\mathit{u}}_h^n,{\mathit{v}}_h)={\int }_{\mathrm{\Omega }}\left[[\frac{({\overline{c}}_1/2)(\nabla \delta {\mathit{u}}_h^n+{(\nabla \delta {\mathit{u}}_h^n)}^T)+{\overline{c}}_2(\nabla ·\delta {\mathit{u}}_h^n)\mathit{I}}{(1+\beta \mathrm{tr}(\mathbf{ϵ}({\mathit{u}}_h^n)))}-\frac{\beta \mathbb{E}[\mathbf{ϵ}({\mathit{u}}_h^n)](\nabla ·\delta {\mathit{u}}_h^n)}{{(1+\beta \mathrm{tr}(\mathbf{ϵ}({\mathit{u}}_h^n))}^2}]\right]:\mathbf{ϵ}({\mathit{v}}_h)]d\mathit{x},$$

(49)

$$L({\mathit{u}}_h^n;{\mathit{v}}_h)=-{\int }_{\mathrm{\Omega }}[\left[\frac{({\overline{c}}_1/2)(\nabla {\mathit{u}}_h^n+{(\nabla {\mathit{u}}_h^n)}^T)+{\overline{c}}_2(\nabla ·{\mathit{u}}_h^n)\mathit{I}}{(1+\beta \mathrm{tr}(\mathbf{ϵ}({\mathit{u}}_h^n)))}\right]:\mathbf{ϵ}({\mathit{v}}_h)]d\mathit{x}+{\int }_{{\mathrm{\Gamma }}_N}(T_{12}{v}_1+T_{22}{v}_2)\mathrm{dx}.$$

(50)

where

$$\begin{array}{c}\multicolumn{1}{c}{\mathbb{E}[\mathbf{ϵ}({\mathit{u}}_h^n)]={\overline{c}}_1\mathrm{e}({\mathit{u}}_h^n)+{\overline{c}}_2\mathrm{tr}(\mathbf{ϵ}({\mathit{u}}_h^n))\mathit{I}}\\ \multicolumn{1}{c}{=({\overline{c}}_1/2)(\nabla {\mathit{u}}_h^n+{(\nabla {\mathit{u}}_h^n)}^{\mathrm{T}})+{\overline{c}}_2\nabla ·{\mathit{u}}_h^n,}\end{array}$$

(51)

and

$$\mathrm{tr}(\mathbf{ϵ}({\mathit{u}}_h^n)):=\nabla ·{\mathit{u}}_h^n.$$

(52)

The solution at the next iteration level is given by ${\mathit{u}}_h^{n+1}={\mathit{u}}_h^n+{\alpha }^n\delta {\mathit{u}}_h^n$ .

Note that we obtain an initial solution for the Newton’s iteration via solving the linear problem (i.e. with $\beta =0$ in (49)), which is an appropriate guess for the solution of a nonlinear problem.

Remark 3.1. We briefly outline the “line search method (algorithm)” used in the overall implementation of the FEM before we present the solution algorithm. Let us define the residual of the overall differential equation as:

$$\begin{array}{c}\multicolumn{1}{c}{f({\mathit{u}}^n,\mathit{v})=(\mathcal{L}(\mathrm{e}({\mathit{u}}^n)),\mathit{v})}\\ \multicolumn{1}{c}{={\int }_{\mathrm{\Omega }}\frac{\mathbb{E}[\mathbf{ϵ}({\mathit{u}}^n)]}{(1+\beta \mathrm{tr}(\mathbf{ϵ}({\mathit{u}}^n)))}:\nabla \mathit{v}d\mathit{x}.}\end{array}$$

(53)

The following algorithm ( Algorithm 1 ) summarizes the overall steps in the line search method. ⁶³

Algorithm 1: Description of backtrack line search algorithm

Input: $\overline{\alpha }>0$  ( $\overline{\alpha }=1$ guarantees the quadratic convergence), $c_1\in (0,1),r_1\mathrm{\&}{r}_2$ satisfying $0<r_1<r_2<1$ Output: ${\alpha }^n$ set $\alpha =\overline{\alpha }$ ; while $f({\mathit{u}}^n+\alpha {\mathit{p}}^n,\mathit{v})>f({\mathit{u}}^n,\mathit{v})+c_1\alpha \nabla f({\mathit{u}}^n,\mathit{v})·{\mathit{p}}^n$ , where ${\mathit{p}}^n$ is the update at $n^{\mathrm{th}}$ -level, replace $\alpha$ by a new value in $[r_1\alpha ,r_2\alpha ]$ return ${\alpha }^n=\alpha$

The new value for $\alpha$ is obtained by minimizing the overall residual given in (53). There are two possibilities for generating the iterates by appropriately choosing the search direction. One is to select relevant values of $r_1$ and $r_2$ so that the interval $[r_1\alpha ,r_2\alpha ]$ always contains the damping parameter $\alpha$ ; Another simple strategy is to iteratively replace $\overline{\alpha }$ with $\overline{\alpha }/2$ , until the value of $f({\mathit{u}}^n,\mathit{v})$ is sufficiently small.

Finally, the overall algorithm for the boundary value problem with the whole nonlinear density-dependent model is presented in Algorithm 2.

Algorithm 2: Algorithm for the nonlinear density-dependent model
Input: Choose the parameters: $\beta ,E,\nu$
Start with a sufficiently refined mesh;
while $[\mathrm{Iteration}\mathrm{Number}<\mathrm{Max}.\mathrm{Number}\mathrm{of}\mathrm{Iterations}]$ .AND. $[\mathrm{Residual}>\mathrm{Tol}.]$ do
Assemble the equations (49) and (50) for the displacements as primary variables using ${\mathbb{Q}}_1$ shape functions;
Use a direct solver to solve for $\delta {\mathit{u}}^n$
Construct the solution at $n$ -th iteration by a direct solver and then update the solution variable using ${\mathit{u}}_h^{n+1}={\mathit{u}}_h^n+\overline{\alpha }\delta {\mathbf{u}}_h^n$
While $[\mathrm{Line}\mathrm{Search}\mathrm{Iteration}\mathrm{Number}<\mathrm{Max}.\mathrm{Number}\mathrm{of}\mathrm{Line}\mathrm{Search}\mathrm{Iterations}]$  do
Calculate Residual using equation (53);
if $\mathrm{Residual}\le \mathrm{Tol}.$  then
Break
end
Do Algorithm 1 to find the optimal value of ${\alpha }^n$
end
end
Write the final converged solution ${\mathit{u}}_h^{n+1}$ to output files for post-processing;
Compute the crack-tip fields (e.g., stress and strain) for visualization

Numerical experiments and discussion

The primary purpose of this section lies in verifying the constitutive relation of density-dependent material moduli for an elastic porous solid via modeling and computational approaches provided in the previous sections. Using the proposed and classical models, we present several numerical tests to compare the stress and strain distributions in the computational domains under different mechanical loadings. Distributions of preferential stiffness from the density-dependent model can be confirmed particularly near the crack-tip preferentially. Ultimately, our goal is to suggest the rationale for this nonlinear model, from which physically meaningful and reliable responses of the damaged pores or the crack-tip can be described under the mechanical loading.

All numerical implementations are done using the deal.II, ⁶⁴ finite element library. Against the nonlinearity, the Newton’s method is employed with the lower bound of the convergence (i.e. the tolerance) as ${10}^{-8}$ and the maximum number of Newton iteration steps is set as $50$ . For the line search algorithm, a constant damping coefficient for the line search is taken as $0.5$ , that is, iteratively replacing $\overline{\alpha }$ with $\overline{\alpha }/2$ , and the maximum number of steps is taken as $10$ (see Algorithms 1 and 2). A direct solver is used as a linear solver to compute the numerical solution of the linearized system of equations.

$h$ -convergence study

In this section, we introduce a sample problem to verify the mathematical description and algorithm of the numerical model proposed in this study. To this end, $h$ -convergence test is performed with a manufactured solution set as

$$\mathit{u}:=(\sin \left(\frac{\pi }{2}x\right),-\cos \left(\frac{\pi }{2}y\right)),$$

for an unit square domain (1 m × 1 m) as Figure 1. The Dirichlet condition that satisfies the exact solution is applied for all the boundaries (from ${\mathrm{\Gamma }}_{D_1}$ to ${\mathrm{\Gamma }}_{D_4}$ ), resulting in ${\mathrm{\Gamma }}_N=Ø$ . For the modeling parameter, $\beta =1.0$ is taken, and for the linearized material moduli, Young’s modulus of $E=100$  Pa and Poisson’s ratio of $\nu =0.1$ are used. We then globally refine the whole domain with a total six cycles. In the sense of $L^2$ -norm, the optimal convergence order of $2$ for the bilinear polynomial is obtained. See the detailed rate values for convergence in Table 1.

Figure 1.

An unit square domain and the Dirichlet boundary conditions for h-convergence study.

Table 1.

The Results of $L^2$ Error Demonstrate the Rate of Optimal Convergence.

Cycle of refinement	$h$	$L^2$ error	Rate
1	0.5	0.053043115884	-
2	0.25	0.013140372001	2.7317
3	0.125	0.003275513882	2.3635
4	0.0625	0.000818290922	2.1809
5	0.03125	0.000204592608	2.0899
6	0.015625	0.000051207706	2.0433

Boundary value problems with mechanical loadings

Premise and setup for numerical experiments

In this study, the elastic regime is of interest, thus we do not consider any dissipation of energy. We also do not consider any fluid saturation in the solid, thus the porous solid of our interest is in the unsaturated condition. Neither is any explicit porosity concept in the model nor any initial porosity value that is directly measurable through experiments, or calculated via formula. However, the implicit concept of porosity is still considered with the volumetric strain, $\mathrm{tr}(\mathbf{ϵ})$ , such as in (14) that relates the density of the material to the mechanical moduli ((32) or (33)).

From Example 1 to Example 4, we assume homogeneous isotropic material under isothermal conditions. For the same unit square domain in these examples (1 m × 1 m), a total of seven global refinements are performed, resulting in the uniformly refined mesh size of $h=0.0078125$ . For its initial linearized elastic moduli, Young’s modulus $\left(E\right)$ of $100$  MPa and Poisson’s ratio $\left(\nu \right)$ of $0.15$ are taken. For the Neumann boundary condition, the same traction value of $f_u=0.01$  MPa is applied to the top boundary (see Figures 2 and 5) with different modes of loading. The rest detailed boundary conditions for each example are described in each subsection. The sign convention for stress follows such that the tensile stress is positive. As the study focuses on the nonlinear effects and different mechanical responses from the parameter of $\beta$ -values, we compare four different cases, that is, $\beta =-200,-50,+50,$ and $+200$ , with the case of $\beta =0$ , that is, the classical linearized elasticity. Comparisons are then highlighted with stress and strain distributions, and strain density energy for each case, focusing on the area near the crack-tip. Furthermore, we also compare the stress intensity factor, volumetric strain, and bulk modulus between the cases.

Figure 2.

Numerical domains (unit square, 1 m × 1 m) for Example 1 with traction boundary conditions in different modes: (a) Example 1a in mode-I and (b) Example 1b in mode-II. The red-dotted line is the reference line for each example.

Example 1: No crack problems

In this example, no crack or slit exists but an intact porous solid is considered. Note that two different traction loads are applied to the top boundary $\left({\mathrm{\Gamma }}_3\right)$ , that is, the tensile loading in mode-I (Example 1a) and in-plane shear loading in mode-II (Example 1b) in these domains of ${\mathrm{\Omega }}_{1a}$ (Figure 2(a)) and ${\mathrm{\Omega }}_{1b}$ (Figure 2(b)), respectively. The detailed boundary conditions of each problem for the field variables are as follows: Example 1a has

$$\mathit{Tn}=0\mathrm{on}{\mathrm{\Gamma }}_2,{\mathrm{\Gamma }}_4,$$

(54a)

$$\mathit{Tn}·\mathit{v}=T_{21}{v}_1+T_{22}{v}_2=f_u{v}_2\mathrm{on}{\mathrm{\Gamma }}_3,$$

(54b)

$$u_2=0\mathrm{on}{\mathrm{\Gamma }}_1,$$

(54c)

while Example 1b has

$$\mathit{Tn}=0\mathrm{on}{\mathrm{\Gamma }}_2,{\mathrm{\Gamma }}_4,$$

(55a)

$$\mathit{Tn}·\mathit{v}=T_{21}{v}_1+T_{22}{v}_2=f_u{v}_1\mathrm{on}{\mathrm{\Gamma }}_3,$$

(55b)

$$u_1=0\mathrm{and}{u}_2=0\mathrm{on}{\mathrm{\Gamma }}_1.$$

(55c)

Note that we have slightly different Dirichlet boundary conditions for the bottom boundary ( ${\mathrm{\Gamma }}_1$ ) for the two problems; the roller boundary conditions are set for Example 1a, thus there is no displacement to its normal direction, that is, $y$ -direction, but it is free in $x$ -direction. Meanwhile, hinge boundary conditions are set for Example 1b, thus there are no displacements in both $x$ - and $y$ -directions. For these problems, we compare the displacements on the reference line (red-dotted line, $L=1$ m) in the $x$ - and $y$ -directions, that is, parallel and perpendicular to the reference line, respectively.

Figures 3 and 4 demonstrate the displacements for Example 1a and Example 1b, respectively. In each figure, the solutions of $x$ -displacement on the left and of $y$ -displacement on the right are exhibited, where we take absolute values with the semi-log scale in the direction perpendicular to the loading, that is, $x$ -direction for Example 1a and $y$ -direction for Example 1b. In the perpendicular direction to the loading, we find that the nonlinear model with different $\beta$ -values yields distinct responses near the center. In the parallel direction to the loading for mode-I and mode-II, we find that positive and negative $\beta$ -values are reflected in the displacement in the opposite way. For example, positive $\beta$ in tension (Figure 3) yields larger $y$ -displacement than that of the linear, while negative $\beta$ has a smaller one. Meanwhile, we can find that positive $\beta$ in shear (Figure 4) yields less rotation angle from the center ( $x=0.5$ ) from the $x$ -displacement. Thus, it implies that the strengths of a material with tensile and shear (in compressive direction) stresses are different under the nonlinear model.

Figure 3.

Displacements (unit: m) on the reference line ( $L=1$  m) for Example 1a: (left) x-displacement in absolute value and (right) y-displacement.

Figure 4.

Displacements (unit: m) on the reference line ( $L=1$  m) for Example 1b: (left) x-displacement and (right) y-displacement in absolute value.

To investigate the preferential stiffness with these positive and negative $\beta$ -values, we focus on the stress and strain distributions depending on the volumetric strain changes, particularly near an edge crack in the following examples; Henceforth, the same numerical domain is addressed, having different boundary conditions depending on the mode of loading in each example but with the same right-edge crack. The crack is expressed with the slit ( ${\mathrm{\Gamma }}_C$ , the blue lines in Figure 5(a)–(c)), where zero traction is applied. The traction value ( $f_u=0.01$  MPa) and other mechanical properties are the same as Example 1. For these examples, the main comparisons lie in the stress and strain distributions including the stress intensity factor on the reference line (the red-dotted lines in Figure 5, $L=0.5$  m), where the volumetric strain and bulk modulus are also compared between the models. These post-processing works for the variables are based on the average values in each grid element of interest using quadrature points inside.

Figure 5.

Numerical domains (unit square, 1 m × 1 m) with cracks for Examples 2–4 with traction boundary conditions in different modes: (a) Example 2 in mode-I, (b) Example 3 in mode-II, and (c) Example 4 in mixed-mode. Crack ( ${\mathrm{\Gamma }}_C$ ) in blue in each domain has the same geometry with the length of 0.5 m and the red-dotted line is the reference line for each example.

Example 2: Tensile loading with crack

For Example 2, a tensile loading in pure mode-I is applied to a numerical domain with a right-edge crack as illustrated in Figure 5(a), which has the following boundary conditions:

$$\mathit{Tn}=0\mathrm{on}{\mathrm{\Gamma }}_2,{\mathrm{\Gamma }}_4,{\mathrm{\Gamma }}_C,$$

(56a)

$$\mathit{Tn}·\mathit{v}=T_{21}{v}_1+T_{22}{v}_2=f_u{v}_2\mathrm{on}{\mathrm{\Gamma }}_3,$$

(56b)

$$u_2=0\mathrm{on}{\mathrm{\Gamma }}_1.$$

(56c)

The roller boundary conditions are applied to ${\mathrm{\Gamma }}_1$ , resulting in the homogeneous Dirichlet boundary in $y$ -direction, that is, $u_y=0$ , without a constraint in $x$ -direction.

Stress and strain

In Figure 6(a)–(c), $T_{22}$ with $\beta =-200$ , $\beta =0$ , that is, the linear model, and $\beta =+200$ , respectively, are illustrated, and the corresponding ${\epsilon }_{22}$ are shown in (d)–(f). For (g)–(i) in Figure 6, the strain energy density ( $\mathrm{SED}$ ) with the unit of $N/m/m$ or Pa illustrated corresponding to each case. For simplicity, the density-dependent nonlinear model with $\beta =+50$ and $\beta =-50$ results are omitted for now. Note that $T_{22}$ and ${\epsilon }_{22}$ are parallel to the mode-I loading (normal to the horizontal line). It is found that the larger the value a variable has, the narrower is its distribution with more localization. For the tensile stress $\left(T_{22}\right)$ of which sign convention is positive, we identify that $\beta =+200$ yields the smallest with the largest of ${\epsilon }_{22}$ , which implies the material strength becomes weaker against the tension. The smallest tensile strain is obtained for $\beta =-200$ . From the $\mathrm{SED}$ results, we confirm that the maximum strain energy density for each case occurs right in front of the crack-tip. We also identify that $\mathrm{SED}$ for the case of $\beta =+200$ shows the largest with the maximum positive axial strain; thus, ${\epsilon }_{22}$ obtained from the positive $\beta$ (e.g. $\beta =+200$ ) is larger than that from the linear and even larger than that from the negative $\beta$ (e.g. $\beta =-200$ ) near the tip. The detailed max and min values for each case are presented in Table 2.

Figure 6.

[Example 2] (a)–(c) stress $T_{22}$ (unit: Pa), (d–f) strain ${\mathbf{ϵ}}_{22}$ , and (g–i) strain energy density ( $\mathrm{SED}$ , unit: Pa) distributions: (left) with $\beta =-200$ , (middle) with $\beta =0$ (i.e. the linear model), and (right) with $\beta =+200$ : (a) $T_{22}\mathrm{with}\beta =-200$ , (b) $T_{22}\mathrm{with}\beta =0(\mathrm{Linear})$ , (c) $T_{22}\mathrm{with}\beta =+200$ , (d) ${\mathbf{ϵ}}_{22}\mathrm{with}\beta =-200$ , (e) ${\mathbf{ϵ}}_{22}\mathrm{with}\beta =0(\mathrm{Linear})$ , (f) ${\mathbf{ϵ}}_{22}\mathrm{with}\beta =+200$ , (g) $\mathrm{SED}\mathrm{with}\beta =-200$ , (h) $\mathrm{SED}\mathrm{with}\beta =0(\mathrm{Linear})$ , and (i) $\mathrm{SED}\mathrm{with}\beta =+200$ .

Table 2.

[Example 2] The maximum and minimum values of the variables in Figure 6 for each case using the linear ( $\beta =0$ ) and nonlinear models with $\beta =-200$ and $\beta =+200$ .

Variable	$\beta =-200$		$\beta =0$ (Linear)		$\beta =+200$
	Max	Min	Max	Min	Max	Min
$T_{22}$ [MPa]	$0.25$	$-0.022$	$0.21$	$-0.025$	$0.17$	$-0.028$
${\mathbf{ϵ}}_{22}$ [-]	$0.0014$	$-0.00023$	$0.0019$	$-0.00025$	$0.0028$	$-0.00027$
$\mathrm{SED}$ [Pa]	290	$3.5\times {10}^{-7}$	330	$3.8\times {10}^{-7}$	370	$4.2\times {10}^{-7}$

We also investigate with the stress intensity factor (SIF) calculated on the reference line in Figure 5(a). The calculation of SIF for the mode-I ( ${\mathrm{K}}_I$ ) is based on the following equation:

$${\mathrm{K}}_I:=\underset{r\to 0}{lim}\sqrt{2\pi r}{T}_{22}(r,\theta =0),$$

(57)

where $r=0.5$ (equivalently, x/L  $=0$ in Figure 7) located on ${\mathrm{\Gamma }}_4$ in Figure 5(a) and $r=0$ (equivalently, x/L  $=1$ in Figure 7) at the tip of crack contacting the reference line. Meanwhile, the stress intensity factor in mode-II loading, that is, ${\mathrm{K}}_{\mathrm{II}}$ is based on the following:

$${\mathrm{K}}_{\mathrm{II}}:=\underset{r\to 0}{lim}\sqrt{2\pi r}{T}_{21}(r,\theta =0),$$

(58)

with the same location for $r$ . In each example henceforth, we plot these two SIFs together, even though the primary mode of loading is different for each example. Note that the negative value implies the compression due to various conditions such as the boundary condition and the Poisson’s ratio. From Figure 7 (left), we confirm that its maximum ${\mathrm{K}}_I$ occurs right in front of the tip for each case. Although the difference between the cases is not much discernible, the calculated SIF of ${\mathrm{K}}_I$ for $\beta =-200$ has the largest of at the tip in the same context as the stress distributions in Figure 6. Interestingly, the case with $\beta =+200$ has the smallest value for the negative region of ${\mathrm{K}}_I$ against the compression, which is analogous to the hardening behavior. In addition, through ${\mathrm{K}}_{\mathrm{II}}$ (Figure 7 (right)), the shear stress with compressive direction in front of the tip is identified with its negative values for each case. It is also found that $\beta =-200$ has the smallest value in the calculated SIF of ${\mathrm{K}}_{\mathrm{II}}$ . Thus from ${\mathrm{K}}_I$ and ${\mathrm{K}}_{\mathrm{II}}$ , the preferential stiffness can be determined based on the parameter, that is, positive and negative $\beta$ -values of the density-dependent material moduli model.

Figure 7.

[Example 2] Stress intensity factor (unit: Pa ${\mathrm{m}}^{1/2}$ ) in mode-I (left, ${\mathrm{K}}_I$ ) and in mode-II (right, ${\mathrm{K}}_{\mathrm{II}}$ ) on the reference line with $L=0.5$  m.

Volumetric strain and bulk modulus

Here, we illustrate the change of the drained bulk modulus $\left(K_{\mathrm{dr}}\right)$ in (33) and the volumetric strain $(\mathrm{tr}(\mathbf{ϵ}))$ on the reference line. Note that the bulk modulus is the inverse of the compressibility of the skeleton of the porous solid; the greater it is, less compressible it is. Figure 8 demonstrates $K_{\mathrm{dr}}$ on the left, $\mathrm{tr}(\mathbf{ϵ})$ in the middle, and their relation on the right. As previously figured in ${\mathrm{K}}_I$ (Figure 7 (left)), we see that a slight compression occurs for each case from the plot for $\mathrm{tr}(\mathrm{e})$ shown in the $x$ -axis upto around x/L  $=0.25$ (Figure 8 (middle)). From their relation plot (Figure 8 (right)), we see that the mechanical property (i.e. $K_{\mathrm{dr}}$ here) and the implicit porosity (i.e. $\mathrm{tr}(\mathbf{ϵ})$ here) is in its reverse relation for the positive $\beta$ model with some nonlinearity. In line with the previous finding that the material property of the negative $\beta$ -value for the density-dependent model becomes stiffer against the dilation from tensile loading, we find from Figure 8 (middle) that the volumetric strains become notably differentiable approaching the vicinity of the tip based on the modeling parameter; $\beta =+200$ has the largest for this pure mode-I loading, while $\beta =-200$ has the smallest. As compared to the linear model, it is phenomenologically similar to the strain hardening or softening in the elastoplasticity regime, although the mechanical response regime remains under the pure elasticity.

Figure 8.

[Example 2] Bulk modulus ( $K_{\mathrm{dr}}$ , unit: Pa) and volumetric strain (tr $\left(\mathbf{ϵ}\right)$ on the reference line ( $L=0.5$  m): (left) $K_{\mathrm{dr}}$ , (middle) tr $\left(\mathbf{ϵ}\right)$ , and (right) $K_{\mathrm{dr}}$ versus tr $\left(\mathbf{ϵ}\right)$ .

Example 3: In-plane shear loading with crack

For Example 3, we have the pure mode-II loading (see Figure 5(b)), the volumetric strain changes under the shear force are addressed. The problem has the boundary conditions as follows:

$$\mathit{Tn}=0\mathrm{on}{\mathrm{\Gamma }}_2,{\mathrm{\Gamma }}_4,{\mathrm{\Gamma }}_C,$$

(59a)

$$\mathit{Tn}·\mathit{v}=T_{21}{v}_1+T_{22}{v}_2=f_u{v}_1\mathrm{on}{\mathrm{\Gamma }}_3,$$

(59b)

$$u_1=0\mathrm{and}{u}_2=0\mathrm{on}{\mathrm{\Gamma }}_1.$$

(59c)

Stress and strain

For this in-plane shear loading, we plot $T_{21}$ , ${\epsilon }_{21}$ , and $\mathrm{SED}$ in Figure 9 for the same $\beta$ -values: $\beta =-200$ , $\beta =0$ (the linear), and $\beta =+200$ , in the same arrangement for the subplots. Due to the in-plane shear loading, $T_{21}$ ((a)–(c) in Figure 9) and ${\epsilon }_{21}$ ((d)–(f) in Figure 9) have their distributions of positive and negative values across, and this traverse pattern exists particularly around the tip. Note that the positive shear strain in the first quadrant of $x$ - $y$ plane implies a decrease of the right angle. Unlike Example 2, the smallest positive $T_{21}$ is for the $\beta =-200$ case, while the largest belongs to the $\beta =+200$ case. We can see more clear distinctions between the cases in the distribution of ${\epsilon }_{21}$ . As positive and negative shear strains intersect around the tip, compressive stress induces the volume decrease and porosity compaction. In the same context, the case of $\beta =+200$ exhibits much smaller positive values of ${\epsilon }_{21}$ compared to the case of $\beta =-200$ , which implies that the positive $\beta$ case has more resistance against the compression (or decreasing of the right angle) along with the sliding force applied. From the $\mathrm{SED}$ results, we also confirm that the maximum strain energy density for each case occurs right in front of the crack-tip. Likewise, the detailed max min values for each case are found in Table 3.

Figure 9.

[Example 3] (a–c) stress $T_{21}$ (unit: Pa), (d–f) strain ${\mathbf{ϵ}}_{21}$ , and (g–i) strain energy density ( $\mathrm{SED}$ , unit: Pa) distributions: (left) with $\beta =-200$ , (middle) with $\beta =0$ (i.e. the linear model), and (right) with $\beta =+200$ : (a) $T_{21}\mathrm{with}\beta =-200$ , (b) $T_{21}\mathrm{with}\beta =0(\mathrm{Linear})$ , (c) $T_{21}\mathrm{with}\beta =+200$ , (d) ${\mathbf{ϵ}}_{21}\mathrm{with}\beta =-200$ , (e) ${\mathbf{ϵ}}_{21}\mathrm{with}\beta =0(\mathrm{Linear})$ , (f) ${\mathbf{ϵ}}_{21}\mathrm{with}\beta =+200$ , (g) $\mathrm{SED}\mathrm{with}\beta =-200$ , (h) $\mathrm{SED}\mathrm{with}\beta =0(\mathrm{Linear})$ , and (i) $\mathrm{SED}\mathrm{with}\beta =+200$ .

Table 3.

[Example 3] The maximum and minimum values of the variables in Figure 9 for each case using the linear $(\beta =0)$ and nonlinear models with $\beta =-200$ and $\beta =+200$ .

Variable	$\beta =-200$		$\beta =0$ (Linear)		$\beta =+200$
	Max	Min	Max	Min	Max	Min
$T_{21}$ [MPa]	$0.12$	$-0.046$	$0.17$	$-0.062$	$0.20$	$-0.071$
${\mathbf{ϵ}}_{21}$ [-]	$0.0037$	$-0.00072$	$0.0019$	$-0.00072$	$0.0013$	$-0.00064$
$\mathrm{SED}$ [Pa]	1500	$8.5\times {10}^{-7}$	1100	$9.5\times {10}^{-7}$	1500	$1.1\times {10}^{-6}$

Two SIFs on the reference line (Figure 5(b)) are illustrated in Figure 10 using (57) and (58). We find that at about x/L  $=0.6$ in the reference line, ${\mathrm{K}}_{\mathrm{II}}$ for the case of $\beta =+200$ surpasses the rest cases. In parallel, the case of $\beta =+200$ has the smallest concentration for ${\mathrm{K}}_I$ (Figure 10 (left)) in front of the tip with the negative tensile stress, that is, under compression. Therefore for the in-plane shear loading, we find the model with negative $\beta$ -values has the weaker material strength against the shear or compression.

Figure 10.

[Example 3] Stress intensity factor (unit: Pa ${\mathrm{m}}^{1/2}$ ) in mode-I (left, ${\mathrm{K}}_I$ ) and in mode-II (right, ${\mathrm{K}}_{\mathrm{II}}$ ) on the reference line with $L=0.5$  m.

Volumetric strain and bulk modulus

In Figure 11, the bulk modulus and volumetric strain for each case are plotted. We confirm that the stiffness of each case is illustrated in the opposite compared to Example 2 (see Figure 8), as the volumetric strains are plotted in the reverse direction. Until about the half of reference line (x/L  $=0.5$ ) (see Figure 11 (middle)), positive tr( $\mathbf{ϵ}$ ), that is, the dilation of implicit porosity, is decreased to zero in each case, which can also be figured in Figure 11 (left) with $K_{\mathrm{dr}}$ . Note that more variation of tr( $\mathbf{ϵ}$ ) with a larger range is observed for the case of $\beta =-200$ , particularly at the tip, thus we confirm that the density-dependent model with negative $\beta$ is relatively weak in the shear (or compressive) loading.

Figure 11.

[Example 3] Bulk modulus ( $K_{\mathrm{dr}}$ , unit: Pa) and volumetric strain (tr( $\mathbf{ϵ}$ )) on the reference line (L = 0.5 m): (left) $K_{\mathrm{dr}}$ , (middle) tr( $\mathbf{ϵ}$ ), and (right) $K_{\mathrm{dr}}$ versus tr( $\mathbf{ϵ}$ ).

Example 4: mixed-mode loading with crack

In this last example, the mixed-mode (i.e. the mode-I and II) of loading is applied. The boundary conditions are as follows:

$$\mathit{Tn}=0\mathrm{on}{\mathrm{\Gamma }}_2,{\mathrm{\Gamma }}_4,{\mathrm{\Gamma }}_C,$$

(60a)

$$\mathit{Tn}·\mathbf{v}=T_{21}{v}_1+T_{22}{v}_2=f_u{v}_1+f_u{v}_2\mathrm{on}{\mathrm{\Gamma }}_3,$$

(60b)

$$u_1=0\mathrm{and}{u}_2=0\mathrm{on}{\mathrm{\Gamma }}_1.$$

(60c)

Note that we have the same bottom boundary $\left({\mathrm{\Gamma }}_1\right)$ condition as Example 3 (see Figure 5(c)) considering the mode-II loading.

Stress and strain

We plot both $T_{22}$ with ${\epsilon }_{22}$ and $T_{21}$ with ${\epsilon }_{21}$ in Figures 12 and 13, respectively. For $T_{22}$ , unlike the pure mode-I loading, we find compressive stress near the tip, which is due to the in-plane shear or mode-II loading imposed simultaneously. As consistent with the fact that the positive $\beta$ case has more resistance against the compression and less against the tension, the smallest negative ${\mathbf{ϵ}}_{22}$ is obtained for the case of $\beta =+200$ as seen in Figure 12. In addition, we figure the smallest positive ${\mathbf{ϵ}}_{22}$ is for $\beta =-200$ as the negative $\beta$ case is relatively stiffer against the tensile loading. The same pattern of distribution is also found; the larger a variable has for the value, the narrower its distribution with more localization. As the positive $\beta$ has more resistance against the compression, compressive stress of the case of $\beta =+200$ is more concentrated with larger value (colored in blue in Figure 12(c)). While for strain ( ${\mathbf{ϵ}}_{22}$ ), the case of $\beta =-200$ has wider tensile strain region with smaller positive values (colored in red in Figure 12(a)), that is, more resistance. About $T_{21}$ with ${\mathbf{ϵ}}_{21}$ in Figure 13, similar patterns to Example 3 with the pure mode-II loading are shown; the positive $\beta$ with the case of $\beta =+200$ is stiffer against the compression with the smallest positive value in ${\mathbf{ϵ}}_{21}$ and the largest negative ${\mathbf{ϵ}}_{21}$ . We see that the maximum strain energy density $\left(\mathrm{SED}\right)$ is shown in front of the tip for each case (Figure 13). See the detailed max min values for each case in Table 4.

Figure 12.

[Example 4] (a–c) stress $T_{22}$ (unit: Pa) and (d–f) strain ${\mathbf{ϵ}}_{22}$ distributions: (left) with $\beta =-200$ , (middle) with $\beta =0$ (i.e. the, $\beta =+200$ : (a) $T_{22}\mathrm{with}\beta =-200$ , (b) $T_{22}\mathrm{with}\beta =0(\mathrm{Linear})$ , (c) $T_{22}\mathrm{with}\beta =+200$ , (d) ${\mathbf{ϵ}}_{22}\mathrm{with}\beta =-200$ , (e) ${\mathbf{ϵ}}_{22}\mathrm{with}\beta =0(\mathrm{Linear})$ , and (f) ${\mathbf{ϵ}}_{22}\mathrm{with}\beta =+200$ .

Figure 13.

[Example 4] ((a–c) stress $T_{21}$ (unit: Pa), (d–f) strain ${\mathbf{ϵ}}_{21}$ , and (g–i) strain energy density ( $\mathrm{SED}$ , unit: Pa) distributions: (left) with $\beta =-200$ , (middle) with $\beta =0$ (i.e. the linear model), and (right) with $\beta =+200$ : (a) $T_{21}\mathrm{with}\beta =-200$ , (b) $T_{21}\mathrm{with}\beta =0(\mathrm{Linear})$ , (c) $T_{21}\mathrm{with}\beta =+200$ , (d) ${\mathbf{ϵ}}_{21}\mathrm{with}\beta =-200$ , (e) ${\mathbf{ϵ}}_{21}\mathrm{with}\beta =0(\mathrm{Linear})$ , (f) ${\mathbf{ϵ}}_{21}\mathrm{with}\beta =+200$ , (g) $\mathrm{SED}\mathrm{with}\beta =-200$ , (h) $\mathrm{SED}\mathrm{with}\beta =0(\mathrm{Linear})$ , and (i) $\mathrm{SED}\mathrm{with}\beta =+200$ .

Table 4.

[Example 4] The maximum and minimum values of the variables in Figures 12 and 13 for each case using the linear $(\beta =0)$ and nonlinear models with $\beta =-200$ and $\beta =+200$ .

Variable	$\beta =-200$		$\beta =0$ (Linear)		$\beta =+200$
	Max	Min	Max	Min	Max	Min
$T_{22}$ [MPa]	$0.14$	$-0.10$	$0.13$	$-0.12$	$0.12$	$-0.13$
${\mathbf{ϵ}}_{22}$ [-]	$0.0011$	$-0.0013$	$0.0013$	$-0.00098$	$0.0014$	$-0.00087$
$T_{21}$ [MPa]	$0.073$	$-0.0041$	$0.092$	$-0.006$	$0.12$	$-0.0077$
${\mathbf{ϵ}}_{21}$ [-]	$0.0015$	$-0.000049$	$0.0011$	$-0.000069$	$0.00078$	$-0.000085$
$\mathrm{SED}$ [Pa]	470	$6.4\times {10}^{-7}$	410	$6.7\times {10}^{-7}$	370	$7.0\times {10}^{-7}$

Using (57) and (58), the SIF of each mode on the reference line is presented in Figure 14. Compared to Example 2 in the pure mode-I loading, the maximum values of ${\mathrm{K}}_I$ of all cases are on the left end boundary with much greater values, which is due to the aforementioned bottom boundary condition, that is, the hinge at the bottom is the same as the one in Example 3. Approaching the tip, ${\mathrm{K}}_I$ for each case is decreased to the negative value, and we confirm the compressive stresses as seen in Figure 12. At the tip, the case of $\beta =+200$ has the minimum as it is stiffer against the compression. Meanwhile for ${\mathrm{K}}_{\mathrm{II}}$ , the case of $\beta =+200$ has its maximum at the tip, whereas the case of $\beta =-200$ shows its maximum about the point of x/L  $=0.5$ . ${\mathrm{K}}_{\mathrm{II}}$ of the case of $\beta =+200$ is found to surpass the rest cases near x/L  $=0.7$ on the reference line, confirming compressive stress applied more approaching the tip.

Figure 14.

[Example 4] Stress intensity factor (unit: Pa ${\mathrm{m}}^{1/2}$ ) in mode-I (left, ${\mathrm{K}}_I$ ) and in mode-II (right, ${\mathrm{K}}_{\mathrm{II}}$ ) on the reference line with $L=0.5$  m.

Volumetric strain and bulk modulus

From Figure 15 (middle) and (right), we identify that the volumetric strains for all cases are about the same. Until about the point of x/L  $=0.6$ on the reference line, positive tr( $\mathbf{ϵ}$ ) implying the dilation of porosity is shown in each case. With about the same change in volumetric strain (or implicit porosity) for each case on the reference line, we note preferentially different and distinctive mechanical responses, thus depending solely on the parameter with $\beta$ -values (see (32) and (33)).

Figure 15.

[Example 4] Bulk modulus ( $K_{\mathrm{dr}}$ , unit: Pa) and volumetric strain (tr( $\mathbf{ϵ}$ )) on the reference line ( $L=0.5$  m): (left) $K_{\mathrm{dr}}$ , (middle) tr( $\mathbf{ϵ}$ ), and (right) $K_{\mathrm{dr}}$ versus tr( $\mathbf{ϵ}$ ).

Conclusion

The purpose of this paper is to study mechanical responses of an elastic porous solid with preferential stiffness of which material moduli are dependent upon the density, and to provide a stable finite element solution of stress and strain fields, particularly around the crack-tip. The proposed model based on the same linearization as the linearized elasticity, that is, the gradient of the displacement is infinitesimal, cannot stem from the conventional framework of Cauchy elasticity. It is structured on a special constitutive relation (27) for elastic porous solids that show significantly distinctive responses from those using classical linearized models. For our numerical approaches, we employ a universal and computationally efficient method of Newton’s method and FEM under the framework developed in Yoon et al. ⁴³ and Lee et al. ⁴⁴ to overcome the nonlinearity of the model with partial differential equation. The proposed algorithm in the study is verified for the optimal convergence rate using the method of manufactured solution. Three different types of loading are considered in this work. Even though the constitutive relation studied in this paper is simplified via a single parameter describing the mechanical response of the material under scrutiny, the distinctive variations in the fields of stress and strain are found around the crack-tip or damaged pores. Some key findings of this work are:

1. In a domain without a crack, the modeling parameter “ $\beta$ ” affects the deformations, which is found in both parallel and perpendicular displacements to the loading directions along the reference line. This parameter controls the strength of a porous material with the change of volumetric strain; the parameter $\beta$ with its sign and magnitude is closely related to the preferential stiffness, and when $\beta =0$ , the classical linearized model can be recovered.

2. For a domain with an edge crack under three different types of loading, both crack-tip stress concentration and strain energy density are found to be the maximum directly in front of the crack-tip, which is consistent with the observation obtained within the classical linearized elasticity model. This implies that the density-dependent moduli model can apply to the crack propagation phenomenon.

3. For a domain with the edge crack under pure tensile loading (or mode-I type of loading), the crack-tip axial strains are larger with smaller bulk modulus for higher positive $\beta$ -values, which implies that the material behaves less stiff (or more compressible). But for lower negative $\beta$ -values, its response is weaker against compression. Thus, under the mode-I tensile loading, the material behavior with positive $\beta$ can be compared to the strain softening, whereas one with the negative $\beta$ to the strain hardening in the elastoplasticity, even though only the pure elastic regime is considered in this study.

4. In the example with in-plane shear loading (or mode-II type of loading), the density-dependent model for positive $\beta$ -values distinctively shows more resistance against compression with the higher stress concentration therein.

5. Finally, for the mixed-mode loading (combination of mode-I and mode-II loadings), under about the same amount of volumetric strain changes on the reference line, the deviations between the cases are concentrated still around the crack-tip, and the positive $\beta$ -values show larger in ${\mathrm{K}}_{\mathrm{II}}$ but less in ${\mathrm{K}}_I$ in front of the tip. Thus, under a similar type or mixed-mode loading, a certain failure may occur preferentially with different material properties; for example, the shear failure may occur first for the brittle material (e.g. rock or ceramic). This anisotropic stiffness can be modeled straightforwardly with the $\beta$ -values in the density-dependent material moduli.

As the fracture toughness or its propagation is known to be related to the newly generated porosity or damaged pores, the model for an elastic porous solid with the preferential material property can be expanded to study a (quasi-static) crack evolution in the porous material via an appropriate numerical method such as the regularized phase-field approach.^43,44 Application of the extended finite element method (similar to the works in Riazi et al. ⁶⁵ and Afrazi et al. ⁶⁶ ) to study quasi-static evolution problems within the context of the proposed model can be another key idea for future investigation. Another important next topic can include bridging the modeling parameter to the material parameter via comparing with some experimental data and identifying its role for the purpose of characterization.