# FEM boundary problem in structural mechanics

I have spent way too much time in the last weeks/months implementing and writing my own finite element solver from scratch. I successfully implemented 2D and 3D elements and wanted to verify a few of my results. As an example, I went ahead and modelled a simple beam with two-dimensional quadratic quad elements.

Its dimensions are $$500 \times 50 \times 1 [mm]$$. I loaded the beam on one end and constrained the other side. To be more exact, I constrained every node at the very left in both (x,y) directions.

Due to $$\nu=0.3$$, the von mises stresses are not meaningful for me. I am particularly interested in the analytical results of $$\sigma_x$$ and the simulation results.

I did the computation by hand myself:

## $$\sigma_{x,max} = \dfrac{M}{W} = \dfrac{F\cdot l}{\frac{I}{e_{max}}}= \dfrac{F\cdot l}{\frac{a^2b}{6}} = 1200 \text{Mpa}$$

Looking at the simulation, it yields the following (displaying $$|\sigma_x|)$$:

If you look closer, you will see the highest values to be found at the boundary:

There seems to be a singularity at the top and bottom node which I constrained. The values next to those nodes have regular stresses around 1200 Mpa which is exactly what I computed.

My question is: How could I reduce this boundary node effect? Maybe my fem-code has some bugs and this is actually not supposed to happen. I am very happy for any kind of feedback.

# Edit 1:

Every node on the left of the beam is constrained the following:

The force on the right end of the beam is applied equally to all nodes (mid and corner-nodes). I know that this is not ideal but since i am mainly focusing on stresses at the boundary, I am fine with that.

• Please add a schematic figure that shows your boundary conditions. Oct 7, 2021 at 20:13
• The stress diagram does not match the type of load (M), how much offset is the axial load (Fx)? What type of support is on the left?
– r13
Oct 7, 2021 at 22:51
• I am loading it with 1000N split over alle the elements on the rhs. I will add an image for the boundary too. Oct 8, 2021 at 4:26
• @r13 Why do you think so? The second image is zoomed into the top left corner. Also note that, for simplicity and color-management I plotted the absolute stress in x direction. So every minus-sign has flipped. Oct 8, 2021 at 4:35
• 1) If the beam is subject to bending where is the compressive stress? The scale has the positive numbers only, which I assume representing the tensile stress. 2) If all supports are pinned, I think you should have received an error message about stability.
– r13
Oct 8, 2021 at 4:51

## 2 Answers

Whatever is causing this is a bug, either in your input, your your FE code, or your graphics display package.

Since we don't have access to any of those, we can't help you find it. The only thing we can do is tell you that it shouldn't happen!

UPDATE: Reading the OP's comments, there seem to be three problems here.

• The constraints are not consistent with the intended stress and displacement field. Specifically, for a 2-dimensional model you only need to fix ONE degree of freedom in the Y direction (e.g. at the middle note of the beam). That will eliminate any issues caused by non-zero Poisson's ratio.
• The loads are also inconsistent. For 4-node elements, the axial load on the two "corner" nodes at the right hand end should be half the load on each of the other nodes. For higher order elements, google for "equivalent nodal loads" to find how to calculate them correctly. In fact this mistake will have minimal effect on the OP's test problem, but that doesn't make it less wrong.
• The algorithm you chose to calculate the element stresses is probably the worst possible one. For displacement-formulation elements it is a terrible idea to calculate surface stresses directly from the element shape functions.
The standard way is to calculate the stresses at the internal points in the element where they are optimally accurate, and then extrapolate to the boundaries. For a 4-node element, calculate the stress at the element centroid and assume it is constant over the whole element.
For higher order elements, the optimal points depend on the integration rule used. The optimal points are sometimes referred to as "Barlow points" after their discover (see https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.1620100202 for the original paper).

To summarize all this: inventing a subject as complex as Finite Element Analysis all by yourself is a fun project, but if you want to produce useful software, you need to learn to avoid the pitfalls that others have discovered during the last 50 years or so.

• So I do know about issue 1), 2) and 3). Extrapolation becomes slightly harder with reduced integration points so I sticked with this for now. It’s inaccurate for sure but does the job for now. I was actually also thinking about only constraining only one DOF in y-direction and it does help :) thanks a lot! Oct 10, 2021 at 8:20

{I am going to assume you are applying a axial loading condition on the right end of the beam).

That seems a singularity and formally, it should happen if you are constraining all the DOF's of all the nodes on the left boundary. But since as it appears that it is only occurring on the top left node and not on bottom left node, I am guessing there exists some bugs within your code. If you want to completely get rid of this problem, then you can use a command within your code where all the nodes on the boundary are able to still deform (to comply with the effect being caused by Poisson's ratio) but the net translation and rotation of this boundary will still be zero, so it won't cause a rigid body motion. What I mean is that you can basically use a deformable fixed support there. This will make the singularities at the top left node (and principally which should also happen at the bottom left node) disappear.

• I am applying a vertical loading condition (in y-direction). I did zoom in the second image. The singularity can be seen at the top and bottom node. What exact reason is there for it to happen? Could this be caused by poissons-ratio? I mean the strain in y direction is very small so I don’t expect any stresses in x direction coming from that Oct 8, 2021 at 4:19
• I had to think about your idea for a while. it does make sense and is probably a lot more realistic than just fixed rigid boundaries in cases where the boundary has equal material properties. Thanks for that idea :) Oct 8, 2021 at 4:39
• @FinnEggers, If you are using the exact and correct gonverning equations and the correct mathematical models and methods to determine the deformations and then strains and then stresses from a set of loading and boundary conditions (as it is done in top FEA softwares like ANSYS or ABAQUS), then it happens because the load ultimately comes to the corner at a single node (rather be at the top corner or bottom corner). This results in a singularity since point load (rather applied or reaction) at a single node result in a singularity. Oct 8, 2021 at 9:07
• Other nodes along that boundary somewhat cannot be thought of as point loads onto the nodes, but the upper and lower corner nodes on the boundary can. Plus, the support that you are providing on the boundary is rigid, i.e. the nodes cannot move along the plane. That also causes a singularity stress problem at the top and bottom corner nodes only (because of poisson's ratio effect). Oct 8, 2021 at 9:08
• That makes sense. I also set poisson = 0 and the entire problem mostly vanished Oct 8, 2021 at 9:10