Hello everyone, I have a rectangular domain where the boundary condition has to be made under the assumption of simply supported beam with bottom left corner point is fixed and the bottom right corner point is constraint only along the y direction.
I have made the boundary condition as per the above condition, but when I visualize the deformation of beam with provided traction, the bending is not at all happening. The displacement magnitude along y direction seems to be zero, My code with imposed boundary condition is,
import os
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import ufl
from dolfinx.fem import Function
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.io import XDMFFile
from dolfinx.nls.petsc import NewtonSolver
from ufl import ds, dx, grad, inner
from dolfinx import fem, default_scalar_type,mesh
#Creating mesh
msh = mesh.create_rectangle (comm=MPI.COMM_WORLD, points=((0.0, 0.0), (6.0,1.0)), n=(96,16), cell_type=mesh.CellType.quadrilateral)
#Create function space
V = ufl.VectorElement("Lagrange", msh.ufl_cell(),1) #For deformation map
P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(),1)
ME = fem.FunctionSpace(msh, ufl.MixedElement([V, P1, P1]))
#Creating Test functions and function
a, q, v= ufl.TestFunctions(ME)
u = Function(ME) # current solution
u0 = Function(ME) # solution from previous converged step
# Split mixed functions
d, p, mu= ufl.split(u)
d0, p0, mu0= ufl.split(u0)
# Zero u
u.x.array[:] = 0.0
# Interpolate initial condition for the concentration
const_val = 0.5
u.sub(1).interpolate(lambda x: np.full((x.shape[1],), const_val))
u.x.scatter_forward()
#Applying boundary condition to beam
def left_corner(x):
return np.logical_and(np.isclose(x[1],0.0),np.isclose(x[0],0.0))
ME0, _ = ME.sub(0).collapse()
bottom_value1 = fem.Function(ME0)
bottom_value1.x.array[:]= 0
left_corner_vertex = mesh.locate_entities_boundary(msh,0,left_corner)
dofs = fem.locate_dofs_topological((ME.sub(0),ME0),0,left_corner_vertex)
bc_left_corner = fem.dirichletbc(bottom_value1,dofs,ME.sub(0))
def right_corner(x):
return np.logical_and(np.isclose(x[1],0.0),np.isclose(x[0],6.0))
ME0y, _ = ME.sub(0).sub(1).collapse()
bottom_value2 = fem.Function(ME0y)
bottom_value2.x.array[:]= 0
right_corner_vertex = mesh.locate_entities_boundary(msh,0,right_corner)
right_dof = fem.locate_dofs_topological((ME.sub(0).sub(1),ME0y),0,right_corner_vertex)
bc_right_corner = fem.dirichletbc(bottom_value2, right_dof, ME.sub(0).sub(1))
def top(x):
return np.isclose(x[1],1)
top_facets = mesh.locate_entities_boundary(msh,msh.topology.dim-1,top)
marked_facets = top_facets
marked_values = np.full_like(top_facets,1)
facet_tag = mesh.meshtags(msh,msh.topology.dim-1,marked_facets, marked_values)
bcs = [bc_left_corner, bc_right_corner]
##### Mechanical solver terms #####
#Defining traction and body forces
Tor=fem.Constant(msh,default_scalar_type((0,-20)))
B=fem.Constant(msh,default_scalar_type((0,0)))
#Defining kinematic quantities
#Spatial dimension
dim=len(d)
#Identity tensor
I=ufl.variable(ufl.Identity(dim))
#Deformation gradient
F=ufl.variable(I+grad(d))
# Right Cauchy-Green tensor
C = ufl.variable(F.T * F)
# Invariants of deformation tensors
Ic = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F))
# Elasticity parameters
E =default_scalar_type(1e4)
nu = default_scalar_type(0.3)
mu_mech = fem.Constant(msh, E / (2 * (1 + nu)))
lmbda = fem.Constant(msh, E * nu / ((1 + nu) * (1 - 2 * nu)))
# Stored strain energy density (compressible neo-Hookean model)
psi = (mu_mech / 2) * (Ic - 3) - mu_mech * ufl.ln(J) + (lmbda / 2) * (ufl.ln(J))**2
# Hyper-elasticity
P = ufl.diff(psi, F)
##### Phase-field solver terms #####
#Defining parameters
p = ufl.variable(p)
epsilon = default_scalar_type(1/64) # interface thickness
areainteng = default_scalar_type(0.3) #Area specific interface energy
dt = 5.0e-06 # time step
k = default_scalar_type(0.05) #Dissipation term
f = areainteng*(((3/2)*epsilon*(grad (p))**2)+(6*(1/epsilon)*p**2*(1-p)**2))
dfdc = ufl.diff(f, p)
# mu_(n+theta)
theta =1
mu_mid = (1.0 - theta) * mu0 + theta * mu
metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=msh,subdomain_data=facet_tag)
dx = ufl.Measure("dx", domain=msh, metadata=metadata)
#Weak form equations
F0 = inner(grad(a), P) * dx - inner(a, B) * dx - inner(a, Tor) * ds(1)
F1 = inner(p, q) * dx - inner(p0, q) * dx + dt * k * inner(grad(mu_mid), grad(q)) * dx
F2 = inner(mu, v) * dx - inner(dfdc, v) * dx - 3*epsilon*areainteng*inner(grad(p), grad(v)) * dx
FF = F0+F1+F2
#Solving Nonlinear problem
problem=NonlinearProblem(FF,u,bcs)
#Solver configurations
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
#solver.convergence_criterion = "residual" #For check
solver.max_it = 10 # Adjust the maximum number of iterations
solver.atol = 1e-8 # Adjust the absolute tolerance
solver.rtol = 1e-6 # Adjust the relative tolerance
# Output file
file = XDMFFile(MPI.COMM_WORLD, "displacement.xdmf", "w")
file.write_mesh(msh)
# We can customize the linear solver used inside the NewtonSolver by
# modifying the PETSc options
ksp = solver.krylov_solver
opts = PETSc.Options() # type: ignore
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
# Step in time
t = 0.0
# Reduce run time if on test (CI) server
if "CI" in os.environ.keys() or "GITHUB_ACTIONS" in os.environ.keys():
T = 1 * dt
else:
T = 50 * dt
V0, dofs1 = ME.sub(1).collapse()
p = u.sub(1)
V1, dofs2 = ME.sub(0).collapse()
disp = u.sub(0)
u0.x.array[:] = u.x.array
while (t < T):
t += dt
r = solver.solve(u)
print(f"Step {int(t/dt)}: num iterations: {r[0]}")
u0.x.array[:] = u.x.array
file.write_function(disp, t)
file.close()
Could anyone confirm me why it is showing zero displacement, even if the traction is given on the top surface. Unfortunately I am unable to locate the beam bending concept regarding my query in previous topics. Thanks for all your time and considerations,