*Article* **A New Material Model for Agglomerated Cork**

**Gabriel Thomaz de Aquino Pereira 1, Ricardo J. Alves de Sousa 2,3, I-Shih Liu 4, Marcello Goulart Teixeira 1,\* and Fábio A. O. Fernandes 2,3**


**Abstract:** It is increasingly necessary to promote means of production that are less polluting and less harmful to the environment following the UN 2030 agenda for sustainable development. Using natural cellular materials in structural applications can be essential for enabling a future in this direction. Cork is a natural cellular material with an excellent energy absorption capacity. Its use in engineering applications and products has grown over time, so predicting its mechanical response through numerical tools is crucial. Classical cork modeling uses a model developed for foam material, including an adjustment function that does not have a clear physical interpretation. This work presents a new material model for an agglomerated cork based solely on well-known hypotheses of continuum mechanics using fewer parameters than the classical model and further a finite element framework to validate the new model against experimental data.

**Keywords:** agglomerated cork; material modeling; successive linear approximation; finite element

#### **1. Introduction**

The use of cellular materials in engineering applications has been established worldwide. This kind of material has excellent crashworthiness and insulation properties. Styrofoam (expanded polystyrene), derived from oil, is widely present in packaging and safety apparel, with excellent cost/performance ratios but with well-known recyclability and biodegradability issues after usage. Cork, the outer bark of *Quercus suber* L. tree, is a natural cellular material by excellence, allying crashworthiness and insulation properties. Compared with its synthetic counterparts, it is a sustainable and eco-friendly solution.

Dart and Eugene [1] made one of the first attempts to characterize cork mechanical behavior under compression loading. They found that each stress–time curve for various compressions can be obtained from each other by scaling. Such scalability was justified by the shape of the stress–strain compression curve. A highly non-linear stress–strain response curve characterizes agglomerated cork under compression. This curve is divided into three parts: a small elastic linear behavior at the beginning, followed by a plateau, and finished with a highly non-linear densification part.

All types of cork, natural, agglomerated, and expanded, are characterized by this compression behavior. The material properties vary with density, cellular dimensions, and porosity, as seen in [2,3]. The Poisson effect in natural cork compression was studied by [4]. They found that for compression in axial and tangential directions, the Poisson's ratio is almost null, while in the radial direction, cork presents a Poisson's ratio of 0.3. The cell geometry of the cork explained this effect.

Unlike natural cork, which presents the natural material anisotropy [5], cork agglomerates are obtained by compressing together randomly oriented cork grains [6]. This manufacturing procedure promotes a fairly regular isotropic mechanical behavior in cork agglomerates.

**Citation:** Pereira, G.T.d.A.; Sousa, R.J.A.d.; Liu, I.-S.; Teixeira, M.G.; Fernandes, F.A.O. A New Material Model for Agglomerated Cork. *Math. Comput. Appl.* **2022**, *27*, 92. https:// doi.org/10.3390/mca27060092

Academic Editor: Gianluigi Rozza

Received: 3 September 2022 Accepted: 7 November 2022 Published: 9 November 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

The authors performed experimental tests to obtain the Young and shear modulus in [7,8]. To determine the Young modulus, a tensile test was conducted, and, in the case of the shear modulus, they performed a torsion test on a cork cylinder. Ref. [8] also concluded that the amount of strain necessary to fracture in the radial direction is much larger than in the other directions. More recently, agglomerated cork is being used as an ideal core material for sandwich components of lightweight structures [9–14]. This structure is interesting because cork is used as an energy dissipator.

Usually, a material characterized by a mechanical behavior similar to cork's under compression is numerically modeled as a foam material model [15]. The Ogden-Hill hyperelastic model, [16,17], is a well-known numerical model for elastomeric polymer foams and hence a good model for agglomerated cork. More recently, in [18], an investigation was carried out to determine the influence of the variability related to the material properties of natural cork based on a numerical homogenization approach to predict the temperature-dependent equivalent elastic properties. These approaches consider only the elastic behavior of the material since the plastic period is only reached through substantial deformation and high strain energies.

This work focuses on the modeling of the mechanical response of cork. We present a new material model for cork agglomerates based on an extension of the Mooney–Rivlin material. The material parameterization will be made through an optimization problem using uniaxial and equibiaxial experimental results and the analytical solutions for these compression tests. To the authors' best knowledge, it is also the first time that experimental data from biaxial compressions of cork has been reported. After the parameterization, a finite element framework is used to validate this new material model against experimental data.

Due to the non-linear nature of both the material model and the large deformation, we use the successive linear approximation method for numerical simulation purposes to deal with the problem. This method calculates the constitutive equations at each state, with the reference configuration updated for each time step. The new reference configuration is the current configuration of the body. Assuming that in each time step that occurs, a small deformation is added, both the constitutive equations and the PDE system are linearized.

In Section 2, we briefly present the Successive Linear Approximation method used in this work, and in Section 3, the hyperelastic material model is considered. Section 4 presents the parameter fitting for the cork agglomerates model, and finally, in Section 5, we show the comparison between numerical results, obtained by a finite element code against analytical and experimental results.

#### **2. Successive Linear Approximation Method**

The Successive Linear Approximation (SLA) Method [19], a relative Lagrangian formulation based on the "small-on-large" idea, allows implementing the solution to large deformation in a successive incremental manner. In other words, at each time step, the constitutive function is calculated at the present state of deformation, which will be regarded as the reference configuration for the next state. From this point of view, it is considered a relative motion description; see, for example [20]. Assuming that the deformation to the next state is small, the constitutive function and the partial differential equation can be linearized. This procedure for large deformation problems was presented in [21].

This approach has significant fields of application, such as salt tectonics [22]. In [23], the influence of temperature in salt domes, a thermoviscoelastic material, was studied.

#### *2.1. Relative Motion Description*

Let *κ*<sup>0</sup> be the preferred reference configuration of an elastic body B, B<sup>0</sup> = *κ*0(B), and **x** = *χ*(*X*, *t*), with *X* ∈ B<sup>0</sup> be the motion of the body, see Figure 1. The configuration with specific material symmetries, such as isotropy, is usually chosen as the reference configuration. Such material symmetries may be lost in the deformed configuration in general.

