Skip to main content

Defining Ansys Superelement SUB File Manually

Photo by  James Owen  on  Unsplash A surprisingly popular blog-post written here is Exporting Stiffness Matrix from Ansys . A sensible follow up question is what can one do with the exported stiffness matrix? In a recent Xansys Forum post, a question was raised on how we can edit the stiffness matrix of a superelement and use it for our model.  An approach presented below is to first create a superelement that has the same number of DOF and nodal location that will serve as a template. An APDL script can then be written to edit the stiffness matrix entries as desired before exporting to a new superelement *.SUB file for use in future models. The self-contained script below demonstrates this.  /prep7 et ,1, 185 mp , ex, 1, 200e3 mp , prxy, 1, 0.33 w = 0.1 ! single element (note nodal locations) n , 1, w, -w, -w n , 2, w, w, -w n , 3, -w, w, -w n , 4, -w, -w, -w n , 5, w, -w, w n , 6, w, w, w n , 7, -w, w, w n , 8, -w, -w, w e , 1, 2, 3, 4, 5, 6, 7, 8 /solu antype , substr     ! analy

Mesh Convergence Study

In a previous post, I mentioned mesh convergence study is not done as often it should. I've been reminding myself to be careful in reporting stress results using sanity checks. A nagging question persists in my mind: "How much should the mesh be refined before calling it good? Is a factor of 2 good enough?"

A presentation here at Ruhr-Universität Bochum is quite interesting. It presented a linear Richardson Extrapolation being:
fex = fn+ k(fm - fn)
fex is the exact solution
fm is the solution at mesh size m
fn is the solution at mesh size n
k is what I'll call Richardson Extrapolation Factor:

k = (mα) / (mα - nα)

where α=2 for linear elements and α=3 for quadratic elements. To plot this factor k relative to the reduction factor (n/m) in Figure 1:

Figure 1: Richardson Extrapolation Factor

It is interesting that for quadratic elements, with a mesh size reduction by a factor of two, the factor is 114.3%. Hypothetically, wouldn't it be great if, for example, one could mesh once with quadratic elements, get 10MPa, mesh again with half the size mesh, get 11MPa, then say with confidence the stress should converge to 11+1.143(11-10) = 12.143MPa?!

To see if this works in practice, a simple model was created. Here's the maximum stress of a simple supported steel beam with concentrated load at the center. The beam is 10mm by 10mm with a length of 100mm. 100N was applied at the mid-span. The model is made of 2D Quadrilateral elements that are square in shape.

Figure 2: Normal Stress (X Axis)

Figure 3: Maximum Bending Stress Vs Element Size

In this simple model, maximum bending stress for both the linear and quadratic elements converges to the hand calculated 15MPa. Note that even with 2 elements through the thickness, the linear element had a maximum stress of 14.025 or an error of 6.5%. That's not too shabby for a quick and dirty model. Real life design margins are practically much larger than this.

Figure 4: Normal Stress X Axis - Linear Elements

Here's the zoomed in view of the stress convergence plot:

Figure 5: Maximum Bending Stress Vs Element Size

Applying the formula in that paper, the "exact solution" was predicted higher for most cases. See the Figure 6. In the following case, the coarsest solution was used as mesh m (size of 5mm), then each consecutive solutions are individual n.

Figure 6: Exact Solution Vs Element Size Reduction Factor

In short, I would still want to do at least 3 different cases for of mesh size density factors (e.g. 1, 2 & 4) to verify the mesh result converges. Secondly, the design margins should be healthy relative to the changes in stresses during mesh convergence. And finally, be mindful of other errors that could be persistent between meshes such as erroneous material properties, boundary condition etc.

Ansys V18.1 Archived Project: link
Octave/Matlab script: link


  1. Hello Jason,

    '11+1.143(11-10) = 12.143MPa?!' mentioned in the third paragraph.

    As per formulae, it should be '11+1.143(10-11)'. Am I right?

    Also please share any practical FEA calculation done (mesh convergence study) using Richard interpolation.


    1. Hi Naveen,

      Applying your equation, the value is 9.857MPa. This goes in the opposite direction of the convergence. I am not sure if this was an error of the source article so I flipped the gradient sign to point in the right direction.

      The conclusion for this article for me is that the Richard interpolation should not solely be relied on and instead, do at least 3 different cases to best approximate converged value.



Post a Comment