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