Results obtained by FEniCS mismatches Abaqus

I/O
1

I performed the following code:

domain = Rectangle(Point(0, 0), Point(3, 3))
mesh = generate_mesh(domain, number_of_element)
Vu = VectorFunctionSpace(mesh, "CG", 1)
u = Function(Vu, name="Displacement")
# Solve displacement u-problem
        problem_u.solve(u.vector())

###define functionspace
V=FunctionSpace(mesh,'P',1)
W=TensorFunctionSpace(mesh,"DG",0)

###get displacement###
u_magnitude = sqrt(dot(u,u))
u_magnitude = project(u_magnitude,V)
u_x, u_y = u.split(deepcopy=True)     
u_x=project(u_x,V)

###get stress component
epsilo=0.5*(nabla_grad(u)+nabla_grad(u).T)    #strain
sigma=E*nu/((1+nu)*(1-2*nu))*tr(epsilo)*Identity(2)+2*(E/(2*(1+nu)))*epsilo
sigma=project(sigma,W)
s11,s12,s21,s22=sigma.split(deepcopy=True)

###get von-mises
s = sigma_ - (1. / 3) * tr(sigma_) * Identity(2)
von_Mises = sqrt(3. / 2 * inner(s, s))
von_Mises = project(von_Mises, V)

the von-mises obtained by FEniCS have 1.3% error Abaqus, but strange, the displacement u_magnitude and stress component s11 both have large error (about 30%) compared with Abaqus. Why this could happen? (I guess the problem could be split command? )

this is a plane strain problem and choose linear triangle element in Abaqus.
the similar number of element is applied in Abaqus and FEniCS

Thanks in advance

2

The code you supplied is not complete, and therefore there is not much one can do to investigate this.
Note that the von Mises stresses should be projected into W (the discontinuous space).

3

Hi dokken, thanks for your note. von_mises should be projected into DG. but interestingly, when i project it into V=FunctionSpace(mesh,'P',1) ,i can still get the same answer with W=FunctionSpace(mesh,"DG",0) . why this can happen?

for whole code, you can download at here. In this code, I want to extract the stress and displacement. then I insert the following code after original code-print("force=",results[i + 1, 0])

###get displacement###
u_magnitude = sqrt(dot(u,u))
u_magnitude = project(u_magnitude,V)
u_x, u_y = u.split(deepcopy=True)     
u_x=project(u_x,V)

###get stress component
epsilo=0.5*(nabla_grad(u)+nabla_grad(u).T)    #strain
sigma=E*nu/((1+nu)*(1-2*nu))*tr(epsilo)*Identity(2)+2*(E/(2*(1+nu)))*epsilo
sigma=project(sigma,W)
s11,s12,s21,s22=sigma.split(deepcopy=True)

###get von-mises
s = sigma_ - (1. / 3) * tr(sigma_) * Identity(2)
von_Mises = sqrt(3. / 2 * inner(s, s))
von_Mises = project(von_Mises, V)

but i can not get the correct u_x results compared with Abaqus. then i want to check the code,for simplicity, i eliminate phase field by inserting the following code:

Vd = FunctionSpace(mesh, "CG", 1)
d = Function(Vd, name="Damage")
dold = Function(Vd, name="Previous damage")
bcd1 = DirichletBC(Vd, Constant(0.),top)
bcd2 = DirichletBC(Vd, Constant(0.), bottom)
bcd=[bcd1,bcd2]
d_=Expression('0.0',degree=1)
d=interpolate(d_,Vd)

then the whole domain will have no damage. problem is elastic. but i can not still get the correct u_x displacement.

4

How do you compare your solution with abaqus? Do you use paraview, or integrate some reference quantities. If you are using paraview, you Need to use XDMFFile.write_checkpoint to save the solution as something other than a CG-1 interpolation. See: Loading xdmf data back in .
I guess that @bleyerj might be able to give you some more pointers given that he is on of the authors of the code.

5

thanks dokken, i will read it late.

6

HI,
sorry but it’s hard to help you like that, the problem could well be from your Abaqus input or outputs. I suggest you check with an independent FEniCS standard elasticity implementation to see if the problem comes from the use of MFront.

1 Like
7

hello again jeremy, I have double check the results and code have no problem. the reason for the error happening is I wrongly write a parameter.
Additionally, regarding to “phase field approach to brittle fracture” code, i have changed it to suit fgm-analysis. but i don’t know how to add a newton-load on boundary like the following picture show:
image

many thanks

8

Hello, jeremy, sorry for asking for help too fast, i have got the answer from your demo “Small-strain von Mises elastoplasticity” .thanks a lot