**Figure 1.** The deformation and deformation gradient diagram.

In this paper, we consider the time *t* as the present time. Thus, *κ<sup>t</sup>* is the deformed configuration at time *t*, B*<sup>t</sup>* = *κt*(B). Now, we can calculate the deformation gradient with respect to configuration *κ*<sup>0</sup> as

$$\mathbf{F}(X,t) = \nabla\_X \chi(\mathbf{X}, \mathbf{t}).\tag{1}$$

Now, let *κτ* be the deformed configuration at time *τ*. The deformation *ξ* = *χ*(*X*, *τ*) from *κ*<sup>0</sup> can be described in the current configuration at the present time *t* by

$$\mathbf{x} = \chi(X, \tau) := \chi\_t(\mathbf{x}, \tau) \quad \text{ for} \quad \mathbf{x} = \chi(X, t), \tag{2}$$

where *χt*(·, *τ*) : B*<sup>t</sup>* −→ B*<sup>τ</sup>* is called relative deformation. With these in mind we can define the relative displacement vector **u** as

$$\mathbf{u} = \boldsymbol{\xi} - \mathbf{x} = \mathbf{u}\_l(\mathbf{x}, \boldsymbol{\tau}) = \chi\_l(\mathbf{x}, \boldsymbol{\tau}) - \mathbf{x}, \quad \mathbf{x} \in \mathcal{B}\_l. \tag{3}$$

Taking the gradient relative to *X* and **x**, respectively, we have

$$\begin{aligned} \mathbf{H}\_{l}(\mathbf{x},\tau)\mathbf{F}(X,t) &= \nabla\_{\mathbf{x}}\mathbf{u}\_{l}(\mathbf{x},\tau)\mathbf{F}(X,t) = \mathbf{F}(X,\tau) - \mathbf{F}(X,t),\\ \mathbf{H}\_{l}(\mathbf{x},\tau) &= \nabla\_{\mathbf{x}}\chi\_{l}(\mathbf{x},\tau) - \mathbf{I} = \mathbf{F}\_{l}(\mathbf{x},\tau) - \mathbf{I},\end{aligned} \tag{4}$$

where **I**, **H***t*(**x**, *τ*) and **F***t*(**x**, *τ*) are the identity tensor, relative displacement gradient and relative deformation gradient, respectively. We can rewrite the above equations as

$$\mathbf{F}\_{l}(\mathbf{x},\tau) = \mathbf{F}(X,\tau)\mathbf{F}(X,t)^{-1}, \quad \mathbf{F}(X,\tau) = (\mathbf{I} - \mathbf{H}\_{l}(\mathbf{x},\tau))\mathbf{F}(X,t). \tag{5}$$

With all these equations, we can define the motions of a body, and it is called relative description formulation.

#### *2.2. Linearized Constitutive Equation*

Consider *τ* = *t* + Δ*t*. By taking Δ*t* small enough, we can assume that the gradient of the displacement is small, and for simplicity, we denote

$$\mathbf{H}(\tau) := \mathbf{H}\_l(\mathbf{x}, \tau)), \text{ with } \parallel \mathbf{H}(\tau) \parallel \ll 1,\tag{6}$$

moreover, from (4), we have

$$\mathbf{F}(\tau) - \mathbf{F}(t) = \mathbf{H}(\tau)\mathbf{F}(t), \quad \mathbf{F}\_l(\tau) = \mathbf{I} + \mathbf{H}(\tau). \tag{7}$$

Without loss of generality, the Cauchy stress tensor of an elastic body in the preferred reference configuration is given by

$$\mathbf{T}(X,t) = -p\mathbf{I} + \mathcal{F}\_{\mathbb{X}\_0}(\mathbf{F}(X,t)),\tag{8}$$

where, for compressible bodies, *p* = *p*(**F**). We shall assume that this dependence is only through the determinant, or equivalently, from the mass balance, depending only on the mass density, i.e., *p* = *p*(*ρ*) where *ρ* = *<sup>ρ</sup>*<sup>0</sup> det**<sup>F</sup>** .

For time *τ* = *t* + Δ*t*, we have

$$\rho(\tau) - \rho(t) = \rho\_0(\det \mathbf{F}(\tau)^{-1} - \det \mathbf{F}(t)^{-1}) = \rho(t)(\det(\mathbf{F}(t)\mathbf{F}(\tau)^{-1}) - 1)$$

$$= \rho(t)(\det(\mathbf{I} + \mathbf{H}(\tau))^{-1} - 1) = -\rho(t)\text{tr}\mathbf{H}(\tau) + o(2),\tag{9}$$

where in the passages, (5) and (7) were used. Using Taylor series expansion in (8) and (9), it follows that

$$\mathbf{T}(\tau) = \mathbf{T}(t) - \left(\frac{\partial p}{\partial \rho}\right)\_t (\rho(\tau) - \rho(t))\mathbf{I} + \nabla\_\mathbf{F} \mathcal{F}(\mathbf{F}(t))(\mathbf{F}(\tau) - \mathbf{F}(t))$$

