Industry Best Practices & Commercial Solvers

CHAPTER 14Advanced FEA

Industry Best Practices & Commercial Solvers

Physical Intuition and the Engineering Reality

As you transition from academic formulations to professional engineering practice, you will rarely write finite element codes from scratch. Instead, you will design models using commercial FEA software suites such as ANSYS, Abaqus, MSC Nastran, Siemens Simcenter, or Altair OptiStruct.

These commercial packages are incredibly powerful: they contain highly optimized solvers, advanced sparse solver algorithms, and robust non-linear contact and material models. However, it is critical to remember that commercial solvers are just calculators: they will happily solve an incorrect model and output beautiful, colorful, but completely wrong stress plots. This is the classic engineering warning: *garbage in, garbage out*.

An analyst must act as a gatekeeper of physical correctness. You must translate a real-world physical assembly (consisting of complex geometries, contact interfaces, bolts, welds, and material behaviors) into a clean, mathematically stable, and computationally efficient finite element model. This requires master-level knowledge of:

Mesh Quality and Convergence: Rushing to view results with a single, unverified mesh is the most common error in industry. You must prove that your results are mathematically independent of the size of the elements in your mesh.

Geometric Singularities: Sharp corners and point loads create infinite stress fields in mathematical models. You must distinguish between physical stress peaks and numerical singularities.

Solver Dynamics: Knowing when to use direct sparse solvers vs. iterative solvers, and when to run an implicit vs. explicit analysis.

Core Concepts

1. Element Quality Criteria (Diagnostics Lab)

When a continuous CAD geometry is meshed, elements are stretched, sheared, and twisted to fit the boundary shapes. If elements are distorted too much from their ideal shapes (equilateral triangles, perfect squares, cubes), the mathematical coordinate mapping (isoparametric transformation) becomes inaccurate. This distortion introduces numerical errors in stress calculations.

Analysts use several key geometric metrics to check mesh quality:

A. Aspect Ratio

The aspect ratio is the ratio of the largest dimension of an element to its smallest dimension. For a perfect square element, the aspect ratio is $1.0$.

•

Calculation: For a quadrilateral, it is the ratio of the maximum node-to-node distance to the minimum distance.

•

Limit: In industry, elements in critical, high-stress regions should have an aspect ratio close to $1.0$. In non-critical regions (e.g., long uniform sections of a beam), aspect ratios up to $5.0$ or $10.0$ are acceptable.

B. Jacobian Ratio (Jacobian Determinant)

As discussed in Chapter 10, the Jacobian matrix $[J]$ maps natural coordinates $(\xi, \eta)$ to physical coordinates $(x, y)$. If an element is distorted, the scaling varies across the element.

•

Jacobian Ratio: The ratio of the maximum determinant $|J|_{max}$ to the minimum determinant $|J|_{min}$ evaluated at the integration points.

•

Limit: A perfect element has a ratio of $1.0$. If $|J|$ becomes zero or negative anywhere within the element, the coordinate transformation is invalid (indicating the element has folded over itself), causing the solver to crash. Industry guidelines recommend keeping the Jacobian ratio below $30$ for general analysis and below $10$ in critical regions.

C. Skewness

Skewness measures how close a triangular or quadrilateral element is to its ideal equilateral shape.

•

Calculation: It is the angular deviation from orthogonal angles (e.g., $90^\circ$ for quadrilaterals, $60^\circ$ for triangles).

•

Limit: Skewness should be kept below $0.5$ in critical stress regions. Skewness above $0.9$ indicates highly distorted elements that will degrade solution accuracy.

D. Warping Factor

Warping occurs when the four nodes of a quadrilateral element do not lie in the same flat plane.

•

Calculation: The normal distance of the nodes from the best-fit flat plane.

•

Limit: A flat quad has a warping factor of $0.0$. Warping factors above $0.2$ can lead to severe stiffness errors, particularly in shell structures.

2. Stress Concentrations vs. Singularities

One of the most common pitfalls in FEA is attempting to converge stress values at a numerical singularity.

A. Physical Stress Concentration

A stress concentration occurs when the flow of stress is disrupted by a change in geometry, such as a circular hole, a shoulder fillet, or a keyway.

•

