MITC4 Plate Bending Elements 3

[DCBII]

I'm a full time engineer, but I've put together a sort of hobby computer program that does basic finite element analysis to help me learn how it works. I've been trying to add MITC4 elements for linear analysis (nothing too fancy) but I'm having some trouble. This is my first time creating isoparametric elements, so I'm learning.

I've modeled a 10' wide x 20' tall x 1' thick concrete rectangular wall panel (E = 3823 ksi, nu=0.17) fixed on on four sides with a uniform surface pressure of 250 psf. I've meshed it incredibly fine with 6"x6" rectangular plates.

The shear stresses my model is calculating are spot on with values published by Timoshenko in "Theory of Plates and Shells", but the bending stresses that are calculated are about 15% smaller. The displacement calculated at the center of the panel is also slightly larger (~10%) than Timoshenko predicts. I've plotted contours and the displaced shape and the shape of the bending moment, and they both look correct. Only the magnitudes appear to be slightly off. The excess deflection and low moments make me think my flexural stiffness is too low, but my stress-strain matrix [Cb] is the easiest matrix to calculate and I've checked it multiple times.

I've calculated my stresses as [Mx, My, Mxy] = [Cb][d], with calculated at the 2x2 gauss integration points (+/-0.5773). I then extrapolated those stresses to the nodes using the shape functions. The extrapolation made very little difference with such a fine mesh, as the stress is nearly constant over the plate's area at locations of maximum moment. I don't think my extrapolation is the problem.

Any finite element geniuses out there with an idea what might be going on here? Do MITC4 bending elements require me to mesh even finer? My slow personal computer might blow up if I do. That seems like extremely poor convergence.

 

[rb1957]

