Compressible Neo Hookean Model Approach

General
1

Dear all,

I have a code which models the Incompressible Neo Hookean Model but now I want to change it to a nearly compressible model and want to monitor the effects of different bulk modulus. I guess I should solve the weak form not with a mixed formulation. What changes should I do, can u help me?

2 

from dolfin import *
from mshr import *
import matplotlib.pyplot as plt
import matplotlib.tri as tri
import numpy as np
import matplotlib.patches as patches
import ufl

# Parameters and settings
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, "eliminate_zeros": True,
               "precompute_basis_const": True, "precompute_ip_const": True}

lpx = 100
lpy = 10
lpz = 10

n_x = 20
n_y = 4
n_z = 4

mesh = BoxMesh(Point(0, 0, 0), Point(lpx, lpy, lpz), n_x, n_y, n_z)

# Define elements for displacement field) and pressure field
V_element = VectorElement("CG", mesh.ufl_cell(), 2)  # Displacement field
P_element = FiniteElement("CG", mesh.ufl_cell(), 1)  # Pressure field

# Mixed Function Space
W = FunctionSpace(mesh, MixedElement([V_element, P_element]))


def boundary_left(x, on_boundary):
    return on_boundary and near(x[0], 0)


def boundary_right(x, on_boundary):
    return on_boundary and near(x[0], lpx)


def boundary_bottom(x, on_boundary):
    return on_boundary and near(x[2], 0)


bc_bottom = DirichletBC(W.sub(0).sub(2), Constant(0.0), boundary_bottom)
bc_left = DirichletBC(W.sub(0).sub(0), Constant(0.0), boundary_left)

# Material parameters
mu = 10
nu = 0.499
Kappa = 2 * mu

# Define trial and test functions
w = Function(W)
(u, p) = split(w)
(v, q) = TestFunctions(W)

# Kinematics
d = len(u)
I = Identity(d)

# Renamed
F_mat = I + grad(u)
J_mat = det(F_mat)

# Incompressible Neo-Hookean energy density function
#psi = (mu / 2) * (tr(F_mat.T * F_mat) - 3) - Kappa * p * (J_mat - 1)
#Pi = psi * dx

hydraulic_pressure = Kappa * (J_mat - 1)


def neo_hookean_PK1_stress(F, p, mu):
    dim = F.ufl_shape[0]
    I = Identity(dim)
    J = det(F)
    return mu * (F - inv(F).T) - p * J * inv(F).T


degree = 2  # Example degree, adjust as necessary
dx = Measure('dx', domain=mesh, metadata={'quadrature_degree': degree})

#pk2 = neo_hookean_PK1_stress(F_mat, p, mu)
pk2 = mu * I - p * inv(J_mat**(-2/3)*F_mat.T * F_mat)
F1 = inner(pk2, dot(F_mat.T, grad(v))) * dx
F2 = hydraulic_pressure * q * dx
F_total = F1 + F2  # Total weak form combining displacement and pressure fields

Jacob = derivative(F_total, w, TrialFunction(W))


# Strain and Stres Tensors Calculations
def green_lagrange_strain(F):
    C = F.T * F
    I = Identity(3)
    return 0.5 * (C - I)


def neo_hookean_cauchy_stress(F, p, mu):
    b = F * F.T  # Left Cauchy-Green deformation tensor (3x3 in 3D)
    I = Identity(3)  # 3x3 Identity matrix
    return mu * (b - I) - p * I



def neo_hookean_PK1_stress(F, p, mu):
    FinvT = inv(F).T
    return mu * (F - FinvT) - p * FinvT


def second_PK_stress(F, P):
    F_inv = inv(F)
    return F_inv * P


# Stretch values from 1.0 to 2.0 in steps of 0.1
T_11_values = []
S_11_values = []
stretch_values = []

for stretch in [1.0 + i * 0.1 for i in range(41)]:
    bc_right = DirichletBC(W.sub(0).sub(0), Constant((stretch - 1.0) * lpx), boundary_right)
    bcs = [bc_left, bc_right, bc_bottom]

    problem = NonlinearVariationalProblem(F_total, w, bcs, J=Jacob)
    solver = NonlinearVariationalSolver(problem)
    solver.solve()

    (u, p) = w.split()

    F_val = grad(u) + Identity(3)
    J_val =det(F_val)
    P_val = neo_hookean_PK1_stress(F_val, p, mu)
    
    S_tensor = mu * I - p * inv(J_mat ** (-2 / 3) * F_val.T * F_val)
    # Compute and project tensors
    #S_tensor = second_PK_stress(F_val, P_val)
    T_tensor = neo_hookean_cauchy_stress(F_val, p, mu)

    S_tensor_field = project(S_tensor, TensorFunctionSpace(mesh, "DG", 0))
    T_tensor_field = project(T_tensor, TensorFunctionSpace(mesh, "DG", 0))

    # Arbitrary point for calculation the stress values

    x_coord = 100  # x coordinate
    y_coord = 5  # y coordinate
    z_coord = 5
    point = Point(x_coord, y_coord, z_coord)

    # Evaluate tensors at the specified point
    S_at_point = S_tensor_field(point)
    T_at_point = T_tensor_field(point)

    S_11_values.append(S_at_point[0])
    T_11_values.append(T_at_point[0])
    stretch_values.append(stretch)

    # Plots

    plt.figure()
    plt.plot(stretch_values, S_11_values, 'o-')
    plt.xlabel('Stretch Ratio (λ)')
    plt.ylabel('S_11')
    plt.title('S_11 vs Stretch Ratio (λ) for Neo-Hookean Material')
    plt.grid(True)
    plt.savefig('NeoHookeanModel_S11_vs_Stretch.png')
    plt.show()
3

Compressible Neo-Hookean model is covered here with two examples of possible models:
https://bleyerj.github.io/comet-fenicsx/intro/hyperelasticity/hyperelasticity.html
and here:
https://jsdokken.com/dolfinx-tutorial/chapter2/hyperelasticity.html