Behavior: As the mesh is refined, the peak stress value increases, but it eventually converges to a finite, physical value. This value matches classical elasticity theory (e.g., $3 \sigma_{\infty}$ for a hole in a plate).

•

Remedy: Mesh refinement is required to capture the steep stress gradients around the fillet or notch.

B. Numerical Stress Singularity

A singularity is a point in the model where the mathematical stress goes to infinity.

•

Behavior: As you refine the mesh (making elements smaller and smaller), the peak stress never converges; it continues to increase without bound:

$$\lim_{h \to 0} S_{peak} = \infty$$

•

Causes:

Sharp Re-entrant Corners: A corner with a fillet radius of zero ($r=0$). According to elasticity theory, stress concentration factor $K_t \propto 1/\sqrt{r}$. If $r \to 0$, stress goes to infinity.

Point Loads and Point Constraints: Applying a force $F$ at a single node. Since a node has zero area ($A=0$), the pressure is $P = F/A = \infty$.

Sharp Re-entrant Corners: A corner with a fillet radius of zero ($r=0$). According to elasticity theory, stress concentration factor $K_t propto 1/sqrt{r}$. If $r o 0$, stress goes to infinity.

Point Loads and Point Constraints: Applying a force $F$ at a single node. Since a node has zero area ($A=0$), the pressure is $P = F/A = infty$.

•

Remedy: Do not attempt to converge stress at a singularity! Instead:

•

Add realistic fillet radii to sharp corners in high-stress zones.

•

Distribute point loads over a small, realistic patch area using rigid connector elements (RBE3/MPC) or pressure bounds.

3. Solver Types and Settings

Commercial packages solve the global system $[K]{u} = {f}$ using two primary classes of solvers:

A. Direct Solvers (e.g., Direct Sparse Solver)

Direct solvers solve the system using mathematical elimination (similar to Cholesky decomposition $LDL^T$).

•

Characteristics: Extremely robust, handles poorly conditioned matrices (e.g., structures with very thin elements next to thick ones).

•

Scaling: Computations scale with the bandwidth of the matrix ($O(N cdot b^2)$).

•

Best For: Small to medium models (under 1 million degrees of freedom), shells, and complex non-linear contact problems.

B. Iterative Solvers (e.g., Preconditioned Conjugate Gradient - PCG)

Iterative solvers start with a guess and refine the solution step-by-step until the residual error is below a tolerance.

•

Characteristics: Requires significantly less memory than direct solvers.

•

Scaling: Computations scale almost linearly with model size ($O(N)$).

•

Best For: Large, bulky 3D solid meshes (over 2 million DOFs) with good element shapes.

Step-by-Step Derivations

Richardson Extrapolation for Mesh Convergence

When performing a mesh convergence study, we solve the model on a sequence of refined meshes to estimate the exact, infinite-mesh-density solution.

Let $S_h$ be a physical quantity (such as peak stress or deflection) calculated on a mesh with characteristic element size $h$. The mathematical relationship between the numerical solution and the exact analytical solution $S_{exact}$ can be expressed using a Taylor series expansion:

$$S_{exact} = S_h + C_1 h^p + C_2 h^{p+1} + \dots$$

For a sufficiently refined mesh (where the element size $h$ is small), we can neglect higher-order terms:

$$S_{exact} \approx S_h + C h^p$$

where:

•

$C$ is a constant related to the geometry and formulation.

•

$p$ is the asymptotic rate of convergence. For a linear element, $p \approx 1$ for stresses; for a quadratic element (e.g., LST, Q8), $p \approx 2$.

1. Determining the Rate of Convergence $p$

To determine $p$, we run three nested meshes with a constant refinement ratio $r$:

$$r = \frac{h_{coarse}}{h_{medium}} = \frac{h_{medium}}{h_{fine}}$$

Let the three numerical solutions be:

•

$S_1$ (on the coarse mesh, $h_1 = r^2 h$)

•

$S_2$ (on the medium mesh, $h_2 = r h$)

•

$S_3$ (on the fine mesh, $h_3 = h$)

Using the convergence relation:

$$S_{exact} \approx S_1 + C (r^2 h)^p$$

$$S_{exact} \approx S_2 + C (r h)^p$$

$$S_{exact} \approx S_3 + C h^p$$

Subtracting the equations to eliminate $S_{exact}$:

$$S_2 - S_1 = C h^p (r^p - r^{2p}) = C h^p r^p (1 - r^p)$$

$$S_3 - S_2 = C h^p (1 - r^p)$$

Taking the ratio of these differences:

$$\frac{S_2 - S_1}{S_3 - S_2} = \frac{C h^p r^p (1 - r^p)}{C h^p (1 - r^p)} = r^p$$

Taking the natural logarithm of both sides allows us to calculate the actual rate of convergence $p$:

$$p = \frac{\ln\left( \frac{S_2 - S_1}{S_3 - S_2} \right)}{\ln(r)}$$

2. Extrapolating the Converged Solution

Once $p$ is determined, we can calculate the extrapolated exact solution $S_{exact}$ by substituting $C h^p = \frac{S_3 - S_2}{r^p - 1}$ into the equation for $S_3$:

$$S_{exact} \approx S_3 + \frac{S_3 - S_2}{r^p - 1}$$

Or, in terms of coarse and fine designations:

$$S_{exact} \approx S_{fine} + \frac{S_{fine} - S_{coarse}}{r^p - 1}$$

This equation is the Richardson Extrapolation formula. It provides an estimate of the converged solution without requiring you to run an extremely dense, computationally expensive mesh.

3. The Grid Convergence Index (GCI)

To quantify the uncertainty in your mesh convergence, ASME standards recommend calculating the Grid Convergence Index (GCI). The GCI is an estimate of how far the current fine mesh solution is from the converged limit:

$$\text{GCI}_{fine} = F_s \frac{|e_{fine}|}{r^p - 1}$$

where:

•

$e_{fine} = \frac{S_{fine} - S_{coarse}}{S_{fine}}$ is the relative change between meshes.

•

$F_s$ is a safety factor. If you only have two meshes, use $F_s = 3.0$. If you have three nested meshes and have verified the rate of convergence $p$, use $F_s = 1.25$.

Solved Example: Richardson Extrapolation and GCI Verification

Problem Statement: An aerospace engineer is designing a turbine shaft shoulder fillet. To verify the safety factor, she runs three models with nested meshes of quadratic elements (using a constant refinement ratio $r = 2.0$). The peak stresses calculated are:

•

Mesh 1 (Coarse, $h = 4.0\text{ mm}$): $S_1 = 285.0\text{ MPa}$

•

Mesh 2 (Medium, $h = 2.0\text{ mm}$): $S_2 = 312.0\text{ MPa}$

•

Mesh 3 (Fine, $h = 1.0\text{ mm}$): $S_3 = 319.5\text{ MPa}$

Perform a mesh convergence check:

Calculate the actual rate of convergence $p$.

Calculate the extrapolated asymptotic peak stress $S_{exact}$.

Calculate the Grid Convergence Index $\text{GCI}_{23}$ for the fine mesh and verify if the solution lies within the asymptotic range.

Step 1: Calculate differences between solutions

$$\Delta_{12} = S_2 - S_1 = 312.0 - 285.0 = 27.0\text{ MPa}$$

$$\Delta_{23} = S_3 - S_2 = 319.5 - 312.0 = 7.5\text{ MPa}$$

Notice that the difference is decreasing ($27.0 \to 7.5$), which indicates that the solution is converging.

Step 2: Calculate the rate of convergence $p$

Using the refinement ratio $r = 2.0$:

$$p = \frac{\ln\left( \frac{S_2 - S_1}{S_3 - S_2} \right)}{\ln(r)}$$

$$p = \frac{\ln\left( \frac{27.0}{7.5} \right)}{\ln(2.0)} = \frac{\ln(3.6)}{\ln(2.0)} \approx \frac{1.2809}{0.6931} \approx 1.848$$

The calculated rate of convergence is $p \approx 1.85$. This is close to the theoretical value of $2.0$ for quadratic elements, confirming that the mesh is converging as expected.

Step 3: Calculate the extrapolated exact stress $S_{exact}$

Using the Richardson Extrapolation formula with the fine mesh ($S_3 = 319.5\text{ MPa}$):

$$S_{exact} \approx S_3 + \frac{S_3 - S_2}{r^p - 1}$$

$$S_{exact} \approx 319.5 + \frac{7.5}{2.0^{1.848} - 1}$$

Let us compute $2.0^{1.848}$:

$$2.0^{1.848} \approx 3.600$$