when you say "6"x6" rectangular plates" … how many elements through the thickness (1') ? minimum 3 for a bending field, the more the better.

btw, your mesh (800 elements) is not, IMHO, "incredibly fine".

you should do a single element check and compare your results with mainstream FEA codes. Then model structures … model simple ones first, ones you can hand calc the exact solution, like round plates.

another day in paradise, or is paradise one day closer ?

 

[DCBII]

These are MITC4 elements, so there's only 1 element through the thickness. They are 4 node quadrilateral bending elements. The reason I say 6"x6" is incredibly fine is because I've run this same problem through other software (STAAD) and gotten the results I expected.

 

[rb1957]

still do a single element test. apply a moment causing bending across the thickness; how does the element deflect ?

another day in paradise, or is paradise one day closer ?

 

[nlgyro]

10'x 20'x 1' plate with pressure load on a fine grid model... sounds like a geometrically non-linear problem to me. how about doing a large disp non-linear analysis instead of a linear analysis and then calculating your responses??

 

[DCBII]

This is a concrete wall. Linear analysis is almost always used for this (see the Portland Cement Association's handbooks on concrete tank design). It's very arduous to model a composite anisotropic material such as cracked reinforced concrete nonlinearly. Usually the linear analysis is sufficient and the complex nonlinear effects are accounted for in the ACI design equations and effective stiffnesses codes prescribe.

 

[nlgyro]

okay I did a quick hand-calc using your panel dimensions/pressure-loading/properties and checked using Timoshenko and must admit the peak deflection values are very small compared against the panel thickness. So a linear analysis would suffice.

One thing that caught my eye in your explanation is you claim the FEA deflections results are higher than what Timoshenko predicts. This could be due to the fact that shell elements use transverse shear flexibility into their stiffness formulation whereas Timoshenko only considers the plate to be a bending panel. The presence of transverse shear flexibility (in addition to membrane and bending stiffness) makes the panel more flexible and thereby could increase the peak panel deflections. Is there a way you can ignore/turn-off the transverse shear flexibility option and run the panel just as a bending panel (meaning membrane plus bending stiffness only included)??

 

[IDS]

I have run the model as described in a commercial FEA package (Strand7), and compared with results from Roarke Table 11.4.8, which should be close to Timoshenko results.

I ran with a Poisson's Ratio of 0.17 and also 0.316, which is the default value given in Roarke. Results were (in SI units, sorry):

Max top face stress at mid-span
Stress, kPa FEA/Roarke
PR=.17 2.46E+03 0.993
PR=.316 2.45E+03 0.992

Deflection at mid-span
Deflection, mm FEA/Roarke
PR=.17 4.43E-04 1.067
PR=.316 4.11E-04 0.989

So stresses were very close to Roarke, regardless of the Poisson's Ratio, and deflections were only 7% greater with the low Poisson's Ratio.

I can only echo rb1957's suggestions.

Are you able to post the code for your plate element (or a link if there is an on-line source)?

Doug Jenkins
Interactive Design Services

http://newtonexcelbach.wordpress.com/

 

[nlgyro]

I finally got around to using MITC4 in ADINA. Modelled the OP's problem with his dimensions/properties/loading and compared the results with Timoshenko, and this is what I get:

dmax = max. deflection at the center of the panel (in inches)
Smax = max. bending stress at x=a/2, y=0 (using Timoshenko's co-ordinates) (in psi)

Timoshenko (dmax/Smax) = 0.00161/86.354
MITC4 (dmax/Smax) = 0.00176/73.473
MITC4/Timoshenko = +9%/-15%

MITC4 over-estimates the deflections by +9% and under-estimates the peak bending stress by -15% essentially matching the OP's claim in one of his earlier posts.

Now for the given dimensions, this plate has a L/t=20 making it borderline between thin and thick. It would be interesting to study the effects of transverse shear flexibility in such a scenario. The MITC4 documentation in ADINA does not give me the option to turn-off transverse shear effects (essentially moving between a thick to a thin plate modelling). I use NX/NASTRAN extensively. The CQUAD4 element present in NASTRAN gives me the option to easily switch b/w thin and thick plate modelling. So I tried to benchmark where the CQUAD4 element stands vis-a-vis standard MITC4 for the given problem. So here are the results:

Timoshenko (dmax/Smax) = 0.00161/86.354
MITC4 (dmax/Smax) = 0.00176/73.473
CQUAD4 (dmax/Smax) = 0.00181/73.413
MITC4/Timoshenko = +9%/-15%
CQUAD4/Timoshenko = +12%/-15%

The CQUAD4 is pretty much equivalent to the MITC4.

The above FEA results were obtained using transverse shear flexibility.

Timoshenko does not take into account the transverse shear effects in his solutions, so I re-ran the CQUAD4 FEA without the transverse shear effects and here are the results:

Timoshenko (dmax/Smax) = 0.00161/86.354
CQUAD4 (dmax/Smax) = 0.00161/73.531
CQUAD4/Timoshenko = 0%/-15%

As can be seen the deflection results perfectly match with that of FEA thereby showing the transverse shear flexibility of the plate has a role to play due to the L/t ratio. Not accounting for this is a limitation in closed form solution which the use of the FEA (in this case) gives the 'true' deflections.

So I would essentially go with the deflection results of the FEA in this case.

Now the underestimation of the stresses is purely due to the slow convergence of the 4-noded shell element. Both MITC4 and CQUAD4 has this issue. This only means that the current mesh size of 6 inches is not sufficient in the regions of high Mc/I stress and further refinement either in terms of a finer 4-noded shell mesh or switching to a 8-noded shell is required for convergence.

One useful tool CQUAD4 offers its users is the option of cubic-bending correction. This uses the grid point displacements and rotations to form a cubic equation for extrapolating the forces/stresses from centroids to the nodes. This works very well when using a coarse grid FEM to get internal forces/stresses of great accuracy.

So using the same mesh size of 6 inches and using this option of recovering nodal stresses on CQUAD4, I get:

Timoshenko (dmax/Smax) = 0.00161/86.354
CQUAD4 (dmax/Smax) = 0.00161/86.232
CQUAD4/Timoshenko = 0%/0%

So the stresses perfectly match. MITC4 does not use this option for force/stress extrapolation, so the only option is use a fine-grid model or a 8-noded shell model while dealing with MITC4 for stress convergence.

 

[rb1957]

does "cubic correction" mean mid-side nodes ? (CQUAD8) ?

transverse shear effects are real (yes?), which Timoshenko neglected in his analysis. Has anyone analyzed a plate in bending including transverse shear effects ?

It would seem like we're (quite correctly) trying to reverse engineer Timoshenko's analysis in FEA, and then roll back there effects (to arrive at a better answer than Timoshenko … wot, better than Timoshenko??!!) ?

another day in paradise, or is paradise one day closer ?

 

[IDS]

I have now found the Timoshenko solution to this problem. For moments in the short span direction Roarke and Timoshenko agree exactly. Roarke doesn't give any results in the long span direction and I'll comment on the Timoshenko long span results later.

For deflections the Timoshenko solution include the effect of Poisson's Ratio, and the results with PR = 0.3 are very close to the Roarke results (about 0.1% difference). Although the main text in Roark says the rectangular plate results used a PR of 0.316, the start of the tables says 0.3, which based on the deflection values seems to be right.

I have also re-run my FEA with imperial units as in the original post, and with a PR of 0.17 and 0.30.

The results with PR = 0.17 agree almost exactly with nlgyro's last run, ie:
dmax = 1.608E-03 in
smax(edge) = 86.17 psi

smax(mid-span) = 42.63 (Timoshenko = 42.92)

With PR = 0.3 the deflection reduces to 1.510E-3 (Timoshenko = 1.512)

So the FEA results with cubic bending correction, and without shear deflections, agree almost exactly with Timoshenko results for deflection and moments/stresses in the short span direction (but see below for the other direction).

In summary, 4 node elements with dimensions as in the OP can give good results for bending, but only if the cubic bending correction is applied.

For the example given including shear deflections increased the maximum deflection by about 10%. I'd suggest that this is a bit academic for a concrete structure, since the combined effect of concrete cracking, shrinkage and creep may well be more than an order of magnitude greater.

Strand7 does not have an option for including shear deflections in the plate results, but using beam elements to model a fixed end span of 120 in, including shear deflections had the following results:
PR = 0.17: deflection increased by 11.2%
PR = 0.30: deflection increased by 12.5%
which is consistent with the increased deflections found with the MITC4 elements.

Regarding the moments and stresses in the long-span direction, I found a paper at:

https://www.nafems.org/downloads/nbr4.pdf

looking at stresses in a rectangular steel plate, simply supported along the short edges, comparing experimental results with Timoshenko and FEA results.

In summary, the paper found a significant difference between the stress ratio for mid-span in the X and Y directions as found from experiment and FEA (about 10) compared with Timoshenko (about 12). It suggests that this was probably a typo, but in my analysis I found a similar discrepancy for mid-span moments in the long direction between FEA and Timoshenko (FEA results were about 10% lower for PR = 0.17 and 15% higher for PR = 0.3). Since the Timoshenko table does not adjust any of the stress results for Poisson's Ratio, whereas the FEA results show a large difference in the long span stresses at midspan, this suggests the difference is a shortcoming in the theoretical values, rather than a typo.

Doug Jenkins
Interactive Design Services

http://newtonexcelbach.wordpress.com/

 

[IDS]

does "cubic correction" mean mid-side nodes ? (CQUAD8) ?

No, my analysis was with the Strand7 default 4 Node elements, and I assume nlgyro's last results were 4 node elements as well.

I don't know the details of how it works, but since it seems to give significantly better results I don't see the reason for not applying this correction

Doug Jenkins
Interactive Design Services

http://newtonexcelbach.wordpress.com/

 

[IDS]

I just tried the same model with 8 node plates:

Moments and stresses: no significant changes from 4 node, in either direction.
Deflections: increased to 1.801E-3 at mid-span.

Doug Jenkins
Interactive Design Services

http://newtonexcelbach.wordpress.com/

 

[rb1957]

yes, I'd've been surprised if mid-side nodes helped on thru-thickness bending.

this must be quite a clever element compared to the CQUAD, which I thought performed badly for thru-thickness bending.

I think it'd still be worth running a single element test on this guy.

If we're tracking down a typo in Timoshenko, can't we (someone) redo the math to find the typo ?

another day in paradise, or is paradise one day closer ?

 

[IDS]

Some more background on the discrepancy between my results and Timoshenko for the long span moments:

Checking the Strand7 Theoretical Manual, it confirms that their 8 node plate element includes shear deflections, and the 4 node doesn't. Since the 4 node results are in good agreement with the Timoshenko deflections, and the 8 node results are in good agreement with nlgyro's results including shear deflection this suggests the Timoshenko numbers are correct for bending deflections only, ignoring shear strain.

For the bending moments, the Timoshenko table is based on a Poisson's Ratio of 0.3, and gives no adjustment for different values. Clearly (from symmetry) square plate would have equal bending moments in both directions, so there would be no adjustment in the moments for Poisson's Ratio. For a rectangle with side ratio of 2 though, the shape of the bending moment diagrams in each direction is very different:

Checking my results for Poisson's ratio = 0.3, the My (long span moment) is only about 3% greater than the Timoshenko value. I had previously compared the maximum moment, which occurs just past 1/4 span, rather than the mid-span value. This suggests that the Timoshenko My is a reasonable approximation for Poisson's ratio = 0.3, but the 0.17 My is significantly lower.

I have also plotted the midspan deflections along the x and y axes to a normalised scale, again showing very different shapes for the two directions, which suggests to me that there is no simple "exact" solution for rectangular plates:

In summary:

1) The proposed element size can give good results for both deflections and bending moments, but it needs a suitable plate element. The MITC4 element doesn't seem to give good results.

2) The Poisson's ratio has a significant effect on the mid-span moment in the long span direction, but only a small effect on the others (at least for the 2:1 shape examined).

3) Allowing for shear deflections has a significant effect, but for a concrete structure this will be minor compared with other approximations, such as the effects of concrete cracking, creep and shrinkage.



