Different results using different solvers in thermal-elasticity problem

Errors
1

Hi everyone,

I am solving a thermal-elasticity problem for a 5510 box, with a uniform temperature difference of 2000K. My code seems to be working. But when I change the solver (Gmres, Mumps, Petsc, Umfpack), my displacement values seem to be different. Anyone has any idea about this issue?

The result comparison is shown in figure.

Thanks

Here is my code:

from dolfin import *

import os

from subprocess import call

Read_Mesh(parameters_dic)

file_name_Basis_Thermal = ‘box’

convert mesh

if os.path.exists(’{}.xml.gz’.format(file_name_Basis_Thermal)):

os.remove(’{}.xml.gz’.format(file_name_Basis_Thermal))

if os.path.exists(’{}_facet_region.xml.gz’.format(file_name_Basis_Thermal)):

os.remove(’{}_facet_region.xml.gz’.format(file_name_Basis_Thermal))

if os.path.exists(’{}_physical_region.xml.gz’.format(file_name_Basis_Thermal)):

os.remove(’{}_physical_region.xml.gz’.format(file_name_Basis_Thermal))

call(‘dolfin-convert {}.msh {}.xml;gzip {}*.xml’.format(file_name_Basis_Thermal,file_name_Basis_Thermal,file_name_Basis_Thermal),shell=True)

mesh = Mesh(’{}.xml.gz’.format(file_name_Basis_Thermal))

boundaries = MeshFunction(“size_t”, mesh, “{}_facet_region.xml.gz”.format(file_name_Basis_Thermal))

subdomains = MeshFunction(“size_t”, mesh, “{}_physical_region.xml.gz”.format(file_name_Basis_Thermal))

#################### Expansion #####################

file_name_Basis_Elasticity = ‘block’

lmbda = 10.0

mu = 5.0

alpha = 1.25e-5

def eps(v):

return sym(grad(v))

def sigma(v, dT):

return (lmbdatr(eps(v))- alpha(3lmbda+2mu)dT)Identity(3) + 2.0mueps(v)

Vu = VectorFunctionSpace(mesh, ‘CG’, 1)

du = TrialFunction(Vu)

u_ = TestFunction(Vu)

bcu = DirichletBC(Vu.sub(2), Constant(0), boundaries, 13) # bottom

#define measures

ds = Measure(‘ds’, domain=mesh, subdomain_data=boundaries)

dx = Measure(‘dx’, domain=mesh, subdomain_data=subdomains)

Define variational form

Delta_T = 2000

Wint = inner(sigma(du, Delta_T), eps(u_))*dx(14)

aM = lhs(Wint)

LM = rhs(Wint)

Solve problem

u = Function(Vu, name=“Displacement”)

solve(aM == LM, u, bcu, solver_parameters={‘linear_solver’: ‘mumps’})

pvd

pvd_file_u = File(“box.pvd”)

pvd_file_u << u

results
results2445×1089 109 KB

2

Please properly format your code using 3x`.

Also note that that the problem that your solving has a null-space, which you are not properly constraining. (You allow rotations in the x-y plane, making the solution non-unique).
See for instance: Bitbucket on how to set the solver options (and further up in the code how to generate the null-space).

4

Thank you for your suggestion. Actually, all the results shown in the picture have been post-processed. I found the displacement (vector) of point(0,0,0), and let it be the reference point. Maybe there is rotation in x-y plane, but I think the magnitude of the displacement can’t be so different.

5

But you do not know what point the cube is rotating around, this is up to the solver. It looks like you can quite clearly see this for the “PETSc” solver. I still recommend fixing the null space.

2 Likes
6

Hi, sorry for not replying to you for such a long time. After fixing the null space, all the results seems perfect. Thank you so much!