Mooney Rivlin Material Model Code

General
1

Dear all, I want to code the Mooney Rivlin Material model in fenics and calculate different stress values for different stretch values. My problem is that my PK2_11 values are going to negative values but I am applying uniaxial tension to my domain. What is my mistake could you please help me ? Attached you can find my code:


from fenics import *
import matplotlib.pyplot as plt
import matplotlib.tri as tri

# Geometry and mesh
lpx, lpy = 100, 10
mesh = RectangleMesh(Point(0, 0), Point(lpx, lpy), 10, 10)

V_element = VectorElement("CG", mesh.ufl_cell(), 2)  # Displacement field

# Function space
V = FunctionSpace(mesh, V_element)

# Boundary conditions
def left(x, on_boundary):
    return on_boundary and near(x[0], 0)

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

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

# Material parameters for Mooney-Rivlin model
C1, C2,  = 1.113, 0.12 # Example values
mu =10
Kappa = 10*mu

# Test and trial functions
u = Function(V)
du = TrialFunction(V)
v = TestFunction(V)

# Kinematics
d = len(u)
I = Identity(d)
F = I + grad(u)
C = F.T * F
Ic = tr(C)
J = det(F)

# Invariants of the deformation
I1 = tr(C)
I2 = 0.5 * (I1**2 - tr(C * C))

# Mooney-Rivlin energy function
Psi_MR = C1 * (I1 - 3) + C2 * (I2 - 3) - Kappa/2 * (ln(J))**2
Pi = Psi_MR * dx

# Compute first variation of Pi (the derivative of Pi with respect to u in the direction of Test_u)
F_variational = derivative(Pi, u, v)

# Compute Jacobian of F
J_variational = derivative(F_variational, u, du)

# Apply boundary conditions
bc_left = DirichletBC(V.sub(0), Constant(0.0), left)
bc_bottom = DirichletBC(V.sub(1), Constant(0.0), bottom)

E_11_values = []
S_11_values = []
T_11_values = []
stretch_values = []

# Loop over stretch values
for stretch in [1.0 + i * 0.1 for i in range(51)]:
    # Update right boundary condition
    bc_right = DirichletBC(V.sub(0), Constant(stretch * lpx - lpx), right)
    bcs = [bc_left, bc_right, bc_bottom]

    # Solve the problem
    problem = NonlinearVariationalProblem(F_variational, u, bcs, J=J_variational)
    solver = NonlinearVariationalSolver(problem)
    solver.solve()


    # Plot the deformed configuration
    p = plot(u, mode="displacement")
    plt.colorbar(p)
    plt.title(f"Deformed configuration at stretch = {stretch:.2f}")
    plt.savefig(f" MooneyRivlinModel_deformed_configuration_{stretch:.2f}.png")
    plt.close()

    def Cauchy_stress_MR(F, C1, C2):
        I = Identity(2)
        B = F * F.T  
        J = det(F)

        stress_tensor = C1 * I + 2 * C2 * B - (C1 + 2 * C2) * J**(-2/3) * B*B

        return stress_tensor


    def PK1_stress_MR(F, C1, C2):
        I = Identity(2)
        J = det(F)
        B = F * F.T  # Left Cauchy-Green deformation tensor
        FinvT = inv(F).T  # F'nin tersinin transpozesi

        # Mooney-Rivlin modeline göre PK1 gerilme tensörünün hesaplanması
        PK1 = F * (C1 * I + 2 * C2 * B - (C1 + 2 * C2) * J**(-2/3) * B*B)

        return PK1
    
    def PK2_stress(F, P):
        F_inv = inv(F)
        return  F_inv * P
    
    def green_lagrange_strain(F):
        C = F.T * F
        I = Identity(2)
        return 0.5 * (C - I)
    
    F_val = grad(u) + Identity(2)
    P_val = PK1_stress_MR(F_val, C1, C2)

    S_tensor = PK2_stress(F_val, P_val)
    E_tensor = green_lagrange_strain(F_val)
    T_tensor = Cauchy_stress_MR(F_val, C1, C2)

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

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

    S_at_point = S_tensor_field(point)
    E_at_point = E_tensor_field(point)
    T_at_point = T_tensor_field(point)

    print(f"#######################################              STRETCH: {stretch}            ###########################################")
   
    print(f"Cauchy Stress Tensor at ({x_coord}, {y_coord}): {T_at_point}")
    print(f"PK2 Stress Tensor at ({x_coord}, {y_coord}): {S_at_point}")
    print(f"GL Strain Tensor at ({x_coord}, {y_coord}): {E_at_point}")
    print("________________________________________________________________________________________________________________________________")
    print("********************************************************************************************************************************")

    E_11_values.append(E_at_point[0])
    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 (λ)')
    plt.ylabel('S_11')
    plt.title('MooneyRivlinMaterialModel S_11 vs Stretch (λ)')
    plt.grid(True)
    plt.savefig('MooneyRivlinMaterialModel S11 vs Stretch.png')
    plt.show()

    plt.figure()
    plt.plot(stretch_values, T_11_values, 'o-')
    plt.xlabel('stretch')
    plt.ylabel('T_11')
    plt.title('MooneyRivlinMaterialModel T_11 vs Stretch (λ)')
    plt.grid(True)
    plt.savefig('MooneyRivlinMaterialModel T_11 vs Stretch.png')
    plt.show()
