Dear all, I want to implement the FEM simulation for Mooney Rivlin Model to fenics and compare the results with the analytical solution. But my results are not converges. Below you can find my code and could you please help me for correcting my fault. Thank you in advance.
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
from math import sqrt
# 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
n_x = 50
n_y = 10
domain = Rectangle(Point(0, 0), Point(lpx, lpy))
mesh = RectangleMesh(Point(0, 0), Point(lpx, lpy), n_x, n_y)
# 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[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)
# Mooney-Rivlin mat pars
c1 = 1.1
c2 = 0.12
mu = 10
Kappa = 10 * mu
# Trial and Test Functions
w = Function(W)
(u, p) = split(w)
(v, q) = TestFunctions(W)
d = len(u)
I = Identity(d)
F_mat = I + grad(u) # Deformasyon gradyanı
J_mat = det(F_mat) # Jakobiyen
def mooney_rivlin_PK1_stress(F_mat, p, c1, c2):
J_mat = det(F_mat)
B = F_mat * F_mat.T
I1 = tr(B)
I2 = 0.5 * (I1**2 - tr(B * B))
J23 = J_mat**(-2.0/3)
J43 = J_mat**(-4.0/3)
return 2 * (c1 * J23 * F_mat + c2 * J43 * (I1 * F_mat - (F_mat * F_mat.T) * F_mat) - p * J_mat * inv(F_mat).T )
def mooney_rivlin_PK2_stress(F_mat, P):
F_inv = inv(F_mat) # F'nin tersi
S = F_inv * P # PK2 gerilme tensörü hesaplama
return S
def mooney_rivlin_energy_density(F_mat, p, c1, c2, Kappa, J_mat):
B = F_mat * F_mat.T
I1 = tr(B)
I2 = 0.5 * (I1**2 - tr(B * B))
energy = c1 * (J_mat**(-2.0/3) * I1 - 3) + c2 * (J_mat**(-4.0/3) * I2 - 3)
compressibility_term = (Kappa/2) * (J_mat - 1)**2
return energy + compressibility_term
hydraulic_pressure = Kappa/2 * (J_mat - 1)
psi = mooney_rivlin_energy_density(F_mat, p, c1, c2, Kappa, J_mat)
Pi = psi * dx - dot(p, J_mat - 1) * dx
pk1 = mooney_rivlin_PK1_stress(F_mat, p, c1, c2)
pk2 = mooney_rivlin_PK2_stress(F_mat, pk1)
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
Jacobian = derivative(F_total, w, TrialFunction(W))
P_11_values = []
S_11_values = []
stretch_values = []
for stretch in [1.0 + i * 0.1 for i in range(51)]:
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, Jacobian)
solver = NonlinearVariationalSolver(problem)
solver.solve()
(u, p) = w.split()
F_val = I + grad(u)
P_val = mooney_rivlin_PK1_stress(F_val, p, c1, c2)
S_val = mooney_rivlin_PK2_stress(F_mat, P_val)
P_tensor_field = project(P_val, TensorFunctionSpace(mesh, "DG", 0))
S_tensor_field = project(S_val, TensorFunctionSpace(mesh, "DG", 0))
x_coord = 100 # x coordinate
y_coord = 5 # y coordinate
point = Point(x_coord, y_coord)
P_at_point = P_tensor_field(point)
S_at_point = S_tensor_field(point)
P_11_values.append(P_at_point[0])
S_11_values.append(S_at_point[0])
stretch_values.append(stretch)
plt.figure()
plt.plot(stretch_values, P_11_values, 'o-')
plt.xlabel('Stretch Ratio (λ)')
plt.ylabel('P_11')
plt.title('P_11 vs Stretch Ratio (λ) for Mooney Rivlin Model Material')
plt.grid(True)
plt.savefig('MooneyRivlin_Model_P11_vs_Stretch.png')
plt.show()
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 Mooney Rivlin Model Material')
plt.grid(True)
plt.savefig('MooneyRivlin_Model_S11_vs_Stretch.png')
plt.show()
end
def mooney_rivlin_analytical_stress_P11(stretch, c1, c2):
# Tek eksenli gerinim durumu için P11 bileşenini hesaplar
lambda_x = stretch # x eksenindeki gerinim oranı
lambda_y = lambda_z = 1 / sqrt(lambda_x) # y ve z eksenlerindeki gerinim oranları
# P11 gerilme bileşenini hesaplayan formül
P11 = 2 * c1 * (lambda_x - lambda_y**2) + 2 * c2 * (1 - lambda_y**4) / lambda_x
return P11
def mooney_rivlin_analytical_stress_S11(stretch, c1, c2):
# Tek eksenli gerinim durumu için S11 bileşenini hesaplar
lambda_x = stretch # x eksenindeki gerinim oranı
lambda_y = lambda_z = 1 / sqrt(lambda_x) # y ve z eksenlerindeki gerinim oranları
# P11 gerilme bileşenini hesaplar
P11 = mooney_rivlin_analytical_stress_P11(stretch, c1, c2)
# S11 gerilme bileşenini hesaplayan formül
S11 = P11 / lambda_x
return S11
# Initialize lists for storing stretch ratios and corresponding S_11 values
analytical_stretch_values = []
analytical_S_11_values = []
analytical_P_11_values = []
# Loop through stretch values from 1 to 6 in increments of 0.1
for lam_fem in [1.0 + i * 0.1 for i in range(51)]:
APK1_11 = mooney_rivlin_analytical_stress_P11(lam_fem, c1, c2)
APK2_11= mooney_rivlin_analytical_stress_S11(lam_fem, c1, c2)
analytical_stretch_values.append(lam_fem)
analytical_S_11_values.append(APK2_11)
analytical_P_11_values.append(APK1_11)
plt.figure()
plt.plot(stretch_values, S_11_values, 'o-', label='FEM Solution')
plt.plot(analytical_stretch_values, analytical_S_11_values, 'x-', label='Analytical Solution')
plt.xlabel('Stretch Ratio (λ)')
plt.ylabel('Second Piola-Kirchhoff Stress (PK2_11)')
plt.title('FEM vs Analytical Results for PK2_11 vs Stretch Comparison')
plt.legend()
plt.grid(True)
plt.savefig('Combined_PK2_vs_Stretch.png')
plt.show()
plt.figure()
plt.plot(stretch_values, P_11_values, 'o-', label='FEM Solution')
plt.plot(analytical_stretch_values, analytical_P_11_values, 'x-', label='Analytical Solution')
plt.xlabel('Stretch Ratio (λ)')
plt.ylabel('PK1 Stress (T_11)')
plt.title('FEM vs Analytical Results for PK_11 vs Stretch Comparison')
plt.legend()
plt.grid(True)
plt.savefig('Combined_PK1_vs_Stretch.png')
plt.show()