Substitute this back:

$$S_{exact} \approx 319.5 + \frac{7.5}{3.6 - 1} = 319.5 + \frac{7.5}{2.6} \approx 319.5 + 2.88 \approx 322.38\text{ MPa}$$

The extrapolated asymptotic peak stress is $322.4\text{ MPa}$.

Step 4: Calculate the Grid Convergence Index (GCI)

First, calculate the relative error between the medium and fine meshes:

$$e_{23} = \left| \frac{S_3 - S_2}{S_3} \right| = \left| \frac{319.5 - 312.0}{319.5} \right| = \frac{7.5}{319.5} \approx 0.02347 \quad (2.35\%)$$

Using a safety factor $F_s = 1.25$ (since we have three nested meshes and have verified $p$):

$$\text{GCI}_{23} = F_s \frac{e_{23}}{r^p - 1} = 1.25 \frac{0.02347}{2.6} \approx 1.25 \times 0.00903 \approx 0.01129 \quad (1.13\%)$$

The GCI of the fine mesh is $1.13\%$. This small value indicates that further mesh refinement will change the solution by less than $1.13\%$, indicating a highly refined, convergent mesh.

Step 5: Verify Asymptotic Range

For three meshes to be in the asymptotic range of convergence, the ratio of GCIs must satisfy:

$$\frac{\text{GCI}_{12}}{r^p \cdot \text{GCI}_{23}} \approx 1.0$$

Let us calculate the relative error between the coarse and medium meshes:

$$e_{12} = \left| \frac{S_2 - S_1}{S_2} \right| = \left| \frac{312.0 - 285.0}{312.0} \right| = \frac{27.0}{312.0} \approx 0.08654 \quad (8.65\%)$$

Now, calculate the GCI for the coarse-medium step:

$$\text{GCI}_{12} = F_s \frac{e_{12}}{r^p - 1} = 1.25 \frac{0.08654}{2.6} \approx 0.0416 \quad (4.16\%)$$

Check the ratio:

$$\frac{\text{GCI}_{12}}{r^p \cdot \text{GCI}_{23}} = \frac{0.0416}{2.0^{1.848} \times 0.01129} \approx \frac{0.0416}{3.6 \times 0.01129} = \frac{0.0416}{0.0406} \approx 1.02$$

Since the ratio ($1.02$) is very close to $1.0$, the mesh convergence study is within the asymptotic range, validating the convergence of the stress calculations.

Extended chapter reference

Model planning, verification, validation, convergence, and decision-ready reporting

Professional FEA is an evidence chain rather than a contour plot. The analyst must document idealization choices, verify numerical implementation, validate physical assumptions, and communicate margins with traceable uncertainty.

Create an analysis plan linking requirements to model outputs.

Separate verification, validation, calibration, and sensitivity studies.

Design mesh-convergence studies around decision quantities.

Report assumptions, uncertainty, margins, and limitations clearly.

Visual map

From formulation to verified result

Each stage produces evidence for the next. A reliable result preserves traceability from the physical idealization to the final engineering check.

Equation sheet

Governing relationships

Safety factor

$$n=\frac{S_{allow}}{S_{demand}}$$

Useful only when demand and allowable definitions match.

Relative mesh change

$$e_h=\frac{|Q_h-Q_{h/2}|}{|Q_{h/2}|}$$

Tracks convergence of the decision quantity Q.

Combined uncertainty

$$u_c=\sqrt{\sum_i u_i^2}$$

Root-sum-square estimate for independent uncertainty sources.

Engineering comparison

Selection and interpretation table

Activity	Question answered	Evidence	Owner
Verification	Did we solve the equations right?	Benchmarks, residuals, convergence	Analyst
Validation	Are these the right equations?	Test correlation	Analysis and test teams
Calibration	Which uncertain inputs fit data?	Parameter estimation	Model owner
Prediction	What happens in the target case?	Validated model plus uncertainty	Design authority

Verification and application

Checks before accepting the result

Loads and constraints have traceable sources
Decision outputs are mesh-converged
Singular peaks are not used directly
Report includes limitations and sign-off criteria

Industry application

A lifting lug assessment should combine global load-path checks, contact and bolt idealization, local mesh convergence away from singular edges, material allowables, and correlation against a hand-calculated bearing/net-section model.