0
$\begingroup$

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.

Force Vector equation

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 ];

fef = H .* dist_load;
% 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?

Thank you for your help!

|cite|improve this question
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Bill Greene Feb 27 '17 at 17:57
  • $\begingroup$ 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) $\endgroup$ – Additional Pylons Mar 1 '17 at 8:22
  • $\begingroup$ 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) $\endgroup$ – Additional Pylons Mar 1 '17 at 8:28
  • $\begingroup$ 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). $\endgroup$ – Additional Pylons Mar 1 '17 at 9:18
  • $\begingroup$ 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. $\endgroup$ – Bill Greene Mar 1 '17 at 16:58
0
$\begingroup$

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.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Additional Pylons Mar 3 '17 at 5:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.