# In the Realm of Hyperelasticity: Constitutive Modeling of Rubber using Lion's Model

M Usman Sikandar, Rameez Ul Haq, Ali Mehrasa, Göktuğ Tekcan

► To cite this version:

M Usman Sikandar, Rameez Ul Haq, Ali Mehrasa, Göktuğ Tekcan. In the Realm of Hyperelasticity: Constitutive Modeling of Rubber using Lion's Model. EKTAM, Gazi University, Ankara, Turkey; METU, Aerospace Engineering, Turkey; METU, Biomedical Engineering, Turkey. 2023. [hal-04328089](https://hal.science/hal-04328089v1)

HAL Id: hal-04328089

https://hal.science/hal-04328089v1

Submitted on 6 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# In the Realm of Hyperelasticity: Constitutive Modeling of Rubber using Lion's Model

M. Usman Sikandar¹, Rameez ul Haq², Ali Mehrasa³, Göktuğ Tekcan²

1. Additive Manufacturing Technologies Application & Research Center, EKTAM, Ankara, Turkey

2. Middle East Technical University, Department of Aerospace Engineering, Ankara, Turkey

3. Middle East Technical University, Department of Biomedical Engineering, Ankara, Turkey

## ABSTRACT

This study presents a comprehensive approach to rubber hyperelasticity modeling using the Lion model. Stress tensors, specifically the first and second Piola-Kirchhoff stress tensors (P and S), were derived from the Lion model's free energy function (Eq. 7). Leveraging the deformation gradient F, Cauchy (σ), and Kirchhoff (τ) stress tensors were subsequently obtained. To bridge theory with experimentation, a MATLAB code was meticulously developed, employing optimization techniques to reconcile experimental and analytical data for the precise determination of material parameters. The culmination of this effort resulted in the acquisition of stress-stretch curves, providing a thorough evaluation of the Lion model's performance in representing rubber hyperelasticity.

**Keywords:** Hyperelasticity Modeling, Lion Model, Constitutive Modeling, Stress-Stretch Curves, Rubber

## 1. Introduction

### 1.1 Continuum Mechanics

The mechanical behavior of materials that are treated as a continuous mass rather than as discrete particles is the focus of the branch of mechanics known as continuum mechanics. The first to develop such models in the 19th century was the French mathematician Augustin-Louis Cauchy. The mathematical foundation for analyzing large-scale forces and deformations in materials is based on the idea of a continuum. Despite the fact that materials are made up of discrete atoms and molecules that are separated from one another by void space, microscopic cracks, and crystallographic flaws, physical phenomena are frequently modeled by taking into account a substance that is dispersed throughout a particular region of space. A continuum is a body that can continuously be divided into infinitesimal elements with locally specified physical characteristics at any given position. Therefore, continuous functions can be used to characterize the bulk material's properties, and calculus can be used to study how these properties change over time
[1]. Two additional independent assumptions are frequently used in the study of continuum mechanics in addition to the continuity assumption. These are isotropy and homogeneity (identical properties at all locations) and isotropy (properties invariant about direction).

## 1.2 Hyperelasticity

A type of constitutive model for optimally elastic material for which the stress–strain relationship is obtained from a strain energy function is known as a hyperelastic material[1]. The observed behavior of many materials is not adequately captured by linear elastic models; rubber, for instance, has a non-linearly elastic, isotropic, and incompressible stress-strain relationship. It is possible to simulate the stress-strain behavior of such materials using hyperelasticity [2].

## 1.3 Hyperelastic Material Models

The techniques to model rubber as a hyperelastic material can be categorized into two primary groups: (i) phenomenological approaches and (ii) micromechanical-based models. The stress-strain curves obtained from uniaxial, equibiaxial, and pure shear experiments serve as the foundation for the phenomenological models and are fit by invariant or principle stretch-based mathematical formulas for the free energy function to describe the mechanical characterization of rubber-like materials. Contrary to phenomenological models, micromechanical models use the physical and geometrical features at the molecular level of the material microstructure [3].

## 1.4 Preliminaries

The following equations from (Dal, 2019) and (Dal et al., 2021) show the mechanical response of an elastic body is characterized by the deformation map $\varphi(X,t)$ which describes the transformation of the material points $X$ onto spatial points $x=\varphi_t(X)$ at time $t$. The deformation gradient $F$ is given by:

$$F: T_{x} \in \mathcal{B}_{0} \rightarrow T_{x} \in \mathcal{B}, \quad F := \frac{\partial \varphi_{t}(X)}{\partial X} \qquad (1)$$

is considered as the mapping between the tangent spaces in Lagrangian $(T_{x}\mathcal{B}_{0})$ and Eulerian $(T_{x}\mathcal{B})$ configurations. Subsequently, $F^{-T}$ characterizes the mapping between the cotangent spaces in Lagrangian $(T_{x}^{*}\mathcal{B}_{0})$ and Eulerian $(T_{x}^{*}\mathcal{B})$ configurations.
$$F^T: T^*_{x_0} \to T^*_{x_0} \mathcal{B}, \quad n = F^{-T} N \qquad (2)$$

Where $N$ and $n$ denote unit normal vectors of Lagrangian and Eulerian configurations, respectively.

The right and inverse of left Cauchy-Green tensors are defined as:

$$C = F^T g F \quad \text{and} \quad F = F^{-T} G F^{-1} \qquad (3)$$

where $G$ and $g$ denote the covariant reference and current metric tensors [3,4].

The energy stored in an isotropic material is governed by principle stretches or invariants by the principles of material objectivity and material frame indifference. The right Cauchy-Green tensor $C$ has the following eigenvalue decompositions:

$$C := \sum_{i=1}^{3} \lambda^2 N^a \otimes N^a \qquad (4)$$

where $\lambda_i$ and $N^a$ represent the principal stretches and directions, respectively [3]. Three isotropic invariants of the right Cauchy Green tensor can then be written as:

$$I_1 = \lambda_1^2 + \lambda_2^2 + \lambda_3^2, \quad I_2 = \lambda_1^2 \lambda_2^2 + \lambda_2^2 \lambda_3^2 + \lambda_3^2 \lambda_1^2, \quad I_3 = \lambda_1^2 \lambda_2^2 \lambda_3^2 \qquad (5)$$

The above-mentioned tensors can also be geometrically interpreted from Figure 1[4].

<figure><img src="image_2.png" /><figcaption>Figure 1. Geometrical representation of metric and stress tensors [4]</figcaption></figure>

Geometrical representation of the Cauchy ($\sigma$), Kirchhoff ($\tau$), first ($P$), and second Piola-Kirchhoff stress tensors ($S$) are also depicted in Figure 2 [4]. Push forward and pull back of these tensors will be used in the derivation section.
<figure><img src="image_3.png" /><figcaption>Figure 2. Geometrical representation of stress tensors [4]</figcaption></figure>

## 1.5 Free Energy Function and Lion Model

The finite elasticity of a particular class of materials is determined by the free energy function $\psi$. A category of phenomenological hyperelastic material models is based on principal invariants. The generalized Rivlin model is the cornerstone of many invariant-based models, the closed-form representation of which is:

$$ \psi(\mathbf{F}) = \sum_{i=0}^{M_1} \sum_{j=0}^{M_2} C_{ij}(I_1 - 3)^i (I_2 - 3)^j \quad (6) $$

Lion model, a special case of the Rivlin model, considers linear and fifth-order terms for $I_1$ and one linear term for $I_2$ such that:

$$ \psi(\mathbf{F}) = C_{10}(I_1 - 3) + C_{01}(I_2 - 3) + C_{03}(I_1 - 3)^5 \quad (7) $$

Material parameters $C_{10}$, $C_{01}$, and $C_{03}$ are determined by optimization using experimental data.

## 1.6 Biaxial Tension Experiments

The homogeneous deformation state for uniaxial tension and biaxial tension of a unit-cube element can be depicted in Figure 2 where the incompressibility assumption is imposed on the deformation gradient.
<figure><img src="image_4.png" /><figcaption>Figure 3. Homogeneous deformation state for uniaxial and biaxial tensions of a unit-cube element [5]</figcaption></figure>

Biaxial tension experiments are commonly used in the mechanical characterization of rubberlike materials and become vital in parameter optimization for constitutive models. In a real-world application, a rubber part is rarely subjected to a unidirectional loading. Thus, the constitutive models are mainly assessed based on the performance of their fits to biaxial tests. Considering more elaborate models, parameters that fit uniaxial tension, equibiaxial tension, and pure shear experiments simultaneously also predict the biaxial test results. The deformation gradient for the biaxial test of a rubber part with incompressibility assumption reads:

$$[F]_{ij} = \begin{bmatrix} \lambda_1 & 0 & 0 \\ 0 & \lambda_2 & 0 \\ 0 & 0 & \frac{1}{\lambda_1 \lambda_2} \end{bmatrix} \qquad (8)$$

where the lambdas ratio defines as the biaxiality ratio. A representative example of the dataset for the work is shown in Figure 4 [5].
<figure><img src="image_5.png" /><figcaption>Figure 4. The biaxial test dataset for different biaxiality ratios (n), [red] n = 1, [green] n = 1.5, [blue] n = 2, and [magenta] n = 2.5 respectively. The uniaxial tension behavior is also shown in the figure (black curve) for comparison with the biaxial test results [5]</figcaption></figure>

## 1.7 Objectives

Stress tensors including **P** and **S**, first and second Piola Kirchhoff stress tensors were first driven from the free energy function of the Lion model (Eq. 7). Then by using the deformation gradient **F**, **σ** and **τ**, Cauchy, and Kirchhoff stress tensors were driven. In the next step, MATLAB code was developed, and optimization was used to fit the experimental to analytical data to determine material parameters. Finally, the stress-stretch curves were obtained for evaluation.

## 2. Derivations

The general equation for the second Piola-Kirchhoff stress can be directly driven by the free energy function:

$$S = 2\partial_c\psi = 2\left(\frac{\partial\psi}{\partial I_1}\frac{\partial I_1}{\partial C} + \frac{\partial\psi}{\partial I_2}\frac{\partial I_2}{\partial C} + \frac{\partial\psi}{\partial J}\frac{\partial J}{\partial C}\right) \quad (9)$$

Where:

$$\frac{\partial I_1}{\partial C} = 1, \qquad \frac{\partial I_2}{\partial C} = I_1 1 - C, \qquad \frac{\partial J}{\partial C} = \frac{1}{2} J C^{-1} \quad (10)$$

Substitution of Eq. 10 into Eq. 9 and rearranging the terms gives the invariant-based representation of the second Piola-Kirchhoff stress for isotropic hyperelastic solid:

$$S = 2(c_1 + I_1 c_2)1 - 2c_2C - pC^{-1} \quad (11)$$
where

$$c_1 = \frac{\partial \psi}{\partial I_1}, \qquad c_2 = \frac{\partial \psi}{\partial I_2}, \qquad p = -J \partial_j \psi \qquad (12)$$

Using the deformation gradient, and **P** = **FS**, the first Piola-Kirchhoff stress tensor components can be driven:

$$\begin{align*}
\text{Uniaxial} \qquad P_1 &= 2 \left(c_1 + \frac{c_2}{\lambda}\right) \left(\lambda - \frac{1}{\lambda^2}\right) \qquad P_2 &= 2 \left(c_1 + c_2 \lambda_1^2\right) \left(\lambda_1 - \frac{1}{\lambda_1^3 \lambda_2^2}\right) \tag{13} \\
\text{Biaxial} \qquad P_1 &= 2 \left(c_1 + c_2 \lambda_2^2\right) \left(\lambda_1 - \frac{1}{\lambda_1^3 \lambda_2^2}\right) \qquad P_2 &= 2 \left(c_1 + c_2 \lambda_1^2\right) \left(\lambda_1 - \frac{1}{\lambda_1^2 \lambda_2^3}\right)
\end{align*}$$

Using Eq. 7 and Eq. 12:

$$c_1 = \frac{\partial \psi}{\partial I_1} = C_{10} + 5C_{03}(I_1 - 3)^4, \qquad c_2 = \frac{\partial \psi}{\partial I_2} = C_{01} \qquad (14)$$

### a) Biaxial Tension

The first step, which was described in the objectives, is to derive an analytical formulation of the first Piola-Kirchhoff stress tensor. The deformation gradient and the first Piola-Kirchhoff stress tensors are given as:

$$F = \begin{bmatrix} \lambda_1 & 0 & 0 \\ 0 & \lambda_2 & 0 \\ 0 & 0 & \frac{1}{\lambda_1 \lambda_2} \end{bmatrix}, \qquad P = \begin{bmatrix} P_1 & 0 & 0 \\ 0 & P_2 & 0 \\ 0 & 0 & 0 \end{bmatrix} \qquad (15)$$