$$= \mathbf{T}(t) - \beta(\text{tr}\mathbf{H}(\tau))\mathbf{I} + \nabla\_\mathbf{F} \mathcal{F}(\mathbf{F}(t))(\mathbf{H}(\tau)\mathbf{F}(t), \tag{10}$$

or

$$\mathbf{T}(\tau) = \mathbf{T}(t) + \mathcal{L}(\mathbf{F}(t))[\mathbf{H}(\tau)],\tag{11}$$

where *β* = *ρ ∂p ∂ρ* is a material parameter and

$$\mathcal{L}(\mathbf{F}(t))[\mathbf{H}(\tau)] := -\beta(\text{tr}\mathbf{H}(\tau))\mathbf{I} + \nabla\_{\mathbf{F}}\mathcal{F}(\mathbf{F}(t))[\mathbf{H}(\tau)\mathbf{F}(t)] \tag{12}$$

defines the fourth-order elasticity tensor, relative to the current configuration *κt*. Furthermore, the first Piola–Kirchhoff stress tensor at time *τ* relative to the current configuration, is given by

$$\begin{aligned} \mathbf{T}\_{\mathbb{K}\_{t}}(\tau) &= \det \mathbf{F}\_{t}(\tau) \mathbf{T}(\tau) \mathbf{F}\_{t}(\tau)^{-T} = \det(\mathbf{I} + \mathbf{H}) \mathbf{T}(\tau) (\mathbf{I} + \mathbf{H})^{-T} \\ &= \det(\mathbf{I} + \mathbf{H}) (\mathbf{T}(t) + \mathcal{L}(\mathbf{F})[\mathbf{H}] + o(2)) (\mathbf{I} + \mathbf{H})^{-T} \\ &= (\mathbf{I} + \mathbf{t} \mathbf{r} \mathbf{H}) (\mathbf{T}(t) + \mathcal{L}(\mathbf{F})[\mathbf{H}]) (\mathbf{I} - \mathbf{H}^{T}) + o(2) \\ &= \mathbf{T}(t) + (\mathbf{t} \mathbf{r} \mathbf{H}) \mathbf{T}(t) - \mathbf{T}(t) \mathbf{H}^{T} + \mathcal{L}(\mathbf{F})[\mathbf{H}] + o(2) . \end{aligned} \tag{13}$$

and we write the linearized first Piola–Kirchhoff as

$$\mathbf{T}\_{\mathbb{K}\_l}(\tau) = \mathbf{T}(t) + \mathcal{K}(\mathbf{F}(t), \mathbf{T}(t))[\mathbf{H}(\tau)],\tag{14}$$

where

$$\mathcal{K}(\mathbf{F}, \mathbf{T})[\mathbf{H}] := (\mathbf{tr}\mathbf{H})\mathbf{T}(t) - \mathbf{T}(t)\mathbf{H}^T + \mathcal{L}(\mathbf{F})[\mathbf{H}] \tag{15}$$

is the fourth-order elasticity tensor for the first Piola–Kirchhoff stress tensor.

To define the fourth-order elasticity tensor we have to define the material model. In Section 3 we will propose one based on Mooney–Rivlin hyperelastic material model.

#### **3. Hyperelastic Material Model**

To obtain a constitutive model for large strains of an isotropic elastic solid, we shall start with the free energy function *ψ* = *ψ*(I**B**, II**B**, III**B**), and the Cauchy stress given by

$$\mathbf{T} = 2\rho \frac{\partial \boldsymbol{\psi}}{\partial \mathbf{B}} \mathbf{B}\_{\prime} \tag{16}$$

where **B** is the left Cauchy–Green strain tensor and I**B**,II**B**,III**<sup>B</sup>** are the first, second and third invariants of **B**. By Taylor series expansion, we can write

$$\begin{split} \psi(\mathrm{I\_{B}}, \mathrm{II\_{B}}, \mathrm{III\_{B}}) &= \psi\_{0}(\mathrm{3,3}, \mathrm{III\_{B}}) + \psi\_{1}(\mathrm{I\_{B}} - \mathrm{3}) + \psi\_{2}(\mathrm{II\_{B}} - \mathrm{3}) \\ &+ \frac{1}{2}\psi\_{3}(\mathrm{I\_{B}} - \mathrm{3})^{2} + \psi\_{4}(\mathrm{I\_{B}} - \mathrm{3})(\mathrm{II\_{B}} - \mathrm{3}) + \frac{1}{2}\psi\_{5}(\mathrm{II\_{B}} - \mathrm{3})^{2}, \end{split} \tag{17}$$

where *ψ<sup>k</sup>* = *ψk*(III**B**) for *k* = 0, 1 ... , 5, which stands for the expansion up to the second order for moderate strains.

From (16), by the use of the relations

$$\frac{\partial \mathbf{I\_B}}{\partial \mathbf{B}} = \mathbf{I}, \qquad \frac{\partial \Pi\_\mathbf{B}}{\partial \mathbf{B}} = \mathbf{I\_B} \mathbf{I} - \mathbf{B}, \qquad \frac{\partial \Pi \Pi\_\mathbf{B}}{\partial \mathbf{B}} = \Pi \Pi\_\mathbf{B} \mathbf{B}^{-1}, \tag{18}$$

we obtain

$$\begin{aligned} \mathbf{T} + p\mathbf{I} &= 2\rho [\psi\_1 + \psi\_3(\mathbf{I\_B} - 3) + \psi\_4(\mathbf{II\_B} - 3)] \mathbf{B} \\ &- 2\rho [\psi\_2 + \psi\_4(\mathbf{I\_B} - 3) + \psi\_5(\mathbf{II\_B} - 3)] (\mathbf{B}^2 - \mathbf{I\_B} \mathbf{B}), \end{aligned} \tag{19}$$

where terms with the identity tensor are lumped into the pressure *p*(I**B**, II**B**, III**B**). By the use of the Cayley–Hamilton theorem,

$$\mathbf{B}^2 - \mathbf{I}\_\mathbf{B}\mathbf{B} = \Pi\mathbf{I}\_\mathbf{B}\mathbf{B}^{-1} - \Pi\_\mathbf{B}\mathbf{I},\tag{20}$$

it becomes

$$\begin{aligned} \mathbf{T} + p\mathbf{I} &= 2\rho[\psi\_1 + \psi\_3(\mathbf{I\_B} - 3) + \psi\_4(\Pi\_\mathbf{B} - 3)]\mathbf{B} \\ &- 2\rho[\psi\_2 + \psi\_4(\mathbf{I\_B} - 3) + \psi\_5(\Pi\_\mathbf{B} - 3)]\Pi\mathbf{I\_B}\mathbf{B}^{-1}, \end{aligned} \tag{21}$$

in which the term II**BI** is absorbed into the indeterminate pressure *p***I**.

We can rewrite the above equation with the definition of the parameters,

$$\begin{aligned} s\_1 &= 2\rho(\psi\_1 - 3\psi\_3 - 3\psi\_4), \\ s\_2 &= 2\rho \Pi \Pi\_\mathbf{B} (\psi\_2 - 3\psi\_4 - 3\psi\_5), \\ s\_3 &= 2\rho \psi\_3, \quad s\_4 = 2\rho \Pi \Pi\_\mathbf{B} \psi\_{4\prime} \\ s\_5 &= 2\rho \psi\_4, \quad s\_6 = 2\rho \Pi \Pi\_\mathbf{B} \psi\_{5\prime} \end{aligned} \tag{22}$$

and propose an isotropic elastic model as:

**An extended Mooney–Rivlin material.** *The constitutive equation*

$$\mathbf{T} = -p\mathbf{I} + (s\_1 + s\_3 \mathbf{I}\_\mathbf{B} + s\_5 \mathbf{II}\_\mathbf{B})(\mathbf{B} - \mathbf{I}) - (s\_2 + s\_4 \mathbf{I}\_\mathbf{B} + s\_6 \mathbf{II}\_\mathbf{B})(\mathbf{B}^{-1} - \mathbf{I}),\tag{23}$$

*is an extended version of Mooney–Rivlin model for isotropic solids at large strain. We shall assume that six parameters s*1, ...,*s*<sup>6</sup> *are material constants, and the pressure is a function of mass density only, p*(*ρ*) *so that*

$$
\beta = \rho \frac{\partial p}{\partial \rho'} \tag{24}
$$

*is a material parameter.*

With the proposed constitutive equation, from (12) we can derive the fourth-order elastic tensor. Therefore, in Einstein notation, it follows that

$$\begin{split} \mathcal{L}\_{ijkl} &= \beta \left( \delta\_{ij} \delta\_{kl} \right) \\ &+ \left( s\_1 + s\_3 \mathbf{I\_B} + s\_5 \mathbf{II\_B} \right) \left( B\_{lj} \delta\_{ki} + B\_{il} \delta\_{kj} \right) \\ &+ \left( s\_2 + s\_4 \mathbf{I\_B} + s\_6 \mathbf{II\_B} \right) \left( B\_{ik}^{-1} \delta\_{jl} + B\_{kj}^{-1} \delta\_{il} \right) \\ &+ 2 \left( s\_3 B\_{kl} + s\_5 \left( \Pi\_\mathbf{B} \delta\_{kl} - \Pi \Pi\_\mathbf{B} B\_{kl}^{-1} \right) \right) \left( B\_{ij} - \delta\_{ij} \right) \\ &- 2 \left( s\_4 B\_{kl} + s\_6 \left( \Pi\_\mathbf{B} \delta\_{kl} - \Pi \Pi\_\mathbf{B} B\_{kl}^{-1} \right) \right) \left( B\_{ij}^{-1} - \delta\_{ij} \right) . \end{split} \tag{25}$$

**Remark on** *β* **value for cork:** Since, in this paper, we are considering the cork agglomerated as a material, from the definition of *β* and the compressibility of cork, it follows that *β* ≈ 0.

**Remarks on linear model:** Let **<sup>B</sup>** <sup>=</sup> **<sup>I</sup>** <sup>+</sup> <sup>2</sup>**E**, hence **<sup>B</sup>**−<sup>1</sup> <sup>=</sup> **<sup>I</sup>** <sup>−</sup> <sup>2</sup>**<sup>E</sup>** and tr**<sup>H</sup>** <sup>=</sup> tr**E**, for small linear strain **E**, then

$$\mathbf{I\_B = tr(I + 2E) = 3 + 2trE},\tag{26}$$

$$\mathbf{I}\_{\mathbf{B}} = (\det \mathbf{B}) \mathbf{I}\_{\mathbf{B}^{-1}} = (1 + \text{tr}\mathbf{E}) \text{tr}(\mathbf{I} - 2\mathbf{E}) \tag{27}$$

$$\mathbf{E} = (\mathbf{1} + \mathbf{tr}\mathbf{E})(\mathbf{3} - 2\mathbf{tr}\mathbf{E}) = \mathbf{3} + \mathbf{tr}\mathbf{E},\tag{28}$$

and, using a Taylor expansion of *p*(*ρ*), the stress becomes

$$\mathbf{T} = -(p\_0 - \beta\_0 \text{tr}\mathbf{E})\mathbf{I} + 2[s\_1 + s\_3(\mathbf{3} + 2\text{tr}\mathbf{E}) + s\_5(\mathbf{3} + \text{tr}\mathbf{E})]\mathbf{E} \tag{29}$$

$$-2\left[\mathbf{s}\_2 + \mathbf{s}\_4(\mathbf{3} + 2\mathbf{tr}\mathbf{E}) + \mathbf{s}\_6(\mathbf{3} + \mathbf{tr}\mathbf{E})\right] \mathbf{E} \tag{30}$$

$$=-p\_0 \mathbf{I} + \lambda(\text{tr}\mathbf{E})\mathbf{I} + 2\mu \mathbf{E} + o(2),\tag{31}$$

where

$$
\lambda = \pounds\_{0\prime} \tag{32}
$$

$$
\mu = s\_1 - s\_2 + \Im(s\_3 + s\_5 - s\_4 - s\_6) \tag{33}
$$

are the elastic Lamé constants. Note that we have the relation from linear elasticity,

$$
\beta\_0 = \lambda = \frac{2\nu\mu}{1 - 2\nu'}\tag{34}
$$

where *ν* is the Poisson ratio.

#### **4. Parameter Fitting for Cork Agglomerates**

Cork agglomerate can be modeled as a hyperelastic material. As any other hyperelastic material, to parametrize an agglomerated cork we need at least two different deformation modes to capture the correct material behavior [24].

In this paper, we will use the uniaxial and the equibiaxial compression tests, schematically represented in Figure 2, to parametrize the extended Mooney–Rivlin model proposed before. For the parametrization procedure adopted we need the experimental data, the analytical solutions for both deformation modes and a good minimization function.

**Figure 2.** Uniaxial and equibiaxial compression tests schemes.

#### *4.1. Uniaxial Compression*

The uniaxial extension is given by

$$x = \lambda X, \quad y = \mathbf{Y}, \quad z = Z,\tag{35}$$

where *λ* is the stretch along the loading. For agglomerated cork, in uniaxial extension, the transversal stretches have the same value and they are near 1, which is why we are not considering them. Therefore, the gradient deformation, the left Cauchy–Green strain tensor and its inverse are given, respectively, by

$$\mathbf{F} = \begin{bmatrix} \lambda & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}, \quad \mathbf{B} = \begin{bmatrix} \lambda^2 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}, \quad \mathbf{B}^{-1} = \begin{bmatrix} \frac{1}{\lambda^2} & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}.$$

and from that, the first, second and third invariants of **B** are

$$\Pi\_{\mathsf{B}} = 2 + \lambda^2, \qquad \Pi\_{\mathsf{B}} = 2\lambda^2 + 1, \qquad \Pi \Pi\_{\mathsf{B}} = \lambda^2. \tag{36}$$

Consequently, the Cauchy stress tensor, *T*(*λ*), can be obtained by inserting (36) into (23):

$$\begin{split} \mathbf{T}(\boldsymbol{\lambda}) &= -p\mathbf{I} + [s\_1 + s\_3(2 + \boldsymbol{\lambda}^2) + s\_5(2\boldsymbol{\lambda}^2 + 1)](\mathbf{B} - \mathbf{I}) \\ &- [s\_2 + s\_4(2 + \boldsymbol{\lambda}^2) + s\_6(2\boldsymbol{\lambda}^2 + 1)](\mathbf{B}^{-1} - \mathbf{I}). \end{split} \tag{37}$$

Since there is no stretch in the second and third direction we have *TU*2(*λ*) = *TU*3(*λ*) = 0. Therefore, *p* = 0 and the principal stress in the first direction is defined by:

$$\begin{split} T\_{lI1}(\lambda) &= [s\_1 + s\_3(2 + \lambda^2) + s\_5(2\lambda^2 + 1)](\lambda^2 - 1) \\ &- [s\_2 + s\_4(2 + \lambda^2) + s\_6(2\lambda^2 + 1)](\lambda^{-2} - 1), \end{split} \tag{38}$$

