Hyperelasticity: uniaxial tension of Gent material

An incompressible Gent model is written as

\psi=-\frac{\mu}2J_m\ln\left(1-\frac{I_1-3}{J_m}\right)

Given a Cartesian coordinate, a uniaxial tension is applied along x_1 direction. The deformation gradient is assumed as

\bold{F}=\begin{bmatrix}
\lambda_1&&\\
&\lambda_0&\\
&&\lambda_0
\end{bmatrix}

Where stretch along x_2,x_3 should be same, \lambda_2=\lambda_3=\lambda_0, due to isotropy. We can therefore calculate C matrix and invariant

\bold{C}=\bold{F}\bold{F}^T=\begin{bmatrix}\lambda^2_1&&\\&\lambda^2_0&\\&&\lambda^2_0\end{bmatrix},\ I_1=tr(\bold{C})=\lambda^2_1+2\lambda^2_0

Stress is calculated by

P_{iA}=\frac{\partial \psi}{\partial F_{iA}}

So we will get

P_1=\frac{\partial\psi}{\partial\lambda_1}=\frac{\mu}2\frac{J^2_m}{J_m-(\lambda^2_1+2\lambda^2_0-3)}\left(2\lambda_1+4\lambda_0\frac{d\lambda_0}{d\lambda_1}\right)

and

P_2=\frac{\partial\psi}{\partial\lambda_0}=\frac{\mu}2\frac{J^2_m}{J_m-(\lambda^2_1+2\lambda^2_0-3)}\left(4\lambda_0+2\lambda_1\frac{d\lambda_1}{d\lambda_0}\right)

Incompressibility applies constraints on deformation gradient

\det(\bold{F})=1\Rightarrow\lambda_1=\frac{1}{\lambda^2_0}

Now we want to replace \lambda_0 in stress with \lambda_1

\begin{align*}
P_1&=\mu\frac{J^2_m}{J_m-(\lambda^2_1+2\lambda^{-1}_1-3)}\left(\lambda_1-\lambda^{-2}_1\right)\\
P_2&=2\mu\frac{J^2_m}{J_m-(\lambda^2_1+2\lambda^{-1}_1-3)}\left(\lambda^{-1/2}_1-\lambda^{-5/2}_1\right)\\
\end{align*}

A stress-strain curve will help us understand role of J_m

clear
clc
close all

Jm = [1.5 2 4 6 8];
lambda = linspace(1,7,2000);

for i = 1:length(Jm)
    [lam, I1, psi, P1, P2] = gent(lambda,Jm(i));

    loglog(lam,P1,'linewidth',2);hold on
end
xlabel('\lambda_1')
ylabel('P_1')
legend('J_m = 1.5','J_m = 2','J_m = 4','J_m = 6','J_m = 8','location','northwest')
grid on

function [lam, I1, psi, P1, P2] = gent(lambda,Jm)
P1 = Jm^2./(Jm-(lambda.^2+2*lambda.^(-1)-3)).*(lambda-lambda.^(-2));
P2 = 2*Jm^2./(Jm-(lambda.^2+2*lambda.^(-1)-3)).*(lambda.^(-1/2)-lambda.^(-5/2));
I1 = lambda.^2+2*lambda.^(-1);
psi = -0.5*Jm*log(1-(I1-3)./Jm);

% x in log(x) must be positive
II1 = (I1-3)/Jm;
index = find(II1<=1);
index(1) = [];
lam = lambda(index);
I1 = I1(index);
psi = psi(index);
P1 = P1(index);
P2 = P2(index);
end