Can't get the exact same Young's Modulus for simple sample

1

Hello, I want to get the exact same Young’s Modulus as my setting through the load-displacement curve using a simple computation. However, the result I got is slightly different and the Young’s Modulus raised from 2.96GPa to 3.5GPa, Why is that?

Thanks in advance!

Here are the codes:

import os
import shutil
import datetime
import dolfin
import mshr
import math
import multiprocessing
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from fenics import *
from scipy import integrate
%matplotlib inline

parameters['form_compiler']['cpp_optimize'] = True
parameters['form_compiler']['optimize'] = True
set_log_active(False)

mesh = UnitSquareMesh(32, 32)
plot(mesh)

# Materials parameters for matrix.
E = 2.96e3  # Young's modulus MPa
nu = 0.35  # Poisson's ratio
lmbda = E * nu / ((1 + nu) * (1 - 2 * nu))  # Lame constants compression ratio MPa
mu = E / (2 * (1 + nu))  # Lame constants shear modulus MPa
G = E / (2 * (1 + nu))  # Lame constants shear modulus MPa

# Define Space
W = VectorFunctionSpace(mesh, 'CG', 1)
WWW = VectorFunctionSpace(mesh, 'DG', 4)
u, v = TrialFunction(W), TestFunction(W)
T = TensorFunctionSpace(mesh, 'DG', 0)

# Constituive functions
def epsilon(u):
    # Computation of strain tensor.
    return sym(grad(u))

def sigma(u):
    # Computation of Cauchy stress.
    return 2.0 * mu * epsilon(u) + lmbda * tr(epsilon(u)) * Identity(len(u))

# Boundary conditions
top = CompiledSubDomain("near(x[1], 1) && on_boundary")
bot = CompiledSubDomain("near(x[1], 0) && on_boundary")
load = Expression("t", t = 0.0, degree=1)
bcbot= DirichletBC(W, Constant((0.0,0.0)), bot)
bctop = DirichletBC(W.sub(1), load, top)
bc_u = [bcbot, bctop]

boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
top.mark(boundaries,1)
ds = Measure("ds")(subdomain_data=boundaries)
n = FacetNormal(mesh)

# Loading Conditions
load_y = Expression("t", t=0.0, degree=1)

# Boundary conditions and loading applied to boundaries
bcbot = DirichletBC(W, Constant((0.0, 0.0)), bot)
bctop_y = DirichletBC(W.sub(1), load_y, top)
bc_u = [bcbot, bctop_y]

# Variational form
unew, uold = Function(W), Function(W)
E_du = inner(grad(v), sigma(u)) * dx
p_disp = LinearVariationalProblem(lhs(E_du), rhs(E_du), unew, bc_u)
solver_disp = LinearVariationalSolver(p_disp)

# Initialization of the iterative procedure and output requests
t = 0
u_r = 0.016
deltaT = 0.01
tol_phi = 0.001
tol_u = 0.001

# Staggered scheme
while t <= 1.0:
    t += deltaT
    t = round(t, 3)
    load_y.t = t * u_r
    iteration = 0
    err = math.inf
    err_u = math.inf

    while err_u > tol_u:
        iteration += 1
        solver_disp.solve()
        u1 = interpolate(uold, WWW)
        u2 = interpolate(unew, WWW)
        err_u = sqrt(assemble((u1-u2)**2 * dx))
        print('version: ', version, 'Error U: ', round(err_u, 5))
        uold.assign(unew)

        if err_u < tol_u:

            print('version: ', version, ', iterationations:', iteration, ', total time', t, 'time: ',
                  datetime.datetime.now().strftime('%Y-%m-%d %H:%M:%S'))
#             logs.write('iterationations: ' + str(iteration) + ', Total time' +
#                        str(t) + 'time: ' +
#                        datetime.datetime.now().strftime('%Y-%m-%d %H:%M:%S') +
#                        '\n')
            
            # For Stress-Strain Plot
            Traction = dot(sigma(unew), n)
            fy = Traction[1] * ds(1)
            print("vertical forces: " + str(assemble(fy)))
#             fname.write(str(t * u_r) + "\t")
#             fname.write(str(assemble(fy)) + "\n")

            
            if round(t * 1e3) % 100 == 0:
                if str_flag == True:
                    signew = project(sigma(unew), T)

fname.close()
logs.close()
print('Simulation completed')
2

Hi,
plane strain vs plane stress : see 2D linear elasticity — Numerical tours of continuum mechanics using FEniCS master documentation