Doug Jenkins
Interactive Design Services

http://newtonexcelbach.wordpress.com/

 

[rscassar]

Are you putting your code to github. I've been trying to get general solution to 4-node quadrilateral plate bending element together but can't get the link between the polynomials matrix to the matrix. Have a closed form solution when the quadrilateral is a rectangle, but can't get it to work for general form.

 

[DCBII]

IDS you are brilliant. Thank you so much. You went above and beyond the call of duty! I wish I could give you 10 stars.

rscassar, I have posted my code to GitHub here: Quad3D.py This element still needs some polishing to finish it and make it user friendly, but it sounds from IDS's posts that my formulation is performing correctly so far. It sounds like it's not the formulation causing my problems, but the assumptions inherent to this element and Timoshenko's solutions.

Polynomial formulations work for rectangular elements, but I found for quadrilateral elements they are impossible to derive. I handed the matrix math to Sympy to crunch a solution and it choked. I tried again and left it solving overnight and it still choked. I've since learned that for this element you must abandon the polynomial displacement function formulation and use an isoparametric formulation and numerical integration at gauss points. That is where the matrix comes from. It's actually easier than it sounds. You can look at my element and see how I did it. While I've used a lot of programs with plates, I've never derived them myself until now. I bought several books on FEA and spent the last few months trying to understand how this works. Authors in this field have a hard time dumbing this down enough for people like me, but I think I finally got it. Hopefully I've provided something that can be educational for you and others. Now that I know what's going on with this element, I'll soon have these MITC4 elements implemented as "Quad3D" elements in my program.

This program is a hobby I've been doing for quite some time in an effort to learn finite element analysis better and to learn Python. The overall project is called "PyNite" (a stupid play on words I thought was clever: Python + Finite) and can be found here here. If you're familiar with pip, you can install easily using "pip install PyNiteFEA". It works well for 3D frame structures and rectangular plates. I've been trying to add general plate elements (quads) lately, hence this thread. Feel free to question anything in this program.

 

[IDS]

DCBII - Thanks. I have now downloaded from your PyNite site, which is worth far more than 10 stars



Doug Jenkins
Interactive Design Services

http://newtonexcelbach.wordpress.com/

 

[IDS]

This paper might be of interest:
A new 4-node MITC element for analysis of two-dimensional solids and its formulation in a shell element
Yeongbin Ko a, Phill-Seung Lee a,⇑, Klaus-Jürgen Bathe

https://www.sciencedirect.com/science/article/abs/pii/S0045794917307757

Doug Jenkins
Interactive Design Services

http://newtonexcelbach.wordpress.com/

 

[FEA way]

Great, this PyNite is one of the very few open-source codes that can solve plates. My friend recommended it to me a few days ago. I must give it a try. Will it support uniformly distributed loads ?