Substitution of Eq. 13 into Eq. 12 gives the components of the first Piola-Kirchhoff stress tensor in terms of material parameters:
$$P_1 = 2 [C_{10} + 5C_{03}(I_1 - 3)^4 + C_{01}\lambda_1^2] \left(\lambda_1 - \frac{1}{\lambda_1^3\lambda_2^2}\right) \qquad (16)$$

$$P_2 = 2 [C_{10} + 5C_{03}(I_1 - 3)^4 + C_{01}\lambda_1^2] \left(\lambda_1 - \frac{1}{\lambda_1^2\lambda_2^3}\right)$$

where the first and the second invariants are:

$$I_1 = \lambda_1^2 + \lambda_2^2 + \frac{1}{\lambda_1^2\lambda_2^2}, \qquad I_2 = \lambda_1^2\lambda_2^2 + \frac{1}{\lambda_1^2} + \frac{1}{\lambda_2^2} \qquad (17)$$

The next step is to determine the analytical formulation of the second Piola Kirchhoff stress tensor. The second Piola-Kirchhoff stress tensor is defined as $S = F^{-1}P$, where:

$$F^{-1} = \begin{bmatrix} \frac{1}{\lambda_1} & 0 & 0 \\ 0 & \frac{1}{\lambda_2} & 0 \\ 0 & 0 & \lambda_1\lambda_2 \end{bmatrix}, \qquad P = \begin{bmatrix} P_1 & 0 & 0 \\ 0 & P_2 & 0 \\ 0 & 0 & 0 \end{bmatrix} \qquad (18)$$

So, the components of the second Piola-Kirchhoff stress tensor will read:

$$S_1 = \frac{P_1}{\lambda_1}, \qquad S_2 = \frac{P_2}{\lambda_2}, \qquad (19)$$

The Kirchhoff stress tensor is given by $\tau = PF^T$ where:

$$F^T = \begin{bmatrix} \lambda_1 & 0 & 0 \\ 0 & \lambda_2 & 0 \\ 0 & 0 & \frac{1}{\lambda_1\lambda_2} \end{bmatrix}, \qquad P = \begin{bmatrix} P_1 & 0 & 0 \\ 0 & P_2 & 0 \\ 0 & 0 & 0 \end{bmatrix} \qquad (20)$$

So, the components of the Kirchhoff stress tensor will read:

$$\tau_1 = P_1\lambda_1, \qquad \tau_2 = P_2\lambda_2 \qquad (21)$$

Since $J=1$ (the jacobian of the deformation gradient = 1 because the material is incompressible and the deformation is volume preserving which is also confirmed from the components of the
deformation gradient) and $\tau = J\sigma$, the Kirchhoff stress and the Cauchy stress have the same value.

## b) Biaxial Tension

Now we will do the same steps for the uniaxial case. The deformation gradient and the first Piola-Kirchhoff stress tensors are given as:

$$ F = \begin{bmatrix} \lambda & 0 & 0 \\ 0 & 1/\sqrt{\lambda} & 0 \\ 0 & 0 & 1/\sqrt{\lambda} \end{bmatrix}, \qquad P = \begin{bmatrix} P_1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \qquad (22) $$

Substitution of Eq. 13 into Eq. 12 gives the component of the first Piola-Kirchhoff stress tensor in terms of material parameters:

$$ P_1 = 2 [C_{10} + 5C_{03}(I_1 - 3)^4 + C_{01}\lambda_2^2] \left(\lambda - \frac{1}{\lambda^2}\right) \qquad (23) $$

where the first and the second invariants are:

$$ I_1 = \lambda^2 + \frac{2}{\lambda}, \qquad I_2 = 2\lambda + \frac{1}{\lambda^2} \qquad (24) $$

Determining the analytical formulation of the second Piola Kirchhoff stress tensor using $S = F^{-1}P$
Where:

$$ F^{-1} = \begin{bmatrix} \frac{1}{\lambda} & 0 & 0 \\ 0 & \sqrt{\lambda} & 0 \\ 0 & 0 & \sqrt{\lambda} \end{bmatrix}, \qquad P = \begin{bmatrix} P_1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \qquad (25) $$

It gives the component of the second Piola Kirchhoff stress tensor:

$$ S_1 = \frac{P_1}{\lambda} \qquad (26) $$

And using $\tau = PF^T$, the component of the Kirchhoff stress is:

$$ \tau_1 = P_1 \lambda \qquad (27) $$
Again, since $J=1$ and $\tau = J\sigma$, the *Kirchhoff* stress and the *Cauchy* stress have the same value.

### 3. MATLAB Code & Optimization

All the equations driven in the previous section were implemented in the MATLAB. Non-linear least squares regression is employed for the optimization. For the case of $n=1$ where $P_{11} = P_{22}$ the code was developed as written below:

```matlab
% For n = 1 case (P_11 = P_22)
function err = n1LSQfunphi)
n = 1;
P_11 = xlsread('C:\Rubber Lion\d10.xlsx', 'n1', 'B1:B100');
lam_1 = xlsread('C:\Rubber Lion\d10.xlsx', 'n1', 'A1:A100'); lam_2
= n * lam_1;
I1 = (lam_1).^2 + (lam_2).^2 + 1./((lam_1).^2.* (lam_2).^2);
I2 = ((lam_1).^2.* (lam_2).^2) + 1./((lam_1).^2) +
1./((lam_2).^2);

% Since n = 1, this means P_11 = P_22, hence only P_11 is
computed below
% Now define the unknowns constants as phi
C_10 = phi(1);
C_03 = phi(2);
C_01 = phi(3);
P_11_model = 2*(C_10 + 5*C_03*((I1 - 3).^4) +
C_01*(lam_2).^2).* (lam_1 - (1./((lam_1).^3.* (lam_2).^2)));
err = abs(P_11 - P_11_model);
end
```

And to determine the material parameters $C_{10}$, $C_{01}$, and $C_{03}$

```matlab
phi0 = [1.75,1.75,1.75]; %Initial guess for C_10, C_03 and C_01
phi = lsqnonlin(@p) n1LSQfun(p), phi0)
```

For the case of $n = 2$, two separate sets of values for $C_{10}$, $C_{01}$, and $C_{03}$ were computed by the 'lsqnonlin' function: one by optimizing for $P_{11}$ and the other for $P_{22}$. The code to optimize for $P_{11}$ was written as:

```matlab
% For n = 2 case, and P_11
function err = n2_1LSQfun(alpha)
n = 2;
P_11_2 = xlsread('C:\Rubber Lion\d20.xlsx', 'n2', 'B1:B100');
lam_1 = xlsread('C:\Rubber Lion\d20.xlsx', 'n2', 'A1:A100'); lam_2
```
```matlab
% Since n = 2, this means P_11 and P_22 are different from each other
% Now define the unknowns constants as alpha
C_10 = alpha(1);
C_03 = alpha(2);
C_01 = alpha(3);
P_11_model2 = 2*(C_10 + 5*C_03*((I1 - 3).^4) + C_01*(lam_2).^2).*
    (lam_1 - (1./((lam_1).^3.* (lam_2).^2)));
err = abs(P_11_2 - P_11_model2);
end
```

And to determine the material parameters C₁₀, C₀₁, and C₀₃

```matlab
alpha0 = [1.9,1.9,1.9]; %Initial guess for C_10, C_03 and C_01
alpha = lsqnonlin(@(a) n2_1LSQfun(a), alpha0)
```

The code to optimize for P₂₂ was written as:

```matlab
% For n = 2 case, and P_22
function err = n2_2LSQfun beta
n = 2;
P_22_2 = xlsread('C:\Rubber Lion\d20.xlsx', 'n2', 'C1:C100');
lam_1 = xlsread('C:\Rubber Lion\d20.xlsx', 'n2', 'A1:A100'); lam_2
= n * lam_1;
I1 = (lam_1).^2 + (lam_2).^2 + 1./((lam_1).^2.* (lam_2).^2);
I2 = ((lam_1).^2.* (lam_2).^2) + 1./((lam_1).^2) +
    1./((lam_2).^2);

% Since n = 2, this means P_11 and P_22 are different from each other
% Now define the unknowns constants as beta
C_10 = beta(1);
C_03 = beta(2);
C_01 = beta(3);
P_22_model2 = 2*(C_10 + 5*C_03*((I1 - 3).^4) +
    C_01*(lam_1).^2).*
    (lam_2 - (1./((lam_1).^2.* (lam_2).^3)));
err = abs(P_22_2 - P_22_model2);
End
```

And to determine the material parameters C₁₀, C₀₁, and C₀₃

