Hello, I am trying to solve a dynamic linear elastic problem. I have gone through numerous tutorials. But the solution is simply diverging. Any help would be appreciated.
import dolfinx as fe
from dolfinx.fem.petsc import assemble_vector, assemble_matrix, create_vector, apply_lifting, set_bc
import numpy as np
from mpi4py import MPI
import ufl
from petsc4py import PETSc
def main():
L = 1
W = 0.2
mesh = fe.mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, W, W])],
[20, 6, 6], cell_type=fe.mesh.CellType.hexahedron)
element = ufl.VectorElement("Lagrange", mesh.ufl_cell(), 1)
V = fe.fem.functionspace(mesh, element)
def boundary_1(x):
return np.isclose(x[0], 0)
fdim = mesh.topology.dim - 1
boundary_1_facets = fe.mesh.locate_entities_boundary(mesh, fdim, boundary_1)
bc = fe.fem.dirichletbc(fe.default_scalar_type((0.0, 0.0, 0.0)), fe.fem.locate_dofs_topological(V, fdim, boundary_1_facets), V)
u = ufl.TrialFunction(V) # Incremental displacement
v = ufl.TestFunction(V) # Test function
u_dot = fe.fem.Function(V) # First time derivative of displacement
u_2dot = fe.fem.Function(V) # Second time derivative of displacement
u_old = fe.fem.Function(V)
u_dot_old = fe.fem.Function(V)
u_2dot_old = fe.fem.Function(V)
rho = 1
d = len(u)
I = ufl.variable(ufl.Identity(d)) # Identity tensor
F = ufl.variable(I + ufl.grad(u)) # Deformation gradient
E = ufl.variable(F.T*F) # Right Cauchy-Green tensor
J = ufl.variable(ufl.det(F))
metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=mesh, metadata=metadata)
dx = ufl.Measure("dx", domain=mesh, metadata=metadata)
maximum_pressure = 15
minimum_pressure = 5
time_step = fe.default_scalar_type(0.0005)
num_steps = int(1/0.0005)
pressure = np.append(np.linspace(minimum_pressure, maximum_pressure, num=10), np.linspace(maximum_pressure, minimum_pressure, num=num_steps))
p = fe.fem.Constant(mesh, pressure[0])
normal = ufl.FacetNormal(mesh)
E = fe.default_scalar_type(1.0e4)
nu = fe.default_scalar_type(0.3)
mu = fe.fem.Constant(mesh, E / (2 * (1 + nu)))
lmbda = fe.fem.Constant(mesh, E * nu / ((1 + nu) * (1 - 2 * nu)))
S = 2.0 * mu * ufl.sym(ufl.grad(u)) + lmbda * ufl.tr(ufl.sym(ufl.grad(u))) * I
F = rho * ufl.inner(u_2dot, v) * dx + ufl.inner(ufl.grad(v), S) * dx - ufl.inner(v, p * normal) * ds
bilinear_form = fe.fem.form(ufl.lhs(F))
linear_form = fe.fem.form(ufl.rhs(F))
A = assemble_matrix(bilinear_form, bcs=[bc])
A.assemble()
b = create_vector(linear_form)
solver = PETSc.KSP().create(mesh.comm)
solver.setInitialGuessNonzero(True)
solver.setOperators(A)
solver.getPC().setType(PETSc.PC.Type.SOR)
solver.view()
displacement = fe.fem.Function(V)
alpha = fe.default_scalar_type(0.5)
betta = fe.default_scalar_type(0.5)
def update_u_2dot(u_prime, u_old_prime, u_dot_old_prime, u_2dot_old_prime):
return 2/(betta * time_step**2) * (u_prime - u_old_prime)\
- 2/(betta * time_step) * u_dot_old_prime\
- (1/betta - 1) * u_2dot_old_prime
def update_u_dot(u_dot_old_prime, u_2dot_old_prime, u_2dot_prime):
return u_dot_old_prime\
+ time_step * (1-alpha) * u_2dot_old_prime\
+ time_step * alpha * u_2dot_prime
t = 0
for i in range(num_steps):
print(i)
t += time_step
p.value = pressure[i]
# Update the right hand side reusing the initial vector
with b.localForm() as loc_b:
loc_b.set(0)
assemble_vector(b, linear_form)
# Apply Dirichlet boundary condition to the vector
apply_lifting(b, [bilinear_form], [[bc]])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, [bc])
# Solve linear problem
solver.solve(b, displacement.vector)
displacement.x.scatter_forward()
print(displacement.x.array)
u_2dot_prime = update_u_2dot(displacement.x.array[:], u_old.x.array[:], u_dot_old.x.array[:], u_2dot_old.x.array[:])
u_dot_prime = update_u_dot(u_dot_old.x.array[:], u_2dot_old.x.array[:], u_2dot_prime)
u_2dot.x.array[:] = u_2dot_prime
u_dot.x.array[:] = u_dot_prime
if __name__ == "__main__":
main()