Hello everyone, I am currently solving the 3D unconstrained thermal expansion problem at a uniform temperature, using the null-space elimination model for translation/rotation.
For a cylinder with an axial mesh of 10 layers, when the length is 10mm, I plotted the displacement results along a diameter on the top/bottom surface (as shown in the figure). Ideally, the displacement in the x-direction should be 0, and the result at this point is on the order of 1e-7, indicating that the null-space is useful. However, when I change the length to 1000mm, this value increases to 1e-3 (indicating that the effect of null-space is minimal), which seems to vary proportionally with the size. Does anyone know what might be causing this, and how I can resolve this issue?
Here is my geo file(The length is 10mm):
r1 = 7;
r2 = 5;
N1 = 16;
N2 = 4;
Point(1) = {0, 0, 0, 1.0};
Point(2) = {-r1, 0, 0, 1.0};
Point(3) = {0, r1, 0, 1.0};
Point(4) = {0, -r1, 0, 1.0};
Point(5) = {r1, 0, 0, 1.0};
Point(6) = {0, -r2, 0, 1.0};
Point(7) = {0, r2, 0, 1.0};
Point(8) = {r2, 0, 0, 1.0};
Point(9) = {-r2, 0, 0, 1.0};
Circle(1) = {3, 1, 5};
Circle(2) = {5, 1, 4};
Circle(3) = {4, 1, 2};
Circle(4) = {2, 1, 3};
Line(5) = {7, 8};
Line(6) = {8, 6};
Line(7) = {6, 9};
Line(8) = {9, 7};
Line(9) = {3, 7};
Line(10) = {5, 8};
Line(11) = {4, 6};
Line(12) = {9, 2};
Transfinite Curve {1:8} = N1 Using Progression 1;
Transfinite Curve {9:12} = N2 Using Progression 1;
Curve Loop(1) = {8, -9, -4, -12};
Plane Surface(1) = {1};
Curve Loop(2) = {1, 10, -5, -9};
Plane Surface(2) = {2};
Curve Loop(3) = {2, 11, -6, -10};
Plane Surface(3) = {3};
Curve Loop(4) = {3, -12, -7, -11};
Plane Surface(4) = {4};
Curve Loop(5) = {8, 5, 6, 7};
Plane Surface(5) = {5};
Transfinite Surface{1:5};
Recombine Surface{1:5};
Extrude {0, 0, 10} {Surface{1:5};Layers {10};Recombine;}
Physical Surface("surface", 123) = {29, 87, 65, 43};
Physical Surface("top", 124) = {34, 100, 78, 56, 122};
Physical Surface("bottom", 125) = {1:5};
Physical Volume("v", 126) = {1:5};
Here is my code:
import dolfinx
from dolfinx import la
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx.io import gmshio
from dolfinx import io
from ufl import (tr, inner, lhs, rhs, Identity, TestFunction, TrialFunction, dx, grad)
from dolfinx.fem import (functionspace, FunctionSpace, Function, form)
from dolfinx.fem.petsc import (assemble_matrix, assemble_vector)
def null_space(V: FunctionSpace):
_x = Function(V)
gdim = V.mesh.geometry.dim
assert gdim == 2 or gdim == 3
dim = 3 if gdim == 2 else 6
nullspace_basis = [
la.vector(V.dofmap.index_map, bs=V.dofmap.index_map_bs, dtype=PETSc.ScalarType)
for i in range(dim)
]
basis = [b.array for b in nullspace_basis]
dofs = [V.sub(i).dofmap.list.reshape(-1) for i in range(gdim)]
for i in range(gdim):
basis[i][dofs[i]] = 1.0
x = V.tabulate_dof_coordinates()
dofs_block = V.dofmap.list.reshape(-1)
x0, x1, x2 = x[dofs_block, 0], x[dofs_block, 1], x[dofs_block, 2]
if gdim == 2:
basis[2][dofs[0]] = -x1
basis[2][dofs[1]] = x0
elif gdim == 3:
basis[3][dofs[0]] = -x1
basis[3][dofs[1]] = x0
basis[4][dofs[0]] = x2
basis[4][dofs[2]] = -x0
basis[5][dofs[2]] = x1
basis[5][dofs[1]] = -x2
for b in nullspace_basis:
b.scatter_forward()
_basis = [x._cpp_object for x in nullspace_basis]
dolfinx.cpp.la.orthonormalize(_basis)
assert dolfinx.cpp.la.is_orthonormal(_basis)
local_size = V.dofmap.index_map.size_local * V.dofmap.index_map_bs
basis_petsc = [
PETSc.Vec().createWithArray(x[:local_size], bsize=gdim, comm=V.mesh.comm)
for x in basis
]
return PETSc.NullSpace().create(comm=V.mesh.comm, vectors=basis_petsc)
# mesh
mesh, cell_markers, facet_markers = gmshio.read_from_msh("test.msh", MPI.COMM_WORLD, gdim=3)
# parameters
E = 2e5
nu = 0.3
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
alpha = 1e-5
def eps(v):
return 0.5 * (grad(v) + grad(v).T)
def sigma(v, dT):
return (lmbda*tr(eps(v))-alpha*(3*lmbda + 2*mu)*dT)*Identity(3)+2.0*mu*eps(v)
# Variational formula
VU = functionspace(mesh, ("CG", 1, (mesh.geometry.dim, )))
dU = TrialFunction(VU)
U_ = TestFunction(VU)
Wint = inner(sigma(dU, 1000), eps(U_)) * dx
aM = lhs(Wint)
LM = rhs(Wint)
A = assemble_matrix(form(aM))
A.assemble()
b = assemble_vector(form(LM))
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
# Create Krylov solver
solver = PETSc.KSP().create(A.getComm())
solver.setOperators(A)
null_space = null_space(VU)
A.setNullSpace(null_space)
null_space.remove(b)
U = Function(VU)
solver.solve(b, U.vector)
# save results
vtkFile = io.VTKFile(MPI.COMM_WORLD, "results/U.pvd", "w")
vtkFile.write_function(U)
vtkFile.close()
and here is the visualization of the deviation results:
