# Keeping the surface flat for finite element analysis

I ran the FEM analysis for a simple supported beam problem using solid elements in plane strain condition, with linear kinematics and linear elastic material. I'm not happy with the solution since there is stress concentration at the point boundary condition, see the figure below. I know this is not avoidable for this type of analysis. I think keeping the left and right boundaries straight would help to remedy this problem. However, is it possible to do so, i.e. using special kinematical constraints? I would appreciate if anybody who has experience with this could give some hints.

The boundary condition is:

• ux=uy=0 at x=0,y=0
• uy=0 at x=L, y=0
• pressure=0.5 on y=0.5 (top edge)
• pressure=0.5 on y=-0.5 (bottom edge)

You should never apply loads or boundary conditions on nodes in a FE model. Points have no physical meaning in the real world, at least when speaking of continua.

Normally you would apply a load or boundary condition on a surface. For example you can require that the average vertical and horizontal displacements on the left and the right edge of the domain to be zero, as in the example below.

#!/usr/bin/env python
# -*- coding: UTF8 -*-
###################################################
import getfem as gf
import numpy as np

gf.util_trace_level(1)
gf.util_warning_level(1)

# Input data
NX = 20         # number of elements in horizontal direction
NY = 5          # number of elements in vertical direction

LX = 20.        # length
LY = 5.         # height

nu = 0.3
E = 1000.

#------------------------------------
geotrans = "GT_QK(2,2)"  # geometric transformation
disp_fem_order = 2       # displacements finite element order

#integration_degree = 3   # 4 gauss points per quad
integration_degree = 5   # 9 gauss points per quad
#integration_degree = 7   # 16 gauss points per quad
#------------------------------------

# auxiliary constants
L_RG = 3
B_RG = 4
R_RG = 5
T_RG = 6

m = gf.Mesh("import", "structured",
f"GT='{geotrans}';ORG=[0,0];SIZES=[{LX},{LY}];NSUBDIV=[{NX},{NY}]")
m.set_region(L_RG, m.outer_faces_in_box([-1e-5,-1e-5],[1e-5,LY+1e-5]))      # Left boundary
m.set_region(B_RG, m.outer_faces_in_box([-1e-5,-1e-5],[LX+1e-5,1e-5]))      # Bottom boundary
m.set_region(R_RG, m.outer_faces_in_box([LX-1e-5,-1e-5],[LX+1e-5,LY+1e-5])) # Right boundary
m.set_region(T_RG, m.outer_faces_in_box([-1e-5,LY-1e-5],[LX+1e-5,LY+1e-5])) # Top boundary

# FEM
mfu = gf.MeshFem(m, 2)
mfu.set_fem(gf.Fem("FEM_Q2_INCOMPLETE(2)"))

mfout = gf.MeshFem(m)
mfout.set_classical_discontinuous_fem(disp_fem_order)

# integration method
mim = gf.MeshIm(m, integration_degree)

# Model
md = gf.Model("real")

# Hyperelasticity

md.add_variable("qL", 2)              # average reaction traction on left boundary
md.add_variable("qR", 2)              # average reaction traction on right boundary

# enforcing zero avg. displacements on left boundary
# enforcing zero avg. displacements on right boundary

print("Displacement dofs: %i\nTotal model dofs: %i" % (mfu.nbdof(),md.nbdof()))

for step,fact in enumerate([0.,0.2,0.4,0.6,0.8,1.]):
md.set_variable("pB", fact*0.05*E)
md.set_variable("pT", fact*0.05*E)
nit, converged = md.solve("noisy", "lsolver", "mumps", "max_iter", 20, "max_res", 1e-8, #)[0]
"lsearch", "simplest", "alpha max ratio", 1.5, "alpha min", 0.2, "alpha mult", 0.6)
VM = md.interpolation("Norm(Deviator(tau3d(u)))", mfout)
output = (mfu, md.variable("u"), "Displacements",
mfout, VM, "Von Mises stress")
for i,j in ((1,1),(2,2),(3,3),(1,2)):
sigma = md.interpolation(f"sigma3d(u)({i},{j})", mfout)
output += (mfout, sigma, f"sigma_{i}{j}")
mfout.export_to_vtu(f"simple_supported_beam_{step}.vtu", *output)


If you really want to keep the edges flat, it is also possible but you need to impose an additional constraint to each edge, you can see how it can be done with Lagrange multipliers by comparing the previous code with the following one:

#!/usr/bin/env python
# -*- coding: UTF8 -*-
############################################################################
import getfem as gf
import numpy as np

gf.util_trace_level(1)
gf.util_warning_level(1)

# Input data
NX = 20         # number of elements in horizontal direction
NY = 5          # number of elements in vertical direction

LX = 20.        # length
LY = 5.         # height

nu = 0.3
E = 1000.

#------------------------------------
geotrans = "GT_QK(2,2)"  # geometric transformation
disp_fem_order = 2       # displacements finite element order

#integration_degree = 3   # 4 gauss points per quad
integration_degree = 5   # 9 gauss points per quad
#integration_degree = 7   # 16 gauss points per quad
#------------------------------------

# auxiliary constants
L_RG = 3
B_RG = 4
R_RG = 5
T_RG = 6

m = gf.Mesh("import", "structured",
f"GT='{geotrans}';ORG=[0,0];SIZES=[{LX},{LY}];NSUBDIV=[{NX},{NY}]")
m.set_region(L_RG, m.outer_faces_in_box([-1e-5,-1e-5],[1e-5,LY+1e-5]))      # Left boundary
m.set_region(B_RG, m.outer_faces_in_box([-1e-5,-1e-5],[LX+1e-5,1e-5]))      # Bottom boundary
m.set_region(R_RG, m.outer_faces_in_box([LX-1e-5,-1e-5],[LX+1e-5,LY+1e-5])) # Right boundary
m.set_region(T_RG, m.outer_faces_in_box([-1e-5,LY-1e-5],[LX+1e-5,LY+1e-5])) # Top boundary

# FEM
mfu = gf.MeshFem(m, 2)
mfu.set_fem(gf.Fem("FEM_Q2_INCOMPLETE(2)"))

mfmult = gf.MeshFem(m, 1)
mfmult.set_classical_fem(disp_fem_order)

mfout = gf.MeshFem(m)
mfout.set_classical_discontinuous_fem(disp_fem_order)

# integration method
mim = gf.MeshIm(m, integration_degree)

# Model
md = gf.Model("real")

# Hyperelasticity

md.add_variable("qL", 1)              # average reaction traction on left boundary
md.add_variable("qR", 1)              # average reaction traction on right boundary

# enforcing zero avg. vertical displacement on left boundary
# enforcing zero avg. vertical displacement on right boundary

# enforcing transverse free rotation to the left and right boundaries

print("Displacement dofs: %i\nTotal model dofs: %i" % (mfu.nbdof(),md.nbdof()))

for step,fact in enumerate([0.,0.2,0.4,0.6,0.8,1.]):
md.set_variable("pB", fact*0.05*E)
md.set_variable("pT", fact*0.05*E)
nit, converged = md.solve("noisy", "lsolver", "mumps", "max_iter", 20, "max_res", 1e-8, #)[0]
"lsearch", "simplest", "alpha max ratio", 1.5, "alpha min", 0.2, "alpha mult", 0.6)
VM = md.interpolation("Norm(Deviator(tau3d(u)))", mfout)
output = (mfu, md.variable("u"), "Displacements",
mfout, VM, "Von Mises stress")
for i,j in ((1,1),(2,2),(3,3),(1,2)):
sigma = md.interpolation(f"sigma3d(u)({i},{j})", mfout)
output += (mfout, sigma, f"sigma_{i}{j}")
mfout.export_to_vtu(f"simple_supported_beam_{step}.vtu", *output)


In the second script, the conditions

$$\text{avg}(u_y)=0$$ at $$x=0$$, and

$$\text{avg}(u_y)=0$$ at $$x=L$$

are imposed with two scalar multipliers $$q_L$$ and $$q_R$$ and the Lagrangians

$$\int_{\Gamma_L} u_y q_L~d\Gamma$$

$$\int_{\Gamma_R} u_y q_R~d\Gamma$$

The straightness of the deformed left edge and right edge, but still allowing the beam to contract or expand in the transverse direction, are implemented with the help of two scalar variables $$s_L$$ and $$s_R$$ (tanthetaL and tanthetaR in the code), and two scalar Lagrange multiplier fields, $$\Lambda_L$$ and $$\Lambda_R$$ (multL and multR in the code). The respective Lagrangians are

$$\int_{\Gamma_L}(u_x-s_L(y+u_y-H/2))\Lambda_L~d\Gamma$$

$$\int_{\Gamma_R}(u_x-s_R(y+u_y-H/2))\Lambda_R~d\Gamma$$

All necessary weak forms are derived by differentiating these Lagrangians with respect to $$u$$, $$s_L$$, $$s_R$$, $$\Lambda_L$$, $$\Lambda_R$$.

The origin of e.g. the first constraint is the line equation

$$\begin{bmatrix}x+u_x-0\\y+u_y-H/2\end{bmatrix}\cdot\begin{bmatrix}1\\s_L\end{bmatrix}=0$$

with $$x=0$$.

This is a general solution for nonlinear kinematics. For linearized kinematics the formulation can obviously be simplified a lot.

• It's actually very cool. Thanks a lot for the example. I have one question though: to enforce the average displacement you are using Lagrange multiplier?
– kstn
Commented Feb 5 at 22:35
• yes, 4 scalar Lagrange multipliers in total, 2 per side, one for each direction. The Lagrangian "u.qL" is equivalent to "u(1)*qL(1)+u(2)*qL(2)". Commented Feb 5 at 22:51
• Do you know if we can use Penalty method to enforce that. The code i'm using is nodal-based and element-based. It's pretty hard to introduce global Lagrange multiplier.
– kstn
Commented Feb 5 at 23:02
• Not with penalization, but most solvers have their way of implementing constraints. If your solver allows it, you could add a constraint that the sum of all nodal displacements on one side, in every direction, is zero. It will be approximately the same as implementing the correct integral. Commented Feb 6 at 9:50
• Added a few clarifying equations. Commented Feb 8 at 15:03