```matlab
beta0 = [1.75,1.75,1.75]; %Initial guess for C_10, C_03 and C_01
beta = lsqnonlin(@b) n2_2LSQfun(b), beta0)
```
The overall constant values for $n = 2$ cases were then computed using the weighted contribution of $P_{11}$ and $P_{22}$. The code below shows the method of approach.

```R
Err_n2_P11 = abs(P_11_2 - P_11_model2);
Err_tot2_P11 = sum(Err_n2_P11);

Err_n2_P22 = abs(P_22_2 - P_22_model2);
Err_tot2_P22 = sum(Err_n2_P22);

% Now sum the errors for n = 2 case, as shown below
Err_tot2 = Err_tot2_P11 + Err_tot2_P22;

% Now determine the ratios of these errors w.r.t. Err_tot2
Perecent_P11_2 = Err_tot2_P11/Err_tot2;
Perecent_P22_2 = Err_tot2_P22/Err_tot2;

% Now find the overall constants for n=2 case by the weighted method below:
C_10_2 = (Perecent_P11_2 * beta(1)) + (Perecent_P22_2 * alpha(1));
C_03_2 = (Perecent_P11_2 * beta(2)) + (Perecent_P22_2 * alpha(2));
C_01_2 = (Perecent_P11_2 * beta(3)) + (Perecent_P22_2 * alpha(3));
```

The constants obtained from here were then used to compute the values of $P_{11}$ and $P_{22}$.

Then for the last step, simultaneous fits were employed by again using the separate weighted contributions between $n = 1$, and $n = 2$ cases. Only the comparisons between the $P_{11}$ for each case were used in the below code.

```R
Err_n1 = abs(P_11 - P_11_model);
Err_tot1 = sum(Err_n1);

Err_n2_P11 = abs(P_11_2 - P_11_model2);
Err_tot2_P11 = sum(Err_n2_P11);

% Now sum the errors for cases n=1 and n=2
Err_tot = Err_tot1 + Err_tot2_P11;
```
<figure><img src="image_6.png" /></figure>

The general simultaneously fitted constants were obtained: $C_{10} = 0.194$, $C_{03} = 0$, $C_{01} = -0.0049$

# 4. Results & Discussion

In this section, results including stress-stretch curves and errors are presented.

## 4.1 Stress-Stretch Curves

### a) Biaxial

Figure 4 shows the comparison between the experimental data and the data obtained from the analytical model for the case of $n = 1$, where $P_{11} = P_{22}$.

<figure><img src="image_7.png" /><figcaption>Figure 5. Comparative stress-stretch curve for n=1 between the experimental and analytical data</figcaption></figure>
Figure 5 shows the comparison between the experimental data and the data obtained from the analytical model for the case of $n = 2$.

<figure><img src="image_8.png" /><figcaption>Figure 6. Comparative stress-stretch curve for $n=2$ between the experimental and analytical data</figcaption></figure>

Figure 6 shows the values of different types of stresses obtained from the analytical model for the case of $n=1$. Notice that $P_{11} = P_{22}$, $S_{11} = S_{22}$, $\tau_{11} = \tau_{22}$, and $\sigma_{11} = \sigma_{22}$.

<figure><img src="image_9.png" /><figcaption>Figure 7. Values of different types of stresses obtained from the analytical model for the case of $n=1$</figcaption></figure>
Figure 7 shows the values of different types of stresses obtained from the analytical model for the case of $n = 2$. Notice that $P_{11} \neq P_{22}$, $S_{11} \neq S_{22}$, $\tau_{11} \neq \tau_{22}$, and $\sigma_{11} \neq \sigma_{22}$.

<figure><img src="image_10.png" /><figcaption>Figure 8. Values of different types of stresses obtained from the analytical model for the case of n=2</figcaption></figure>

a) Uniaxial

Figure 8 shows the comparison between the experimental data and the data obtained from the analytical model.

<figure><img src="image_11.png" /><figcaption>Figure 9. Comparative stress-stretch curve between the experimental and analytical data</figcaption></figure>
Figure 9 shows the values of different types of stresses obtained from the analytical model.

<figure><img src="image_12.png" /><figcaption>Figure 10. Values of different types of stresses obtained from the analytical model</figcaption></figure>

## 4.2 Errors