#### *4.2. Equibiaxial Compression*

The equibiaxial extension is given by

$$\mathbf{x} = \lambda X, \quad \mathbf{y} = \lambda Y, \quad z = Z,\tag{39}$$

Therefore, the deformation gradient, the left Cauchy–Green strain tensor and its inverse are given, respectively, by

$$\mathbf{F} = \begin{bmatrix} \lambda & 0 & 0 \\ 0 & \lambda & 0 \\ 0 & 0 & 1 \end{bmatrix}, \quad \mathbf{B} = \begin{bmatrix} \lambda^2 & 0 & 0 \\ 0 & \lambda^2 & 0 \\ 0 & 0 & 1 \end{bmatrix}, \quad \mathbf{B}^{-1} = \begin{bmatrix} \frac{1}{\lambda^2} & 0 & 0 \\ 0 & \frac{1}{\lambda^2} & 0 \\ 0 & 0 & 1 \end{bmatrix}.$$

and from that, the first, second and third invariants of **B** are

$$\mathrm{II}\_{\mathsf{B}} = 1 + 2\lambda^2, \qquad \mathrm{II}\_{\mathsf{B}} = 2\lambda^2 + \lambda^4, \qquad \mathrm{III}\_{\mathsf{B}} = \lambda^4. \tag{40}$$

Consequently, the Cauchy stress tensor, *T*(*λ*), can be obtained by inserting (40) into (23):

$$\begin{split} \mathbf{T}(\boldsymbol{\lambda}) &= -p\mathbf{I} + [s\_1 + s\_3(1 + 2\lambda^2) + s\_5(2\lambda^2 + \lambda^4)](\mathbf{B} - \mathbf{I}) \\ &- [s\_2 + s\_4(1 + 2\lambda^2) + s\_6(2\lambda^2 + \lambda^4)](\mathbf{B}^{-1} - \mathbf{I}). \end{split} \tag{41}$$

Since there is no stretch in the third direction, we have *TB*3(*λ*) = 0. Therefore, *p* = 0 and the principal stress in the first and second directions are defined by:

$$T\_{B1}(\lambda) = T\_{B2}(\lambda) = [s\_1 + s\_3(1 + 2\lambda^2) + s\_5(2\lambda^2 + \lambda^4)](\lambda^2 - 1)$$

$$- [s\_2 + s\_4(1 + 2\lambda^2) + s\_6(2\lambda^2 + \lambda^4)](\lambda^{-2} - 1). \tag{42}$$

Equations (38) and (42) are the analytical solutions for the uniaxial and equibiaxial compression problems, respectively.

#### *4.3. Experimental Tests*

To parameterize our model, we have to collect some experimental data for both tests, uniaxial, and equibiaxial compression. The quasi-static uniaxial compression tests were performed in a universal testing machine, the Shimadzu AGS-X-10 kN. The equibiaxial compression tests were performed in an in-house built biaxial machine developed and properly validated in [25]. Both types of experiments were carried out at room temperature.

For uniaxial compression tests, 60 mm agglomerated cork cubes were manufactured, as shown in Figure 3. For the equibiaxial tests, 30 mm thick octagon-shaped samples were produced. The other dimensions are shown in Figure 4.

**Figure 3.** Agglomerated cork sample positioned for uniaxial compression testing at quasi-static strain rates.

**Figure 4.** Dimensions of the octagon-shaped sample and the experimental setup for the equibiaxial compression tests.

In the equibiaxial compression tests, we consider the values of stretch and stress at the central part of sample, a square with a 70 mm side. The values for stress and stretch for both tests are presented in Table 1.


**Table 1.** Values of experimental tests of uniaxial and equibiaxial compression.

#### *4.4. Curve Fitting*

To parametrize the material we use the Wolfram Mathematica software through the function Nminimize. This function implements four methods that don't need derivatives for minimization: Nelder–Mead, Simulated Annealing, Differential Evolution and Random Search. The results presented by the four were very similar with a slight advantage for the Nelder–Mead, used in this work.

This algorithm tries to find a maximum or minimum of an objective function by doing a direct search. This search uses simplexes as input data for the objective function and evaluates how close the minimum of this function is. One of its features is that the derivatives of the function may not be known.

