PDE Toolbox FEA Units

Perhaps I am not understanding the units analysis correctly, but what is the default unit behavior for the PDE toolbox? From my research it should be unitless/based upon your input values, correct?

My max stress seems to be off by a multiplier of 100 when compared to some other FEA tools... am I just missing something with my unit conversions or is there some underlying plot and VonMises stress behavior I don't understand?

Secondly, since my input geometry is mm and my input force is N I should get stress in N/mm^2, or MPa.... however, this can't be possible because these stress values would be astronomical! If they are in Pa then I am off by a factor of 100... if they are in MPa then it is FAR above what should be expected....

I looked at a list of the stress values, to ensure there are no nodal stress spikes from bad mesh, and it all looks smooth and continuous.

Answer 1

[Ravi Kumar]

Hello Colton,

You are right, PDE Toolbox does not have any units. Two things you can check again:

1. You save E as 10300000 psi in the comment, which converts to about 7.1E10 Pa, this not the value you use for E.
2. Note that surface traction also has the units of Pa, if your model is in mm, then you are applying 4448 N/mm^2, not 4448 N.

If you are exporting the STL from SolidWorks, you can set the unts to be in m to avoid all these conversion from m to mm.

Regards,

Ravi

[colton campbell]

Ravi, thank you, I think my problem is the surface traction! I thought it was a force that the software interpolated across the surface, rather than me defining the pressure.

The old E is from something else and is, as you mentioned, incorrect.

But I think the surface traction will fix my issue, I will give it a shot thank you!!

[colton campbell]

Ravi, perhaps you can help me with the next problem... I updated my surface traction value by dividing it by the area (in meters, 0.000629 m^2) and then applying it in the same fashion as before. My displacement values are extremely close to what a NASTRAN solver gets me, but my stress values are still off by a factor of x10 for some reason.... I verified this against two NASTRAN variations and got the same result.

What is going on here?

[Ravi Kumar]

Hi Colton,

Few more things to check:

1. Are the stress differece uniform or the max values?
2. How different are the mesh between the two models you are comparing?
3. Are the element formlation (quadratic tet) same?

I do see minimum quality is very bad, is the difference you are observing corresponding to those elements?

Regards,

Ravi

[colton campbell]

Ravi, the differences are pretty uniform amongst the stress distributiont throughout the part. Mesh size is almost the same, and refining both meshes didn't help. The max stress values don't appear to be coming from the poor aspect ratio elements.

I do believe the NASTRAN solvers are also using the Quad Tet mesh elements as well.

I am inclined to believe it is a unit error, except that the displacement values seem spot on!

[Mario Malic]

There is a thing that NASTRAN might assume differently than PDE toolbox, which may account for closeness of displacement results, but I doubt that relates to such significant difference in stress values.

NASTRAN by default uses incompatible mode elements, at least for 4-node shells, but might also be the case with your type of mesh 8-node solid?

I see also that you have a triangular mesh in MATLAB, which "makes" your model "stiffer"? Quotas, since I am not sure about this statement.

[colton campbell]

To eliminate the NASTRAN questions, I just did the same setup in Solidworks (COMSOL), with similar mesh size, and I also get max stress on in the order of 10^7 Pa, similar to the NASTRAN solvers.

I am really not sure what is going on in the Matlab solver to make the stress be 10x higher. And it isn't just a nodal spike either, the average stress over the max 100 nodes is 3.5e+08, 10x higher than the other 3 solvers

[Mario Malic]

Type of mesh is same/similar in all your analyses?

[Ravi Kumar]

Hi Colton,

When you say displacements agree, how close are they? Derivative could get amplified.

Regards,

Ravi

[colton campbell]

Yes, quad-tet mesh in all programs

Displacements (x max)

COMSOL: -0.1736mm

NASTRAN: -0.23mm

Matlab: -0.27mm

Stress Max (Von Mises, no stress spikes)

COMSOL: 2.28e+07 Pa

NASTRAN: 2.48e+07 Pa

Matlab: 4.98e+08 Pa

SIMPLE BEAM MODEL (10x10x100mm)

Deflection, NASTRAN: 13.35mm

Deflection, MATLAB: 13.54mm

Stress, NASTRAN: ~1.0e+08 Pa

Stress, MATLAB: Max (spike) 1.9e+09, average (100 nodes) 3.0e+08

So obviously we can rule out the stress spikes, but the average is still 3x higher than other solvers...

BY Hand Calc for Beam

Max stress = moment*deflection(max)/inertia

Max Stress = 200N*100mm*13.5mm/833.3mm^4

Max Stress = 324 MPa = 3.24e+08 Pa

Very very close to what Matlab gives me when averaging VonMises stress across max 100 nodes....

So maybe the problem with the other analysis was the STL conversion of the geometry...?

[Ravi Kumar]

Hi Colton,

The beam calculation is a good idea. However, I think that might not exactly pinpoint what is the source of difference. Can you try this script and tell me what is the results your obtain for max VM stress:

model = createpde('structural');
model.Geometry = multicuboid(0.5,0.1,0.1);
figure
pdegplot(model,'FaceLabels','on')
generateMesh(model);
figure
pdemesh(model)
structuralBC(model,'Face',5,'Constraint','fixed');
structuralBoundaryLoad(model,'Face',2,'Pressure',1E6);
structuralProperties(model,'YoungsModulus',210E9,'PoissonsRatio',0.28);
R = solve(model);
figure
pdeplot3D(model,'ColorMapData',R.VonMisesStress)
disp('Max von Mises stress')
max(R.VonMisesStress)

Also, is it possible for you to compare this beam results to other tools you are using? Please also tell which version of MATLAB you are using. In many tools, reduced integration is enabled by default that could have a significant impact on the stress, so does any any special averaging mechanisms.

I looked at your original model, I notice that both your loaded and constrained faces are cylindrical. In this thread we also discussed dividing the load by area. Note that STL is a tessellated geometry, if you calculated area using perfect cylinder kind of face, then it brings in some discretization error. And second descretization error is about the100 or so bad elements in the mesh with mesh quality < 0.3. Would it be possible to share the original CAD model. If not, can you try to export the finer STL and see how big an impact does STL refinement (not mesh refinement) on the values.

[colton campbell]

In response to the STL quality, running the simulation with a higher tesselated STL didn't make any change for me unfortunately. On a side note, it is a big pain to select individual faces from a plot to choose constraint locations in complex geometry... is there no better way to do this?

I am using 2020b

When I run your code the stress output is: ans = 8.1630e+07