2

Shouldn’t the last term be + instead of -?

1 Like
3

You are right but when I change it to + , nothing changes

4

I don’t have time to actually run the script, but I don’t think you’re updating J_variational inside the loop. I’m pretty sure that J_variational stays the same as you defined outside the loop.

5

Hello,

One thing that seems to be wrong in your implementation is the way you have applied the uniaxial tension. Note that the stretch should be applied uniformly on your domain (NOT just on the the right edge!). In addition, the edge on the left, should NOT be fixed in both x and y direction. You need to fix it only in y direction. If you apply the stretch uniformly on the boundaries of the domain, the left edge is automatically fixed in x direction as well. Lets take a closer look:

First, you can store all values for the boundary conditions in an array:

# Defining Stretches
stretch = [0.15,0.2,0.25,...]
BC = []

for x in range(len(stretch)):
              N = stretch[x] - 1.0
              BC.append(N)

Lets say the domain is a unit cube:

# Defining the mesh which is a single hexahedron element
mesh = UnitCubeMesh.create(1,1,1,CellType.Type.hexahedron)

You can define the upper face of the cube as:

tol = 1E-14
# Define boundaries
def FRONT(x, on_boundary):
              return on_boundary and abs(x[2] - 1.0) < tol

This is how you can define the uniform uniaxial tension in x direction:

  def border(x, on_boundary):
                return on_boundary
  
  bound_x =  Expression(("t*x[0]"), degree=1, t=0)

This is the loop for solving the problem:

for i in range(len(BC)):

          bound_x.t = BC[i]

          bc_x = DirichletBC(W.sub(1).sub(0), bound_x, border)
          bc_front = DirichletBC(W.sub(1).sub(2), Constant((0)), FRONT)

          bc_all = [bc_x,bc_front]
          .
          .
          .

Not that how the boundary conditions should be applied to simulate the uniaxial tension.
I would suggest to have a look at my implementation and validation of the uniaxial tension on a unit cube for Neo-Hookean material model HERE

2 Likes
6

Dear @Leo ,

Could you please elaborate why you need to apply stretch uniformly on the domain?
If you were to perform uniaxial tension test in real life, you basically stretch from left and right edge. How is it different in FEM setting?

Thanks!

Best,

7

Hello,

The analytical solution for uniaxial tension is based on homogeneous deformation. With that being said, it is assumed that every two lines and planes which are parallel before deformation, remain parallel in deformed configuration (Rreference)

Here, we need to avoid applying any restriction that violates this assumption. It should be noted that by fixing the left face of the cube in all directions and applying the uniaxial stretch on the right face, this assumption does not hold true anymore (If you try this boundary conditions, you will see that the deformed cube is NOT a cube anymore and instead it will end up a pyramid shape).

In general, the domain, should have enough freedom to deform in a way that the deformation is homogenous. When it comes to validation of the FEM results with the analytical solution, this is very important to keep this in mind to apply the correct boundary conditions.

1 Like
8

That makes sense. Thank you so much for the explanation!

1 Like