As we have to consider both, uni and equibiaxial tests, our objective function has to take both experiments into account. Let us define *TExp <sup>U</sup>* , <sup>Λ</sup>*Exp <sup>U</sup>* , *<sup>T</sup>Exp <sup>B</sup>* , <sup>Λ</sup>*Exp <sup>B</sup>* as stress and stretch in uniaxial test and stress and stretch in equibiaxial test, respectively. Therefore, our objective function is given by

$$
\Theta(\mathbf{s}\_1, \dots, \mathbf{s}\_6) = Q\_1(\mathbf{s}\_1, \dots, \mathbf{s}\_6) + Q\_2(\mathbf{s}\_1, \dots, \mathbf{s}\_6), \tag{43}
$$

with,

$$Q\_1(s\_1, \ldots, s\_6) = \sum\_{r=1}^{R} \left[ \frac{T\_{\mathcal{U}\_r}^{Exp} - T\_{\mathcal{U}1}(\Lambda\_{\mathcal{U}\_r}^{Exp})}{T\_{\mathcal{U}\_r}^{Exp}} \right]^2,$$

$$Q\_2(s\_1, \ldots, s\_6) = \sum\_{s=1}^{S} \left[ \frac{T\_{B\_s}^{Exp} - T\_{B1}(\Lambda\_{B\_s}^{Exp})}{T\_{B\_s}^{Exp}} \right]^2. \tag{44}$$

where *<sup>R</sup>*, *<sup>S</sup>* are the number of experimental data for each case and *TU*1(Λ*Exp Ur* ) and *TB*1(Λ*Exp Bs* ) are the uniaxial and biaxial stress calculated, respectively, by (38) and (42), considering experimental deformations Λ*Exp Ur* and <sup>Λ</sup>*Exp Bs* . *Q*<sup>1</sup> and *Q*<sup>2</sup> are the normalized error function

for each of the tests and both depend on the parameters *s*1, ... ,*s*6. This objective function was based on that used in the [24].

In Table 2, we present the obtained values for each parameter. The experimental and analytical curves, using the obtained parameters, are presented in Figure 5.



**Figure 5.** Comparison between analytical and experimental results for uniaxial and equibiaxial tests, respectively.

#### **5. Linearized Partial Differential Equations**

Let <sup>Ω</sup> <sup>=</sup> {*<sup>x</sup>* <sup>∈</sup> *<sup>κ</sup>t*(B)} ⊂ <sup>R</sup><sup>3</sup> be the region occupied by the body at the present configuration *κt*, and let *∂*Ω = Γ<sup>1</sup> ∪ Γ<sup>2</sup> be the disjoint unions of its boundary. Let *n*(*x*, *t*) be the exterior unit normal to *∂*Ω at the present time.

At time *τ* > *t*, we shall consider an initial boundary value problem in Lagrangian formulation, with the present state at time *t* as the reference configuration, given by

$$
\rho(t)\,\ddot{u}\_t(\tau) - \text{div}\_x T\_t(\tau) = \rho(t)\mathbf{b}(\tau), \quad \text{in } \Omega,\tag{45}
$$

$$T\_t(\tau)\,\mathfrak{u}(t) = f(\tau), \quad \text{on } \Gamma\_{1\prime} \tag{46}$$

$$\mathfrak{u}\_{l}(\tau) = \mathfrak{d}(\tau), \quad \text{on } \Gamma\_{2\prime} \tag{47}$$

$$
\dot{\mathfrak{u}}\_l(t) = 0, \quad \dot{\mathfrak{u}}\_l(t) = \mathfrak{v}(t), \quad \text{in } \Omega. \tag{48}
$$

The body is subjected to the surface traction *f*(*x*, *τ*), the boundary displacement *d*(*x*, *τ*) on the respective parts of *∂*Ω at time *τ* > *t*, and the initial velocity *v*(*x*, *t*) in Ω at the present time *t*. Note that unlike the explicit time dependence in the above expressions, the spatial dependence is implicitly understood and is not explicitly indicated for simplicity.

In these relations, for simplicity, div*x* stands for the divergence operator relative to the coordinate *x* ∈ *κt*(B), which is the same as the operator div*κ<sup>t</sup>* for the reference configuration *κ<sup>t</sup>* in this case.

Together with the constitutive Equation (23), the mechanical problem is to be solved for the relative displacement vector *ut*(*τ*). Since the constitutive function **T** in (23) is generally nonlinear for finite deformations, the partial differential equation of this problem is genuinely nonlinear. However, in the relative Lagrangian formulation, for small enough incremental time Δ*t* = *τ* − *t*, we can easily linearize the constitutive equations relative to the present state at time *t*, so that the boundary value problem becomes linear.

By regarding the present state as the reference state, it is assumed that the state variables are all known functions at the present time *t*. Those include the deformation gradient *F*(*t*) and the stress *T*(*t*).

By use of the linearization (14), the partial differential equation of the problem become [26]

$$\rho(t)\,\ddot{u}\_t(\tau) - \text{div}\_x\Big(\mathcal{K}(t)[\nabla\_x u\_t(\tau)]\Big) = \text{div}\_x T(t) + \rho(t)b(\tau),\tag{49}$$

where the relevant material function K is defined in (15) and (25).

The Equation (49) is a linear partial differential equation for the relative displacement vector *ut*(*x*, *τ*) to be solved with the corresponding initial boundary conditions in the problems (45), for which the state variables of the body at time *t* are known and the external supplies *b*(*τ*) is given so that the right-hand side of the Equation (49) is known quantity.

In this linearization, we do not assume the deformation is small, rather only the relative displacement gradient with respect to the present state is assumed to be small. This is the idea of "small-on-large", the same as the well-known problem of small deformations superposed on finite deformation in the literature [27]. Therefore, the overall deformation may be of finite values, in contrast to the usual theory of linear elasticity which linearizes the problem with respect to the fixed reference configuration assuming a small deformation gradient.

#### **6. Variational Formulation**

Consider

$$V = \{ \mathfrak{v} \in \left( H^1(\Omega) \right)^n ; \mathfrak{v} = 0, \quad \text{on} \quad \Gamma\_2 \}.$$

