Hello all, I am trying to implement the Incompressible Neo Hookean model into fenics and try to compare it with the analytical results and the results for different bulk modulus. Attached you can see the code that I wrote. My problem is that although my code works, It gives me the same result for different Kappa values. What is wrong with me ?
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
Parameters and settings
parameters[“form_compiler”][“cpp_optimize”] = True
ffc_options = {“optimize”: True, “eliminate_zeros”: True,
“precompute_basis_const”: True, “precompute_ip_const”: True}
Geometry and mesh
lpx = 50
lpy = 10
n_x = 1
n_y = 1
domain = Rectangle(Point(0, 0), Point(lpx, lpy))
mesh = RectangleMesh(Point(0, 0), Point(lpx, lpy), n_x, n_y)
Define elements for displacement (vector field) and pressure (scalar field)
V_element = VectorElement(“CG”, mesh.ufl_cell(), 2) # Displacement field
P_element = FiniteElement(“CG”, mesh.ufl_cell(), 1) # Pressure field
Create a 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):
# Check if a point is at the bottom boundary (y=0)
return on_boundary and near(x[1], 0)
bc_bottom = DirichletBC(W.sub(0).sub(1), Constant(0.0), boundary_bottom)
bc_left = DirichletBC(W.sub(0).sub(0), Constant(0.0), boundary_left)
Sol kenarda hem x hem y için sıfır deplasman
Material parameters
mu = Constant(10.0) # Shear Modulus
nu = 0.499
Kappa = Constant(1000.0) * mu
Define trial and test functions
w = Function(W)
(u, p) = split(w)
(v, q) = TestFunctions(W)
Kinematics
d = len(u)
I = Identity(d)
F = I + grad(u)
J = det(F)
C = F.T * F
Ic = tr(C)
Incompressible Neo-Hookean energy density function
psi = (mu / 2) * (Ic - 3) - Kappa * p * (J - 1)
Pi = psi * dx
Compute the first variation (derivative) of Pi
F = derivative(Pi, w, TestFunction(W))
Compute Jacobian of F
J = derivative(F, w, TrialFunction(W))
def green_lagrange_strain(F):
C = F.T * F
I = Identity(2)
return 0.5 * (C - I)
Stress and strain tensor functions
def neo_hookean_cauchy_stress(F, p, mu):
b = F * F.T
I = Identity(2)
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]
# Define and solve the problem
problem = NonlinearVariationalProblem(F, w, bcs, J=J)
solver = NonlinearVariationalSolver(problem)
solver.solve()
# Extract solutions
(u, p) = w.split()
F_val = grad(u) + Identity(2)
P_val = neo_hookean_PK1_stress(F_val, p, mu)
# 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))
# Example Values
x_coord = 50 # x coordinate
y_coord = 5 # y coordinate
point = Point(x_coord, y_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)
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()
plt.figure()
plt.plot(stretch_values, T_11_values, 'o-')
plt.xlabel('Stretch Ratio (λ)')
plt.ylabel('T_11')
plt.title('T_11 vs Stretch Ratio (λ) for Neo-Hookean Material')
plt.grid(True)
plt.savefig('NeoHookeanModel_T11_vs_Stretch.png')
plt.show()