Bending of the simply supported beam

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,

Hello @dokken. Hope you are doing well. As per the last instruction I have added respective boundary condition for my beam problem. Even after corresponding bcs applied, while visualising it is not showing any displacement on y direction. Could you help me in confirming whether my implementation is correct or not?

Also one more question. Could you please provide some information on applying point load to a non linear problem like the above. As i am unable to locate how to proceed for such a case. Everytime I am getting converged to same topic, How to construct a sparse matrix used for calculating function values at multiple points. Could you please help for my case of problem to apply a point load.

Thanks for your time and considerations. Any advice on the above queries is greatly appreciated.

1 Like

Hello @kumar
Have you finally been able to apply a point load with dolfinx ?
I am stocked to that point too like you were a few month ago…

I’ve shown how to apply a point load in

Hey sorry for the late response. I tried but not succeed. So I applied it as traction over finely discretized element