Individual fitting for the case of $n = 1$ yielded the material parameters as $C_{10}= 0.1942$, $C_{03}= 0.0000$, $C_{01} = -0.0046$. Figure 10 depicts the comparison of analytical and experimental data in this step.

<figure><img src="image_13.png" /><figcaption>Figure 11. Individual fitting for the case of $n = 1$</figcaption></figure>
Individual fitting for $P_{11}$ for the case of $n = 2$ yielded the material parameters as $C_{10} = 0.1955$, $C_{03} = 0.0000$, $C_{01} = -0.0028$. Individual fitting for $P_{22}$ for the case of $n = 2$ yielded the material parameters as $C_{10} = 0.1900$, $C_{03} = 0.0000$, $C_{01} = -0.0108$. Figure 11 depicts the comparison of analytical and experimental data in this step.

<figure><img src="image_14.png" /><figcaption>Figure 12. Individual fitting for the case of $n = 2$</figcaption></figure>

Errors associated with the individual fitting for the cases of $n = 1$ and $n = 2$ are shown in Figure 12.

<figure><img src="image_15.png" /><figcaption>Figure 13. Errors associated with the individual fitting</figcaption></figure>
Simultaneous fitting for *n* = 2 by separate weighted contributions from *P*<sub>11</sub> and *P*<sub>22</sub> yielded the material parameters as *C*<sub>10</sub> = 0.1935, *C*<sub>03</sub> = 0.0000, and *C*<sub>01</sub> = -0.0058. Figure 13 depicts the comparison of analytical and experimental data in this step.

<figure><img src="image_16.png" /><figcaption>Figure 14. Simultaneous fitting for n=2</figcaption></figure>

Finally, simultaneous fitting by separate weighted contributions from *n* = 1 and *n* = 2 yielded the material parameters as *C*<sub>10</sub> = 0.194, *C*<sub>03</sub> = 0.0000, and *C*<sub>01</sub> = -0.0049. Errors associated with the simultaneous fitting are shown in Figure 14.

<figure><img src="image_17.png" /><figcaption>Figure 15. Errors associated with the simultaneous fitting</figcaption></figure>
Table 1 summarizes the errors associated with the fitting of experimental data and analytical data using optimization.

<table><caption>Table 1. Biaxial tension results for the Lion model</caption><thead><tr><th>Biaxial Tension</th><th>Error</th></tr></thead><tbody><tr><td>Simultaneous fit for n=2</td><td></td></tr><tr><td>P<sub>11</sub></td><td>0.8338</td></tr><tr><td>P<sub>22</sub></td><td>0.2444</td></tr><tr><td>Simultaneous fit for n=1 and n=2</td><td></td></tr><tr><td>P<sub>11</sub></td><td>3.333</td></tr></tbody></table>

## 5. Concluding Remarks

Biaxial experimental data were used to mechanically characterize hyperelastic material using the Lion model. A non-linear least squares regression approach was taken to optimize and fit the experimental and analytically obtained data to determine material parameters. According to the results obtained, although the Lion model accurately models the hyperelastic material behavior in non-Gaussian regime, there are errors in the small-stretch region.

## 6. Funding Statement

A part of this independent work is conducted under two funding sources: the European Union's Horizon 2020 program with Marie Skłodowska-Curie grant agreement no. 10103442 and the Scientific and Technological Research Council of Türkiye (TUBITAK) with grant no. 120C158 for the A2M2TECH project under TUBITAK's 2236/B program.

## 7. References

[1] Wikipedia contributors. (2023, November 6). Continuum mechanics. Wikipedia. https://en.wikipedia.org/wiki/Continuum_mechanics

[2] Muhr, A. H. (2005). Modeling the stress-strain behavior of rubber. Rubber chemistry and technology, 78(3), 391-425.

[3] Dal, H., Açıkgöz, K., & Badienia, Y. (2021). On the performance of isotropic hyperelastic constitutive models for rubber-like materials: A state of the art review. *Applied Mechanics Reviews*, 73(2), 020802.

[4] Dal, H., 2019, "A Quasi-Incompressible and Quasi-Inextensible Element Formulation for Transversely Isotropic Materials: A Mixed Element Formulation for Anisotropic Materials," *Int J Numer Methods Eng*, **117**(1), pp. 118-140.

[5] Badienia, Y., 2019, "A comparative study of the fitting performance of hyperelastic constitutive models" (Master's thesis, Middle East Technical University).