Lets *u*, *v* ∈ *V*. By multiplying Equation (49) for *v*, integrating it by parts over Ω and using the definition of the inner product of second-order tensors (if *A* and *B* are second-order tensors, *<sup>A</sup>*.*<sup>B</sup>* <sup>=</sup> tr(*ABT*)), we have, <sup>∀</sup>*<sup>v</sup>* <sup>∈</sup> *<sup>V</sup>*,

$$\begin{split} \int\_{\Omega} \rho(t) \, \ddot{\mathfrak{u}}\_{t}(\tau) \cdot \mathfrak{v} d\Omega &+ \int\_{\Omega} \text{tr} \Big[ \Big( \mathcal{K}(t) [\nabla\_{x} \mathfrak{u}\_{t}(\tau)] \Big) \nabla\_{x} \mathfrak{v}^{T} \Big] d\Omega \\ &- \int\_{\partial\Omega} \underline{\mathfrak{v}} \cdot \Big( \mathcal{K}(t) [\nabla\_{x} \mathfrak{u}\_{t}(\tau)] \Big) \mathfrak{u}\_{t} d\Gamma = \int\_{\Omega} \text{div}\_{x} \, T\_{\varepsilon}(t) \cdot \mathfrak{v} d\Omega \\ &+ \int\_{\Omega} \rho(t) \, \mathfrak{b}(\tau) \cdot \mathfrak{v} d\Omega. \end{split} \tag{50}$$

By the boundary condition,

$$\begin{split} -\int\_{\Gamma\_{1}} \boldsymbol{\mathfrak{v}} \cdot \left( \mathcal{K}(t) [\nabla\_{\boldsymbol{x}} \boldsymbol{\mathfrak{u}}\_{t}(\boldsymbol{\tau})] \right) \boldsymbol{\mathfrak{u}}\_{k} d\boldsymbol{\Gamma} &= -\int\_{\Gamma\_{1}} \boldsymbol{\mathfrak{v}} \cdot \boldsymbol{f}(\boldsymbol{\tau}) d\boldsymbol{\Gamma} \\ -\int\_{\Gamma\_{1}} \boldsymbol{\mathfrak{v}} \cdot \boldsymbol{T}\_{t}(t) d\boldsymbol{\Gamma}\_{\prime} &\quad \forall \boldsymbol{\mathfrak{v}} \in V. \end{split} \tag{51}$$

Replacing (51) in (50) and using, again, integration by parts, we have the weak formulation,

$$\begin{aligned} \int\_{\Omega} \rho(t) \, \ddot{\boldsymbol{u}}\_{t}(\tau) \cdot \boldsymbol{\sigma} d\Omega + \int\_{\Omega} \text{tr}\left[ \left( \boldsymbol{\mathcal{K}}(t) [\nabla\_{\boldsymbol{x}} \boldsymbol{u}\_{t}(\tau)] \right) \nabla\_{\boldsymbol{x}} \boldsymbol{\sigma}^{T} \right] d\Omega &= \\ \int\_{\Omega} \rho(t) \, \boldsymbol{b}(\tau) \cdot \boldsymbol{\sigma} d\Omega - \int\_{\Omega} \boldsymbol{T}\_{\boldsymbol{\varepsilon}}(t) \cdot \nabla\_{\boldsymbol{x}} \boldsymbol{\sigma} d\Omega & + \int\_{\Gamma\_{1}} \boldsymbol{\sigma} \cdot \boldsymbol{f}(\tau) d\Gamma, \quad \forall \boldsymbol{\sigma} \in V, \end{aligned} \tag{52}$$

or we can rewrite this in terms of the bilinear forms

$$\begin{split} \left(\rho(t)\,\ddot{u}\_{t}(\tau),\boldsymbol{\upsilon}\right) + a\Big(\boldsymbol{\kappa}(t)[\nabla\_{\boldsymbol{x}}u\_{t}(\tau)],\boldsymbol{\upsilon}\Big) &= \left(\rho(t)\,b(\tau),\boldsymbol{\upsilon}\right) - \\ \left(T\_{\boldsymbol{\varepsilon}}(t),\nabla\_{\boldsymbol{x}}\boldsymbol{\upsilon}\right) + \left(\boldsymbol{\upsilon},f(\tau)\right)\_{\Gamma\_{1}} \qquad \forall \boldsymbol{\upsilon} \in V. \end{split} \tag{53}$$

Let us consider *V<sup>h</sup>* a finite subspace of *V*. By restricting the formulation (53) to the space *V<sup>h</sup>* we have

$$\begin{split} \left(\rho(t)\,\dot{\mathfrak{u}}\_{t}^{h}(\tau),\boldsymbol{\mathfrak{v}}^{h}\right) + a\Big(\mathcal{K}(t)[\nabla\_{\boldsymbol{x}}\boldsymbol{\mathfrak{u}}\_{t}^{h}(\tau)],\boldsymbol{\mathfrak{v}}^{h}\Big) &= \left(\rho(t)\,\mathfrak{b}(\tau),\boldsymbol{\mathfrak{v}}^{h}\right) - \\ \left(T\_{\varepsilon}(t),\nabla\_{\boldsymbol{x}}\boldsymbol{\mathfrak{v}}^{h}\right) + \left(\boldsymbol{\mathfrak{v}}^{h},f(\boldsymbol{\tau})\right)\_{\Gamma\_{1}} \qquad \forall \boldsymbol{\mathfrak{v}}^{h} \in \boldsymbol{V}^{h}. \end{split} \tag{54}$$

Now consider {*ϕ*1, *<sup>ϕ</sup>*2, ··· , *<sup>ϕ</sup>n*} a basis of the subspace *<sup>V</sup>h*, that is, all elements *<sup>u</sup><sup>h</sup>* <sup>∈</sup> *<sup>V</sup><sup>h</sup>* can be expressed as

$$u^h = \sum\_{i=1}^n b\_i \varphi\_i. \tag{55}$$

