# 3D Solid 8 Node FEM Matlab Code

So this semester, I'm taking a Finite Element Method course at my graduate school. We started out making codes for 1D bars and came all the way to 8 node solid elements. However, I seem to have run into a wall, as I have made my code and combed through it for the past week, making sure all the shape functions and mathematics were correct. In the end, everything looks correct. However, while I am supposed to get deflections in the z direction around the ballpark of 0.3 m, I keep getting deflections in the magnitude of E-04. I am very much at my wit's end. Also, because I genuinely don't know which StackExchange community to post this question to, I have tried Mathematics and Overflow. I would like to emphasize that I'm not trying to just spam my question. If you kindly explain to me why my question doesn't belong here, I will delete it quickly. But I have done my research and I have made my own code, so I'm not trying to something out of nothing. Any and all help is appreciated. I will detail the problem below:

• There is a bar that is fully clamped at one end.
• It is 1[m] long in the x-dir and 0.1[m] long in the y-dir and z-dir.
• A point force of 100[kN] is applied in an upward direction at the free end of the bar. The number of elements in the y-dir and z-dir is to always be 1 element only, but the I'm supposed to increase the number of x-dir elements to show that the displacement values converge. Currently, not only do they not converge, but the values are much too small.

Once again, thank you for any and all help.

• There are some low-level tests you should do when writing code for finite element matrices. First, you should print det(Jacob) at every integration point; you will immediately see a problem in your code. Then print Jacob at every integration point. Does what you see make sense? When you sort out those problems, simply calculate the volume of the element by summing over your eight integration points. If you don't get the correct volume (,01 in your case), obviously you won't get the correct stiffness matrix. – Bill Greene Nov 24 '16 at 16:12
• Hi! Thanks for responding. Sorry if I'm being incredibly thick, but I can't really see a problem with my det(Jacob) and my Jacobian Matrix in general. when I display it, it rounds it up from 0.00125 to 0.0013, but changing the format to long doesn't change anything. Also, I believe that the Jacobian is meant to be diagonal, especially since I programmed it to follow the math shown here. Could you perhaps be a little more specific on what is wrong with my Jacobian? – Additional Pylons Nov 25 '16 at 4:17
• Sorry, my mistake! When I rerun your code more carefully, the jacobian and it's determinant look OK to me. The calculated volume is also correct. As a guess, the results from your 1-element case also look OK to me. Haven't tried more than one element. Why do you expect a tip deflection of around .3m? – Bill Greene Nov 25 '16 at 13:36
• I can't comment on the correctness of your model, but I do have a question about how you apply your boundary conditions. I always zero-out the rows of the stiffness matrix and set the boundary node value to 1, as you do, but I leave for the columns alone. Why do you set these to zero as well? – cbcoutinho Nov 25 '16 at 20:11
• Using the Euler-Bernoulli equations I obtained a deflection in the tip of $1/30\text{ m}$. Using, elasticity for your problem you would end up with a deflection that depends on the third power of your coordinates... probably a single element is not enough to capture this behavior. – nicoguaro Nov 27 '16 at 20:30

After some investigation, I believe your FE code is correct. Congratulations!

From beam theory, deflection of a cantilever beam with a point load at the tip is $$u_{tip} = PL^3/(3EI)$$ For the parameters of your model, $u_{tip}=.03333m$.

With a single element, your code gives $u_{tip}=8.4697 \times 10^{-4}m$. For 30 elements, $u_{tip}=0.027862m$

So why isn't the FE solution closer to the beam solution?

These fully-integrated eight-node hex elements (as well as the corresponding 2D, four-node quadrilateral) are notoriously bad when modeling bending behavior. Remedies for this problem have been studied for more than 30 years. One crude solution is to simply under-integrate the terms when forming the stiffness matrix (e.g. a single gauss point). But this causes new problems. A better approach is to use a more sophisticated approach than the simple displacement formulation. You can find a wealth of information if you search for "shear locking in finite elements."

• Thank you so much for dedicating so much time to my problem! I know how tedious debugging code can be, let alone for someone else's code. The reason why I was expecting a deflection value of around 0.3[m] was because the professor gave us that value to help us debug. I realize now that he made a mistake and told us 0.3 instead of 0.03, but I never questioned it because my values didn't converge and I thought I was wrong. I will definitely look more into the locking effect. Your analysis has been of tremendous help. So thank you once more for you help! – Additional Pylons Nov 26 '16 at 3:03
• Second order elements would yield better results. They cause a lot less trouble when it comes to locking. – P. G. Nov 27 '16 at 9:45