# How to Apply a Uniform Load to a Laminate Plate in FEM

As the title says, I'm trying to make my code include a uniform pressure/force to the surface of a plate for analysis.

Structural/Element Modeling

I am modeling a 2D 9 node Mindlin plate that has a clamped boundary condition on one side. It is 0.2 [m] by 0.2 [m]. Currently, I am using 6x6 elements. The plate has 6 layers, consisting of, in order from top to bottom,

• 0.1 [mm] thick piezoelectric layer at a 0 degree angle
• 0.25 [mm] thick composite material layer at a -45 degree angle
• 0.25 [mm] thick composite material layer at a 45 degree angle
• 0.25 [mm] thick composite material layer at a 45 degree angle
• 0.25 [mm] thick composite material layer at a -45 degree angle
• 0.1 [mm] thick piezoelectric layer at a 0 degree angle

The Problem

As far as calculating the electric load, I am 100% sure that it is correct. Same with calculating the various effects of the different orientations of the plates. The code works perfectly fine when I only include the electrical load, so that's not the problem.

After some tinkering and some astute remarks, I have realized that the issue is not with the construction of the Force Matrix, but rather, the way I implemented the uniform load. I had originally just made a Force Matrix and multiplied the entire thing with the load I wanted to apply. This was a fundamental flaw. So I started my studies again and found that the force matrix must be calculated as shown in the image.

Since in my case, the fz was a constant, I integrated the nine-node shape functions w.r.t to the area and got the following:

H1 =  (1/4) * (1 - xi) * (1 - eta) * xi * eta --> 1/9
H2 = -(1/4) * (1 + xi) * (1 - eta) * xi * eta --> 1/9
H3 =  (1/4) * (1 + xi) * (1 + eta) * xi * eta --> 1/9
H4 = -(1/4) * (1 - xi) * (1 + eta) * xi * eta --> 1/9
H5 = -(1/2) * (1 - xi^2) * (1 - eta) * eta    --> 4/9
H6 =  (1/2) * (1 + xi) * (1 - eta^2) * xi     --> 4/9
H7 =  (1/2) * (1 - xi^2) * (1 + eta) * eta    --> 4/9
H8 = -(1/2) * (1 - xi) * (1 - eta^2) * xi     --> 4/9
H9 =  (1 - xi^2) * (1 - eta^2)                --> 16/9


The Code I Implemented

So, instead of making my Force Vector with the code:

fef = [ H1 H2 H3 H4 H5 H6 H7 H8 H9 ] * P;

for i = 1:length(nodes)

Fe(i,1) = 0;
Fe(i+num_node_ele,1) = 0;
Fe(i+num_node_ele*2,1) = 0;
Fe(i+num_node_ele*3,1) = 0;
Fe(i+(num_node_ele*4),1) = fef(i);

end

f = f + Fe * (detjacob);


I changed it by making P into a matrix like so:

a = 1/9;
b = 4/9;
c = 16/9;

dist_load = P * [ a a a a b b b b c ];

% the rest is identical


Doing so gives me deflections of reasonable ranges and the deflection converges as I increase the number of elements. So now, I must ask, did I do this correctly? Is there something I missed?

• If the displacement vector is all zero, my guess is that this is just a basic programming error. What do you get when you print Fg for a one-element model? The nd vector of indices looks odd to me. – Bill Greene Feb 27 '17 at 17:57
• Hi! You always seem to help me out with my FEM debugging and I truly appreciate it! Thank you. So first, your advice was well taken, and the loading of the initial vector was a bit wonky. So I fixed that and added a rather redundant and length section at the end for global assembly. Also, instead of doing 2 nested for loops for the Gauss Points, I just did a single for loop of 4, thus simplifying things. However, now the problem is the opposite: the smallest force deforms the plate too much! And to top it off, increasing the number of elements increases the load! (1/2) – Additional Pylons Mar 1 '17 at 8:22
• I have checked the Jacobians and they are diagonal and everything else seems to be in place. I just don't know what the issue is. Also, for the sake of clarity, the matrix is ordered so that it first lists all the x-directional force for the all nodes, then the y-directional, then the x-directional angular, then the y-directional angular, then finally all the z-directional load. Thus, the Fg matrix is all 0 until the end. Any advice is appreciated! (2/2) – Additional Pylons Mar 1 '17 at 8:28
• Sorry, You haven't replied yet, but I just had a quick question: Is there some factor that I should multiply the force P with? It is a uniform load and in one of my previous questions, that was a problem I had with a 2D plate and a uniform force was applied along the length of the plate and I had to multiply it with Lx/(2*num_ele_x). – Additional Pylons Mar 1 '17 at 9:18
• The sum of the nine nodal forces obviously should equal P * Area. And the nodal forces should be symmetric wrt to axes through the element center. When debugging, I always add an area( or volume) calculation to my integration loop as a simple test. Those tests should be enough to find the problem. There is no need to consider more than a single element while debugging. – Bill Greene Mar 1 '17 at 16:58

For a plate element lying in the x-y plane, the consistent nodal loads corresponding to a pressure load on the face are given by

$$f_i = \int_{A_e} N_i P(x,y) dx dy$$

where $N_i$ is the shape function at the ith node and $P$ is the pressure. Of course, if $P$ is constant, as it is in your case, it can be moved outside the integral so the nodal loads are simply the integral of the shape functions multiplied by the pressure

$$f_i = P\int_{A_e} N_i dx dy$$

The only remaining step is to sum these elemental load components into the appropriate locations in the global load vector.

That's all there is to it.

• Hello! Sorry, but I don't think I was clear enough. I apologize. In my new implementation of code, the load_dist 'scaling' is simply stating the integrals of the shape functions in a matrix and multiplying it all with the applied load P. – Additional Pylons Mar 3 '17 at 5:13