Replacing in (54) *uh* by (55) and taking *<sup>v</sup><sup>h</sup>* <sup>=</sup> *<sup>ϕ</sup>j*, 1 <sup>≤</sup> *<sup>j</sup>* <sup>≤</sup> *<sup>n</sup>*, we have

$$
\mathcal{A}\ddot{\mathbf{b}} + \mathcal{L}\mathbf{b} = \mathcal{N},
\tag{56}
$$

where, for 1 ≤ *i*, *j* ≤ *n*

$$\begin{aligned} \mathcal{A}\_{ij} &= \left(\rho(t)\,\,\varrho\_{i\cdot}\,\,\rho\_{\dot{j}}\right), \\ \mathcal{L}\_{ij} &= a\Big(\mathcal{K}(t)[\nabla\_{x}\varrho\_{\dot{i}}], \,\rho\_{\dot{j}}\big), \\ \mathcal{N}\_{\dot{j}} &= \left(\rho(t)\,\,b(\tau), \,\rho\_{\dot{j}}\right) - \left(T\_{\epsilon}(t), \nabla\_{x}\varrho\_{\dot{j}}\right) + \left(\varrho\_{\dot{j}\cdot}f(\tau)\right)\_{\Gamma\_{1}}. \end{aligned} \tag{57}$$

#### **7. Numerical Validation**

For the numerical simulation, as the two cases studied are examples of large deformations, the SLA method that was presented in Section 2 was used. We consider the material model developed in Section 3 with the parameters defined in Section 4, Table 2. In both cases, we are interested in the deformation after each loading step. Thereby, we are considering quasi-static, time-independent problems.

Regarding the implementation, all the finite element codes were developed by the author. The solver of the linear algebra problem that naturally arises in the finite element context was also implemented by the author. This decision was taken, principally, to solve the linear system with a good algorithm that takes into account the finite element matrix. The Gaussian quadrature was the integration method adopted. All the code is developed in an object-oriented programming language, C++.

For uniaxial compression, we consider a mesh of 50 mm × 100 mm, with 10 elements in the x-direction and 10 elements in the y-direction. The boundary condition applied considered that a displacement was prescribed on the upper surface and the lateral surface is free, while the base is free horizontally but not vertically. It applied a prescribed displacement equal to 7 mm in each load step.

In the case of the equibiaxial test, we consider a mesh of 200 mm × 200 mm, with 10 elements in the x-direction and 10 elements in the y-direction. The symmetry of the problem on the x and y axes was considered, and with only a quarter of the model was modeled. Still considering the symmetry, equal displacements were prescribed on both x and y surfaces on the top and the right side, while the left was free to move vertically and the bottom was free to move horizontally. In this case, the prescribed displacement in each load step is equal to 10 mm both in *x* and *y* directions.

For both tests, we use 100 steps of SLA and Q4 element, a bilinear quadrilateral element which combines two sets of Lagrange polynomials, i.e., linear isoparametric quadrilateral elements with four nodes. Figure 6 shows the finite element mesh and the boundary conditions for both cases.

**Figure 6.** Finite element mesh and boundary condition for both tests, uniaxial (**left**) and equibiaxial (**right**).

In Figure 7, we can see the results obtained for the two cases studied. From left to right, we can see the comparison between the numerical and analytical results of the equibiaxial and uniaxial tests. As expected, the equibiaxial test has a smaller plateau area and densification occurs earlier than in the uniaxial test.

**Figure 7.** Comparison between numerical, analytical and experimental solutions for uniaxial and equibiaxial compression tests, respectively.

These two simulations show that the agglomerated cork model presented represents with good accuracy the real behavior of the material, and therefore, it is a good option to be studied from now on.

#### **8. Conclusions and Future Work**

In this paper, we present a new material model for agglomerated cork. This model is based on the Mooney–Rivlin hyperelastic model and, as the cork agglomerate has a Poisson ratio near zero, has six material parameters. Analytical solutions for uniaxial and equibiaxial tests were developed and used for parameterization. To parameterize, the Nelder–Mead algorithm was used.

We also describe and use the successive linear approximation method (SLA) for the numerical simulation. The description of the SLA was quite general, and can be used for any material.

Our results show that the presented model is a good model and was validated through numerical simulation. One of the advantages of our model is the mathematical simplicity in relation to the classical model [16,17], being a direct extension of a Mooney–Rivlin model.

For future work, it will be interesting to try to relate the parameters *s*1, *s*2, *s*3, *s*4, *s*<sup>5</sup> and *s*<sup>6</sup> with physical parameters of the material, as was performed directly with the *β* and the Poisson ratio. Another interesting approach to test the efficiency of the SLA will be the study of a contact-impact problem.

**Author Contributions:** Conceptualization, G.T.d.A.P., R.J.A.d.S., I.-S.L., M.G.T. and F.A.O.F.; methodology, G.T.d.A.P., R.J.A.d.S., I.-S.L., M.G.T. and F.A.O.F.; formal analysis, G.T.d.A.P., I.-S.L. and M.G.T.; investigation, G.T.d.A.P. and M.G.T.; resources, R.J.A.d.S., M.G.T. and F.A.O.F.; data curation, G.T.d.A.P.; writing—original draft preparation, G.T.d.A.P.; writing—review and editing, G.T.d.A.P., R.J.A.d.S., I.-S.L., M.G.T. and F.A.O.F.; visualization, I.-S.L. and F.A.O.F.; supervision, R.J.A.d.S., I.-S.L. and M.G.T.; project administration, M.G.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** The author GTAP would like to thank Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) for their support (grant number 88881.188598/2018-01). The authors FAOF and RJAS acknowledge the support given by Fundação para a Ciência e a Tecnologia (FCT), through the projects UIDB/00481/2020 and UIDP/00481/2020. In addition, CENTRO-01-0145- FEDER-022083—Centro Portugal Regional Operational Programme (Centro2020) under the PORTU-GAL 2020 Partnership Agreement through the European Regional Development Fund.

**Data Availability Statement:** All data is reported within the article.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Mathematical and Computational Applications* Editorial Office E-mail: mca@mdpi.com www.mdpi.com/journal/mca

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel: +41 61 683 77 34

www.mdpi.com ISBN 978-3-0365-6757-0