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

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?

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.

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)
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)

bcbot = DirichletBC(W, Constant((0.0, 0.0)), bot)
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)
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')
``````

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