# 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. 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? 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? 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? 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. Nov 27 '16 at 20:30

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$