WebCAE Reference Manual

1.1What WebCAE is

WebCAE is a full-cycle finite element analysis package. It covers geometry preparation (CAD import and parametric modelling), meshing, assignment of materials, boundary conditions and loads, solution, and result post-processing. Its main area is structural strength and dynamics: solid bodies, shells, beams, and combinations of them. Heat transfer, coupled thermomechanical problems, acoustics, and fluid flow are also supported.

The product runs on two platforms with the same computational core:

• desktop application (Windows, macOS): native solver, multithreading, problem size limited only by the machine's memory;
• browser: the same core compiled to WebAssembly; convenient for getting acquainted with the product and for small models, limited by the memory of the tab and by a single thread.

The native desktop path is the reference. Verification numbers are taken on it. If the results of the two platforms differ on some model, the native result is considered correct, and the discrepancy is described explicitly (see Platforms, performance and quality control).

1.2How to read this reference

The reference describes exactly what the program computes: in what formulation, with which elements and methods, under which assumptions. It does not replace a finite element textbook, but the formulas in it are taken from the implementation, not from a textbook. Where the implementation is simpler than the classical formulation, this is stated directly.

The sections follow the order of the analysis cycle:

1. Geometry and CAD and Meshes: where the analysis model comes from.
2. Material models, Boundary conditions, loads and constraints, Contact and connections: what the user specifies.
3. Finite element library: how the model is discretized.
4. Linear analysis types, Nonlinear analysis and dynamics, Heat transfer and multiphysics: which problems are solved and in what formulation.
5. Solvers: how the resulting systems of equations are solved.
6. Results, post-processing and strength assessment: what the engineer receives and how it is computed.
7. Projects, scripts and automation, Platforms, performance and quality control: how to work with all this and why the results can be trusted.

Typical structure of a subsection: purpose, formulation and equations, inputs and outputs, implementation, limitations, verification (links to verification cases), example script command.

1.3Readiness statuses

Each subsection is marked with one of three statuses:

A status is assigned conservatively: if there was doubt between two, the stricter one was chosen. Statuses are reviewed with each release.

1.4Units of measurement

Inside the solver, all quantities are in a consistent SI system: length in metres, force in newtons, stress, pressure and elastic modulus in pascals, density in kg/m³, temperature in kelvins. The solver is dimension-agnostic: it solves a system of equations, and the answer is correct only if the input units are consistent. Therefore unit conversion is performed in one place, at the boundary between the model and the solver, and nowhere else.

The user independently chooses two things:

• the input system for loads and properties: SI, metric engineering (mm, N, MPa), or imperial;
• the scene length unit in which the geometry is stored; usually millimetres, as in the source CAD.

Results are shown with labelled units. The formulas in this reference are written for any consistent system: unless stated otherwise, SI is assumed.

1.5Testing and verification

Accuracy claims rest on an open verification programme: more than two hundred problems with analytical solutions, NAFEMS and MacNeal–Harder problems, and nonlinear and dynamic benchmarks. For each problem, the formulation, the reference value, the obtained value, the deviation, and the solver version on which the run was made are published. Throughout the reference, links lead to specific cases; the full list is on the Verification page.

The tolerance for comparison with the reference in the verification programme is 5 %. It must not be loosened to make a case pass: a discrepancy larger than the tolerance is considered a defect and is fixed in the computation, not in the criterion.

1.6Related pages

• Element library: summary cards of the elements.
• DSL scripting language: complete command reference.
• MCP server: controlling the model from AI assistants.
• Verification: all cases with numbers.

An analysis in WebCAE starts with geometry: solid bodies, sheets, and lines in boundary representation (B-Rep). A body comes from a CAD file, is built from primitives and sketches, or is assembled from bars and shell patches. Everything that follows (mesh, boundary conditions, results) is attached to the faces, edges, and vertices of this representation. For this reason the section describes not only ways to obtain geometry, but also how its entities can be addressed and what happens to the attachment when the geometry changes.

Every operation goes through a single command layer. A ribbon button, a context menu item, and a DSL line execute the same command, so anything built with the mouse can be reproduced by a script (syntax on the page /docs/dsl). The exceptions are named in the subsections: midsurfaces and Onshape synchronization are available from the interface only, and file import by path from the DSL works in batch mode.

2.1Geometry kernel and unitsstable

Geometry is built and checked by the open-source Open CASCADE Technology kernel (OCCT 7.9). A body in a document can be a solid, a sheet (shell), a wire, or a point. Boolean operations, fillets, sewing, and validity checking are performed by OCCT. The kernel is built for the desktop application and for the browser (WebAssembly) from the same source code.

Units are organized in three layers.

• Geometry lives in millimetres. The coordinates of primitives, sketches, and imported models are scene units (mm).
• The solver receives data in SI (metres, newtons, pascals, kilograms). Conversion is performed at the solver boundary in a single place, and the user does not see it.
• The script unit is set by the command units m|mm|in. It applies to dimensions that are not scene coordinates: element size size=, radii and sizes of refinement zones, procedural beams, shells, and point masses. The coordinates of body-building commands (box, primitive, move) are not converted and remain scene millimetres.

units mm
box dx=20 dy=20 dz=100
list faces of body1

2.2STEP importstable

STEP (AP203, AP214, AP242) is the main exchange format. The file is read by the OCCT STEP reader; file units are converted to millimetres. The content is split into bodies: solid bodies, or shells if there are none, otherwise the whole result as a single body. Body names are taken from the STEP data where possible. Large assemblies are displayed progressively, as tessellation proceeds.

After reading, each body undergoes a validity check (analogous to BRepCheck), and the model receives one of the readiness statuses:

• valid solid: a volume mesh can be built;
• valid sheet: only a surface mesh is available;
• needs review: some bodies of the assembly have issues, the rest are usable;
• unsupported or read error: the analysis does not start, and the reason is reported as a code.

If the check was not performed for some reason, the status is marked as unchecked and is not passed off as "usable". A volume mesh of a sheet body is blocked explicitly, with a reason code, rather than silently replaced by another one.

In the interface, a file is loaded with the import button or by dragging it onto the canvas. In the DSL, the import command with a file path works in batch (node) mode; inside the application, use loading through the interface.

import "bracket.step"
list faces of body1

2.3IGES and STL importbeta

IGES is read by the OCCT kernel; after reading, surfaces are sewn into shells and, if closed, into a solid. In the DSL the format is set by the parameter format=iges or by the file extension. As with STEP, import by path from the DSL is performed in batch mode.

STL is a mesh of triangles, not a B-Rep. WebCAE reads binary and text STL, merges coincident vertices, checks watertightness (each edge belongs to exactly two triangles), and on that basis classifies the body as solid or sheet. So that faces can be addressed by number, the triangles are segmented into "faces" by the dihedral angle between neighbouring triangles (default threshold about 35°); segment boundaries become edges. For such a faceted body the volume mesh is built from the triangles themselves.

2.4Topology numberingstable

Every face, edge, and vertex receives a number. This is the same number that is shown in the model tree next to the node, and the same address that the DSL accepts: surface number N, edge number N, vertex number N. For bodies, numbering goes in the tree by body.

Numbering rules:

• a number is issued once per entity and is tied to its stable identifier in the topology registry, not to a position in an array;
• a number is not reused after deletion: the counter only grows, so that a reference in a report or script never starts pointing at a different entity;
• if a face survived an operation with the same stable identifier, its number is unchanged; new faces receive new numbers;
• the kind word (surface, edge, vertex) is mandatory and is checked: surface number 27 is rejected if 27 is a vertex;
• addresses by extent (face z-, vertex "x-,y-,z-") and bare indices are not supported and produce an explicit error. In an assembly there is no such thing as "the z- face", and such an address would silently land on someone else's entity.

Numbers can be looked up with the commands list faces of <body> and list edges of <body>: they print the number, area, centroid, and normal of each face, or the length, midpoint, and direction of each edge.

list faces of body1
fix surface number 3
load surface number 5 pressure=-1000000

2.5Primitives, Boolean operations and transformationsstable

The basic set of operations on bodies:

• primitives: box, cylinder, sphere, cone, torus, wedge, pipe, hollowBox, hemispheres (primitive <shape> parameters);
• Boolean operations boolean union|cut|intersect target= tool=, sewing and merging (sew, union);
• fillets and chamfers of edges (fillet, chamfer), a hollow body of given thickness (hollow);
• transformations: move, rotate, mirror, scale, linear and circular patterns (pattern);
• construction from a profile: sketch (contour by points), extrude, revolve, sweep, loft;
• direct face editing: offset (offset_face, pull), move (move_face), and draft (draft_face).

A Boolean operation is built by OCCT with history: the kernel reports which faces were modified, generated, or deleted, and numbers and attachments are carried over according to this history.

units mm
box dx=60 dy=40 dz=10
primitive cylinder r=8 h=10
boolean cut target=body1 tool=body2

2.6Parametric CAD modelbeta

A parametric model stores not a body but a recipe: a sketch, then operations on it. Editing a sketch or a dimension recomputes the chain and replaces the body in place; it keeps its name, colour, and numbers.

A sketch lies on a plane (xy, xz, yz, a reference plane, or a face) in local coordinates (u, v) in millimetres. Sketch entities: polyline and rectangle, circle, arc. Operations on a sketch: extrude (add, cut, or new body), revolve, sweep along a path sketch, loft through sections, thread (cosmetic or geometric thread). Sketch editing is available: move, rotate, mirror, contour offset, trim, and extend.

Sketch constraints include coincident points, horizontal, vertical, parallel, perpendicular, collinear, equal, concentric, midpoint, symmetry, tangent, fix, and dimensions (distances, angle, radius, diameter). Dimensions can be named and changed by number.

The constraint solver reduces the constraints to residuals \mathbf{r}(\mathbf{x}) in the sketch coordinates and minimizes \tfrac{1}{2}\lVert \mathbf{r}\rVert^{2} by the Levenberg–Marquardt method with a numerical Jacobian \mathbf{J}:

The number of degrees of freedom of a sketch is n_{\text{dof}} = n_{\text{var}} - \operatorname{rank}(\mathbf{J}). Determination by rank is correct even with duplicate constraints, where the difference "variables minus equations" would give a negative number. A sketch receives the status of fully constrained, under-constrained (with the number of degrees of freedom), or over-constrained (with a list of conflicting constraints). If the system cannot be solved, the sketch geometry is not partially changed.

Interface: the parametric operations panel (FeaturePanel), the "Operations" group on the ribbon, branches in the model tree. DSL: the cad command with subcommands sketch, extrude, revolve, sweep, loft, thread, constraint, dimension, solve, recompute, move, rotate, mirror, offset, trim, extend.

cad sketch plane=xy points=0,0,40,0,40,20,0,20 id=sketch1
cad extrude sketch=sketch1 distance=10
cad recompute

2.7Geometry healing and auto-fixbeta

An imported model is often unfit for meshing: gaps between faces, degenerate edges, small chamfers, and holes smaller than the element. WebCAE distinguishes three actions.

1. Healing on import. A budget-limited set of stages: shell sewing, shape fixing by OCCT (ShapeFix), merging of coplanar and coaxial faces (UnifySameDomain), assembly of a closed shell into a solid. The result is checked for validity, and a repeated run on the same input gives the same shape.
2. Auto-fix before meshing. Works on a copy. The source body, the model tree, and the raw STEP are not changed under any outcome. Small features (short edges, small faces, narrow slivers, holes and fillets smaller than the element) are patched over so that the mesher does not refine the mesh for them. Thresholds are tied to the element size set by the user, not to absolute millimetres. Each stage must prove that things did not get worse: the result is valid, the number of bodies is unchanged, the volume has not drifted more than allowed. Otherwise the stage is rolled back and the reason goes into the report.
3. Explicit commands. defeature <body> maxHoleRadius= maxFaceArea= removes holes and faces by your thresholds; repair <body> [merge] [heal] tolerance= merges seams and fixes gaps with a given tolerance; simplify discards small bodies of an assembly by volume and diagonal thresholds; shrinkwrap builds the outer shell of an assembly. Automatic mode is switched on and off by the command autofix on|off; autofix element=<mm> prints a report for a given element size.

Boundary conditions are not lost silently. If an auto-fix stage would remove a face to which a boundary condition is attached, the stage is rolled back to a less aggressive configuration. If the attachment is lost anyway (geometry edit, repeated import), the boundary condition record stays in the model and is flagged as orphaned, and the report lists such records by name. An analysis with an orphaned condition does not start until the attachment is restored: a mesh without the required face would look like a successful analysis of a different problem.

autofix on element=2
defeature body1 maxHoleRadius=1.5
repair body1 merge heal tolerance=0.01

2.8Imprint and splittingstable

An imprint puts a line or a plane on a face and splits the face along it; internal edges and new faces with their own numbers appear in the B-Rep. This is the main way to obtain stringer and frame lines, load application zones, or constraint regions on a sheet without changing the volume of the body. Imprint edges become mesh lines, so bars seated on them share nodes with the shell.

Command variants:

• imprint <body> points=x,y,z,…: a polyline projected onto the sheet along a given direction or along the normal (project=);
• imprint <body> plane point=x,y,z normal=x,y,z: section of a sheet or body by a plane;
• imprint <body A> <body B>: imprint of two bodies onto each other (common traces on touching faces);
• split target=<body> plane=<name> or point= normal=: splitting a body into independent parts.

The parameter fuzzy= sets the tolerance of the operation in script units.

units mm
imprint body1 points=-10,100,0, 2010,100,0
list edges of body1

2.9Assemblies of multiple bodiesbeta

An assembly is a set of bodies in one document. A body can be hidden, suppressed (excluded from the analysis), moved, and rotated. Intersection of bodies is determined by measuring the volume of their common part. Touching bodies can be glued into a single model or left separate with contact. When gluing, WebCAE aligns the topology of shared faces with a Boolean operation (fragment), and for a real volumetric overlap of an interference fit it offers three modes: cut (the volume is given to the smaller body), split into parts, or leave as is (then the material in the interference is counted twice, and a warning is shown). The pair selection tolerance and the tolerance of the Boolean operation itself are set separately.

Before the analysis, the connectivity of the assembly mesh is checked: bodies that "glue all" could not join (no contact) form separate pieces without constraints, and a linear analysis on them diverges. Such cases are caught before the solver. Details of mesh gluing are in the section "Meshes".

2.10Midsurfacesexperimental

For thin-walled parts a volume model is excessive: a 2 mm sheet is more conveniently analysed as a shell. The "midsurface" operation (interface only, the geometry operations panel) looks for pairs of nearly parallel faces whose separation lies in a given thickness range and builds a sheet body between them. Pairs of planar and cylindrical faces are found; regions meeting at corners are closed by pulling free edges to the junction planes. The thickness found for a pair is shown in the sheet name; the shell section thickness (section shell t=) must be set separately.

2.11Frames, beams and shell patchesstable

For trusses, frames, and stiffened shells, WebCAE builds the analysis model directly, bypassing the solid body. There are two ways.

Parametric frame frame. Nodes and bars are numbered from 1 in order of creation. A bar is a Timoshenko beam with a section from the library or given by dimensions (rectangle, circle, tube, I-beam, channel, angle, tee, box); a bar can be given a different section at its second end (taper) and hinges at the ends. A shell is defined by patches of quadrilateral cells: a planar quadrilateral (geom=quad), a cylinder or cone (geom=cyl), a ruled surface between two polylines (geom=ruled), or an explicit grid of points (geom=points). Patch nodes coincide with existing ones within tolerance, so the shell and bars share common nodes. Stiffeners are added by the command frame stiffen along the grid lines of a patch.

Beams on B-Rep edges. assign beam edges=… section=rect w= h= flush=outer|inner [ecc=] seats Timoshenko bars on selected edges of a sheet body (including imprint edges). The parameter flush presses the edge to the outer or inner surface of the skin, and ecc offsets the bar axis along the normal. The offset is accounted for kinematically, so the effective stiffness includes the term of the parallel axis theorem I = I_{0} + A\,e^{2} (here e is the distance from the bar axis to the midsurface of the skin), and not only the intrinsic moment of inertia of the section.

What the user gets: a model in which bars and shells are connected by shared nodes, rather than by approximate rigid constraints. The analyses are static, modal, and buckling; for modal analysis of a free body (ship hull, aircraft) there is a mode without the mechanism check, in which modes one through six are rigid-body modes.

units mm
frame node at=0,0,0
frame node at=1000,0,0
frame member 1 2 section=rect w=20 h=40
frame support node=1 dofs=all
frame load node=2 force=0,0,-500
frame solve static

2.12Onshape synchronizationbeta

The application connects to an Onshape document (OAuth 2.0 authorization), receives the model in STEP format with the version indicated (microversion), and loads it as an ordinary import. A repeated synchronization replaces the imported bodies with the new version and tries to carry boundary conditions, loads, constraints, and contacts over to the changed geometry.

The transfer is based on geometric matching: for each face, edge, and vertex of the old model, the entity of the new model closest in position, orientation, and size is sought. Records for which no reliable match was found are not deleted: they are returned in a "needs review" list, and the final message states the number of transferred and untransferred conditions. The new geometry is addressed only by document numbers.

2.13FreeCAD connectorexperimental

The FreeCAD add-on adds the command "Analyze in WebCAE". It exports the selected bodies, or all visible bodies of the document, to STEP and opens them in WebCAE: in the browser or in the desktop application via its own link. The transfer is one-way and one-time.

2.14Mesh import and exportexperimental

A ready mesh can be taken from an external program and handed out. Supported are Abaqus (.inp), Gmsh (.msh; export is format 2.2), and Nastran (.bdf, .nas). Besides meshes, the deformed surface and the topology optimization isosurface are exported (STL, STEP).

• The "Project" menu exports a mesh with bodies (SOLID_<k>) and surfaces (SURF_<N>, that is, surface number N) and imports them back. Coordinates in the file are in SI (metres).
• DSL: import mesh file="model.inp" and export mesh file="out.msh" [deformed] [scale=k]. Mesh import in the DSL reads nodes and tetrahedra, without bodies and surfaces; export writes coordinates in scene units (mm) without conversion.

import mesh file="part.inp"
export mesh file="part_out.bdf"

A mesh turns B-Rep geometry into finite elements. WebCAE builds tetrahedra (tet4, tet10) for arbitrary bodies, hexahedra, prisms, and pyramids where the geometry allows, quadrilateral and triangular shell elements on sheets, and beam elements along lines. Element types and their formulations are described on the page /docs/elements; here we cover how a mesh is built, what controls it, and how its quality and sufficiency are assessed.

Meshing works in two environments: the desktop application (a native mesher process) and the browser (WebAssembly). Both are built from the same source code, but on the same geometry the result may differ in the number of elements and in quality. The native path is the reference; the web is for reference only. Some capabilities exist only on the desktop (boundary layers, meshing of large assemblies without limits); this is stated in the subsections.

3.1Tetrahedral mesherstable

Tetrahedra are built by the open-source Netgen 6.2 mesher on top of OCCT geometry. The surface mesh is built on the B-Rep faces; the volume mesh uses the Delaunay method with advancing-front completion and subsequent optimization of element shape (swaps, smoothing); the method is described in detail in Schöberl J., NETGEN: An advancing front 2D/3D-mesh generator based on abstract rules, Computing and Visualization in Science, 1997. The mesh is conforming: nodes coincide on the common boundary of two faces.

Result of the build:

• nodes and tetrahedra of first (tet4) or second (tet10) order;
• a region (body number) for each element, so that materials are assigned by body;
• boundary triangles, each carrying the number of the B-Rep face it was generated from. For this reason conditions specified as surface number N survive mesh rebuilding: they are attached to the face, not to nodes.

The user sets the size, the transition smoothness (grading), the curvature-based density (safety), the minimum size, and the element order (size control). There are three quality profiles: "Draft" (fast coarse mesh without optimization, first order), "Standard", and "Final" (full optimization for a production analysis).

If the mesh could not be built or turned out to be of unusable quality, the mesher retries with modified parameters. This never happens silently: the application shows a banner stating which bodies were meshed with changed settings, which are coarser than the requested size, and why (profile, density, or a browser limitation). Bodies that could not be meshed are listed separately and are excluded from the analysis only with your consent.

The build can be interrupted with the "Cancel" button or the mesh cancel command. For an assembly, the partial result (the bodies already meshed) is kept.

units mm
box dx=60 dy=30 dz=10
mesh tet size=3 order=1

3.2Element order P1 and P2stable

A linear tetrahedron tet4 gives constant strains within the element and "locks" in bending and in nearly incompressible materials. A quadratic tet10 has ten nodes (four corners and six edge midpoints), the strains are linear, and the edge and face can be curved. On the same mesh, tet10 is usually noticeably more accurate than tet4; the price is the number of nodes, which for a typical tetrahedral mesh grows by roughly a factor of 7 to 8.

A P2 mesh is built in two steps. The mesher builds a P1 mesh, then adds edge midpoint nodes and moves those lying on the boundary onto the CAD surface: the edge and face of the element acquire the curvature of the source geometry. After that, each element is checked for a positive Jacobian. The determinant of the tet10 mapping is exactly a cubic polynomial, so it is expanded in the Bernstein basis of degree 3 on the tetrahedron: if all 20 coefficients have the same sign, the Jacobian is positive throughout the element (Johnen A., Remacle J.-F., Geuzaine C., Geometrical validity of curvilinear finite elements, J. Comput. Phys., 2013). In elements where the condition is violated, the midside nodes are straightened locally, and the mesh remains usable for analysis. The criterion is sufficient, not necessary: straightening may be triggered even where the element was admissible.

The order is set globally by the parameter order=1|2 of the mesh command or by the "Linear/Quadratic" switch in the mesh panel. For hexahedra, the order is set by the same switch: after a hex8 is built, midside nodes are placed at the middle of a straight edge, and on a curved edge on the surface of the primitive (hex20). Prisms and pyramids of hybrid meshes remain first order; quadratic hybrids are disabled (an explicit refusal, see hybrid zones).

units mm
mesh tet size=4 order=2

3.3Local element orderbeta

A part can be analysed on tet4, keeping quadratic elements only in the required region. The size h and the order p remain independent settings: the zone refine zone … order=2 keeps tet10 in the elements intersecting the sphere, plus one layer of transition elements, while tet4 remain outside the zone.

To keep the displacement field continuous, on the common edge of a P2 element and a neighbouring tet4 the edge midpoint is tied to its end points by a constraint equation:

where u_{a}, u_{b} are the displacements of the corner nodes of the edge, and u_{\text{m}} is the displacement of its midside node. The constraint removes the quadratic increment exactly at the boundary between the regions and is passed to the solver as an ordinary constraint equation.

units mm
mesh tet size=6 order=1
refine zone at=30,15,5 radius=8 size=1.5 order=2

3.4Element size controlstable

Element size is set from several sources simultaneously, and the outcome is determined not by priority but by the minimum. At a point \mathbf{p} the mesher obtains

where h_{\text{glob}} is the global size, h_{\text{face}} is the size set on a face, R is the local radius of curvature, s is the curvature density parameter (safety), h_{\text{zone},k} are the sizes of the refinement zones covering the point, and h_{\text{min}} and h_{\text{max}} are the lower and upper size limits. Consequences that are easy to get wrong: a size on a face is not "more important" than a zone (inside a zone the smaller one wins), and no input can make the result coarser than the global size.

Parameters:

• global size: set in millimetres or by a preset based on the characteristic dimension L: "Draft" L/10, "Standard" L/20, "Fine" L/40;
• safety s: curvature density h \le R/s, range 0.3 to 6, default 1.5. A smaller s gives a coarser mesh on fillets and holes. This is the only lever when the size is determined by curvature rather than by the slider: the interface shows how many faces are limited by curvature and predicts the number of elements;
• grading: smoothness of the transition from small elements to large ones; smaller means smoother;
• minimum size: the lower limit; 0 means automatic (a fraction of the maximum size);
• size on a body, face, edge, vertex: the "Element size" and "Local size" tabs; in the DSL, zones attached to surface number N;
• refinement zones: a sphere (centre, radius, size), refine zone at=x,y,z radius=r size=h. For a local hexahedral layer there are a cylinder, a ring around a hole, a box, a fillet, and toroidal zones (see hybrid zones);
• zone reset: refine clear.

The "Size field" mode colours bodies by the effective element size and shows which source determined it and whether it hits the lower limit. From the quality histogram, the worst elements can be selected and a local size placed on the CAD faces they touch ("Refine zone").

units mm
mesh tet size=5 safety=1.2
refine zone at=50,20,0 radius=15 size=1.5
refine clear

3.5Boundary layersexperimental

Near-wall zones (fluid dynamics, thin contact layers) need prismatic layers along the wall. Layers are defined by the number of layers n, the growth ratio q (range 1 to 5), and the thickness of the first layer h_{1}. The total thickness is the sum of a geometric progression:

The interface (the "Boundary layers" panel) shows this total thickness. The layers are built by Netgen on the selected faces; in the analysis mesh, prisms are split into tetrahedra. Layers are set in the desktop application interface; the same parameters are available from the built-in JS console.

3.6Mesh quality metricsbeta

For tetrahedra, the application computes three groups of quantities.

Element shape ("Mesh quality" panel). The edge ratio \mathrm{ER} = \ell_{\max}/\ell_{\min} over the six edges of the tetrahedron. The histogram is built over q = 1/\mathrm{ER} in ten bins from 0 to 1; hovering over a bar highlights the elements on the canvas, and a click leaves only them.

Normalized Jacobian (SICN). In the mesh passport the field is called SICN (signed inverse condition number), and the tolerance thresholds are compared against it. In fact, for a tetrahedron the mesher computes the normalized determinant of the edges emanating from its first vertex:

A regular tetrahedron gives q_{J} = 1/\sqrt{2} \approx 0.71; a right-angle corner of mutually orthogonal edges gives 1. This is not the full definition of SICN from the literature (by which a regular tetrahedron gives 1): the quantity is computed from one corner and is suitable for comparing meshes within WebCAE, but not for comparison with numbers from other packages. An element is considered inverted if the sign of its determinant differs from the one prevailing in the mesh; a zero determinant means degenerate.

Tolerance. The acceptance threshold q_{J,\min} \ge 0.1 is applied as follows:

• accepted: there are no inverted or degenerate elements, \min q_{J} \ge 0.1;
• warning: the Jacobian is positive everywhere, but \min q_{J} < 0.1 (slivers): the analysis is allowed, and a hint about possible unreliability is shown;
• rejected: there are inverted or degenerate elements, or an empty mesh: the analysis is blocked;
• disconnected mesh pieces (the assembly fell apart) block the analysis unconditionally: the stiffness matrix is singular.

For hexahedra the threshold is the same, and the measure is the minimum scaled Jacobian over all eight corners of the element (Verdict definition); non-positive values are not allowed. For tet10, validity is determined by the Bernstein coefficients (see element order).

A linear hex8 locks in bending: with a small number of elements through the thickness of a layer the deflection comes out underestimated, and the interface shows such a warning next to the mesh.

3.7Volume hexahedral meshbeta

Structured hexahedra give a regular mesh with a predictable number of elements through the thickness, but they cannot be built automatically for every body. The application tries routes in decreasing order of simplicity and, on failure, honestly falls back to tetrahedra.

1. Structured primitives. Box, cylinder (O-grid topology: a central square block and four polar blocks), sphere, truncated cone, hemispheres: transfinite interpolation, mesh hex nx=… ny=… nz=….
2. Sweeping. The body is recognized as prismatic ("2.5D"): a cross-section with holes is swept along a direction; bodies of revolution are allowed.
3. Decomposition. The body is divided into prisms, the cross-section of each being a quadrilateral block; interfaces are matched, then a structured hex mesh is built in each block and the nodes are stitched.
4. Failure. If no route has produced an accepted mesh, a tetrahedral mesh is built, and the reason is reported as a code: there is no silent substitution.

A built hex mesh passes acceptance before use: the mesh volume is compared with the exact volume of the body, free faces must lie on the B-Rep boundary, the mesh must be one connected component with no dangling nodes and no faces shared by more than two elements, and the minimum scaled Jacobian over the eight corners must not be below the quality threshold (0.1). A mesh that fails acceptance does not go into the analysis.

In the desktop application, external generators are additionally available if installed: HOHQMesh (planar parts of constant thickness) and the HexaBlock block method from SALOME. In "Auto" mode they are chained as "structured sweep → HexaBlock → tetrahedra".

units mm
box dx=60 dy=30 dz=10
mesh hex size=5

3.8Local hex zones and hybrid meshesbeta

When a hexahedral layer is needed in one place (a hole edge, a shaft neck, a contact pad), the rest of the body remains tetrahedral. A zone is defined relative to a selected face or volume; the other parameters are the layers, the progression, and the coupling method.

Zone types:

• refine tube surface number N thickness=t: a layer along a cylindrical face (hole edge, shaft neck);
• refine inflate surface number N thickness=t: a layer along a planar or spherical face (contact pad);
• refine box at=… dx= dy= dz=: a box inside the body, filled with structured hexahedra (clipped to the boundary with clip=on);
• refine sweep surface number A to surface number B: a prism between two parallel planar faces of one body.

Layer parameters: layers (1 to 12), bias, the progression ratio of thicknesses from the wall (1 to 3), elements=wedge|hex, the layer element type, conform=weld|mpc, the coupling method.

Interface consistency. The bilinear trace of a quadrilateral hex face does not coincide with the piecewise-linear trace of two tetrahedra over the same face, so simply gluing hex8 to tet4 does not pass the patch test. By default the layer is built from prisms (wedge6): the prism over each interface triangle is consistent with the tetrahedra by construction. With elements=hex the transition is made by pyramids (pyr5): the interface hex8 cell is divided toward its own centre into pyramids, and on the faces on the tetrahedra side the pyramids are replaced by two tetrahedra with the same diagonal. The method conform=mpc builds an independent layer and ties it to the tetrahedra with nodal constraint equations.

Measured agreement of mixed meshes in the patch test (maximum relative displacement error): prisms glued to tetrahedra, about 10^{-14} (web platform); direct gluing of hex8 to two tet4, 3.1 \times 10^{-3} (native); transition hex8 → pyr5 → tet4, 9 \times 10^{-15} (native). Numbers of different platforms are not mixed.

For the Kirsch problem (plate with a hole), a hybrid mesh with prisms at the edge gives an error of the maximum \sigma_{\theta} within 0.1 to 0.9 % and a root-mean-square error along the hole arc of 0.7 to 1.2 % at 1635 to 4473 degrees of freedom. A purely tetrahedral mesh at 1374 to 12828 degrees of freedom fluctuates within −17 to +4 % on the peak, and the root-mean-square error along the arc does not drop below 5.3 %.

units mm
refine tube surface number 12 thickness=3 layers=3 bias=1.2
refine sweep surface number 1 to surface number 2 layers=4

3.9Meshes of assemblies and contact pairsbeta

An assembly has two tasks, and they are opposite. Bodies connected rigidly (welding, adhesive, interference fit) must share nodes at the interface, otherwise the stiffness matrix falls apart into disconnected blocks. Bodies that may separate or slide (contact) are meshed independently, and the interface is handled by the solver's contact algorithm.

Gluing (bonding=1). Touching bodies are meshed as one conforming model: common faces receive common nodes by construction. Modes: glue all, glue selected groups, do not glue (bonding=0). For contacts, gaps within tolerance, and overlaps, a dedicated Boolean operation is performed; the pair selection tolerance (bondtol=) and the tolerance of the operation itself (fuzzy=) are set separately. The joint mode glue=auto|contact|full: fast gluing by contact, a full intersection search, or automatic selection. Overlap is handled by the parameter interference=cut|split|ignore: give the volume to the smaller body (default), split the intersection (the volume is exact, but the interference fragment may be a thin film), or leave as is (material in the interference is counted twice).

Contact pairs. Bodies not glued to each other are each meshed separately; nodes on contact surfaces need not coincide. Before the analysis, interface readiness is checked: contact regions may have collapsed into one region, no pairs may have been found, the gap may be too large, or a glued pair may be non-conforming. Each such outcome is reported explicitly.

Parallel meshing. With separate meshing, bodies are built independently and in parallel, progress is shown per body, and the result is assembled into a single mesh. If a body could not be meshed in native mode, it may be completed by the built-in web mesher; the application warns that such meshes may differ slightly. With gluing, the assembly is meshed as a single model.

Welding of a ready mesh. The command weld [tol=] and the "Weld" panel merge coincident nodes of an already built mesh within a tolerance (0 means automatic, a small fraction of the bounding box diagonal). This is welding of nodes, not of geometry: useful for shell patches that are meshed separately on purpose.

units mm
mesh tet size=8 bonding=1 glue=auto interference=cut

3.10Shell and beam meshesbeta

Shells. A sheet (B-Rep shell) is meshed with a surface mesh: mesh shell size=h element=tri|quad. Triangles are built on each face in parametric space: a constrained Delaunay triangulation along the face contour with projection onto the surface. Quadrilaterals are obtained from triangles by the "1 → 3" scheme: each triangle is divided into three quadrilaterals through the centroid and the edge midpoints. The edge midpoint is common to neighbours, so the mesh is conforming at face seams, and the elements are entirely quadrilateral for any number of triangles. The quadrilateral edge size is from 0.29 to 0.5 of the edge length of the source triangle, so the source triangulation is built noticeably coarser.

For quadrilateral faces with straight sides (boxes, plates), the parameter pattern=structured (together with element=quad) builds a structured mesh by Coons transfinite blending with consistent edge divisions. If the shape contains an unsuitable face (a hole, an arc), the command refuses with the face number instead of mixing incompatible meshes.

Beams. A beam mesh is built along selected lines (B-Rep edges): the number of elements or the size is set; common nodes with the adjacent shell arise when a beam is seated on an edge (Frames and beams). Elements for the assembly are compiled at the analysis step.

units mm
mesh shell size=10 element=quad
assign beam edges=3,7 section=rect w=10 h=30 flush=outer

3.11Zienkiewicz–Zhu error indicator and adaptive refinementbeta

The maximum stress does not distinguish a real concentration (hole, fillet) from an artifact of a coarse mesh. For this, the Zienkiewicz–Zhu error indicator is used. The finite element stress field \boldsymbol{\sigma}_{h} is discontinuous between elements. A continuous field \boldsymbol{\sigma}^{*} is recovered by local smoothing over the patch of elements around a node, and the error on an element is estimated by the energy norm of the residual:

where \mathbf{C} is the elasticity matrix. The relative model error is \eta_{\text{rel}} = \eta \,/\, \sqrt{\lVert\boldsymbol{\sigma}^{*}\rVert_{E}^{2}+\eta^{2}}. For tet10, recovery is performed at the superconvergent Gauss points by the SPR method (Zienkiewicz O.C., Zhu J.Z., The superconvergent patch recovery and a posteriori error estimates, Int. J. Numer. Meth. Eng., 1992); the integral over an element is taken exactly, without quadrature error. For tet4 and mixed meshes, nodal averaging is used, which has lower accuracy.

The command hotspots finds hotspots and their sizes from the indicator and the stress field, and the command refine adaptive closes the loop: analysis → hotspots → refinement zones in the hotspots → mesh rebuild → analysis. The loop stops when one of the conditions is met: the grid convergence index (GCI, see the convergence study) is not above the target, a singularity is diagnosed, the relative error \eta_{\text{rel}} is not above the given value, or the number of passes is exhausted.

Parameters of refine adaptive: passes, the number of passes from 1 to 8 (default 3); gci, the target GCI band in percent (default 5); tol, the target relative ZZ error in percent; refine, the factor by which the element size in a hotspot is reduced per pass (default 1.8); minsize, the lower size limit; grading, max_hotspots, at=, radius=, the transition smoothness, the number of hotspots, and an explicit measurement region. After execution, the measurements adaptive_passes, adaptive_verdict, adaptive_gci_percent, adaptive_richardson_peak, adaptive_zz_error_percent are available.

units mm
solve
refine adaptive passes=3 gci=5 tol=2 refine=1.6 minsize=0.5
probe adaptive_verdict

Submodelling (submodel at=… box=… size=… order=2) continues the same idea: a subregion is cut out around a concentrator and meshed finely (including with quadratic tet10), and the displacements of the global coarse model are transferred to the cutting boundaries by the Saint-Venant principle (the boundary displacement interpolation method). The displacements of the nodes of the cutting boundaries, including midside nodes, are prescribed by a nonhomogeneous Dirichlet condition through constraint equations.

3.12Mesh convergence studybeta

Mesh convergence is not solver convergence. Here a series of independent analyses on an increasingly fine mesh answers the question: what does the chosen quantity tend to (most often the stress peak at a concentrator)? If it does not settle, any number in the results table is a mesh artifact.

The series metric is the maximum of the field inside a sphere around the hotspot of the previous iteration, not the global maximum: on different meshes the global maximum may sit in different places of the model (at the hole on one, at the clamp on another), and then different points are being compared. The characteristic element size is taken from the mesh itself in the sphere as h = \overline{\Delta V_{e}}^{\,1/3} (Celik I.B. et al., Procedure for estimation and reporting of uncertainty due to discretization in CFD applications, J. Fluids Eng., 2008), not from the requested size: the mesher delivers not the size that was asked for.

For the last three meshes with sizes h_{1} > h_{2} > h_{3} and values f_{1}, f_{2}, f_{3}:

Here p is the observed order of convergence, f_{\text{ext}} is the Richardson extrapolation to an infinitely fine mesh, GCI is the grid convergence index after Roache (ASME V&V 20), and the uncertainty band equals \mathrm{GCI}\cdot\lvert f_{3}\rvert. The order p is computed as signed (in Celik the absolute value of the whole expression is taken): the sign of p is itself the diagnosis of divergence. For unequal r_{21}, r_{32} the geometric mean r_{\text{eff}} is used, not the iterative formula of Celik.

Verdict for the series:

• converged: GCI is not above the target (default 5 %), the result is recorded as f_{\text{ext}} \pm band;
• keep refining: the series converges but the band is still wide; the same for alternating-sign increments (the series oscillates around the limit, extrapolation is unreliable);
• probable singularity: p \le 0, the increments do not decrease: the peak grows with refinement. This is typical of a sharp re-entrant corner (stress grows as h^{-\lambda}, for a 270° corner \lambda \approx 0.456), a point force, or a rigid clamp along an edge. The mesh will not help; add a fillet or change the formulation;
• insufficient data: fewer than three meshes; invalid: the size does not decrease, the peak is zero or non-numeric, or the mesh did not change between iterations.

There is no silent stop at the iteration limit: the verdict is always one of those listed.

This section describes the constitutive relations that WebCAE computes: what the law contains, which parameters the user sets, in which analysis types the law works, and how it has been verified. Some models are implemented in the computational core but are not brought to the point of analysing a body; such models are marked as experimental, and their subsections say so directly.

Conventions. Tension is positive. Stresses and strains in the laws are written in Voigt order [xx,\,yy,\,zz,\,xy,\,yz,\,zx]; shear strains are engineering (\gamma_{ij} = 2\varepsilon_{ij}), shear stresses are tensorial (\tau_{ij} = \sigma_{ij}). At the solver boundary everything is in SI (Pa, K, kg/m³, s); in the interface, input fields accept MPa, GPa, psi, °C and convert to SI on entry. In the DSL, stresses and moduli are given in pascals, plasticity temperatures in kelvins.

4.1Linear elastic isotropic materialstable

The basic law of all linear analyses and the elastic part of all nonlinear ones. It is sufficient as long as the stresses do not reach yield and the strains are small compared to unity.

Hooke's law in Voigt form:

The block \mu on the diagonal corresponds to engineering shears: \tau_{xy} = \mu\gamma_{xy}. The admissible range is -1 < \nu < 0.5. As \nu \to 0.5 the matrix \mathbf{C} degenerates in its volumetric part, and linear tetrahedra lock; for nearly incompressible materials use second-order elements (see elements).

The user sets Young's modulus E, Poisson's ratio \nu, and density \rho (see the next subsection). The result is displacements, strains, stresses, and derived quantities (von Mises equivalent, principal), which are covered in the section on linear analyses.

One material is assigned to a whole body. For an assembly of several bodies each body has its own material, including different E, \nu, \rho at the interface between bodies.

material E=200e9 nu=0.3 density=7850

4.2Density, thermal expansion and temperature dependencesstable

Density \rho (kg/m³) is needed by modal analysis, dynamics, gravity, and acceleration; without it a request involving mass is refused rather than using a default. For a shell and a laminate, the mass is taken from the density of the material or of the plies.

Thermal expansion is specified by the coefficient of linear expansion \alpha (1/K). The thermal strain is isotropic, and the stress is determined by the mechanical part:

Thermal loads and coupled problems are described in the section on heat and multiphysics. What matters here is that \alpha is constant: there is no dependence \alpha(T). Temperature dependences in material models exist in three places: the thermal conductivity table \kappa(T) (kappa_table="T:k,…", temperature in °C, linear between nodes, constant at the end value outside the table), the temperature factor of the yield strength in the Johnson–Cook model, and the thermal expansion of laminate plies. Young's modulus does not depend on temperature.

4.3Material librarystable

The built-in library is a single table of presets; it is used by the material card in the interface, by the DSL, and by the programming interface. A preset carries E, \nu, \rho, thermal conductivity, heat capacity, \alpha, and the yield strength for the safety factor.

• Steels: general structural (200 GPa, 0.3, 7850 kg/m³), 20, 40Kh, 09G2S, stainless 304.
• Aluminium: general (70 GPa, 0.33), D16T, 6061.
• Titanium: general (116 GPa, 0.32) and VT6; copper (120 GPa, 0.34); plastics ABS and PA6.
• Rubber as a neo-Hookean material: \mu = 3.36 MPa, K = 166.7 MPa.
• Plastic presets: steel S355 (J2, \sigma_y = 355 MPa), aluminium 6061-T6 (J2), copper OFHC and steel 4340 (Johnson–Cook).
• Media for CFD: air, water, engine oil, kerosene.

The values are reference and generic: take heat capacity and limits for a specific grade from your own certificate. The yield strength in an ordinary preset is needed only for the safety factor and does not enable plasticity; plasticity appears in presets with an explicit law and in materials where you have specified it.

The user adds their own materials to the project library, edits them, and assigns them to bodies. For laminated shells there is a separate library of orthotropic ply materials (a built-in unidirectional plastic 140/10 GPa).

In the DSL a preset supplies the base, and explicit parameters override it:

material steel E=210e9

4.43D orthotropic material and orientationbeta

For materials whose properties differ along three axes: unidirectional and woven composites, wood, rolled products. The law is specified by nine independent constants E_1, E_2, E_3, \nu_{12}, \nu_{23}, \nu_{13}, G_{12}, G_{23}, G_{13} in material axes and by the orientation of the axes in each element.

The compliance matrix in Voigt order, engineering shears:

Poisson's ratios are defined as \nu_{ij} = -\varepsilon_j/\varepsilon_i under uniaxial \sigma_i; symmetry gives \nu_{ji} = \nu_{ij}E_j/E_i. Positive definiteness is checked before the analysis: E_i > 0, G_{ij} > 0, 1 - \nu_{ij}\nu_{ji} > 0 for the three pairs, and

A set that fails is refused with the name of the violated condition.

Orientation. At the centroid of each element a right-handed triad (\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3) is built, and the rows of the matrix \mathbf{Q} are its vectors in global axes. The element stiffness matrix is obtained by Bond rotation, and the stresses are returned in material axes:

Triad definition modes:

• surface: the axis \mathbf{e}_3 equals the outward normal of the triangle of the reference face nearest to the centroid; \mathbf{e}_1 is the projection of the given direction axis onto the tangent plane, \mathbf{e}_2 = \mathbf{e}_3\times\mathbf{e}_1. Suitable for laminated parts along a curved surface.
• cartesian: a single triad: \mathbf{e}_3 along normal, \mathbf{e}_1 is the projection of axis.
• cylindrical: \mathbf{e}_3 radially from the cylinder axis through origin.
• spherical: \mathbf{e}_3 radially from the point origin.

The angle angle rotates \mathbf{e}_1, \mathbf{e}_2 about \mathbf{e}_3: \mathbf{e}_1' = \cos\theta\,\mathbf{e}_1 + \sin\theta\,\mathbf{e}_2. If the projection of the given axis is degenerate, a fallback direction is taken; the number of such elements is recorded in the result. The triads can be displayed on the model.

The solution is separate from the general isotropic assembly: for TET4 and TET10 tetrahedra the stiffness is assembled from rotated matrices, the system is solved by the conjugate gradient method with a nodal block preconditioner, and constraints (CE, rigid regions) are introduced by a penalty. Results: displacements, stresses in global axes and in material axes (localStressXX…XZ, probe ortho_local_s11_max and analogues), the triad field.

material E1=140e9 E2=10e9 E3=10e9 nu12=0.3 nu23=0.4 nu13=0.3 G12=5e9 G23=3.5e9 G13=5e9 density=1600
orientation surface number 3 axis=1,0,0 angle=30

In the interface, the orthotropic law is selected in the material card ("Orthotropic 3D", nine constants, reference surface, direction, angle); the cylindrical and spherical modes are available from the DSL.

4.5Laminated shells: laminate theorystable

A laminated shell replaces a set of plies with a single through-thickness averaged law, the stiffness package. It is suitable for composite panels and sandwiches in which the thickness is small compared to the dimensions. Plies are specified from bottom to top (from -z_e to +z_e), each with its own orthotropic material, thickness, and angle.

A ply in material axes, plane stress state:

Rotation by an angle \theta from the element axis x_e counterclockwise about +z_e (c = \cos\theta, s = \sin\theta):

The terms \bar Q_{16}, \bar Q_{26} are proportional to \sin 2\theta and change sign together with the angle; at 0^\circ and 90^\circ they are zero, so checks at only these angles do not catch sign errors and shear terms. The product's check suite runs \pm15^\circ, \pm30^\circ, \pm45^\circ, \pm60^\circ, \pm75^\circ and unsymmetric stacks.

For ply k between z_{k-1} and z_k (measured from the reference surface), the stiffness package per unit width:

Here \boldsymbol\varepsilon_0 = [\varepsilon_{xx}, \varepsilon_{yy}, \gamma_{xy}], \boldsymbol\kappa = [\kappa_{xx}, \kappa_{yy}, \kappa_{xy}] with engineering twist curvature, \varepsilon(z) = \varepsilon_0 + z\kappa. The coupling matrix \mathbf{B} is nonzero for unsymmetric stacks; tension then causes bending. The transverse shear of a ply \bar{\mathbf{Q}}_s is built from G_{13}, G_{23} by the same rotation:

Shear correction for non-uniform shear: by default k_s = 5/6 for all plies. The variant shear=energetic takes k_s^{(k)} = \min\bigl(1,\ \tfrac56\,G_{\max}/G_k\bigr) per ply (a soft sandwich core brings k_s closer to 1); this is an engineering heuristic, not a derivation from equilibrium conditions.

Shifting the reference surface by z_0 from the midplane (offset=bottom|middle|top|number) recomputes the package without reassembly: \mathbf{B}' = \mathbf{B} + z_0\mathbf{A}, \mathbf{D}' = \mathbf{D} + 2z_0\mathbf{B} + z_0^2\mathbf{A}.

The element is the four-node MITC4 (Dvorkin, Bathe 1984; Bathe, Dvorkin 1985): shear from covariant strains at tying points, membrane and bending on 2\times 2 Gauss points, and rotation about the normal stabilized by a Hughes–Brezzi penalty with membrane shear stiffness A_{66}. The package (\mathbf{A},\mathbf{B},\mathbf{D},\mathbf{A}_s) enters the element directly, so an isotropic shell is a special case of the package (\mathbf{A} = Et/(1-\nu^2)\,\mathbf{C}, \mathbf{D} = Et^3/(12(1-\nu^2))\,\mathbf{C}, \mathbf{B}=0).

The axis x_e from which ply angles are measured is set by the projection of the reference axis onto the shell (axis=X|Y|Z), by a coordinate system (cs=), and by a rotation phi=. If the reference axis is nearly normal to the element, a fallback axis is taken, and the number of such elements is reported.

The thermal load of a shell enters through the forces \mathbf{N}_T = \int \bar{\mathbf{Q}}\,\bar{\boldsymbol\alpha}\,\Delta T\,dz and moments \mathbf{M}_T; ply stresses are computed from the mechanical strain \varepsilon - \bar\alpha\Delta T(z). A ply needs \alpha_1, \alpha_2 (along and across the fibre).

What the user specifies: the orthotropic ply material (E1 E2 nu12 G12 G13 G23, density, strengths, \alpha_1, \alpha_2) and the layup (plies="name:thickness:angle,…", symmetric, repeat, offset, shear). Ply thicknesses are given in script length units, moduli in pascals. In the interface there is a ply layup editor in shell sections. The formulation is the same on both platforms: on the desktop (reference) shells are computed by the native path, in the browser by the web build.

ortho name=ud E1=140e9 E2=10e9 nu12=0.3 G12=5e9 G13=5e9 G23=3.5e9
section shell laminate plies="ud:0.25:45,ud:0.25:-45,ud:0.25:0" symmetric=true

4.6Ply failure criteria and interlaminar stressesbeta

A per-ply postprocessor: from the generalized shell strains in an element (or their nodal SPR recovery), ply stresses in material axes (\sigma_1, \sigma_2, \tau_{12}) are computed at the bottom, middle, and top of the ply, then failure indices. The criteria require the ply strengths X_t, X_c, Y_t, Y_c, S in ortho; without any of them the criterion is refused with the field name.

• Maximum stress: FI = \max\bigl(\lvert\sigma_1\rvert/X,\ \lvert\sigma_2\rvert/Y,\ \lvert\tau_{12}\rvert/S\bigr), with X and Y selected by the sign of the stress.
• Tsai–Hill: FI = (\sigma_1/X)^2 - \sigma_1\sigma_2/X^2 + (\sigma_2/Y)^2 + (\tau_{12}/S)^2.
• Tsai–Wu: FI = F_1\sigma_1 + F_2\sigma_2 + F_{11}\sigma_1^2 + F_{22}\sigma_2^2 + F_{66}\tau_{12}^2 + 2F_{12}\sigma_1\sigma_2, F_1 = 1/X_t - 1/X_c, F_{11} = 1/(X_tX_c), F_{12} = -\tfrac12\sqrt{F_{11}F_{22}}.
• Hashin (planar, four modes): fibre rupture in tension (\sigma_1/X_t)^2 + (\tau_{12}/S)^2, fibre buckling (\sigma_1/X_c)^2, matrix cracking in tension (\sigma_2/Y_t)^2 + (\tau_{12}/S)^2, in compression (\sigma_2/2S_{23})^2 + \bigl[(Y_c/2S_{23})^2 - 1\bigr]\sigma_2/Y_c + (\tau_{12}/S)^2; the ply index is the maximum over modes.
• Delamination: FI = (\tau_{13}/S_{13})^2 + (\tau_{23}/S_{23})^2 from interlaminar shear.

The reserve factor RF is the load multiplier at which FI = 1: for a linear criterion 1/FI, for quadratic ones 1/\sqrt{FI}, for Tsai–Wu and Hashin "matrix compression" the positive root of the equation a\lambda^2 + b\lambda = 1. DSL probes: laminate_fi_<criterion>, laminate_rf_<criterion>, laminate_s1, laminate_s2, laminate_t12, laminate_t13, laminate_t23 with the suffix _plyN for ply N (1 is the bottom).

Interlaminar stresses by first-order theory: \tau_{13}, \tau_{23} are constant through the ply thickness. The option interlaminar=equilibrium recovers \tau_{xz}, \tau_{yz} by integrating the equilibrium equations through the thickness from the free bottom face (gradients are taken from neighbouring mesh elements), and \sigma_{zz} as an estimate by the same principle.

4.7J2 plasticitystable

The Huber–von Mises model for metals: elastic behaviour up to the yield strength, then plastic flow with hardening. It is applied in nonlinear analyses of solid bodies (small strains or geometric nonlinearity) and in explicit dynamics. The general scheme of nonlinear iterations is described in the section on material nonlinearity.

The yield surface with back-stress tensor \boldsymbol\alpha:

where \mathbf{s} is the stress deviator. Integration is by radial return (Simo, Hughes, Computational Inelasticity, 1998, ch. 3; Wilkins, 1964). Trial elastic step \mathbf{s}^{\mathrm{tr}} = 2\mu\,\mathrm{dev}(\boldsymbol\varepsilon - \boldsymbol\varepsilon_p); if f^{\mathrm{tr}} \le 0, the step is elastic. Otherwise, with the normal \mathbf{n} = \boldsymbol\xi^{\mathrm{tr}}/\lVert\boldsymbol\xi^{\mathrm{tr}}\rVert:

The scalar equation for the increment of equivalent plastic strain \Delta p is strictly monotonic and has a unique root: for linear hardening \Delta p = f^{\mathrm{tr}}/(3\mu + H), for a curve a safeguarded Newton method with bisection on the bracket [0,\ f^{\mathrm{tr}}/3\mu]. The tangent is algorithmic (consistent with the integrator) and symmetric:

Hardening. The curve \sigma_c(p) is the diagram of uniaxial monotonic loading:

• perfect plasticity: \sigma_c = \sigma_{y0};
• linear by tangent modulus E_t: \sigma_c = \sigma_{y0} + Hp, H = E E_t/(E - E_t);
• multilinear: a table (p_i, \sigma_i), the first point p_0 = 0 sets \sigma_{y0}, p is strictly increasing, \sigma is non-decreasing; beyond the last point the stress is constant (perfect plasticity) rather than continuing along the slope.

The kinematic hardening fraction \beta \in [0, 1] (0 is isotropic, 1 is Prager linear kinematic) mixes the mechanisms: radius R = \sigma_{y0} + (1-\beta)\bigl(\sigma_c(p) - \sigma_{y0}\bigr), centre shift \Delta\boldsymbol\alpha = \sqrt{2/3}\,\beta\,\bigl[\sigma_c(p_{n+1}) - \sigma_c(p_n)\bigr]\,\mathbf{n}. The uniaxial diagram does not depend on \beta, while reverse yielding (the Bauschinger effect) does. Softening is prohibited (H \ge 0).

Kinematics of nonlinear statics: in geometric nonlinearity mode the law accepts the Green–Lagrange strain and returns the second Piola–Kirchhoff stress (plasticity is additive in this pair: small strains, large rotations); the output is the Cauchy stress \sigma = \mathbf{F}\mathbf{S}\mathbf{F}^T/J. In the "material only" mode \varepsilon = \mathrm{sym}\nabla\mathbf{u}. A load history (loadpath="1,0") gives residual strains after unloading and cycles.

What the user specifies: E, \nu (elastic part), \sigma_{y0} (yield=, Pa), hardening (hardening=linear|multilinear, et=, curve="p:σ,p:σ,…"), kinematic=. Result: equivalent plastic strain (nodal field, probe plastic_strain), the fraction of yielded Gauss points, and for perfect plasticity the limit load fraction.

material E=210e9 nu=0.3 plasticity=j2 yield=250e6 hardening=linear et=1e9
nlgeom material loadpath="1,0"
solve

4.8Johnson–Cookstable

A law for metals at high strain rates and heating: yield strength grows with rate and falls with temperature. The main area is explicit dynamics: impact, penetration, acceleration.

The rate factor is switched on only for \dot\varepsilon^* > 1 (below that it equals one: quasistatics serves as the lower bound). The return is the same as for J2, but the equation for \Delta p includes the step rate \dot{\bar\varepsilon}_p = \Delta p/\Delta t and the current temperature factor. Adiabatic heating is switched on by specifying the heat capacity c_p of the law's material: \Delta T = \eta\,q\,\Delta p/(\rho c_p), where \eta is the Taylor–Quinney fraction (0.9 by default); the temperature changes stepwise, and the step's yield is computed at the temperature at the start of the step.

Parameters: A, B (Pa), n, C, m, \dot\varepsilon_0 (epsdot0, 1/s), T_{\mathrm{room}} (troom, K), T_{\mathrm{melt}} (tmelt, K), taylor, cp. Kinematic hardening is not combined with Johnson–Cook (refused).

material E=120e9 nu=0.34 density=8900 plasticity=jc A=90e6 B=292e6 n=0.31 C=0.025 m=1.09 tmelt=1356 cp=385

4.9Material failure criteriabeta

The criteria work inside the J2 or Johnson–Cook plastic law and are needed for element erosion in explicit dynamics; statics and modal analysis ignore them. Triggering of any criterion is failure of the Gauss point. What the product does with the element is described in the section on erosion.

• \varepsilon_p \ge \varepsilon_{\mathrm{fail}}: limit equivalent plastic strain (ep:).
• Accumulated Johnson–Cook damage D = \sum \Delta\varepsilon_p/\varepsilon_f \ge 1 (jc:d1/d2/d3/d4/d5):

The sign of d_3 is as in Johnson, Cook (Eng. Fract. Mech. 21(1), 1985, 31–48): for metals d_3 < 0; \eta is limited to the interval [-3, 3].

• \sigma_1 \ge \sigma_{\max}: largest principal stress (smax:, Pa).
• \sigma_m \ge p_{\mathrm{spall}}: hydrostatic tension, spall (spall:, Pa).

Damage D does not reduce stiffness: it only triggers failure.

material steel plasticity=j2 yield=355e6 fail="ep:0.25,smax:1.2e9"

4.10Neo-Hookean materialstable

A hyperelastic material for rubber and elastomers at finite strains. It works in nonlinear statics of solids in the Total Lagrangian formulation (geometric nonlinearity).

The second Piola–Kirchhoff stress \mathbf{S} = 2\,\partial W/\partial\mathbf{C} and the tangent \mathbb{D} = 4\,\partial^2 W/\partial\mathbf{C}\partial\mathbf{C} are written out analytically (comparison with finite differences: up to 10^{-8}). The output stress is Cauchy, \boldsymbol\sigma = \mathbf{F}\mathbf{S}\mathbf{F}^T/J. The constants \mu, \kappa are taken from E, \nu of the material by the formulas of the subsection on isotropic elasticity (in the rubber card, \mu and K are specified). Elements: HEX8, TET4, WEDGE6, PYR5, HEX20, TET10, HEX27.

material E=1e7 nu=0.45 model="neo-hookean" density=1100
nlgeom on

4.11Mooney–Rivlin, Yeoh, Arruda–Boyce, Ogdenexperimental

Potentials for rubber stretch curves that the neo-Hookean law does not describe. The core contains all four with an analytical tangent; the full chain "material card — body analysis" is absent in the current version, so the section is experimental.

A volumetric part \tfrac\kappa2(J-1)^2 is added to each potential (for Arruda–Boyce \tfrac\kappa2\bigl[\tfrac12(J^2-1)-\ln J\bigr]); \kappa = 2/D_1 if the compressibility parameter D_1 is given. Small-strain moduli: Mooney–Rivlin \mu_0 = 2(C_{10}+C_{01}), Yeoh \mu_0 = 2C_{10}, Arruda–Boyce \mu_0 = \mu, Ogden \mu_0 = \sum\mu_p (Abaqus convention; the classical c_p = \mu_p/\alpha_p is selected by the flag ogden_abaqus=false). For Ogden, the principal stretches are obtained by spectral decomposition of \mathbf{C}.

4.12Finite elastoplastic strains F = Fe·Fpexperimental

Multiplicative plasticity is needed where the additive J2 scheme is wrong: strains of tens of percent. A prototype of the law has been implemented, not connected to the interface or the DSL.

The gradient is decomposed as \mathbf{F} = \mathbf{F}_e\mathbf{F}_p, \det\mathbf{F}_p = 1. Elasticity is Hencky in the elastic configuration: from the spectral decomposition of \mathbf{C}_e = \mathbf{F}_e^T\mathbf{F}_e the logarithmic strains h_a = \tfrac12\ln\lambda_a are obtained, and the Kirchhoff stress is \boldsymbol\tau = K\,\mathrm{tr}(\mathbf{h})\,\mathbf{1} + 2\mu\,\mathrm{dev}(\mathbf{h}). The flow is associated von Mises with linear hardening \sigma_y = \sigma_{y0} + Hp; the increment \Delta p = f^{\mathrm{tr}}/(3\mu + H) is found in closed form in the space of logarithmic strains, and the plastic part is updated by the exponential \mathbf{F}_p^{n+1} = \exp(\Delta p\,\mathbf{N})\,\mathbf{F}_p^{n}, \mathbf{N} = \tfrac32\,\mathrm{dev}(\boldsymbol\tau)/q. The scheme is of the same type as in Eterovic, Bathe (1990), and Simo (1992). For passing to the Newton iteration, the Cauchy stress \boldsymbol\sigma = \boldsymbol\tau/J is converted to the second Piola–Kirchhoff stress, and the tangent \partial\mathbf{S}/\partial\mathbf{E} is taken by fourth-order central differences; on HEX8, HEX20, HEX27, and TET10 elements an \bar F correction against volumetric locking is applied.

4.13Norton–Bailey creepexperimental

Power-law creep of metals at high temperatures and long-duration loads: strain grows at constant stress, and stress relaxes at constant strain.

q is the von Mises equivalent stress; the core has time and strain hardening and the integral form \bar\varepsilon_{cr} = A q^n t^m (NAFEMS Test 10B). The step is implicit Euler with radial return and a scalar equation for q_{n+1}, with an analytical tangent. Parameters in the DSL: creep_a, creep_n, creep_m.

4.14Perzyna viscoplasticityexperimental

Plastic flow with rate dependence: the stress may exceed the static limit, and the excess causes viscoplastic flow.

Implicit Euler, the scalar equation \Delta p - \Delta t\bigl\langle (f^{\mathrm{tr}} - (3\mu + H)\Delta p)/\eta\bigr\rangle^m = 0 on a bracket, with an analytical tangent. As \Delta t\to\infty the law tends to static J2, and for small \eta to rigid plasticity.

4.15Viscoelasticity: Prony series and rheological chainsexperimental

Linear viscoelasticity of polymers: relaxation of the modulus in time. The shear modulus in the form of a Prony series, the bulk modulus analogously:

The internal variables of the branches are updated by an exact exponential recurrence (Simo, Hughes, 1998, §10.2): \mathbf{h}_i^{n+1} = e^{-\Delta t/\tau_i}\mathbf{h}_i^n + 2G_0 g_i\,\beta_i\,\Delta\mathbf{e}, \beta_i = (1 - e^{-\Delta t/\tau_i})\,\tau_i/\Delta t; up to six branches, and if the k_i are absent, k_i = g_i is taken. Separately there are one-dimensional Maxwell, Kelvin–Voigt, Zener, and Burgers chains with an exact exponential step, and a thermomechanical damper with an energy balance.

4.16Drucker–Prager for soils and rocksexperimental

Pressure-dependent plasticity for soil, concrete, and rock: shear strength grows with confinement, and flow is accompanied by dilatancy.

The cone is fitted to the triaxial compression meridian of the Mohr–Coulomb surface, so the limit axial stress at a confining pressure P_c equals

The potential is non-associated: g = q + \beta p, \beta = 6\sin\psi/(3-\sin\psi), dilatancy angle 0 \le \psi \le \varphi. Return to the smooth wall of the cone; when the projection falls beyond the apex, return to the apex (de Souza Neto, Perić, Owen, 2008, §8.3). The tangent is nonsymmetric for \psi \ne \varphi. Parameters in the DSL: cohesion (Pa), friction_angle, dilatancy_angle (degrees).

4.17GTN porous plasticityexperimental

The core contains only the Gurson–Tvergaard–Needleman yield function \Phi = (q/\sigma_y)^2 + 2q_1 f\cosh\bigl(3q_2 p/2\sigma_y\bigr) - (1 + q_3 f^2) = 0 and the curve of limit equivalent stress against pressure. There is no porosity evolution, stress integrator, or connection to analysis.

4.18Strength criteria on the elastic solutionbeta

A safety factor postprocessor on top of nodal stress fields: changing the criterion or the limit does not require a new analysis. It is not part of the material law and not plasticity. For each node the utilization u is computed, and the safety factor n = 1/u (capped from above for zero stress).

A boundary condition or load in WebCAE is a record. A record has four parts: a kind (support, load, rigid link, constraint equation), an element (displacement, force, pressure, acceleration, etc.), a component with a magnitude, and a list of geometry objects to which the record is applied: vertices, edges, faces or a body. Objects are addressed by the numbers shown in the model tree (surface number 3, edge number 7, vertex number 2), not by a direction in the bounding box. The form in the tree and the script command create the same record and share the same logic; the command grammar is described on the page /docs/dsl.

The record compiler gathers all records into a single problem statement and passes it to the solver. Records are additive: a pressure and a force on one face give a single total traction, and supports on different components of one face accumulate. Values are entered in the project unit system (the unit is shown next to the field), geometry is stored in mm, and at the solver boundary everything is converted to SI: metres, newtons, pascals, kilograms. The conversion is performed in one place; there is no repeated conversion further down the pipeline.

Thermal actions (prescribed temperature, heat flux, convection, radiation, volumetric heat generation) and the temperature change for thermoelasticity are described in the thermal analysis section: Heat transfer and multiphysics. Contact and tied bodies are covered in the section Contact and connections. This section collects the mechanical actions and constraints.

5.1Supports and prescribed displacementsstable

A support sets the displacement of the nodes of the selected geometry to zero along one, two or three axes. A prescribed displacement gives these nodes a nonzero value d. Both cases are Dirichlet conditions u_i = d for a set of degrees of freedom S.

The conditions are enforced exactly, without a penalty. The degrees of freedom in S are eliminated from the system: the solution is written as \mathbf{u} = \mathbf{T}\,\mathbf{u}_{\mathrm{f}} + \mathbf{u}_0, where \mathbf{u}_{\mathrm{f}} are the free degrees of freedom and \mathbf{u}_0 carries the nonzero values d, and the reduced system is solved

Lagrange multipliers are not introduced. The same elimination mechanism serves rigid links and MPC equations; details are in the section Solvers.

What the user specifies: the target (vertex, edge, face or body), the component all, x, y or z, and the magnitude. A zero magnitude is a support, a nonzero one is a prescribed displacement in the project length unit. The component all with a nonzero magnitude sets the same value along all three axes. Result: the displacement field and the support reactions (the reactions command, section Results).

Rules worth knowing:

• supports in x, y and z on one face accumulate; three axes make a full fixity;
• one dependent constraint is allowed per degree of freedom. Therefore a support (d=0) takes precedence over a prescribed displacement on the same face and the same axis; of two prescribed displacements on one face and axis, the first in order is applied, and a difference in values produces a warning;
• a prescribed displacement is allowed only on a face. A record on a vertex, edge or body with a nonzero value does not vanish silently: the analysis reports that a face is required for it;
• a model with a prescribed displacement is considered constrained: the "no support" warning is not issued;
• a support on a single vertex of a solid produces a warning: the reaction is collected at one node, the stress there grows as the mesh is refined, and the stress maximum becomes a property of the mesh rather than of the model. Three vertices that suppress the six rigid-body modes are accepted by the formulation;
• for solids, only translational degrees of freedom are constrained. A "rotation" record on a solid face is not accepted by the solver; supports of beams and shells are specified by the frame commands (frame support).

fix surface number 1
fix surface number 2 component=z value=0.5
fix edge number 4

5.2Normal constraint: roller and symmetry planestable

A normal constraint prevents the nodes of a face from moving along the surface normal and leaves them free in the tangent plane. With a zero magnitude this is a roller or a symmetry plane. With a nonzero magnitude it is a prescribed displacement along the outward normal (inflation, compression).

For a node a of the face the equation is written

where \mathbf{n}_a is the nodal normal. It is obtained by averaging the normals of the adjacent boundary triangles with a weight equal to their interior angle at the node (Thürmer G., Wüthrich C., Computing vertex normals from polygonal facets, Journal of Graphics Tools, 1998). Area weighting was tested and rejected: on a flat quadrilateral cut by a diagonal it biases the normal, whereas angle weighting does not depend on the cut.

If a node lies at the junction of two faces with normal constraints, several equations on three degrees of freedom appear at it. They are reduced by full Gauss–Jordan elimination with pivoting, and the normals are not averaged. Two orthogonal symmetry planes give two independent rows, and the node remains free to slide along the line of their intersection. Nearly parallel normals (a smooth surface cut into two faces) give a redundant row, which is discarded. Contradictory constraints produce a warning, and the row is not applied. At a corner of three independent planes the node is fixed.

There is no special "symmetry" type. A symmetry plane is a normal constraint with d=0, and on a flat face perpendicular to an axis it is equivalent to fixing one component. Loads on the symmetric part of a model are specified as a fraction of the total: if a quarter is cut out, a quarter of the force is applied. Antisymmetry (the conditions u_t = 0 and u_n free) is not available in the product.

fix surface number 2 component=normal value=0
fix surface number 3 component=normal value=0.05

5.3Concentrated forces, line forces and body forcesstable

A force is specified by a vector (N) and applied to a vertex, edge, face or body. What exactly the solver does with the resultant depends on the kind of target.

• Vertex. A concentrated force goes to a single mesh node. The vertex must resolve to exactly one node; otherwise the analysis refuses. Under a point force the stress is singular, values near the point depend on the mesh, and the result should be judged from sections far from the point (Saint-Venant's principle).
• Edge. The resultant is divided equally among the nodes lying on the edge. The shares are the same for all nodes of the edge; the lengths of the segments belonging to each node are not taken into account. Therefore on an edge with a nonuniform mesh the load is not a uniform line load but "equal per node".
• Face. The resultant \mathbf{F} is divided by the face area in m²: this gives a uniform traction \mathbf{t}=\mathbf{F}/A, which is integrated consistently (next subsection). Therefore the sum of the nodal forces equals \mathbf{F}.
• Body. The resultant is distributed over the volume: \mathbf{f}_a = (\mathbf{F}/V)\int_{\Omega} N_a\,dV, where V is the volume of the body. For an assembly the body is addressed; without a body specified, the force acts on the whole model.

A concentrated moment on a solid cannot be specified: a node of a volume element has no rotational stiffness. A moment applied to a region is specified by a remote load. For beams and shells, nodal forces and moments are available through the command frame load node=… force=… moment=… (see Frame, beams and shell patches).

load vertex number 3 force=0,0,-100
load edge number 5 force=0,-2000,0
load surface number 2 force=5000 component=z

5.4Pressure and distributed forces on a facestable

Pressure p acts along the normal to the face and is positive when it presses toward the body. A distributed force is a traction \mathbf{t}_d in Pa given as a vector. The traction acting on the face is

where \mathbf{n}_{\mathrm{out}} is the outward normal. The nodal forces are consistent:

The outward normal of each face is determined from the adjacent element: the normal points away from the centroid of the owner element. This gives the correct pressure sign on interior surfaces too (a pipe bore, a drilled hole), where the centroid of the whole model lies outside the material. For flat faces of linear elements (tri3, quad4) the integral is exact. For faces of quadratic elements, tri6 (three-point quadrature) and quad8 (3\times3 Gauss quadrature) are used. On a flat tri6 face the corner nodes receive 0 and the midside nodes one third of the resultant each; on a flat quad8 rectangle the corner nodes receive -1/12 of the resultant and the midside nodes +1/3. These are the classical nodal shares of quadratic elements.

A force on a face (in N) is converted to a traction by dividing by the face area in m²; the distributed_force record is given directly in Pa. Pressure and forces on one face are summed into a single load.

units mm
load surface number 1 pressure=10
load surface number 4 distributed_force=2 component=y

5.5Follower pressurebeta

Follower pressure rotates together with the face and always acts along the normal of the current, deformed configuration. For large displacements and rotations this is the only correct model of gas or fluid pressure. The face is parametrized by coordinates (\xi,\eta), the current coordinates are \mathbf{x}=\sum N_a(\xi,\eta)\,\mathbf{x}_a, where \mathbf{x}_a=\mathbf{X}_a+\mathbf{u}_a. The nodal force at load fraction \lambda is:

This is the Hibbitt scheme (1979); see also Wriggers, Nonlinear Finite Element Methods, 2008, §4.2.5. The faces are tri3, quad4, tri6 and quad8, with the quadrature chosen so that the integral is exact and the sum of forces and moments over a closed surface vanishes to machine precision. The tangent stiffness of the load is analytic and in general unsymmetric (the load on a single face is nonconservative); for a closed surface under constant pressure the sum over faces is symmetric. Details of the formulation and the contribution to the Newton iterations are in the section Nonlinear analysis and dynamics.

The pressure mode is selected in the nonlinear solid analysis settings: "follower" or "dead". When geometric nonlinearity is enabled, follower is the default. In the "material nonlinearity only" mode (small strains) the pressure is dead. For an arbitrary face there is an explicit command nlgeom follower nodes="…" p=… with the face nodes (3, 4, 6 or 8 of them, p>0 presses toward the body).

5.6Acceleration, gravity and rotationstable

Acceleration is a body load on the whole model. A vector \mathbf{a} in m/s² is specified; the body force acting on the material is

Gravity is a special case. For gravity along -z, a_z=-9.81\,\mathrm{m/s^2} is specified: the sign in front of \rho\mathbf{a} is positive, and the vector \mathbf{a} is the direction in which the body force acts. For components along x, y and z the records are summed into a single vector. Nodal forces are computed by consistent integration over the volume, so that the resultant equals \rho V\mathbf{a} to within the integration accuracy.

The material density is required. For a shell without density the acceleration would silently vanish, so the analysis refuses and asks for a density. For assemblies the density is taken per body.

Centrifugal load as a separate static action \omega^2\rho\,r is absent from the list of loads. Rotating structures are treated differently: gyroscopic terms and spin softening are included in the rotor analysis (Rotor dynamics), the initial angular velocity is specified in implicit dynamics, and seismic excitation of the base is specified by the base acceleration in dynamics and spectral analysis (the relative formulation \mathbf{F}=-\mathbf{M}\boldsymbol{\iota}\,\mathbf{a}_0\,k(t)).

load accel=-9.81 component=z

5.7Remote force, moment and displacementbeta

A remote load applies a force \mathbf{F} and a moment \mathbf{M} at an arbitrary point \mathbf{x}_m (which may lie outside the body), transferring them to the selected face, edge, vertex or body. A virtual master node is created at the point. The nodes of the target become dependent according to the rigid-body rule:

and the force and moment are applied to the master node. After the dependent nodes are eliminated, they are distributed over the target as an equivalent load. A remote displacement specifies the degrees of freedom of the master node: three translations and three rotations, each axis can be enabled or disabled. The point, forces and moments are specified in the selected local coordinate system.

remote load surface number 2 point=150,5,5 force=1000,0,0 moment=0,0,5
remote motion surface number 2 point=1,2,3 translate=0.01,0,0 uaxes=1,0,0

5.8Rigid links RBE2 and constraint equationsstable

A rigid region (RBE2) forces the dependent nodes to move together with the master as a single rigid body; the formulas are given in the section Finite element library. The record specifies the master target and the slave targets, any of vertex, edge, face or body.

How it enters the analysis. The nodes of the master and slave targets form a single rigid body; the master is the node closest to the centroid of the master target. For quadratic meshes the midside nodes also become dependent. The dependent degrees of freedom are eliminated by the transformation \mathbf{u}=\mathbf{T}\mathbf{u}_{\mathrm{f}}+\mathbf{u}_0: Lagrange multipliers and penalties are not used, chains of dependencies are collapsed, and the reduced system is solved by the same solver (see Solvers). A conflict of two constraints on one degree of freedom, or a cycle, stops the analysis with a reason code.

So that the link does not contradict other conditions, the following are excluded from it before the analysis:

• nodes fixed by boundary conditions (the number of excluded nodes is reported by a warning);
• nodes already dependent in another rigid link (the intersection is removed);
• nodes of the master target that are already dependent in another link: they are not chosen as master, and if there are no free nodes, the link is skipped.

If no slave node remains, the link is skipped with a warning.

A general constraint equation \sum_i c_i\,u_i = r (MPC) for vertex degrees of freedom can be specified in the project. The first term of the equation is the dependent degree of freedom, the rest are the master ones. Free-form equations are created only in the interface form; there is no script command for them (for frames there is frame constraint). Kinematic joints (revolute, cylindrical, spherical) are assembled from rigid regions and linear equations around virtual master nodes; there is no script command for them and no verification case.

rigid link master vertex number 6 slaves surface number 2

5.9Local coordinate systemsbeta

A local coordinate system (CS) is needed for loads directed along inclined axes. A Cartesian CS is defined by an origin, the z axis and, optionally, the direction of the x axis. The triad is built as follows: \mathbf{e}_z is the normalized z axis vector, \mathbf{e}_x is the x axis orthogonalized by Gram–Schmidt, \mathbf{e}_y=\mathbf{e}_z\times\mathbf{e}_x. If the x axis is not given or is parallel to z, an arbitrary perpendicular is constructed, and in that case the CS can be rotated about z only by specifying x.

The component x of a record bound to a CS, with a given magnitude q, means the vector q\,\mathbf{e}_x (likewise y, z). The rotation is done once, in the load compiler; all quantities further down the pipeline are already global. The CS origin is used only for remote points; a concentrated force as a free vector does not read it. The CS binding acts on force, distributed force, acceleration, remote force and remote displacement.

cs turned origin=0,0,0 axisZ=0,0,1 axisX=1,1,0
load surface number 2 fx=1000 cs=turned

5.10Bolt preloadexperimental

A virtual bolt models joint tightening without building the thread: two flanges are clamped by an axial force F_0, and the bolt carries it as an elastic bar. It is suitable for estimating joint opening and bolt load in assemblies.

How the model is built. Two master nodes are created at the centroids of the bearing faces of the head and the nut; each is tied to the nodes of its face by a rigid region. A Timoshenko beam with cross-section area A and axial stiffness k_b=EA/L is introduced between the master nodes, where L is the distance between the centroids. A self-equilibrated pair of forces \pm F_0 along the bolt axis is applied to the master nodes. The force can be specified either as the force F_0 or as a tightening torque: F_0 = T/(K\,d), where K is the torque coefficient (0.2 by default) and d is the diameter. Cross-section area: nominal \pi d^2/4, the tensile stress area A_s=\tfrac{\pi}{4}(d-0.9382\,P)^2 with the coarse thread pitch (ISO 898-1), or specified manually. The property class (8.8, 10.9, 12.9) determines the default yield strength: 640, 900 and 1080 MPa; the value can be specified manually.

Postprocessing for each bolt: working force F_{\mathrm{work}}=F_0+k_b\,\Delta u, where \Delta u is the relative axial displacement of the master nodes in the solution; stress \sigma=F_{\mathrm{work}}/A; safety factor \mathrm{SF}=\sigma_y/\sigma. Flags are raised for bolt yielding, low safety factor (below 1.25) and joint opening, when the working force drops almost to zero.

bolt surface number 10 surface number 20 preload=20000 diameter=12 class=10.9

5.11Springs and point massesbeta

A spring is an elastic connection between two regions of the model or between a region and a fixed ground; a point mass is a single-node element with no stiffness. Both elements are described in the section Finite element library; here we cover how they enter the list of conditions.

A spring is created by a record with its own target (face, edge or vertex) and a stiffness k in N/m; for the ground, the word ground is given instead of the second target. A transverse stiffness, a torsional stiffness and a preload force F_0 can also be specified. Postprocessing gives the elongation \Delta L and the working force F=F_0+k\,\Delta L. A point mass is the mass command; in the product it affects the modal and dynamic analysis of frames made of beams and shells.

spring surface number 3 ground k=2e6 preload=100
mass node=1000,0,0 m=100

5.12Cyclic symmetryexperimental

Cyclic symmetry makes it possible to analyse one sector of a disk, wheel or bladed ring instead of the entire structure. Formulation: the sector is repeated N times about an axis, the sector angle is \Delta\theta=2\pi/N. The nodes of the lower boundary of the sector are rotated by +\Delta\theta about the axis (Rodrigues' formula) and matched to the nodes of the upper boundary within a tolerance. For each pair three constraint equations are imposed

This is the harmonic with index k=0: the deformed state is identical in all sectors. After the analysis the result can be expanded to the full 360^\circ: displacements are rotated, and the stress tensor is rotated as a tensor (the von Mises equivalent stress does not change).

For harmonics k\ge1 (Floquet–Bloch waves) the code forms a constraint with the complex phase factor \mathrm{e}^{\mathrm{i}k\Delta\theta} in terms of real and imaginary parts. The static solver is real and does not carry such a field, so harmonics k\ge1 are not used in the analysis.

The user specifies the lower and upper boundaries of the sector (faces), the number of sectors, the axis and the tolerance: in the interface, the "Cyclic symmetry" section of the boundary conditions step; in a script, the cyclic command. The sector analysis is linear statics with constraint equations.

cyclic surface number 1 surface number 2 sectors=12 axis=0,0,1 tol=0.05

5.13Steps, phases and time laws of loadingbeta

Load records have no concept of a "stage" of their own, but several analyses use them in different ways.

• Linear statics. The load-step slider shows the scaling of the single solution: the fields are multiplied by the fraction s/N. There is no real step-by-step loading here.
• Nonlinear statics. The load and prescribed displacements grow in proportion to the fraction \lambda; the load path loadpath="1,0" includes unloading, "1,-1,1" a cycle. Details are in the section Nonlinear analysis and dynamics.
• Harmonic response. Force, pressure and distributed force can carry a phase angle \Delta\varphi: F(t)=F\cos(\omega t+\Delta\varphi) (the phase parameter, in degrees). Acceleration carries no phase. A static analysis ignores the phase and reports it.
• Dynamics. The load is multiplied by an amplitude law k(t) given by the table amplitude="t:k,…".
• Birth and death of bodies. Construction and excavation stages switch bodies on and off as a whole: Nonlinear analysis and dynamics.

5.14What happens to records when the geometry is rebuiltbeta

Automatic geometry healing, sewing, body suppression and Boolean operations can remove a face or vertex to which a record was attached. The rule is: a record is never deleted silently.

• The record stays in the project state and is marked as orphaned ("target lost"); the healing report states how many records lost their attachment and which ones.
• An orphaned record is not transferred to other objects. If several faces have merged into one, none of the original records moves to the merged face: the transfer would change the physical meaning until the user selects the face themselves.
• The analysis does not pretend the record is absent. In scripts the analysis refuses with a list of the lost records and their targets. In the interface the record is excluded from the analysis with a warning by name; if no support remains after that, the analysis is blocked.
• A record that has no receiver in the solver (for example, an automatic connection of a midsurface junction while the junction nodes have not been built) is not turned into placeholder equations: it is marked as unimplemented, and the analysis refuses with its name.
• Automatic geometry healing receives the list of faces on which records are placed and does not remove them.

An assembly of several bodies requires a decision on how the bodies are connected at the interface. For each pair of bodies, WebCAE selects one mode: tied, spring connection, frictionless contact, frictional contact, or no connection. The mode is set in the "Connections" panel (a matrix of pairs) or by the contact command; a default mode can be set for the whole assembly and overridden for an individual pair.

The modes differ substantially in formulation. Bonding (bonded, tied) is linear: the interface transmits any forces and does not open. Frictionless contact is unilateral: the bodies transmit only compression and can separate, so the problem is nonlinear and is solved in load steps. Friction adds a dependence on history. Not every mode is available in every type of analysis. If a pair is assigned a mode that the selected analysis does not support, the analysis refuses explicitly; it does not replace the mode with another one and does not run without the pair.

Modal analysis, linear buckling and harmonic response with frictionless contact are linearized: the normal stiffness is taken from the gap at the moment of analysis, and the set of active pairs is not updated in the eigenproblem. The desktop (native) application is the reference; the browser supports tied and spring connections.

6.1Contact pair detectionbeta

The command contact pairs all|visible <mode> and the button in the panel create pairs for all touching bodies of the assembly. The criterion is the distance between the body surfaces, not the intersection of bounding boxes: a box knows nothing about shape, and a part sitting in a cutout of another is indistinguishable to it from a part that merely neighbours the cutout. The distance is computed from the tessellation of the boundary that the OpenCascade kernel builds on import. For each pair the type of interface is determined: touching, a measurable gap (a design clearance), or an overlap with an intersection volume (a press fit or an assembly error). The tolerance is relative: it is tied to the size of the assembly. Pairs near the threshold can be refined by the exact distance between bodies in the B-Rep.

The "Check connectivity" button finds groups of bodies and "floating" components that are connected neither by a support nor by a connection mode. In modal analysis, buckling and harmonic response the analysis refuses with a list of unconnected bodies until the groups are connected.

contact pairs all tied
contact body1 body2 frictionless

6.2Bonding by shared nodes (bonded)stable

The most reliable connection: at the meshing stage the touching faces of the bodies receive shared nodes, and the displacement field is continuous across the interface by construction. The solver has no additional conditions; the bodies simply share degrees of freedom. This bonding introduces no error other than the geometric one and works in all types of analysis.

It is controlled by mesh parameters: bonding=0|1 enables bonding, glue=auto|contact|full selects a fast search for contacts or a full search for intersections, interference=cut|split|ignore determines what to do with a real overlap. In the cut mode the disputed volume is given to the smaller body (a pin displaces the plate material), the interface becomes an exact contact, and the volume equals the union. In the split mode the overlap is fragmented: the volume is exact, but the fit sliver may turn out to be a thin film that breaks the mesh. In the ignore mode the overlap is kept and the material is counted twice.

mesh tet size=125 bonding=1 glue=full

6.3Bonding on nonconforming meshes (tied)beta

Tied bonding is needed when the meshes of two bodies at the interface differ: different density, separate nodes. For each node of the slave side, a projection point on a face of the master side is found, and the displacement of the slave node is expressed through the face nodes with barycentric weights w_k:

This is three constraint equations per node; for quadratic meshes the midside nodes enter the pairs. The constraint is enforced exactly: the slave degrees of freedom are eliminated from the system, as for a rigid link (Solvers), and Lagrange multipliers are not introduced. The method is node-to-surface (NTS); it is based on the general contact formulation in Wriggers P., Computational Contact Mechanics, 2006.

Settings of the "Connections" panel:

• Formulation. "Auto" selects the kinematic constraint for the direct solver and the spring one for the iterative solver with an algebraic multigrid preconditioner on large assemblies. "Kinematic (Exact MPC)" and "Penalty" set the choice explicitly.
• MPC gap limit. The ratio of the slave–master distance to the size of the master face; a pair beyond the limit is rejected (0.8 by default).

Before the analysis an interface check is performed: conforming faces with shared nodes are skipped, and an interface with no pair found, a gap above the limit, or an analysis with a single region refuse with a stated reason.

contact body1 body2 tied

6.4Penalty connection (penalty)beta

A penalty pair uses the same pair geometry as tied, but the connection is compliant: a spring of stiffness k is placed between the slave node and the point on the master face. As k\to\infty the solution tends to tied. The mode is needed when an exact constraint worsens the conditioning on large assemblies or when the interface must yield slightly.

The user specifies a multiplier of the computed stiffness (1.0 by default, 10.0 for bonding) or their own value of k in N/m.

6.5Frictionless contact (frictionless)stable

The bodies transmit only compression along the surface normal. The Signorini conditions for the gap g and the normal pressure p:

with the gap g = g_0 + \mathbf{n}\cdot\bigl(\mathbf{u}_s - \sum_k w_k\,\mathbf{u}_{m_k}\bigr) for a node s of the slave side and a point on the master face. The conditions are regularized by a penalty with the potential \Pi_c=\tfrac12\,k\sum\langle -g\rangle^2 (where \langle x\rangle=\max(x,0)): the pressure is p=k\langle -g\rangle, and the tangent stiffness of an active pair k\,\mathbf{n}\mathbf{n}^{\mathsf T} is symmetric and positive definite. A pair is active upon penetration g<0; so that the set of active pairs does not chatter, the switching is supplemented by hysteresis.

The analysis proceeds in load steps. At each step the nonlinear system is solved by the Newton method with damping and a nonmonotone line search (Armijo condition); the penalty grows to the specified value over the steps. The number of steps is set on the analysis tab. The pairs are searched again after a step on the deformed configuration, and the gap value is carried over so that the reopening of a pair does not add a spurious offset.

The name mortar in the product refers to a formulation with a weighted gap. The constraint is built on the nodes of the slave face: each slave node is matched with up to three nodes of the master face with barycentric weights, and the gap and the normal are taken at the projection point. Full segment-to-segment integration with polygon clipping (Popp A., Gee M., Wall W., A finite deformation mortar contact formulation using a primal–dual active set strategy, 2009) is not available in the product. The augmented Lagrangian (Uzawa method) is implemented in the kernel as a research option; it is not included in the product path, and the solution relies on the penalty. Therefore a residual penetration of the order of p/k remains.

The user specifies the contact stiffness k (N/m) or a multiplier of the computed one, the number of load steps and the bodies of the pair. Results: contact pressure over the faces, status (open, contact), contact force, history over the load steps.

contact frictionless
mesh tet size=0.25 order=1 bonding=0
load surface number 8 force=-25 component=z
solve
extract contact_pressure target="surface number 7"

6.6Coulomb friction contactbeta

Friction limits the tangential force in contact by Coulomb's law: \lVert\mathbf{t}_t\rVert\le\mu\,p_n, with sticking below the limit and sliding at the limit. The law is irreversible, so the answer depends on the loading history: for each active contact point the accepted state is stored.

In implicit transient dynamics, for each slave node the position on the master face and the accumulated elastic slip \mathbf{e} are stored. The tangential force is regularized by a length s_{\mathrm{tol}}, the elastic slip at the Coulomb limit:

where the trial elastic slip \mathbf{e}_{\mathrm{tr}}=\mathbf{e}_{\mathrm{old}}+\Delta\mathbf{s}_t is the sum of the previous value and the tangential displacement of the node relative to the material point of the master face, and is then smoothly limited by the length s_{\mathrm{tol}}. The limitation is continuously differentiable, so the tangent stiffness is defined in the transition zone as well. The law gives genuine static friction: a block under a horizontal force below the limit stays in place. The tangent is found by central differences of the pair force and is consistent with it, including the motion of the projection upon a change of facet.

The user specifies the coefficient \mu in the range 0<\mu\le10 (0.3 by default in the form) and optionally the length s_{\mathrm{tol}} (the slip parameter, m; by default a fraction of the edge of the contact face). Results: contact pressure, status (open, stick, slip), stick and slip counters, work of friction forces, number of facet changes.

The nonlinear contact kernel also has a quasi-static variant with a return to the Coulomb cone (a stick trial and a radial return at \Phi=\lVert\mathbf{t}_{\mathrm{tr}}\rVert-\mu p_n>0, a consistent tangent for sliding); it has been verified on analytical problems and is not connected in the product.

contact body1 body2 frictional mu=0.3
solve transient dt=0.001 tend=0.3

6.7Contact at large displacements and self-contactbeta

In geometrically nonlinear statics of solids, contact is built on the deformed configuration: pairs are searched again at every assembly, that is, at every Newton iteration, and not only between load steps. The surface consists of the triangles of the body faces (hexahedron faces are split into two triangles). For a slave node the closest point on the master triangles is found (Ericson projection), and the gap is

where \mathbf{n} is the contact normal, an interpolation of the nodal normals of the master face: it is continuous across edges and vertices, otherwise the force direction jumps over an edge between noncoplanar facets. Contact is active at g<0; a penalty \varepsilon\,g\,\mathbf{n} acts on the slave node in the kernel's internal-force convention, and an equal and opposite force with barycentric weights acts on the nodes of the master triangle, so that Newton's third law holds exactly. The tangent is consistent with the force and is found by central differences over the twelve coordinates of the pair (it is unsymmetric). Opposing normals cut off spurious pairs: a node of the bottom face of a body does not "penetrate" the top face of the same body through its thickness.

The penalty is chosen automatically: \varepsilon=\kappa\,E\,h, where E is the smaller of the moduli of the two bodies, h is the mean edge length of the faces, and \kappa is a multiplier set by the user (10 by default). Self-contact is enabled by a checkbox on the analysis tab (where \kappa is also set): the surface of a body contacts itself, and the triangles in the topological neighbourhood of the node are excluded from the search. For nodal lists of beams and shells there are the commands nlgeom contact and nlgeom beamcontact.

nlgeom on method=load steps=10
contact body1 body2 frictionless
solve

6.8Contact of beams and shellsexperimental

For beams, line-to-line contact is used (Wriggers P., Zavarise G., On contact between three-dimensional beams undergoing large deflections, 1997): for a pair of segments the closest points of the axes \mathbf{x}_1(\xi), \mathbf{x}_2(\zeta) are taken, and the gap is

where r are the section radii, and the force is a penalty along the normal, transferred to four nodes with linear weights (the sum of forces and the sum of moments are zero). For parallel segments the midpoint of the overlap is taken and the force is multiplied by the overlap fraction so that the result does not depend on the shift of the meshes. Self-contact of a beam is possible.

For shells, contact is built on the midsurface with an offset of half the thickness: g = s\,(\mathbf{x}_s-\mathbf{p})\cdot\mathbf{n}-(o_s+o_m), where o are the offsets and s=\pm1 is the side. The side of the master two-sided shell is remembered per node–triangle pair; otherwise a node that has passed through the midsurface within a single iteration would be pushed out on the opposite side.

Pairs are specified by the commands nlgeom contact and nlgeom beamcontact with lists of nodes in the model numbering, section radii and a penalty.

6.9Contact in explicit dynamicsbeta

The explicit scheme uses eroding contact: a unilateral penalty node-to-segment over the whole outer surface of the live elements. The surface is updated as erosion proceeds, and the impactor meets the newly opened cavity. The penalty stiffness is set as a fraction \alpha of the stiffness of the master face (the contactStiffness parameter, 0.25 by default): too small a value gives deep penetration, too large a value gives an energy gain visible in the energy balance. The history outputs the contact energy and the overall balance.

The scheme, the erosion and the criteria are described in the section Nonlinear analysis and dynamics.

explicit erodingContact=on contactStiffness=0.25

6.10Thermal contactbeta

Thermal contact links the temperatures of two separate bodies through the interface conductance h_c: q_c=h_c\,(T_A-T_B). It does not depend on the mechanical contact and does not affect the mechanics. The formulation, the gap tolerance and the results are described in the section Heat transfer and multiphysics; the comparison with the reference is there as well.

The kernel also has a conductance model dependent on pressure and gap, h_c(p_n,g_n), of the Cooper–Mikic–Yovanovich type (Cooper M., Mikic B., Yovanovich M., Thermal contact conductance, 1969): for a closed interface h_c=h_0+h_p\,(p_n/p_{\mathrm{ref}})^{\alpha}, for an open gap h_c=k_{\mathrm{gas}}/(g_n+g_0) with a cutoff at the limiting gap. This model has been verified only by an internal kernel test; it is not connected in the interface or in scripts: the user specifies a constant h_c, and the mechanical contact pressure does not affect the conductance.

contact body1 body2 thermal hc=1000

6.11Cohesive zones (CZM)experimental

A cohesive zone describes delamination of adhesive and composite joints by a traction–separation law. In the product the law is bilinear: an elastic branch with stiffness K up to the peak traction T_{\max}, then linear softening to complete failure at the opening \delta_f:

where D\in[0,1] is a scalar damage variable that does not decrease on unloading; G_c is the fracture energy (the area under the diagram). For mixed mode the effective opening is \delta_m=\sqrt{\delta_n^2+\delta_s^2}, and the failure criterion is the power law of Benzeggagh and Kenane (Benzeggagh M., Kenane M., 1996):

where G_I=\tfrac12 K_n\delta_n^2, G_{II}=\tfrac12 K_t\delta_s^2, and the exponent \eta is 1.45 by default. In mixed mode \delta_0 and \delta_f are computed for the current mode ratio; damage initiation is equivalent to the quadratic nominal stress criterion under proportional loading (Camanho P., Dávila C., 2002).

The user specifies a pair of surfaces (master and slave) and the parameters; in the "Connections" panel a pair has a cohesive contact checkbox with fields for the strength T_n and the energy G_{Ic}, and the other parameters are available through the script command: stiffnesses K_n, K_t (N/m³), strengths T_n, T_s (Pa), fracture energies G_{Ic}, G_{IIc} (J/m²), the exponent \eta. Result: the maximum and mean damage D over the interface, the peak normal and shear traction, the maximum opening, the number of node pairs processed.

It is important how the law enters the analysis. In the product the interface is evaluated as a postprocessor: from the displacements of the completed solution, the jump at each slave node is compared with the closest master node, and the damage and tractions are computed. Equilibrium with softening is not recomputed inside the global iteration, so the result does not show how the load-carrying capacity of the joint drops as the crack grows. Integration of the law into the nonlinear solution (contact pairs with damage) exists in the kernel at the level of a node-to-surface pair, but is not exposed in the interface or in scripts.

contact cohesive master="surface number 3" slave="surface number 4" kn=1e12 kt=1e12 tn_max=30e6 gc_n=450

6.12The "converged" flag for contactbeta

In a contact problem, Newton convergence does not guarantee physical equilibrium: the set of active pairs may chatter, and a penalty problem may "converge" to a solution with a body floating without support. Therefore the "converged" flag is composed of several conditions, not only of Newton iterations.

• For frictionless contact in step-by-step mode the full load must be reached: if the steps are cut off or cancelled, the flag is false.
• The physical equilibrium residual, computed over the free degrees of freedom (without stabilizing stiffness), must not exceed the tolerance in percent of the load. Only this condition does not depend on what Newton did.
• If Newton has not converged, a fallback criterion of small penetration is allowed, but only if there are pairs with penetration in the final solution. An analysis without contact under load does not receive "converged" because of an empty set.
• The message names the number of active pairs, the equilibrium residual in percent of the load and the maximum displacement; for a body without support ("mechanism") or a contact that has not closed this is stated explicitly.

For nonlinear statics the flag is built on the Newton–Raphson criterion (see Nonlinear analysis and dynamics), and the result is not accepted if the problem has not converged. A mode that the analysis does not support refuses: bodies with shared nodes are not taken into contact, and friction in statics does not fall back to a frictionless pair.

Structural elements are assembled into a single global system: each node has six degrees of freedom [u_x, u_y, u_z, \theta_x, \theta_y, \theta_z]. Solid elements carry only three displacements; the rotational degrees of freedom of nodes to which no beams, shells or rigid links connect have no stiffness and are fixed automatically. In this way solids, shells, beams, springs and point masses coexist in one model.

This section describes which elements exist, in what formulation they are built, and where their limits lie. A short catalogue with type codes and node counts is on the page /docs/elements; the summary table at the end of the section has been checked against the kernel code.

7.1General conventionsstable

All solid elements are isoparametric (except Tet10, see below). The gradients of the shape functions with respect to the physical coordinates are taken through the Jacobian of the mapping:

Strains and stresses are written in Voigt order [\varepsilon_{xx},\varepsilon_{yy},\varepsilon_{zz},\gamma_{xy},\gamma_{yz},\gamma_{zx}] with engineering shear strains \gamma = 2\varepsilon. The matrix \mathbf{B} is assembled from the shape function gradients: the block of node i is \partial N_i/\partial x in the row \varepsilon_{xx} and column u_x, \partial N_i/\partial y in the row \varepsilon_{yy} and column u_y, \partial N_i/\partial z in the row \varepsilon_{zz} and column u_z, and the shear rows carry two gradients each (for example, \gamma_{xy}: \partial N_i/\partial y in the column u_x and \partial N_i/\partial x in the column u_y). Isotropic elasticity:

The element stiffness matrix is Gauss quadrature:

The geometrically nonlinear path (Total Lagrangian) uses the same shape functions and the same quadratures: there the tangent consists of a material and a geometric part, \mathbf{K}_T = \mathbf{K}_{\text{mat}} + \mathbf{K}_{\text{geo}}, and the default material is Saint Venant–Kirchhoff (more in the section Nonlinear analysis and dynamics). Materials with anisotropy, plasticity and hyperelasticity are described in the section Material models; here only the element itself is considered.

Boundary conditions and constraints (rigid regions, constraint equations) are eliminated from the system exactly, by the master–slave scheme: the system \mathbf{T}^{\mathsf T}\mathbf{K}\,\mathbf{T}\,\mathbf{u}_r = \mathbf{T}^{\mathsf T}\mathbf{f} is solved, not a penalty one. This does not spoil the conditioning and does not give approximate enforcement of the constraints.

Units. Geometry is specified in millimetres (or in units declared by the units command), and at the solver boundary everything is converted to SI: metres, newtons, pascals, kilograms. Section dimensions, shell thicknesses and masses of point elements are entered by the user in the same units as the geometry (masses in kilograms).

7.2Linear tetrahedron Tet4stable

A four-node simplex with constant strain. This is the main element of automatic mesh generation for arbitrary CAD geometry: it is robust with respect to body shape but transmits bending poorly.

The shape functions coincide with the barycentric coordinates, N_i = L_i, i = 1,\dots,4. The gradients \nabla L_i are constant, so \mathbf{B} is also constant, and exact integration gives

that is, one Gauss point at the barycentre. The stress is constant over the element, \boldsymbol{\sigma} = \mathbf{C}\mathbf{B}\mathbf{u}_e; the nodal field is built by averaging over adjacent elements.

The element passes the patch test on an arbitrarily distorted mesh: an affine displacement field is reproduced to machine precision. But constant strain means parasitic stiffness in bending. On a cantilever bar of \approx 1400 tetrahedra the first natural frequency is overestimated by 22\,\% (the same element size with Tet10 gives -0.16\,\%). For bending bodies use Tet10.

mesh tet size=10 order=1

7.3Quadratic tetrahedron Tet10stable

A ten-node element: four vertices and six edge midpoints. The quadratic displacement field gives a linearly varying strain within the element and removes the main part of the bending locking of Tet4 on the same mesh.

Shape functions through the barycentric coordinates (N_i is a vertex, N_{ij} is the midpoint of the edge between vertices i and j):

Gradients: \nabla N_i = (4L_i - 1)\nabla L_i, \nabla N_{ij} = 4(L_j\nabla L_i + L_i\nabla L_j). Integration is the four-point rule of degree 2 with points (a,b,b,b) and permutations, a = 0.5854101966, b = 0.1381966011, weights V_e/4. For an element with straight edges the integrand \mathbf{B}^{\mathsf T}\mathbf{C}\mathbf{B} is a polynomial of degree two, and the rule is exact.

An important implementation feature: in the main assembler Tet10 is subparametric. The geometry is defined by the four vertices (the mapping is linear, the Jacobian is constant), while the displacement field is defined by ten nodes. The positions of midside nodes displaced onto the curve of a CAD surface are not used in the stiffness matrix: the element is integrated as a straight-edged one. This preserves nondegeneracy on curved boundaries, but the volume and stiffness of elements at a curved boundary are computed along the chord rather than the arc. For a smooth boundary the error decreases with element size, but there is no separate accuracy gain from curved edges. This is how the main assembler is built; in the multi-body analysis branch with tying and contact, other element implementations are used, and their mapping may differ.

Surface loads on Tet10 faces are integrated consistently over six-node triangles (Tri6), and are not split equally among the vertices.

The kernel computes the stress at the element centroid, where the contribution of the corner nodes vanishes (4L_i - 1 = 0). For the nodal picture the postprocessor computes the stress at each node of the element and averages over the number of adjacent elements, so at material boundaries and near stress concentrators the field is smoothed.

mesh tet size=10 order=2

7.4Hexahedron Hex8stable

An eight-node trilinear element for regular volumes: extended bodies, bars, plates and layers at the wall of a hole. On a structured mesh it gives a smooth stress field and fewer elements than tetrahedra.

Shape functions in natural coordinates (\xi,\eta,\zeta) \in [-1,1]^3:

Integration is full, 2\times2\times2 Gauss points at \pm 1/\sqrt{3}, with unit weights. The matrix \mathbf{K}_e is assembled by the formula from the "General conventions" section without any modifications.

It is important to state plainly here: in the linear path Hex8 has no B-bar (selective volumetric projection), no incompatible modes and no reduced integration. The element suffers shear locking in bending when there are few elements through the beam height, and volumetric locking at \nu \to 0.5. The remedy is mesh refinement through the height or a switch to Hex20. Examples: a cantilever beam on a structured 128\times8\times8 mesh gives a 1.03\,\% deviation from Euler–Bernoulli; a cantilever on a 96\times6\times6 mesh (six elements through the height) gives about 2\,\% because of shear locking. In explicit dynamics Hex8 is a separate one-point element with hourglass control after Flanagan–Belytschko; it is unrelated to the matrix of this subsection.

The stress is computed at the element centre, \xi = \eta = \zeta = 0.

mesh hex size=5 order=1

7.5Hexahedron Hex20 (serendipity)stable

A twenty-node element: eight vertices and twelve edge midpoints, with no nodes at the face centres or at the body centre. The quadratic interpolation makes it substantially softer than Hex8 in bending: a single layer of elements through the thickness is enough for a reasonable estimate. It is used for rings around holes, fatigue stress raisers and contact problems where a smooth boundary geometry is needed.

Shape functions of the serendipity family. A vertex with coordinates (\xi_i,\eta_i,\zeta_i) = (\pm1,\pm1,\pm1):

The midpoint of an edge parallel to the \xi axis (\xi_i = 0):

and analogously for edges along \eta and \zeta. Integration is 3\times3\times3 (27 Gauss points, \pm\sqrt{0.6} and 0 with weights 5/9, 8/9, 5/9). Unlike Tet10, Hex20 is fully isoparametric: the coordinates of all twenty nodes are used, so the curvature of edges on a CAD boundary is accounted for in \mathbf{J}.

The stress is computed at the element centre. The element passes the patch test on a distorted mesh to machine precision; on a single cantilever the relief of locking relative to Hex8 is substantial.

mesh hex size=5 order=2

7.6Hexahedron Hex27 (triquadratic)experimental

A twenty-seven-node element of the full quadratic Lagrange basis: the centres of the six faces and the body centre are added to the nodes of Hex20. The shape functions are products of one-dimensional quadratic Lagrange polynomials:

where (a,b,c) \in \{-1,0,1\}^3 is the position of node n. Integration is 3\times3\times3. The element is more complete than the serendipity one: more accurate under strong distortion and bending, but more expensive.

Hex27 is not part of the product mesh pipeline: the mesher does not build it, so the element is available only when connectivity is supplied directly. The Hex27 mass matrix is not implemented in implicit dynamics, and such an analysis refuses with a stated reason; modal and thermal analysis work.

7.7Wedge Wedge6stable

A six-node linear prism: a triangle at the base and linear interpolation through the height. It arises in prismatic layers at the wall (inflation) and in transitions between hexahedra and tetrahedra.

The shape functions are the product of the barycentric coordinates of the base and a linear factor in height:

where k is the number of the base vertex, with the minus sign for the lower triangle (\zeta=-1) and plus for the upper. Integration is the product of the three-point rule on the triangle (points (1/6,1/6), (2/3,1/6), (1/6,2/3), weights 1/6) and the two-point Gauss rule through the height: six points in all. The sum of the weights equals the volume of the reference prism, 1.

The element is linearly complete, so it passes the patch test at a sewn interface with Tet4. A degenerate (non-positive) Jacobian stops the analysis with a message.

7.8Pyramid Pyr5beta

A five-node linear pyramid: a quadrilateral base and an apex. It is needed as a transition element in hybrid meshes, where a quadrilateral face of a hexahedron meets two triangular faces of tetrahedra.

A direct "degeneration" of a hexahedron into a pyramid (merging the four top nodes) is not linearly complete and does not pass the patch test, so the rational basis of the pyramid element is used (Zgainski et al.). In the implementation it is written in "collapsed" coordinates, \hat a, \hat b \in [-1,1], w \in [0,1), where the physical coordinates of the section shrink as (1-w):

In the physical coordinates \xi = \hat a(1-w), \eta = \hat b(1-w) the basis contains the rational term \xi_i\eta_i\,\xi\eta/(1-w), hence the name. Integration: 2\times2 Gauss points over (\hat a, \hat b) and two Gauss nodes over w on [0,1], eight points in all. The points do not fall on the apex w=1, where the mapping is degenerate. The volume equals one third of the product of the base area and the height, as for a true pyramid. The stress is computed on the pyramid axis at w = 1/2.

7.9Hybrid meshes: Hex8, Wedge6, Pyr5, Tet4beta

A hybrid mesh joins a layer of hexahedra or prisms (at the wall of a hole, in a concentration zone, in an extended block) with a tetrahedral remainder. The interface is made in two ways: by sewing at shared nodes (conformal; transition through Pyr5 and Wedge6) or by node-to-surface constraint equations. For a conformal interface the Hex8 face is divided toward the centre into pyramids, and the pyramid on the side of two Tet4 is divided into two Tet4 along the same diagonal, so the triangular faces coincide.

The correctness of the interface is checked by the patch test: on a sewn "prisms + tetrahedra" mesh a uniform stress field is reproduced to 5\cdot10^{-13}. The stress in a hybrid mesh is computed at the node point of each element and averaged over adjacent elements: the centre of the first layer of hexahedra at the wall of a hole is half a layer away, and the value at the centre underestimates the stress peak by 12–15 %.

7.10Quadratic Hex20–Pyr13–Tet10 transitionexperimental

An interface of quadratic elements requires the traces on the common face to match: the eight-node face of Hex20 must meet two six-node faces of Tet10. There are two prototypes for this; they are not accessible from the user pipeline (script, mesher).

The first is the "common trace" for flat Hex20–Tet10 interfaces. The face carries nine nodes: four corners, four edge midpoints and the midpoint of the triangulation diagonal. The trace-matching condition is three linear relations on the nodal values (the nine values must be values of a single complete quadratic polynomial on a parallelogram). They are applied as constraint equations, and the master degrees of freedom are chosen by a global Gauss–Jordan reduction in exact rational arithmetic, so that the rigid modes and all quadratic fields are preserved. The mesh is assembled by a separate adapter (cutting a rectangular cavity out of a Tet10 body and inserting a Hex20). On a mesh of 96 Tet10 and one Hex20 with an exact quadratic field prescribed on the boundary, the displacement deviation was of the order of 10^{-8}.

The second is a converter of a linear hybrid mesh into a conformal quadratic Hex20–Pyr13–Tet10 mesh with shared Quad8 and Tri6 faces. The thirteen-node pyramid Pyr13 is a second-order pyramid element (5 vertices and 8 edge midpoints). Pyr13 as a stiffness element is absent in the solver. Only the shape functions and a check of the Jacobian at Gauss points exist on the mesh preparation side; this converter was not included in the release, and analysis on such a mesh is unavailable.

7.11Locking: what exists and what does notstable

Locking is parasitic stiffness of an element when its interpolation cannot reproduce the required deformation. Below is an honest list of what is implemented.

• Shear locking of bending elements. The MITC4 shell removes it by assumed transverse shear strains (subsection MITC4 shell). The DKT triangular shell imposes Kirchhoff constraints at the midpoints of the sides and contains no transverse shear. The Timoshenko beam is built by the standard \Phi-modification of the Hermitian matrix; for small \Phi it passes into the Euler–Bernoulli element without locking.
• Shear locking of solid elements. It is treated by order: Tet10, Hex20, Hex27 instead of Tet4, Hex8, Wedge6, Pyr5. For Hex8 there is no B-bar, no incompatible modes and no reduced integration (in the linear path).
• Volumetric locking at \nu \to 0.5. There is no mixed formulation: as \nu approaches 0.5, linear elements become stiffer than they should be. Quadratic elements reduce the locking but do not remove it; for nearly incompressible materials use hyperelastic models (see Material models).
• F-bar. For multiplicative finite-strain plasticity \mathbf{F} = \mathbf{F}_e\mathbf{F}_p on Hex8, Hex20, Hex27 and Tet10 elements, a volumetric averaging of the deformation gradient (F-bar) is introduced: \bar{\mathbf{F}} = (\bar J / J)^{1/3}\mathbf{F}, where \bar J is the element average. This code is present only in the working version and was not included in the release; consider it experimental.
• Explicit dynamics. Hex8 is integrated at one point with an exact volume-averaged gradient and hourglass control after Flanagan–Belytschko (viscous and stiffness forms).
• Membrane locking of shells. The membrane of the quadrilateral shell is an ordinary bilinear Q4 with full integration, without incompatible modes; on coarse meshes in in-plane bending it is stiffer than the exact solution. For the triangular shell the membrane is a constant-strain CST.

7.12Timoshenko beamstable

A two-node spatial bar element for frames, trusses and stiffening ribs. Six degrees of freedom per node, twelve per element. The material and section are specified by six numbers: E, G, A, I_y, I_z, J. The shear modulus is taken as G = E/(2(1+\nu)).

The stiffness in local axes (x along the axis from node 1 to node 2) is in closed form, with no integration over points. The axial and torsional parts:

Bending in the xy plane (degrees of freedom v_1,\theta_{z1},v_2,\theta_{z2}):

In the xz plane the matrix is the same with the replacement I_z \to I_y, \kappa_y \to \kappa_z and opposite signs in the displacement–rotation couplings, since \theta_y = -dw/dx. As \Phi \to 0 the Euler–Bernoulli element (a long slender beam) is obtained; for large \Phi the shear compliance grows; there is no locking. The shear coefficient in the product is fixed: \kappa_y = \kappa_z = 5/6 (rectangle, Timoshenko classic); the value does not depend on the section shape.

Orientation. The local axis x runs from node 1 to node 2. The reference vector is the global Z (for a nearly vertical beam, \lvert x_z\rvert > 0.9999, the global Y). Then y = v_{\text{ref}} \times x, normalized, z = x \times y. The matrix is transformed to global axes by \mathbf{K}_g = \mathbf{T}^{\mathsf T}\mathbf{K}_l\,\mathbf{T}, \mathbf{T} = \operatorname{blockdiag}(\mathbf{R},\mathbf{R},\mathbf{R},\mathbf{R}).

The internal forces at the ends are recovered as \mathbf{f}_l = \mathbf{k}_l\,\mathbf{u}_l with the convention: N > 0 is tension, T > 0 is a moment by the right-hand rule along the axis from node 1 to node 2; shear forces and bending moments are given as internal section forces.

frame node at=0,0,0
frame node at=1000,0,0
frame member 1 2 section=rect w=40 h=80
frame support node=1 dofs="all"
frame load node=2 force=0,0,-5000

7.13Extended beam: shear, offsets and nonstructural massstable

The extended beam type extends the Timoshenko element to twenty-two numbers and removes three assumptions of the basic one: the axis through the centroid, the section in principal axes, and zero roll.

• Roll and reference vector. Rotation of the section about the beam axis by a given angle and an arbitrary reference vector that defines the xy plane (instead of the automatic one).
• Product of inertia I_{yz}. If the section axes are not principal, the element is rotated by the principal axes angle \alpha, \tan 2\alpha = -2I_{yz}/(I_y - I_z), where bending in the two planes is decoupled.
• Offset of the centreline. The section centre is offset from the node at each end independently (a rigid arm): \mathbf{u}_c = \mathbf{u}_n + \boldsymbol{\theta}\times\mathbf{r}, \boldsymbol{\theta}_c = \boldsymbol{\theta}_n. The flexible length is the distance between the section centres at the ends. This models a stiffening rib lying on the surface of a sheet, with the parallel axis theorem taken into account: instead of artificially increasing the moments of inertia, the eccentricity is transferred geometrically.
• Shear centre. The element is built at the shear centre, where bending and torsion are decoupled; it is moved to the section centre by the same rigid arm \mathbf{s}: \mathbf{K}_c = \mathbf{T}_s^{\mathsf T}\mathbf{K}_{sc}\mathbf{T}_s. A bending–torsion coupling appears when the line of force does not pass through the shear centre.
• Shear coefficients \kappa_y, \kappa_z are specified separately in the element data (by default 5/6).
• Nonstructural mass per unit length (kg/m): added as a consistent translational mass.

The node stiffness is obtained as \mathbf{K}_n = \mathbf{A}^{\mathsf T}\mathbf{K}_c\mathbf{A}, where \mathbf{A} maps the node degrees of freedom to the degrees of freedom of the section centres. The same \mathbf{A} are applied to the mass, the geometric stiffness and the thermal force. With zero additional data the matrices coincide with the basic beam. The type is selected automatically: the extended type is enabled for nonzero roll, offsets, shear centre or nonstructural mass.

assign beam edges=12,14 section=rect w=10 h=30 flush=outer

7.14Beam hingesstable

A hinge at a beam end disables the transfer of one of the three moments — M_y, M_z or the torque T — in the beam's local axes. A hinge at either end and a hinge within the span are allowed.

The implementation uses no penalty springs. The beam end with a hinge receives its own twin node at the same coordinates: the three displacements of the twin equal the displacements of the main node, and the rotations are tied only along the unreleased axes by the equation \mathbf{a}\cdot(\boldsymbol{\theta}_d - \boldsymbol{\theta}_n) = 0. The released rotation remains a free degree of freedom of the twin, and its stiffness is provided entirely by the element — this is static condensation performed by the solver. A node at which all beams are released about one axis is fixed in that rotation; a moment applied to such a rotation has nothing to balance it, and the analysis refuses, naming the node.

7.15Beam section librarystable

A section is selected by type and dimensions; from them A, I_y, I_z, J and the extreme fibres c_y, c_z are computed (section moduli W_y = I_y/c_z, W_z = I_z/c_y). Types: rectangle, box section, solid circle, tube, I-beam, T-section, channel, angle. The kernel also provides direct input of A, I_y, I_z, J, c_y, c_z as a set of numbers.

Composite profiles are computed by summing rectangles with the parallel axis transfer, cutouts are rectangles with a negative sign:

The Saint-Venant torsion constant:

• circle and tube — the polar moment, J = \pi (D^4 - d^4)/32;
• solid rectangle with sides a \ge b — Roark's formula, J = a\,b^3\left[\tfrac13 - 0.21\tfrac{b}{a}\left(1 - \tfrac{b^4}{12a^4}\right)\right];
• box — a thin-walled closed contour after Bredt, J = 4A_m^2\,t / s, where A_m is the area enclosed by the midline, s is its perimeter, and t is the smaller wall thickness;
• I-beam, T-section, channel, angle — an open thin-walled profile, J = \tfrac13\sum b_i t_i^3.

For open profiles this is a first-order estimate: fillets and flange–web junctions are not accounted for in A, I and J.

From the internal forces the stress distribution over the section is built: \sigma = N/A \pm M_y z/I_y \pm M_z y/I_z; the shear stresses from shear forces by the Zhuravskii formula (the static moments are computed numerically on a grid of section points), from torsion by the formulas for a circle and a tube, by Bredt for a box and by the theory of thin-walled open profiles for the rest. The equivalent stress is the von Mises one.

7.16Mesh-based sections and the profile catalogueexperimental

A separate module for the properties of an arbitrary section, a polygon with cutouts, has been added to the kernel. The area, the centroid, I_y, I_z, I_{yz} and the principal axes are computed exactly from the contour (Green's theorem). Torsion and shear are found by a finite element analysis on a Tri6 mesh from two Poisson problems for the warping function \omega and the transverse shear function \Psi:

Here I_p is the polar moment, k_s is the energy shear coefficient (the shear energy of the Timoshenko beam is equated to the energy of the elastic solution), and \mathbf{g} is the shear field. The shear centre is also determined. The module returns the shear area k_s A: 5/6 for a narrow rectangle, exactly 1/2 for a thin circular tube. A catalogue of standard profiles has been added to the module (I-beams IPE, HEA, HEB, HEM per EN and I-beams and channels per GOST 8239 and GOST 8240).

7.17Corotational beamstable

For problems with large displacements and rotations the beam is recomputed in the corotational formulation of Battini–Pacoste (Battini J.-M., Pacoste C., CMAME 191, 2002). The configuration is defined by the node positions and the full node rotations \mathbf{R}_g; the rotation increment is the spatial spin \delta\mathbf{w}, and the update is multiplicative, \mathbf{R}_g \leftarrow \exp(\delta\mathbf{w})\,\mathbf{R}_g (addition of angles for finite rotations is incorrect).

The rigid frame of the element is built from the chord and the averaged y axis of the nodes. The deformational quantities in this frame are small for small strain: the elongation \bar u = l - l_0 and the rotations \bar{\boldsymbol\theta}_i = \log\left(\mathbf{R}_r^{\mathsf T}\mathbf{R}_{g,i}\mathbf{R}_0\right). The linear Timoshenko element from the previous subsections is applied to them unchanged, \mathbf{f}_l = \mathbf{K}_l\,[\bar u, \bar{\boldsymbol\theta}_1, \bar{\boldsymbol\theta}_2]. The internal forces in global axes are obtained by exact virtual work; the tangent matrix by numerical differentiation of the analytical internal forces with the same spin update as in Newton's method, and is therefore consistent with the iteration.

The basic beam and the extended one without offsets and shear centre are supported; otherwise the analysis refuses explicitly instead of silently substituting the physics.

7.18MITC4 shellstable

A four-node flat Reissner–Mindlin shell element for sheets, plates and thin-walled structures. Six degrees of freedom per node; the element is built in a local frame and consists of three parts: membrane, bending with transverse shear, and a penalty on rotation about the normal.

The local frame is built from the nodes: the normal is the cross product of the diagonals, the axis e_1 is directed along the edge 1\to2 after orthogonalization to the normal. The four nodes are projected onto the plane of the frame, warping is not modelled. The element matrix in global axes: \mathbf{K}_g = \mathbf{T}^{\mathsf T}\mathbf{K}_l\,\mathbf{T}. For the bending part, instead of the rotation vectors \theta_x, \theta_y, the normal tilt angles \beta_x = -\theta_y, \beta_y = \theta_x are used, so the rotational block of \mathbf{T} contains an additional permutation.

The membrane is a bilinear isoparametric Q4 in plane stress, full 2\times2 integration:

Bending. Curvatures and transverse shear strains:

The bending part is integrated by 2\times2 Gauss. For shear MITC4 is applied (Dvorkin E. N., Bathe K.-J., 1984; Bathe K.-J., Dvorkin E. N., 1985): the covariant shear components are computed only at four tying points, at the midpoints of the sides, and interpolated linearly:

where the tying points A,C lie at \xi=0, \eta=\mp1, and D,B at \xi=\mp1, \eta=0. The covariant component at a tying point, \gamma_\xi = \partial_\xi w - (\partial_\xi x\,\beta_x + \partial_\xi y\,\beta_y), is computed with the Jacobian at that very point. They are converted to the Cartesian [\gamma_{xz},\gamma_{yz}] at each integration point through \mathbf{J}^{-1}. The assumed shear strains remove shear locking as t/L \to 0. On a rectangular element the formulas coincide with the rectangular version to machine precision.

The stresses are computed at the element centre, \xi=\eta=0: the membrane \boldsymbol{\sigma}_m = \mathbf{C}_{ps}\boldsymbol{\varepsilon}_m plus the bending one at the top and bottom fibres, \sigma_b = \pm 6m/t^2, m = \mathbf{D}_b\boldsymbol{\kappa}. The output is \sigma_{xx}, \sigma_{yy}, \tau_{xy} and the plane-stress von Mises equivalent \sqrt{\sigma_{xx}^2 - \sigma_{xx}\sigma_{yy} + \sigma_{yy}^2 + 3\tau_{xy}^2} at each fibre; transverse shear stresses are not recovered by this route. The nodal load from pressure is consistent (\int N_i\,p\,dA with the full Jacobian); the pressure is directed against the outward normal of the element.

The offset of the reference surface z_0 (the parameter offset=bottom|middle|top or a number in the section shell command) is accounted for in the stiffness package: \mathbf{B}' = \mathbf{B} + z_0\mathbf{A}, \mathbf{D}' = \mathbf{D} + 2z_0\mathbf{B} + z_0^2\mathbf{A}; the element is then equivalent to a rigidly offset isotropic one (verified by a test).

units mm
section shell t=4
mesh shell size=40 element=quad pattern=structured

7.19Rotation about the normal (drilling)stable

The Q4 membrane and CST contain no rotational degree of freedom \theta_z about the normal, whereas a shell node has it and it enters the global system. Without an addition the matrix is singular in this degree of freedom. Tying \theta_z to zero (a grounding spring) would overestimate the bending stiffness of curved shells, so the Hughes–Brezzi penalty (Hughes T. J. R., Brezzi F., 1989) is applied: it ties \theta_z not to zero but to the rotation of the element's own displacement field:

Here \gamma is a dimensionless penalty parameter, chosen by tests and small compared with the membrane stiffness; it is not accessible in the interface or in scripts. The penalty introduces no noticeable energy in a state of pure membrane or pure bending, and ties the rotations of adjacent noncoplanar shells at a fold. For an element with a laminate, the membrane shear stiffness of the package is used.

7.20Triangular shell CST and DKTstable

A three-node flat element for unstructured surface meshes: triangles are built by any surface mesher without recombination. It consists of the CST membrane, DKT bending and the same penalty on rotation about the normal.

The membrane is a constant-strain simplex: \mathbf{B} is constant, \mathbf{K}_m = t\,A\,\mathbf{B}^{\mathsf T}\mathbf{C}_{ps}\mathbf{B}, one point. Bending is the Discrete Kirchhoff Triangle (Batoz J.-L., Bathe K.-J., Ho L.-W., 1980): the curvatures are expressed through the nine degrees of freedom [w,\theta_x,\theta_y] of the three nodes with the Kirchhoff hypothesis satisfied discretely at the midpoints of the sides, \boldsymbol{\kappa} = \mathbf{B}_{DKT}(\xi,\eta)\,\mathbf{u}. Integration is three points at the midpoints of the sides with weights A/3. There is no transverse shear in DKT, so there is no shear locking by construction.

Stress: the membrane stress is constant over the element, the bending stress is taken at the centroid, \xi=\eta=1/3, \sigma_b = \pm 6m/t^2 at the top and bottom fibre.

units mm
mesh shell size=20 surface=tri

7.21Laminated shell (ABD package)beta

A four-node MITC4 shell into which, instead of the isotropic E, \nu, t, a constitutive package per unit width is supplied:

The package is assembled ply by ply by classical laminate theory with first-order transverse shear (FSDT). For ply k with boundaries z_{k-1}, z_k (measured from the reference surface):

where \bar{\mathbf{Q}}_k is the reduced stiffness of the ply rotated by the layup angle (the angle is measured from the element axis x_e about the normal), and \bar{\mathbf{Q}}^s_k is the transverse shear stiffness of the ply. The shear coefficient is either uniform 5/6 or "energy-based" per ply, \kappa_k = \min\{1, \tfrac56\,G_{\max}/G_k\} (for sandwiches with a soft core \kappa \to 1). In all other respects the element coincides with MITC4: the shear with assumed strains uses the full matrix \mathbf{A}_s (2\times2), the membrane–bending coupling \mathbf{B} and the rotation penalty are added. The isotropic case is reproduced identically (regression at the 10^{-12} level).

In addition: ply-by-ply stresses in the material axes and failure criteria — maximum stress, Tsai–Hill, Tsai–Wu, Hashin (two-dimensional) and the delamination index from interlaminar shear stresses; the interlaminar stresses are recovered either from FSDT or from the equilibrium equations. The thermal forces \mathbf{N}_T, \mathbf{M}_T are computed from the CTE of the plies.

units mm
ortho name=ud E1=140e9 E2=10e9 nu12=0.3 G12=5e9 G13=5e9 G23=3.5e9
section shell laminate plies="ud:0.25:0,ud:0.25:45" repeat=2 symmetric=true shear=energetic
mesh shell size=40 element=quad pattern=structured

7.22Corotational shellsstable

For large displacements and rotations of shells the EICR corotational formulation (element-independent corotational; Rankin–Nour-Omid; Felippa C. A., Haugen B., CMAME 194, 2005) is applied on top of the linear MITC4 and CST+DKT elements. As with the corotational beam, the configuration is the node positions and full rotations, and the update is multiplicative.

The element frame is determined only by the node positions: for a quadrilateral the normal is from the diagonals, the axis e_1 is the projection of the vector between the midpoints of opposite sides; for a triangle it is along the edge 1\to2. The deformational degrees of freedom are the node displacements in the frame \bar{\mathbf{u}}_i = \mathbf{R}_r^{\mathsf T}(\mathbf{x}_i - \mathbf{x}_c) - \mathbf{R}_0^{\mathsf T}(\mathbf{X}_i - \mathbf{X}_c) and the rotations \bar{\boldsymbol\theta}_i = \log(\mathbf{R}_r^{\mathsf T}\mathbf{R}_{g,i}\mathbf{R}_0). The linear stiffness of the element is applied to them unchanged: \bar{\mathbf{f}} = \mathbf{K}_l\,\bar{\mathbf{d}}. The internal forces are the exact gradient of the energy \tfrac12\bar{\mathbf{d}}^{\mathsf T}\mathbf{K}_l\bar{\mathbf{d}} with respect to the global degrees of freedom (the projector term removes the work on rigid frame modes and centroid translation); the tangent is by differences of the analytical force, consistent with the spin update.

The code also supports the laminate type, but there is no separate verification case for it.

7.23Point massbeta

A single-node element with six degrees of freedom and no stiffness (\mathbf{K} = 0). Parameters: the mass m (kg) and three moments of inertia I_{xx}, I_{yy}, I_{zz} in global axes (in SI units). The mass matrix is diagonal:

It is used for attached equipment, engines and ballast in modal and dynamic analysis of frames and shell models. The point can be specified by a node or by coordinates.

units mm
mass node=1000,0,0 m=100 Ixx=0.5 Iyy=0.5 Izz=0.8

7.24Springbeta

An elastic connection between two regions of the model or between a region and a fixed ground. A region is selected by a face, an edge or a vertex.

Implementation: a master node is created at the centroid of each selected region; it is tied to the nodes of the region by a rigid region (RBE2); a two-node Timoshenko element with equivalent properties is placed between the master nodes. The axial stiffness k (N/m) is specified; E_{\text{eq}} = kL/A_{\text{eq}} is chosen at a fixed area, so that EA/L = k. The transverse stiffness k_t and the torsional k_\theta are optional: I_{\text{eq}} = k_t L^3 / (12 E_{\text{eq}}), J_{\text{eq}} = k_\theta L / G_{\text{eq}}; without them the bending and torsional stiffness are negligibly small. The "three axes" mode places three independent orthogonal springs k_x, k_y, k_z to the ground. A preload F_0 is possible, a self-equilibrated pair of forces along the spring axis. The postprocessor outputs the elongation \Delta L = (\mathbf{u}_B - \mathbf{u}_A)\cdot\mathbf{e} and the working force F = F_0 + k\,\Delta L.

spring surface number 3 ground k=2e6 preload=100

7.25Rigid elements and constraint equationsstable

A rigid region (RBE2) forces the dependent nodes to move together with the master as a single rigid body:

Both forces and moment are thereby transmitted, and the joining of nodes with three and six degrees of freedom (a solid with a beam) is solved without special devices.

A general constraint equation \sum_i c_i u_i = r (in the terms of other systems, CE/MPC) is also eliminated exactly: the first term is the slave degree of freedom, the rest are masters. A conflict of slave degrees of freedom and a cycle of constraints are detected and stop the analysis.

rigid link master vertex number 6 slaves surface number 2

7.26Mass matrices, geometric stiffness, stressesstable

Mass matrix. For solid elements in modal analysis the element mass \rho V_e is divided equally among its nodes (lumped, diagonal), translations only. For a regular Hex8 and Tet4 this coincides with HRZ lumping; for Tet10 and Hex20 the equal division gives the correct total mass, but the frequencies come out inexact (on bending modes of a P2 cantilever it nevertheless gives a good result). There is no consistent mass of solids in modal analysis. For beams and shells both matrices, lumped and consistent, are implemented in the kernel.

• Timoshenko beam: the consistent Przemieniecki mass with the bending rotary inertia \rho I and the torsional inertia \rho J; as \Phi \to 0 it is the Hermitian mass. This inertia is critical for deep beams and high modes (NAFEMS FV5).
• Shell (MITC4, CST+DKT): the consistent translational mass M_{ab} = \rho t\int N_aN_b\,dA, for a triangle in closed form \rho t A/12\cdot[2\ 1\ 1; 1\ 2\ 1; 1\ 1\ 2]; the bending rotary inertia is not accounted for in the modal mass.
• Implicit dynamics of solids: a choice of the consistent \rho\int N^{\mathsf T}N\,dV, the lumped HRZ (the diagonal of the consistent one, scaled to the total mass — always positive, unlike the row sum for serendipity elements) and "equally among nodes". For Hex27 the mass matrix is not implemented in dynamics.
• Beams and shells in implicit dynamics carry the translational mass and the rotary inertia of the node: for a beam \tfrac12\rho L\,\operatorname{diag}(I_y+I_z, I_y, I_z) per end, for a shell \rho t^3/12\cdot A_i, where A_i is the fraction of area belonging to the node. This formulation appeared after the last release.
• Explicit dynamics: a diagonal mass (HRZ lumping); for a regular Hex8 it is \rho V/8 per node.

Geometric stiffness \mathbf{K}_\sigma (for buckling and prestressed frequencies). For solids, from the actual stress tensor, a block for the three displacements:

For a beam, from the axial force N (the Przemieniecki matrix, the contribution N/L in the bending blocks). For shells, only from the membrane forces \mathbf{N} on all three translations, \int \nabla u_k^{\mathsf T}\mathbf{N}\nabla u_k\,dA; the contribution of the rotations is not accounted for. The sign is the same for all types: tension stabilizes.

Stresses at the element level. The basic recovery is by the formula \boldsymbol{\sigma} = \mathbf{C}\mathbf{B}\mathbf{u}_e at a characteristic point of the element (the centre for hexahedra and prisms, the centroid for tetrahedra, the axis at mid-height for a pyramid). The nodal field is averaging over adjacent elements; for Tet10 and hybrid meshes the stress is first computed at the nodes of each element and then averaged. For a beam, the end forces \mathbf{f} = \mathbf{k}\mathbf{u} and an estimate at the extreme fibres; for a shell, at the top and bottom fibres at the element centre. More accurate methods (smoothing, an error indicator) are described in the section Results, postprocessing and strength assessment.

7.27Thermal elementsbeta

The thermal problem is a scalar field T (one degree of freedom per node) on the same geometries: Tet4, Hex8, Wedge6, Pyr5 and the quadratic Tet10, Hex20, Hex27 with the same quadratures and shape functions. The element: K^e_{ij} = \int\kappa\,\nabla N_i\cdot\nabla N_j\,dV, the load from volumetric heat generation F^e_i = \int Q\,N_i\,dV. The faces of quadratic elements, Tri6 and Quad8, are not independent elements but faces for integrating boundary conditions: a six-node triangle with the six-point Dunavant rule of degree 4 and an eight-node quadrilateral with the 3\times3 Gauss rule; flux, convection, radiation and thermal contact are integrated over them. The formulations (steady-state, transient, nonlinear with radiation, thermal contact, thermoelasticity) are described in the section Heat transfer and multiphysics.

7.28Summary tablestable

In the table n is the number of nodes, DOF is the degrees of freedom of the element, "Mass" is which matrix is implemented in modal analysis. Details are in the subsections above.

What the library does not have: Pyr13 and Wedge15 as stiffness elements; quadratic shells; plane (strain, stress) and axisymmetric elements in the structural assembler — such problems are modelled by shells and a narrow wedge of solid elements with symmetry conditions, and plane elements exist only in explicit dynamics; beams with a warping degree of freedom; nonlinear and unilateral springs.

A short catalogue of elements with type codes is on the page /docs/elements.

This section covers analyses in which the structure remains linearly elastic, displacements and rotations are small, and the loading does not depend on them. They include statics, natural frequencies and mode shapes, eigenvalue buckling, steady-state response to harmonic excitation, response to a spectrum and to random vibration, as well as auxiliary techniques: submodeling, superelement reduction and rotor dynamics. Contact, plasticity, large displacements and time-domain transient dynamics are described in the section Nonlinear analysis and dynamics. Elements and their formulations are in the section Finite element library; methods for solving linear systems and eigenproblems are in the section Solvers.

General provisions for all analysis types:

• beams and shells have six degrees of freedom per node (three displacements and three rotations), solid elements have three displacements;
• constraints, constraint equations and rigid regions are eliminated from the system by the master–slave method, so eigenproblems and harmonic analysis use the same constraints as statics;
• the solver boundary uses SI units (m, N, Pa, kg), geometry is stored in millimetres, and the conversion is performed in one place (see Units of measurement);
• the mass matrix of solid elements is lumped: the element mass is divided equally between its nodes, including mid-side nodes for quadratic elements. For beams and shells, a consistent mass matrix is available in the modal analysis of frames;
• the native (desktop) and browser paths are built from the same formulation; the native one is considered the reference.

8.1Linear staticsstable

Statics determines the displacements of a structure under a constant load. Use it when the stresses lie in the elastic range, the displacements are small compared with the dimensions, and inertia can be neglected.

Formulation

The principle of virtual displacements for a linearly elastic body: for any admissible \delta\mathbf{u}

where \mathbf{C} is the elasticity tensor, \mathbf{b} is the body force (including that due to acceleration), \mathbf{t} is the surface load (pressure, distributed force), and \mathbf{F}_i are concentrated forces. After discretisation with the element shape functions, the system is

where \mathbf{B} is the strain–displacement matrix (shape function gradients). For beams and shells, the local stiffness matrices are rotated to the global axes, and all element types are assembled into a single system with six degrees of freedom per node.

Constraints and restraints

Restraints (including those with a prescribed nonzero displacement), constraint equations \sum_i c_i u_i = r and rigid regions are eliminated from the system by a linear transformation:

A restraint is a constraint with no master degrees of freedom (u = \text{const}). A rigid region drives the dependent node by rigid-body kinematics: \mathbf{u}_d = \mathbf{u}_m + \boldsymbol{\theta}_m \times (\mathbf{x}_d - \mathbf{x}_m). It transmits both forces and moment, so it is suitable for joining solid elements to beam and shell elements. Conflicting or cyclic constraints cause the analysis to fail with an error code rather than silently selecting one of them. Degrees of freedom with zero stiffness (for example, rotations at nodes that belong only to solid elements) are fixed automatically so that the matrix is not singular.

The reduced system is symmetric and positive definite if the structure is restrained against rigid-body motion. The solution method (direct factorisation or the preconditioned conjugate gradient method) is chosen according to the problem size; details are in the section Solvers.

What the user specifies and what is obtained

The user specifies materials, restraints, loads (forces, pressure, acceleration, moments, remote loads), constraints and contacts in the linear formulation. The output is nodal displacements, strains and stresses (components, von Mises equivalent, principal), support reactions, beam forces, and membrane and bending resultants in shells. The equivalent stress is computed as

For more on stress recovery and strength assessment, see the section Results, postprocessing and strength assessment.

In a script, statics is the default analysis type:

material E=200000000000 nu=0.3 density=7850
fix surface number 1
load surface number 2 force=5000 component=z
solve static

In the interface, the analysis is started at the "Solve" step, and the results are viewed at the results step.

8.2Modal analysisstable

Modal analysis finds the natural frequencies and mode shapes of free vibration of a restrained structure. It is used to avoid resonance, to assess the stiffness and mass of a structure, and as a basis for harmonic, spectrum and random analyses.

Formulation

The generalised eigenvalue problem for the reduced matrices:

The mode shapes \boldsymbol{\varphi}_r are mass-orthogonal: \boldsymbol{\varphi}_i^{\mathsf{T}} \mathbf{M} \boldsymbol{\varphi}_j = m_i\,\delta_{ij}. In the output, the mode shapes are scaled by the maximum displacement magnitude: \max \lvert \mathbf{u} \rvert = 1. This is only a convenient normalisation for display; the generalised mass m_r is not equal to one under it and is taken into account when computing participation (see response spectrum).

The mass matrix is assembled from the material density (solid elements, shells, beams), the non-structural mass of beams and laminate shells, and user-defined lumped masses. A lumped mass carries no rotational inertia, except where it is explicitly specified for the lumped mass.

Participation factors and effective masses

For a direction d (unit vector \mathbf{t}_d in the degrees of freedom), the following are computed:

The sum of effective masses over all modes tends to the mass of the free translational degrees of freedom (the mass at restrained nodes goes into the base). These quantities are output in spectrum analysis, where they determine the mass coverage; the panel of an ordinary modal analysis shows frequencies and mode shapes.

How it is implemented

The reduced matrices \mathbf{K}_r = \mathbf{T}^{\mathsf{T}}\mathbf{K}\mathbf{T} and \mathbf{M}_r = \mathbf{T}^{\mathsf{T}}\mathbf{M}\mathbf{T} are brought to good conditioning by diagonal scaling; the lowest modes are found by subspace iteration with a single stiffness factorisation (for frames, by the block LOBPCG method). Free structures (beams, shells): the lowest zero modes are returned as zero frequencies, and a spectral shift is applied for the factorisation. Details and the choice of method by problem size are in the section Solvers.

What the user specifies and what is obtained

The number of modes, restraints, materials with density, lumped masses, constraints. The output is frequencies in hertz and mode shapes (for viewing and animation). Density is mandatory: without mass there are no natural frequencies, and the analysis fails.

fix surface number 1
solve modal nmodes=6
probe frequency mode=1

8.3Modal analysis with prestressbeta

Load changes the natural frequencies: tension stiffens (frequencies rise), compression softens (frequencies fall and tend to zero near the stability limit). This analysis is needed for strings, blades, compressed columns and, in general, for structures operating under a large constant load.

Formulation

First, linear statics is solved under the applied load and the stresses \boldsymbol{\sigma}_0 are found. Then the geometric stiffness matrix (stress stiffness matrix) \mathbf{K}_\sigma is assembled from them and the problem

is solved. For a solid element, the block of the matrix \mathbf{K}_\sigma for a pair of nodes a, b is a scalar integral repeated for the three displacements:

The sign is such that tension gives a positive contribution and compression a negative one.

Why exactly this problem is solved rather than the formula f = f_0\sqrt{1 - P/P_{\text{cr}}}. The formula is a single-mode Rayleigh–Ritz reduction, exact only when the vibration mode coincides with the buckling mode. For a pinned beam this holds (0 % deviation); for a cantilever the error is up to 4 %; but for a pinned 2:1 plate at half the critical load the formula underestimates the frequency by 14 %, and at the critical load it reduces it to zero instead of the true 60 % of f_0. Therefore the product solves the full problem.

What the user specifies and what is obtained

Restraints, loads (surface, concentrated, body), the number of modes; in the interface, the "account for the loaded state" option, in a script, prestress=on. The output is frequencies and mode shapes under load.

The analysis fails loudly rather than returning a plausible number in the following cases: no density; nothing is restrained; the load is zero (there is no prestress, this is an ordinary modal analysis); the matrix \mathbf{K} + \mathbf{K}_\sigma has lost positive definiteness (the structure has already buckled, there are no natural frequencies); the eigenproblem has not converged.

fix surface number 1
load surface number 2 force=-1000 component=z
solve modal nmodes=3 prestress=on

8.4Linear bucklingstable

The eigenvalue analysis finds the load multiplier \lambda_{\text{cr}} at which an ideal structure loses stability (bifurcation), and the buckling mode shape. It is applied to columns, frames, shell panels and plates.

Formulation

Linear statics is solved under the reference load \mathbf{f}_{\text{ref}}, the internal forces are found from the displacements, and the geometric stiffness matrix \mathbf{K}_G = \mathbf{K}_\sigma(\mathbf{u}_{\text{ref}}) is assembled. The critical multiplier is the smallest positive \lambda of the problem

The critical load is \lambda_{\text{cr}}\,\mathbf{f}_{\text{ref}}. The entire model load is taken as the reference: surface, concentrated, body and acceleration loads. The multiplier applies to the whole load at once.

\mathbf{K}_G is assembled by element type:

• beam: through the axial force N (averaged over the element ends), using the consistent Przemieniecki geometric stiffness (block below);
• shell: through the membrane forces N = \sigma_m t (quadrilateral, triangular and laminate shells);
• solid element: through the stress tensor \boldsymbol{\sigma}, as in the previous subsection.

For a beam, the deflection–rotation block in the v–\theta_z plane (in the w–\theta_y plane the block is the mirror image with opposite coupling signs):

The sign convention is the same for all types (tension stabilises, compression destabilises), so mixed models are assembled into a single matrix.

Mode ordering

The requested number of smallest positive multipliers is returned in ascending order. Negative \lambda (buckling under reversed load direction) are not output: for them, reverse the sign of the load. Repeated values (for example, symmetric modes of a square plate) are output as separate modes. The mode shape is scaled by the maximum displacement.

How it is implemented

The matrix -\mathbf{K}_G is indefinite and of low rank, so direct inverse iteration on it diverges. A spectral shift and K-orthonormalisation of the subspace (for large problems) or a dense symmetric-pair solver (for small ones) are used. A check against an independent dense solver agreed to six digits. Details are in the section Solvers.

fix surface number 1
load surface number 6 force=-100 component=z
solve buckling nmodes=2
probe load_factor mode=1

8.5Harmonic responsestable

Harmonic analysis determines the steady-state response of a structure to a load varying harmonically over a given frequency range. The output is amplitude and phase frequency response functions and fields at a selected frequency. It is used to assess resonances and vibration amplitudes under vibratory loads.

Formulation

The load F\cos(\omega t + \varphi) is represented by the complex amplitude \hat{\mathbf{F}} = \mathbf{F}\,e^{i\varphi}, and the response \mathbf{u}(t) = \operatorname{Re}\left(\hat{\mathbf{u}}\,e^{i\omega t}\right) satisfies

The phase angle \varphi is set separately for each load (positive leads), so out-of-phase forces can be specified. Loads are grouped by phase, and the modal forces are complex. The output includes the amplitude \lvert \hat{\mathbf{u}} \rvert, the phase, the velocity amplitude \omega\lvert \hat{\mathbf{u}} \rvert and the acceleration amplitude \omega^{2}\lvert \hat{\mathbf{u}} \rvert, the model maximum envelope and the peak frequencies; the displacement field is available at a selected frequency and phase.

Damping models

One of three models is selected; they are mutually exclusive.

Modal. A single fraction of critical damping \zeta for all retained modes: \mathbf{C}_m = \operatorname{diag}(2\zeta\, m_j\,\omega_j).

Rayleigh. \mathbf{C} = \alpha\mathbf{M} + \beta\mathbf{K}. In the modal basis the matrix is exactly diagonal, and the fraction of critical damping of a mode is

The coefficients are specified directly (\alpha in 1/s, \beta in s) or from two "frequency — fraction of critical damping" pairs (f_1, \zeta_1), (f_2, \zeta_2):

A pair that gives a negative \alpha or \beta (that is, negative damping over part of the spectrum) is rejected. In Rayleigh mode, \zeta is not added separately.

Rayleigh per body. Each body of the assembly has its own pair (\alpha_r, \beta_r), \mathbf{C} = \sum_r (\alpha_r\mathbf{M}_r + \beta_r\mathbf{K}_r). The damping is non-proportional, and the modal matrix is full:

At each frequency a complex system whose size equals the number of modes is solved. The diagonal approximation is not used: it loses the coupling of modes across the body boundary. The results panel warns about zero damping, about a fraction of critical damping above one, and about strong non-proportionality.

There is no structural (hysteretic) damping.

Two solution methods

Modal superposition. The response is expanded in the retained mode shapes \hat{\mathbf{u}} = \sum_j \boldsymbol{\varphi}_j q_j, and at each frequency the system

is solved. For a constant \zeta this gives q_j = \hat{f}_j / \left[m_j(\omega_j^{2} - \omega^{2} + 2i\zeta\omega_j\omega)\right]. The method is fast but depends on the number of retained modes: the contribution of the truncated higher modes (the static correction) is not added. The modal force accounts only for the translational components of the load and mode shapes. The frequency grid is linear or logarithmic, with optional refinement around resonances.

Direct method (Direct FRF). Without truncating the basis, the full sparse complex system \left(\mathbf{K} - \omega^{2}\mathbf{M} + i\omega\mathbf{C}\right)\hat{\mathbf{u}} = \hat{\mathbf{F}} is solved at each frequency (sparse LU factorisation with COLAMD ordering; the structure is factorised once for the whole range). Rayleigh damping is specified per element: \mathbf{C} = \sum_e(\alpha_e\mathbf{M}_e + \beta_e\mathbf{K}_e), with a single common pair as a special case. All element types are supported, including beams and shells, as well as homogeneous constraint equations and rigid regions. At each frequency, the residual of the complex system is checked; if it is exceeded, or in case of singularity (undamped resonance), the analysis stops with an error and is not replaced by the modal one. Complex fields are additionally output: support reactions, stresses (linear components as complex amplitudes), beam forces, shell stresses. Nonlinear invariants such as the von Mises equivalent stress are not computed for complex amplitudes. The direct method is enabled explicitly (method=direct) or selected automatically for Rayleigh damping per body if the native solver is available. It requires Rayleigh damping (a constant \zeta is defined only in the modal basis) and is available only in the desktop application.

What the user specifies and what is obtained

The frequency range (f_{\min}, f_{\max}), the number of points, the number of modes for the modal method, the damping model, loads with phases. The output is the amplitude and phase frequency responses at selected nodes, the model envelope, peak frequencies, fields at a frequency, and for the direct method also reactions and stresses.

fix surface number 1
load surface number 2 force=-1000 component=z phase=0
solve harmonic fmin=380 fmax=440 npoints=61 nmodes=6 zeta=0.05
probe harmonic_amplitude x=200 freq=406

Direct method with Rayleigh damping specified by two pairs:

solve harmonic method=direct rayleigh f1=50 z1=0.02 f2=500 z2=0.05 fmin=10 fmax=1000 npoints=40

8.6Response spectrum (RSA)stable

Spectrum analysis estimates the maximum response of a structure to base motion specified by a design acceleration spectrum (seismic, shock and vibration spectra). The calculation is performed on top of a completed modal analysis and requires no time integration.

Formulation

The base moves with acceleration a_g(t) in direction d (vector \mathbf{t}_d):

For mode m with period T_m and the pseudo-acceleration from the spectrum S_a(T_m), the peak values (the participation \Gamma and the effective mass M_{\text{eff}} are as in modal analysis, formulas above) are:

where V_m is the modal base force. The mass is computed from the lumped mass matrix of the eigenproblem; the mode shape is multiplied by \Gamma, so the result does not depend on the normalisation and sign of the mode shape.

Modal combination

The peak contributions r_m of a single quantity (with signs!) are combined by one of the following rules:

• SRSS: R = \sqrt{\sum_m r_m^{2}}, for well-separated frequencies;
• CQC: R = \sqrt{\sum_m\sum_n \rho_{mn}\, r_m r_n}, for closely spaced frequencies, with the Der Kiureghian correlation coefficient (Der Kiureghian, 1981):

• ABS: R = \sum_m \lvert r_m \rvert, a conservative upper bound.

The fraction of critical damping \zeta in CQC sets the correlation; the ordinates of code spectra are specified for 5 % and are not rescaled with \zeta. Excitation directions are combined after the modes: SRSS over directions, R = \sqrt{\sum_d R_d^{2}}, or the 100–30–30 rule (ASCE 7-16, cl. 12.5.3): R = \max_d\left(R_d + 0.3\sum_{e\neq d}R_e\right).

Residual mass

The truncated modes give a rigid response to the zero-period acceleration (ZPA). The residual mass is M_{\text{res},d} = M_{\text{ref},d} - \sum_m M_{\text{eff},m,d}, where M_{\text{ref},d} is the mass of the unrestrained translational degrees of freedom. The correction to the base force V_{\text{res}} = M_{\text{res},d}\cdot \text{ZPA}_d is combined with the modal response by SRSS (NRC Regulatory Guide 1.92, ANSYS: missing mass). For displacement and stress fields, the residual part is found by a separate static solution \mathbf{K}\mathbf{u}_{\text{res}} = \text{ZPA}\,\mathbf{M}\,(\mathbf{t}_d - \sum_m\Gamma_{m,d}\boldsymbol{\varphi}_m) with the same constraints. The codes require at least 90 % mass coverage in each direction (SP 14.13330, ASCE 7-16, cl. 12.9.1.1): with lower coverage, the result is flagged with a warning that it is unsuitable for design, and the number of modes must be increased.

Design spectra

• SP 14.13330.2018: design seismic intensity of the site 7, 8, 9 on the intensity scale, soil category I–IV, importance factor K_0, permissible damage factor K_1, energy dissipation factor K_\psi; for the vertical component, a separate spectrum using the same formulas with K_\psi = 1 and a multiplier of 0.75 (cl. 5.12).
• ASCE 7-16 and 7-22: parameters S_s, S_1 (7-16, factors F_a, F_v from tables) or S_{DS}, S_{D1} (7-22), site class, period T_L, R, I_e.
• User-defined table of pairs (period or frequency, S_a) with linear or logarithmic interpolation, in fractions of g or in m/s².

Code parameters are not substituted silently: anything missing is named in an error before the modal analysis is run.

What is obtained

A table of modes (frequency, period, generalised mass, \Gamma, effective masses and fractions, cumulative fractions, S_a, modal base force), mass coverage by direction, base force by direction and the resultant, fields. Displacements are combined component-wise and are non-negative (an envelope); the sign of the displacement is lost in the combination, as is customary in RSA. Stress and strain components are combined using the sign of the modal contributions before SRSS or CQC is applied. The von Mises equivalent stress, the principal stress and the equivalent strain are not computed from the envelopes of components: they are built as maxima of the quadratic invariant over the ellipsoid of SRSS/CQC combinations ("response-spectrum ellipsoid peak"). This is the peak of the envelope, not the time maximum of a given accelerogram; such fields are available for the combined result with SRSS over directions and SRSS or CQC rules over modes.

fix surface number 1
solve spectrum code=sp14 intensity=8 soil=II dir=x,y combine=cqc dircombine=srss zeta=0.05 nmodes=20
probe spectrum_base_shear dir=x
probe spectrum_cumulative_mass dir=x

8.7Random vibration (PSD)experimental

Random vibration analysis estimates the root-mean-square (RMS) response to stationary Gaussian base excitation specified by the acceleration power spectral density (PSD, ASD), typical for transport, aviation and space engineering.

Formulation

The excitation is a one-sided density S_a(f) in units of g^2/Hz (or in SI units, m²/s⁴ per Hz) given by a table of "frequency — density" points with piecewise logarithmic (log–log) or linear interpolation. Between points (f_1, S_1), (f_2, S_2) with log–log interpolation, S(f) = S_1(f/f_1)^{k}, k = \ln(S_2/S_1)/\ln(f_2/f_1). The input RMS a_{\text{rms}} = \sqrt{\int S_a\,df} is integrated analytically over the segments.

The response is built by modal superposition. The transfer function of mode j for relative displacement and for absolute acceleration:

The covariance of the generalised coordinates is computed by direct integration over frequency (trapezoidal rule on a grid refined around resonances), without assuming well-separated frequencies:

The response variance at a node is \sigma^{2} = \sum_j\sum_k \phi_j\,\phi_k\,V_{jk}, where \phi_j is the value of mode shape j at the node. The CQC correlation matrix \rho_{jk} (the same formula as in the response spectrum) is computed and output as a diagnostic, but the variance itself uses V_{jk}. The output includes the RMS 1σ and the 3σ value (three times the RMS, a non-exceedance probability of 99.73 % for a Gaussian process) for relative displacements, absolute accelerations and von Mises equivalent stresses.

For a single-degree-of-freedom system with white noise, Miles' equation (Miles, 1954) is obtained:

(the second formula is for light damping).

How it is implemented and what is simplified

The frequencies, mode shapes and direction participation factors \Gamma (mass-based, as in the response spectrum) are taken from modal analysis; modal stresses are recovered from the mode shapes. The von Mises stress is estimated as \sigma_{\text{vM}}^{2} \approx \sum_j\sum_k s_j s_k V_{jk}, where s_j is the von Mises equivalent stress of the j-th mode shape. This is a simplification: the von Mises scalar has no sign and no cross-correlations of the tensor components, so the estimate is accurate for one dominant mode and approximate when several modes contribute noticeably. For displacements, the magnitude of the mode shape vector at the node is used, and the result is written to the component along the excitation direction. Fatigue and the distribution of maxima are not computed: 3σ is a linear scaling of the RMS, not a time peak.

fix surface number 1
solve psd psd=[(20,0.01),(80,0.04),(350,0.04),(2000,0.005)] dir=z zeta=0.02 modes=10
probe psd_disp_3sigma
probe psd_stress_1sigma

8.8Submodelingbeta

The cut boundary displacement interpolation method (Cut Boundary Displacement Method) makes it possible to obtain accurate stresses near a stress concentrator (hole, fillet) without refining the mesh of the entire part. It relies on Saint-Venant's principle: far from the concentrator, the displacement field of a coarse mesh is close to the exact one, even if the stress peak at the edge is underestimated by a factor of 1.5–2.

Formulation

1. The global model is solved on a moderate mesh: the field \mathbf{u}_g(\mathbf{x}).
2. A local submodel is cut out around the region of interest: the intersection of the body with a cutting box by a Boolean operation, followed by a fine mesh (linear or quadratic).
3. The submodel faces are divided into cut faces (passing inside the original body) and original free faces (hole, fillet, outer faces).
4. At each node of the cut faces (including the mid-side nodes of quadratic elements), a Dirichlet condition with the interpolated global displacement is imposed:

The interpolation is barycentric over the tetrahedra of the global mesh (linear or quadratic according to the order of the global elements); points outside the global mesh on a curved boundary are projected onto the nearest element.

1. The local system is solved by the same structural solver; free faces remain unloaded.

What the user specifies and what is obtained

In the interface, a submodel is launched in the "Submodeling" section of the results step (it calls the same command as the script). After the global solve, the following are specified: the centre at=x,y,z (by default, the stress maximum hotspot or the detected concentration hotspot), the size of the region radius= or box=dx,dy,dz, the submodel element size size=, and the order order=1|2. The output is the peak von Mises stress in the submodel outside the near-wall layer at the cut boundary, the number of nodes on the cut boundary and the number of elements (submodel_peak_vm, submodel_cut_node_count, submodel_element_count).

mesh tet size=12 order=1
fix surface number 1
load surface number 3 pressure=-100000000
solve static
submodel at=60,30,4 box=36,36,8 size=2.4 order=2
probe submodel_peak_vm

8.9Craig–Bampton reduction (CMS)experimental

The component mode synthesis method (Craig, Bampton, 1968) replaces a large subassembly with a superelement: boundary degrees of freedom and a small number of fixed-interface normal modes. It is applied to assemblies in which a subsystem is used repeatedly in dynamics.

Formulation

The degrees of freedom are divided into boundary b and interior i. Static (constraint) modes and fixed-interface normal modes:

The transformation and the reduced matrices:

The mass matrix \hat{\mathbf{M}} has an identity block for the modal coordinates. The full field is recovered by the inverse transformation.

What actually exists

Dense matrix operations are implemented: static and normal modes, reduced matrices, recovery of displacements and stresses for linear tetrahedra, and the reduced load vector. The module has been verified by automated tests on small matrices (Guyan limit, block-diagonality of \hat{\mathbf{K}}, agreement of the lowest eigenvalues). The cms command and the "Craig–Bampton superelement" section in the results accept a boundary surface and a number of modes, but in the current build the user's finite element mesh is not passed to them: the calculation runs on a built-in test chain of masses and springs. Therefore these commands cannot be used to reduce a real subassembly.

cms boundary="surface number 1" modes=6

8.10Rotor dynamicsexperimental

Rotor dynamics studies how rotation changes natural frequencies: the gyroscopic effect splits the bending frequencies into forward (FW) and backward (BW) whirl. The result is presented as a Campbell diagram (frequencies as a function of rotational speed) with synchronous excitation lines 1× and 2× and critical speeds at their intersections.

Formulation

In a reference frame rotating with angular velocity \Omega about the axis:

Here \mathbf{G} is the skew-symmetric matrix of Coriolis terms (\mathbf{G}^{\mathsf{T}} = -\mathbf{G}), and \mathbf{M}_{\text{spin}} is the softening due to centrifugal forces in the plane perpendicular to the axis. The problem is linearised in state space with the matrix

and solved as a full nonsymmetric eigenvalue problem (Householder reduction to Hessenberg form and the double-shift QR algorithm). From the pairs of complex conjugate eigenvalues, the frequency \operatorname{Im}\lambda and the fraction of critical damping -\operatorname{Re}\lambda/\lvert\lambda\rvert are taken. Critical speeds are the intersections \omega_i(\Omega) = k\,\Omega for k = 1, 2, found by linear interpolation over the speed grid.

What actually exists

This is not a full analysis of the user's rotor but a prototype:

• the matrix \mathbf{G} is built from the Coriolis terms of lumped masses; the gyroscopic moment of discs via the polar moment of inertia, bearings with stiffness and damping, and cross-coupling coefficients of the oil film do not participate in the calculation path;
• the whirl direction is determined by the sign of the frequency shift relative to the non-rotating rotor, and modes are matched by their order number in the frequency list, without tracking by mode shape, so at mode crossings the classification may be wrong;
• the solve rotor command and the "Rotor dynamics" section in the results compute the diagram for a built-in two-degree-of-freedom model (mass 2 kg, stiffness 2·10⁴ N/m, gyroscopic coefficient 0.4), not for the user's mesh; the axis, speed range, number of steps and number of modes parameters control only the sweep.

solve rotor axis=z speed_min=0 speed_max=1000 steps=20 modes=4
probe rotor_crit_speed_1

This section describes analyses in which the stiffness depends on the solution or in which inertia cannot be neglected: large displacements and rotations, loss of stability with passage through limit points, construction and excavation stages, implicit and explicit dynamics. Linear statics, modal analysis and linear buckling are described in the section Linear analysis types. Material models (plasticity, hyperelasticity, creep) are in the section Material models, contact is in the section Contact and connections, and linear solvers are in the section Solvers.

The units are the same as everywhere: geometry in millimetres, SI at the solver boundary (metres, newtons, pascals, kilograms, seconds). Results are returned in interface units. Time in dynamics is specified in seconds.

What is computed where:

An analysis that the platform does not support fails with a named reason. There is no silent substitution of a linear solution: in the nonlinear path, any unsupported element type, constraint or material raises an error with a message rather than a simplified answer.

9.1Geometric nonlinearity of solidsstable

It is enabled when the linear stiffness no longer describes the structure: displacements are comparable to the cross-section dimensions, an element rotates noticeably, and the pressure remains normal to the deforming surface. For solid bodies, the Total Lagrangian formulation is used: all integrals are taken over the initial configuration V_0.

Formulation

The displacement gradient and the deformation gradient at a Gauss point:

The gradient is built from displacements rather than from \nabla(\mathbf{X}+\mathbf{u}): when \lvert\mathbf{u}\rvert \ll \lvert\mathbf{X}\rvert, subtracting nearly equal numbers loses significant digits. For the same reason, the Green–Lagrange tensor is computed through \mathbf{H}:

The second Piola–Kirchhoff stress tensor \mathbf{S} is computed from the material law. By default this is the Saint Venant–Kirchhoff model:

For large deformations of rubber and elastomers, hyperelastic models are selected: the compressible neo-Hookean model with energy W = \tfrac{\mu}{2}(\bar I_1 - 3) + \tfrac{\kappa}{2}(J-1)^2, \bar I_1 = J^{-2/3} I_1, as well as Mooney–Rivlin, Yeoh, Arruda–Boyce and Ogden (parameters are in the section Material models). The formulation is purely displacement-based, without a separate pressure field.

The virtual work of internal forces \int_{V_0} \mathbf{S}:\delta\mathbf{E}\,dV_0 gives the internal force vector and the tangent matrix. In Voigt notation with engineering shear (order 11, 22, 33, 12, 23, 13):

The matrix \mathbf{B}_{NL} depends on \mathbf{F}: for node i and component m, the E_{11} row equals F_{m1}\,\partial N_i/\partial X, the \gamma_{12} row equals F_{m1}\,\partial N_i/\partial Y + F_{m2}\,\partial N_i/\partial X, and so on. The tangent is the sum of the material and geometric parts:

The tangent is analytical and symmetric (except for the contributions of follower loads and contact). For the isotropic Saint Venant–Kirchhoff model, \mathbb{D} is constant. For hyperelastic and elastoplastic laws, \mathbb{D} is the consistent tangent of the law, so Newton convergence is quadratic.

The quadratures are the same as in the linear assembly (section Finite element library): TET4 uses one point, TET10 four, HEX8 2\times2\times2 (full integration), HEX20 3\times3\times3. The kernel also supports prisms, pyramids and HEX27. The Solve button and DSL routes connect first- and second-order tetrahedra and hexahedra and the hybrid Hex8–Pyr5–Tet4 mesh; other combinations of types fail loudly.

What the user obtains

Displacements, support reactions (\mathbf{f}_{\mathrm{int}} - \lambda\,\mathbf{f}_{\mathrm{ext}} in the final configuration), Green–Lagrange strains, true Cauchy stresses and, from them, von Mises equivalent stresses and principal stresses. The Cauchy stresses are obtained by a push-forward:

The element value is the average over the volume V_0 across the Gauss points; the nodal value is the average over the adjacent elements weighted by their volume V_0. This is not a superconvergent recovery: in gradient regions the nodal peaks are smoothed. For the accepted load steps, complete fields are stored, and the postprocessor can switch between steps. The snapshot memory is limited by a budget; when it is exceeded, steps are thinned uniformly, keeping the first and the last.

The settings are shared by the Solution panel and DSL: the number of load steps, the iteration cap per step, the residual tolerance, the pressure mode (follower or dead), and the contact penalty multiplier. Follower pressure acts along the normal of the current configuration: \mathbf{f}_a = -\lambda p \int N_a\,\tilde{\mathbf{n}}\,d\xi\,d\eta, where \tilde{\mathbf{n}} = \mathbf{x}_{,\xi}\times\mathbf{x}_{,\eta}. The tangent of this load is accounted for and is in general nonsymmetric (the load on a single face is nonconservative).

nlgeom on pressure=dead steps=10 maxiter=25 tol=1e-6
solve static

Low-order locking at large displacements was measured on the bending of a cantilever with L/h = 20 at PL^2/EI = 10 (40 elements along the length, one through the thickness). The reference is the Euler elastica, integrated independently. The maximum error of the tip displacement:

This is locking on a specific mesh, not a guarantee for arbitrary geometry. For bending at large displacements, use second-order elements.

9.2Corotational beams and shellsstable

For beams and shells, the full Lagrangian formulation is excessive: strains remain small, while rotations may be large. A corotational approach is used here: the element frame follows the nodes, a linear element is applied in it, and the geometric nonlinearity comes entirely from the rotation of the frame.

The nodal configuration is the position \mathbf{x} = \mathbf{X} + \mathbf{u} and the total rotation \mathbf{R}_g\in SO(3). The rotation increment is the spatial spin \delta\mathbf{w}, and the update is multiplicative:

Angles cannot be added: at large rotations this is incorrect. As a result, the nodal rotation is output as the vector \log\mathbf{R}_g with an angle in the range [0, \pi].

Beam

The formulation of Battini and Pacoste (2002). The element frame \mathbf{R}_r is built from the chord and the mean position of the y axis of the nodal rotations. The deformational quantities in the frame:

They are small at small deformation, so the linear Timoshenko element is applied to them unchanged: \bar{\mathbf{f}}_l = \mathbf{K}_l\,[\bar u,\bar{\boldsymbol\theta}_1,\bar{\boldsymbol\theta}_2]^{T}. The internal forces are obtained by exact virtual work with account of the derivative of the rotation; the tangent is analytical (agreement with finite differences to 10^{-9}).

Shell

The element-independent corotational formulation EICR (Rankin and Nour-Omid; Felippa and Haugen, CMAME 194, 2005) on top of linear shells: a quadrilateral with a Q4 membrane, MITC4 bending and a fictitious drilling stiffness about the normal; a triangle is a CST with DKT bending. The element frame depends only on the nodal positions. The forces are the exact gradient of the energy \tfrac12\bar{\mathbf{d}}^{T}\mathbf{K}_l\bar{\mathbf{d}} with respect to the global degrees of freedom; the projector term removes the work on the rigid modes of the frame. The tangent is analytical.

What the user obtains

Nodal displacements and rotations (rotations as the total rotation vector), reactions, beam forces N, V_y, V_z, T, M_y, M_z at the ends in the current element frame. For shells, the von Mises equivalent stress is computed at the top and bottom fibres at the element centre, from linear relations with the deformational quantities of the frame; the larger of the two is output. There are no through-thickness stress fields in the nonlinear analysis.

Nodal moments are specified by axial vectors fixed in space and conjugate to the spatial spin. In planar problems such a load is conservative; follower moments are not available in three-dimensional problems.

frame solve static nlgeom=on steps=10 maxiter=30 tol=1e-6

9.3Newton–Raphson iterations and load stepsstable

The task is to find \mathbf{u} at which the internal forces balance those applied at the load fraction \lambda. The residual on the free degrees of freedom:

The load increases from \lambda = 0 to \lambda = 1 in steps. At each step the Newton–Raphson method is applied, with the tangent fully recomputed at every iteration:

Rotations are updated multiplicatively (see the previous subsection). The tangent is factorised by sparse LU (the Eigen library), and the symbolic analysis of the sparsity pattern is reused between iterations. LU is used rather than Cholesky because the tangent of corotational elements is nonsymmetric away from equilibrium; for more on linear solvers, see the section Solvers.

Convergence criterion

A step is accepted when the relative residual does not exceed the tolerance:

Here f_{\mathrm{ref}} is the reference load norm (below). The criterion is based on the residual only; there is no criterion on the displacement increment or on energy. The reference norm accounts for follower loads (otherwise a purely follower problem with zero external force in equilibrium would give a relative residual of 1) and, when loading by prescribed displacements or unloading to \lambda = 0, a lower bound from the full model load. An iteration whose residual is already small before the first solve is not counted: at least one increment is performed.

Automatic step cutting

The initial step \Delta\lambda is the fraction 1/N of the specified number of steps. If Newton has not converged within the allotted number of iterations (or the residual is not a number), the step is discarded, the state is rolled back to the last accepted one and \Delta\lambda is halved. If Newton converged quickly, the step grows. Growth and cutting are bounded below and above. Cutting below the minimum is a failure, "auto-step exhausted". All accepted steps and rollbacks are counted and shown in the result.

Load path

For elastoplastic problems and cycles, the load can follow an arbitrary sequence of target fractions \lambda: \{1, 0\} is loading and full unloading (residual deformation), \{1, -1, 1\} is a cycle. The auto-step works on each segment. The load path is supported only under load control.

What the user specifies

The failure reasons visible to the user are fixed: auto-step exhausted (equilibrium not found); step cap reached; singular tangent (mechanism or limit point under load control); landing on \lambda=1 did not converge; invalid input; the path ran onto a nonphysical branch (arc-length). Together with the reason, the achieved \lambda and the number of steps are reported.

9.4Crisfield arc-length methodstable

When the load path contains load limit points (snap-through) or snap-back, load control stops at them: the tangent is singular. The arc-length method advances along the equilibrium curve itself rather than along \lambda.

The unknowns are the increments \Delta\mathbf{u} and \Delta\lambda with an additional constraint equation. Crisfield's cylindrical variant is implemented: the arc-length constraint is imposed only on the displacements of the free degrees of freedom,

The iterative corrector is decomposed as \delta\mathbf{u} = \delta\mathbf{u}_r + \delta\lambda\,\delta\mathbf{u}_t, where \delta\mathbf{u}_r = \mathbf{K}_T^{-1}\mathbf{r}, \delta\mathbf{u}_t = \mathbf{K}_T^{-1}\mathbf{f}_{\mathrm{ref}}, and \delta\lambda is a root of the quadratic equation

Of the two roots, the one is taken for which the updated increment makes the smaller angle with the increment before correction. Complex roots cause the step to be cut.

The predictor is tangent; the sign of the load increment is determined by Feng's rule: by the sign of the scalar product of the previous increment and the tangent vector. This is insufficient for snap-back in many degrees of freedom, so Crisfield's reversal check is added: the accepted step increment must form an acute angle with the previous one, otherwise the step is rejected and the arc is halved. The arc length is adjusted according to the number of corrector iterations. If an arc step has crossed \lambda = 1, it is rolled back and "lands" on \lambda = 1 by load control in small substeps. If landing is impossible (a branch with negative stiffness), the last converged arc step is accepted.

As a result, the user obtains the full history \lambda(u) along the path: with segments where \lambda decreases as u grows and with \lambda<0 after snap-back. The method is selected in the Solve step (the "Method" list) and in DSL (method=arc).

frame solve static nlgeom=on method=arc steps=40 maxiter=40 tol=1e-6

9.5Initial imperfections and postbuckling behaviourbeta

A perfect structure under compression loses stability by bifurcation; a real one, by a limit load that depends on the shape and amplitude of the initial deviations (imperfection sensitivity, Koiter). A typical chain: linear buckling analysis \to introducing an imperfection along the computed mode shape \to nonlinear arc-length analysis.

The imperfection command takes the k-th eigenmode of linear buckling \boldsymbol\phi_k, normalises it by the largest translational displacement and shifts the mesh node coordinates:

Here w_0 is the largest shift in model length units (the scale parameter). Several consecutive commands add up the mode shapes (a multimode imperfection). It is applicable to shell meshes, solid meshes and frames. The command requires a preceding solve buckling on the same model and changes the geometry of the working model: subsequent analyses run on the perturbed mesh.

solve buckling nmodes=3
imperfection mode=1 scale=0.5
solve static nlgeom=on method=arc steps=20 maxiter=30 tol=1e-5

9.6Material nonlinearity in the iterationsstable

Material nonlinearity enters the same Newton loop. At each Gauss point the history of the law (the accepted state) is stored; at an iteration, the law is evaluated from the accepted state rather than accumulated over iterations. The accepted state is updated after a converged step; when a step is rolled back, the state remains unchanged. Therefore repeated assembly at the same displacements is idempotent, and step cutting does not corrupt the history.

The law returns the stress and the consistent tangent \partial\mathbf{S}/\partial\mathbf{E}. For von Mises plasticity (isotropic, kinematic, multilinear hardening) and Johnson–Cook, the Green–Lagrange tensor is the input of the law and the second Piola–Kirchhoff tensor is the output: plasticity is additive in this pair. The description of the laws and their parameters is in the section Material models.

The "material nonlinearity only" mode (nlgeom=material) disables the geometric nonlinearity of solids: \boldsymbol\varepsilon = \mathrm{sym}\,\nabla\mathbf{u}, \mathbf{K}_{\mathrm{geo}} is not assembled, and the pressure is dead.

For perfect plasticity, exhaustion of the auto-step at the last load target means exhaustion of the load-carrying capacity. In this case the product returns the last converged equilibrium as the result and reports the limit load fraction (\lambda_{\mathrm{lim}}); with hardening, such a failure remains an analysis error.

material steel plasticity=j2 yield=250e6
nlgeom material steps=20 loadpath="1,0"
solve static

9.7Prescribed displacements in the nonlinear formulationbeta

Loading can be applied not only by force but also by a prescribed displacement (an actuator, an interference fit, specimen stretching). In nonlinear statics, a nonzero prescribed displacement is introduced by a linear constraint equation on the translational degrees of freedom with a right-hand side proportional to \lambda:

Here c_i are the constraint coefficients (for a prescribed displacement, one degree of freedom with coefficient 1), \bar r is the prescribed value (in metres at the solver boundary); the first term of the equation is the dependent degree of freedom. The displacement grows together with the load fraction and returns to zero in the load paths \{1, 0\}. A constraint on translations remains exact at large displacements. Prescribed displacements are compatible with forces, follower pressure and the load path.

When loading by displacements alone, there is no external force, and the usual reference residual norm is zero: the round-off residual would be divided by zero and the analysis would falsely fail to converge. Therefore the reference norm is taken from the energy scale of the prescribed displacements. This change does not affect analyses loaded by forces.

9.8Implicit dynamics of solidsstable

Transient dynamics with inertia and damping, when the load varies in time quickly compared with the period of the fundamental mode but impacts and stress waves are not needed (for those, use the explicit scheme). The semi-discrete equation of motion is solved:

with the load \mathbf{F}(t) = k(t)\,\mathbf{F}_0, where k(t) is a piecewise linear amplitude from a table of "time — multiplier" pairs.

Schemes

Newmark-\beta and HHT-\alpha (Hilber–Hughes–Taylor) are implemented. The generalised-\alpha scheme (Chung–Hulbert) is not available. Both schemes are written in one form:

Newmark corresponds to \alpha = 0; by default \beta = 1/4, \gamma = 1/2 (the trapezoidal rule: unconditionally stable, without numerical damping). For HHT-\alpha the parameter is \alpha\in[-1/3, 0], with \beta = (1-\alpha)^2/4, \gamma = (1-2\alpha)/2; the numerical damping of high frequencies is stronger the smaller \alpha is: the limit of the spectral radius as \Omega\to\infty is (1+\alpha)/(1-\alpha). General \beta, \gamma are accepted for \gamma \ge 1/2 and \beta \ge (\gamma + 1/2)^2/4, that is, only unconditionally stable ones; the others (including "linear acceleration") are refused to avoid silent instability.

The unknown of the scheme is the step increment \Delta\mathbf{u}. Newton uses the effective tangent

Rayleigh damping is \mathbf{C} = a_R\mathbf{M} + b_R\mathbf{K}_0; \mathbf{K}_0 is the tangent of the initial configuration and is not updated during the analysis (like material damping in commercial codes). The coefficients are specified as the pair (a_R, b_R) or as two "frequency — fraction of critical damping" pairs. If a step has not converged, it is halved recursively to a given depth (substeps within a single output step). For small strains without plasticity (\mathbf{K} constant), a fast linear path with a single factorisation for the whole analysis is used.

What the user specifies

The step h and the end time are mandatory. Mass: consistent (default), HRZ lumped, or "equally among nodes". Kinematics: small strains or large displacements (Total Lagrangian). Initial conditions: statics under the full load ("release") or zeros; an initial velocity for all nodes or for a selected body. Base acceleration is specified in the relative formulation \mathbf{F} = -\mathbf{M}\boldsymbol\iota\,\mathbf{a}_0\,k(t). The output is displacements, velocities and accelerations at probe points versus time, energy histories (kinetic, internal work, dissipation by damping, work of external forces, balance), the maximum von Mises stress and plastic strain, field frames and the von Mises envelope over the nodes.

solve transient dt=1e-4 tend=0.05 scheme=hht alpha=-0.05 amplitude="0:0,0.001:1"
solve transient dt=1e-4 tend=0.05 rayleigh f1=50 z1=0.02 f2=500 z2=0.05

The keys alpha= and beta= here are scheme parameters, not Rayleigh coefficients. An alternative to direct integration is modal superposition (method=modal) with the exact Duhamel integral and a static correction for the residual modes; it is linear and does not replace nonlinear dynamics.

9.9Implicit dynamics of beams and shells, constraints and basebeta

For beams and shells, dynamics is carried out on the rotation group SO(3). The node carries rotational inertia; the Euler equation for rotation in spatial form is

The step rotation increment \Delta\boldsymbol\theta = \log(\mathbf{R}\,\mathbf{R}_n^{T}) enters Newmark/HHT, and the tangent of the rotational block includes the derivative of \mathbf{T}^{-1}. Mass: consistent or lumped; an initial angular velocity \mathbf{w}_0 is supported (rigid rotation about the centre of mass).

For solids, the following are additionally available: CE/MPC and RBE2 constraints, normal constraints, remote forces and moments, prescribed support motion (u = k(t)\bar u in the absolute formulation), and base acceleration (seismic excitation). Dynamics runs in the independent degrees of freedom with the reduced mass \mathbf{M}_r = \mathbf{T}^{T}\mathbf{M}\mathbf{T}. Contact and friction (node-to-surface penalty, regularised Coulomb friction) are described in the section Contact and connections.

9.10Explicit dynamicsstable

The explicit scheme is needed for fast processes: impact, penetration, blast loading, wave propagation. The time step is limited by the Courant–Friedrichs–Lewy condition, but there is no solution of a system of equations and no iterations: one evaluation of internal forces per step. The equations of motion are integrated by the central difference scheme (velocity Verlet form) with a diagonal mass matrix:

HRZ lumped mass (Hinton–Rock–Zienkiewicz): the diagonal of the consistent mass, scaled to the full element mass. It is positive also for second-order elements, where row summation gives negative masses at corner nodes.

Critical time step

The step does not exceed the estimate of the wave travel time across the smallest element:

Here L_e is the characteristic element size (for a hexahedron V/A_{\max}, for a tetrahedron 3V/A_{\max}, where A_{\max} is the largest face), and \alpha_{\Delta t} is the safety factor set by the user (dtScale). For HEX8, the strict element limit 2/\omega_{\max} from the maximum eigenvalue of \mathbf{M}^{-1}\mathbf{K}_e is used. For a cube at \nu=0.3 it equals 0.735\,L_e/c (volumetric "breathing"), so the estimate L_e/c by itself is not safe. The step is recomputed during the analysis: over the live elements and over the current size at finite strains. If the mesh contains degenerate elements with an anomalously small characteristic size, their mass is scaled so that they do not reduce the step of the entire analysis to zero. The scaling distorts the inertia of these elements.

Elements, materials, loads

Three-dimensional TET4, TET10, HEX8; two-dimensional QUAD4 and QUAD8 (plane strain or plane stress). HEX8 is a one-point element with an exact volume-averaged gradient and hourglass control after Flanagan–Belytschko (viscous and stiffness forms, with the coefficient set by the user; the analysis reports the ratio of hourglass energy to internal energy and warns if it is excessive). Materials: linear elastic, J2 with perfect, linear or multilinear hardening (isotropic and kinematic fractions), Johnson–Cook (rate, temperature, adiabatic heating), the Mie–Grüneisen equation of state for shock waves (see the section Material models). The bulk viscosity q_1 (linear) and q_2 (quadratic) suppresses oscillations at fronts.

Loads and conditions: initial velocity of a body or surface, restraint, symmetry, prescribed surface motion (by velocity or displacement from a table), rigid wall, body force and gravity, absorbing boundary. The analysis maintains an energy balance: kinetic and internal energy, hourglass energy, erosion energy, work of external actions.

Finite strains

For TET4 and HEX8, the finite strain mode is enabled (largeStrain=on): updated Lagrangian formulation, gradients on the current mesh, Jaumann hypoelasticity with rotation of accumulated strains after Hughes–Winget, and the current volume and element size in the step estimate. Without it the formulation is geometrically linear (small strains and rotations). For impacts with strains of order one, this mode is necessary: the Taylor impact in the geometrically linear formulation deviates from theory by 1.6 % and 5.1 % (at a dimensionless load \rho U^2/\sigma_y of 0.25 and 0.5).

material E=124e9 nu=0.34 density=8960 plasticity=j2 yield=300e6
explicit endTime=0.00017 dtScale=0.8 largeStrain=on
bc explicit_velocity body vx=-91.5
bc rigid_wall plane x 0 friction=0
solve explicit

9.11Element erosion and failurebeta

For penetration and spall analyses, an element that reaches the failure criterion is removed from the internal forces. Criteria can be combined; failure occurs when any is reached:

• limit equivalent plastic strain \varepsilon_p \ge \varepsilon_{\mathrm{fail}};
• accumulated Johnson–Cook damage D = \sum \Delta\varepsilon_p/\varepsilon_f, where \varepsilon_f = [d_1 + d_2 e^{d_3\eta}]\,[1 + d_4\ln\dot\varepsilon^*]\,[1 + d_5 T^*];
• maximum principal stress \sigma_1 \ge \sigma_{\max};
• spall pressure \sigma_m \ge p_{\mathrm{spall}}.

A removed element stops carrying forces, but the mass of its nodes remains (momentum is conserved), the step is computed over the live elements, the internal energy of the element is recorded as erosion energy, and the energy balance remains closed. The eroding contact is a one-sided node-to-segment penalty over the entire outer surface of the live elements; as erosion proceeds, new faces open up and the projectile meets the cavity that has formed.

material steel plasticity=j2 yield=355e6 fail="ep:0.25,smax:1.2e9"
explicit endTime=0.03 dtScale=0.8 largeStrain=on erodingContact=on

9.12Wave problems in explicit dynamicsbeta

Propagation of elastic waves in soil and solids. Sources: pressure or shear load on a face, a force at a node, a blast cavity with pressure as a function of time, a smooth expansion centre (irrotational body force, a pure P-wave). Time laws: Ricker pulse, exponential, frequency sweep, harmonic, step, table. Elastic layers define a layered half-space. Virtual geophones record displacement, velocity and acceleration at a node or at the node nearest to a given point; PPV, arrival time and wave speed are computed from them. The fields \nabla\cdot\mathbf{u} and \lvert\nabla\times\mathbf{u}\rvert are output to the grid and to the nodes.

The Lysmer–Kuhlemeyer non-reflecting boundary imposes stresses proportional to velocity:

The density and wave speeds are taken from the material of the adjacent element, the normal from the face itself, and the damper enters the step implicitly.

explicit geophone at=0,0,-100 name=g1
explicit layer axis=z at=-300 E=2e9 nu=0.3 density=1800
bc explicit_absorbing surface number 4

9.13Element birth and death, construction and excavation stagesstable

For tier-by-tier construction of a structure, soil excavation, tunnelling and layer-by-layer deposition, the set of active elements changes from stage to stage, loads and restraints are cumulative, and displacements and stresses accumulate. A list of stages is specified: each stage activates bodies (birth) and deactivates bodies (death) at its start.

The "soft element" method: the mesh remains whole, and the size of the system does not change between stages. An inactive element is not deleted; its internal force and tangent are multiplied by a small softness factor s (the stagesoft parameter), the law is elastic and without history, and density and weight are not assigned to it. The error of the algorithm is proportional to the factor: the relative error is of the order of 0.04\,s. Nodes belonging only to inactive elements follow the active part: at the start of a stage their position is set from the nearest boundary node of the active body.

Birth. At the start of a stage, the reference configuration of a new element is taken equal to the current one: \mathbf{H} = \nabla(\mathbf{u} - \mathbf{u}_{\mathrm{ref}}) = \mathbf{0}, \mathbf{S} = \mathbf{0}, and the plasticity history is clean. The element is born stress-free on the deformed configuration rather than on the original mesh. Death. The element becomes soft and its residual force is set to zero; the load it carried is taken over by the structure gradually together with the load increment of the stage.

Stage loading proceeds along the parameter \lambda\in[0,1]:

where \mathbf{b} is the internal force for the new set of active elements on the configuration \mathbf{u}_0 at the start of the stage (in equilibrium by construction), and \mathbf{f}_{\mathrm{stage}} is the full load of the stage. The difference \mathbf{f}_{\mathrm{stage}} - \mathbf{b} carries both the increment of the external load and the "excavation force" of the dead elements. Newton and the auto-step are shared with nonlinear statics. The stage tolerance is measured against the full model load. Weight and body forces are assigned per body: an inactive body has no weight.

Results are stored per stage; the stage view command switches the view and measurements.

nlgeom on pressure=dead
stage add name=L1
stage add name=L2 activate=body2
solve static
stage view 2

9.14Coupling T(t) → u(t)beta

Transient temperature can load a structure in dynamics. The coupling is one-way: the heat conduction step gives T_{n+1}, and the thermal load \mathbf{f}_{\mathrm{th}}(T_{n+1}) enters the mechanics; there is no reverse effect of deformation on temperature (the Gough–Joule term). There are two modes: quasi-statics (at each thermal step \mathbf{K}\mathbf{u} = \mathbf{f}_{\mathrm{ext}} + \mathbf{f}_{\mathrm{th}}(T), without inertia) and dynamics, with the same Newmark/HHT integrator as in implicit dynamics and the thermal load on the right-hand side; HHT takes (1+\alpha)\mathbf{f}_{\mathrm{th}}(T_{n+1}) - \alpha\mathbf{f}_{\mathrm{th}}(T_n). One thermal step corresponds to a specified number of mechanical steps. The formulation of transient heat conduction is in the section Heat transfer and multiphysics.

9.15Cancelling an analysisbeta

The "Cancel" button in the Solve step sets the cancellation flag of the native solver and interrupts waiting for the result. What exactly stops depends on the analysis.

The flag is polled: in the explicit dynamics loop at every step (the analysis returns a partial result with the status "cancelled"), in transient heat conduction after every step (partial frames), and in the load steps of mortar contact. In the Newton loop of nonlinear statics, in implicit dynamics and in the staged analysis there is no polling point: the interface shows the cancellation immediately, but the kernel may run to the end of the current call.

9.16Honesty of the "converged" flagstable

A nonlinear analysis that shows a result disguised as a solution is worse than one that fails. The product follows these rules.

For load control, an analysis is considered converged only if \lambda = 1 has been reached (or the last target of the load path) and the last step has converged on the residual. An exhausted auto-step, the step cap, a singular tangent, or non-numeric values result in a refusal with a reason code, the achieved \lambda and the number of steps, rather than the last step disguised as a result. In the interface, the status reads "Converged: \lambda = 1 reached" or "Not converged: reason". Check: a problem without equilibrium returns converged = false and \lambda < 1.

Partial progress is not masked: the user sees the number of accepted steps, rollbacks and iterations. For perfect plasticity, the limit load is a separate explicit result, marked as the limit load fraction. An inverted element (with a non-positive Jacobian) is considered non-numeric and cuts the step: an analysis under a pressure exceeding the limit one does not return a "converged" solution with negative volume. Native results are checked for completeness: an old binary that does not know plasticity or the load path fails with a request to update the application rather than returning an elastic answer.

This section describes thermal problems (heat conduction, radiation, thermal contact), thermoelasticity, flow of a liquid or gas around bodies (CFD) and cavity acoustics. The modules differ in maturity: heat conduction and thermoelasticity run on both desktop and browser, CFD on desktop only; acoustics and the Gough–Joule effect are so far experimental. The status of each subsection is given in its heading, and limitations are set out in separate callouts.

Units at input and output: geometry in mm, time in seconds, temperature in °C (in the radiation formulas the kernel converts it to kelvin itself), thermal conductivity \kappa in W/(m·K), specific heat c_p in J/(kg·K), density in kg/m³, heat flux density in W/m², heat transfer coefficient h in W/(m²·K), volumetric heat generation Q in W/m³, coefficient of thermal expansion \alpha in 1/K. SI is used at the solver boundary.

The desktop build computes with the native kernel, the browser build with the same code in WebAssembly. The desktop is considered the reference. On the verified problems (nonlinear heat conduction, coupled thermoelastic dynamics) the numbers of the two platforms agree to the printed digits. CFD is available only on desktop.

10.1Steady-state heat conductionbeta

The steady-state problem gives the established temperature field for prescribed temperatures, fluxes, convection and volumetric heat generation. It is used to estimate temperature differences in parts, heat sink power, and as the first step of coupled thermoelasticity.

Formulation

The Fourier equation for an isotropic material in the domain \Omega:

Boundary conditions on non-overlapping parts of the boundary (the normal \mathbf n is outward):

• temperature: T=T_D;
• flux, positive into the body: \kappa\,\partial_n T=q;
• convection: -\kappa\,\partial_n T=h\,(T-T_\infty);
• radiation into the environment: -\kappa\,\partial_n T=\varepsilon\sigma\,(T_K^4-T_{\infty K}^4) — see nonlinear heat conduction.

The weak form for the linear case (without radiation): find T equal to T_D on the prescribed-temperature boundary such that for any test function w

Discretisation gives a symmetric positive definite system \mathbf K\,\mathbf T=\mathbf F with the matrix K_{ij}=\int_\Omega\kappa\,\nabla N_i\cdot\nabla N_j\,dV+\int_{\Gamma_h}h\,N_iN_j\,dA and the vector F_i=\int_\Omega Q\,N_i\,dV+\int_{\Gamma_q}q\,N_i\,dA+\int_{\Gamma_h}h\,T_\infty N_i\,dA.

What the user specifies and what is obtained

The user specifies the material \kappa (separately for each body), temperatures, fluxes and convection coefficients on faces by surface numbers, and the volumetric heat generation Q for the whole body or a selected body. The output is the nodal temperature field, the minimum and maximum, and the power balance in watts:

On the left are heat generation and incoming flux, on the right the heat removal through faces with prescribed temperature (from reactions), convection and radiation. The imbalance shows whether the boundary conditions are consistent with the matrix.

How it is implemented

The elements are linear and quadratic tetrahedra (4- and 10-node types). The kernel also supports 8/20/27-node hexahedra, the wedge and the pyramid. The boundary faces are first- and second-order triangles and quadrilaterals. The matrix is sparse: small systems are solved by direct LDL^T factorisation, large ones by the conjugate gradient method with an incomplete Cholesky preconditioner. Without a prescribed temperature or convection on at least one face, the problem is singular (the temperature is defined up to a constant), and the analysis fails with an explanation.

Interface: thermal analysis in the ribbon, thermal boundary conditions ("Heat flux", "Volumetric heat generation") at the boundary conditions step, and the result is the "Temperature" field and the power balance.

material preset=steel kappa=40
thermal surface number 1 temperature=100
thermal surface number 2 convection=25 ambient=20
load heat_generation=2e6
solve thermal

10.2Transient heat conductionbeta

A transient analysis gives the evolution of temperature under time-varying loads: heating and cooling, thermal shock, cyclic regimes. It is needed when the thermal inertia of a part matters, not only the steady state.

Formulation

The heat capacity matrix C_{ij}=\int_\Omega\rho c_p\,N_iN_j\,dV is either consistent or HRZ diagonal (Hinton, Rock, Zienkiewicz, 1976): the diagonal of the consistent matrix, scaled per element so that the sum equals \int\rho c_p\,dV. The diagonal matrix removes oscillations at a sharp front but is supported only for linear elements.

Time integration is the \theta-method:

At \theta=1 the scheme is backward Euler: first order, L-stable, monotone. At \theta=\tfrac12 it is Crank–Nicolson: second order, A-stable, but it gives decaying oscillations at a load jump — this is a property of the scheme. The step matrix is constant and is therefore factorised once. The step is fitted to the end time: \Delta t_{\text{used}}=t_{\text{end}}/\lceil t_{\text{end}}/\Delta t\rceil. A user-defined non-uniform step schedule is also possible.

What the user specifies

The step \Delta t, the end time and the initial temperature are mandatory: the product does not guess the initial state. In addition: the material density and c_p, the \theta scheme and the heat capacity type, probes (points in mm, snapped to the nearest node), and the number of saved field frames. Each load channel is multiplied by a piecewise linear time table k(t): temperature T_D(t)=T_D\,k_T(t), flux, convection coefficient, environment temperature, heat generation. Outside the table the multiplier is constant. The temperature at nodes with a prescribed value at t=0 equals T_D(0), so a wall jump at t=0 creates no impulse in the first step.

Energy balance

Summing the rows of the step equation gives the exact identity

where P is the total power of heat generation, flux and reactions at prescribed temperature, less convection and radiation. The product outputs the residual of this identity and an independent integral of the accumulated energy \int\rho c_p\,\Delta T\,dV. A residual of the order of machine precision means that the reactions, convection and sources are assembled consistently with the matrix.

Interface: at the solve step, a "Steady-state / transient" switch and fields for the step, time, scheme and k(t) tables.

material preset=steel kappa=40 specific_heat=460
thermal surface number 1 temperature=120
solve thermal transient dt=0.1 tend=16 t0=20 theta=0.5 probes="10,0,0;20,0,0"

10.3Nonlinear heat conduction: κ(T) and radiationbeta

The nonlinear mode is enabled automatically if the model contains a radiating boundary or a \kappa(T) table. It is used for high-temperature regimes: cooling of a heated part, heat removal by radiation, materials with a strong temperature dependence of conductivity.

Formulation

Radiation from a surface into the environment according to the Stefan–Boltzmann law:

here \sigma is in W/(m²·K⁴), and \varepsilon is the emissivity.

The dependence \kappa(T) is a piecewise linear table of "temperature : conductivity" pairs; outside the table the value of the end point is taken. The residual and tangent matrix of the steady-state problem:

The tangent is complete and nonsymmetric (because of \kappa'), so the Newton–Raphson method converges quadratically. In a transient analysis, at each step of the \theta-method the nonlinear equation

is solved with the Jacobian \tfrac{1}{\Delta t}\mathbf C+\theta\,\mathbf J_{\text{int}}. The Newton step is damped by the residual norm; the convergence criterion is the relative residual with allowance for the round-off level.

What the user specifies and what is obtained

Radiation: the emissivity \varepsilon\in(0,1] and the environment temperature in °C. The \kappa(T) table is pairs of "°C : W/(m·K)", temperatures strictly increasing, \kappa>0. The result is supplemented with the number of Newton iterations, the residual history and the radiation power in the balance. If the method has not converged, the temperature is not returned: the analysis fails rather than producing an approximation.

material preset=steel kappa_table="0:20,400:60" specific_heat=460
thermal surface number 2 radiation=0.85 ambient=20
solve thermal transient dt=5 tend=600 t0=900 theta=0.5

10.4Thermal contact at an interfacebeta

Thermal contact models the thermal resistance between two separate bodies of an assembly: the flux through the interface is proportional to the temperature jump. It is needed for bolted and press-fit joints, gaskets, contact pads. Mechanical contact is described in the section Contact; thermal contact does not depend on it and does not affect the mechanics.

Formulation

One body of the pair is the slave (the first in the command), the other is the master. At the Gauss points of the slave faces the projection onto the master surface is found; a point counts as an interface if the distance to the master face does not exceed the tolerance, the normals are nearly antiparallel, and the gap is directed along the normal. The remaining faces of the bodies stay adiabatic. The contribution to the matrix \mathbf K_c=\sum h_c\,w\,\boldsymbol\varphi\boldsymbol\varphi^T, where \boldsymbol\varphi=[N_s;\,-N_m], is symmetric and semidefinite: the interface neither creates nor removes heat by itself. The method is collocation over slave points, without mortar multipliers. Since the contribution is linear, contact by itself does not activate Newton iterations.

What the user specifies and what is obtained

The interface conductance h_c in W/(m²·K) and the gap tolerance in mm (by default, a quarter of the mean face size). The output is the power through the interface, the area of the bonded surfaces, the number of bonded points and the mean temperature jump. If h_c\le0, the interface is not found within the tolerance, or the bodies are glued by common nodes, the analysis fails: the interface is neither silently insulated nor silently perfect.

contact body1 body2 thermal hc=1000
thermal surface number 5 temperature=100
thermal surface number 32 temperature=0
solve thermal

10.5Thermoelasticity (one-way, static)beta

One-way thermoelasticity: temperature affects the stress–strain state, and there is no reverse effect of deformation on temperature. It is used for thermal stresses in parts with a known or computed temperature field.

Formulation

The thermal strain is isotropic: \boldsymbol\varepsilon_{th}=\alpha\,\Delta T\,\mathbf I, \Delta T=T-T_{\text{ref}}, where T_{\text{ref}} is the stress-free temperature. The stress is determined by the mechanical part of the strain:

Solid bodies

For an isotropic material \mathbb C\,\boldsymbol\varepsilon_{th}=c_0\,(1,1,1,0,0,0)^T, where

Nodal forces are assembled for all seven types of solid elements. For the quadratic tetrahedron, the load is integrated with a more accurate quadrature than the stiffness (27 points): the integrand contains the product of the shape gradient and the quadratic field \Delta T. Stresses are recovered from the total strain, so in the postprocessor the volumetric part c_0\,\mathbf I is subtracted from them. The source of \Delta T is either a uniform value set by the user or a nodal field from a steady-state thermal analysis (coupled thermoelasticity).

Beams and frames

A uniform \Delta T and temperature differences between the section faces \Delta T_y, \Delta T_z are specified. The section temperature is T(y,z)=T_{\text{ref}}+\Delta T+g_y\,y+g_z\,z, g_y=\Delta T_y/h_y. The equivalent nodal forces: axial N_T=E\,A\,\alpha\,\Delta T, bending M_{T,z}=E\,I_z\,\alpha\,g_y and M_{T,y}=E\,I_y\,\alpha\,g_z. The same quantities are subtracted from the internal forces. There are no shear forces or torsion from heating.

Shells

For an isotropic shell in plane stress, c_0=E\,\alpha\,\Delta T/(1-\nu) (not 1-2\nu: the third direction is unrestrained). For a laminated shell (MITC4), the thermal forces and moments are computed per layer:

Here \bar{\mathbb Q}_k is the transformed stiffness of the layer in the element axes, \bar{\boldsymbol\alpha}_k=(\alpha_1c^2+\alpha_2s^2,\ \alpha_1s^2+\alpha_2c^2,\ 2(\alpha_1-\alpha_2)\,sc)^T are the layer expansion coefficients rotated by the ply angle (engineering shear), z is measured from the reference surface, and the sign of the bending part is set by the adopted direction of curvature. The temperature through the thickness is linear: the mid-plane \Delta T and the "top minus bottom" difference are specified. Layer stresses are computed from the mechanical strain \boldsymbol\varepsilon-\bar{\boldsymbol\alpha}_k\,\Delta T(z).

What the user specifies

The material coefficient \alpha and the heating: a uniform \Delta T (in DSL deltaT=) or a temperature field from a thermal analysis (solid bodies only). For a beam, dt, dty, dtz; for a shell, shell_temperature dt=… gradient=… or temperature=… reference=… (stress-free temperature). A zero value is not counted as a load, and \alpha is then not required; with a nonzero \Delta T and no \alpha the analysis fails loudly.

Interface: the "Thermoelasticity" analysis type (\Delta T and \alpha fields), "Coupled thermoelasticity" at the solve step (heat → stresses), and for beams and shells, temperature fields in their panels.

material preset=steel
fix surface number 1
solve thermoelastic deltaT=80

10.6Coupled thermoelastic dynamicsbeta

A transient "temperature → displacements" analysis: the thermal field changes in time, and the mechanics responds with inertia taken into account. It is used for thermal shock, rapid heating, thermoelastic waves, where the quasi-static approximation is unacceptable.

What "coupled" means here

The coupling is one-way: the temperature T(t) affects the displacements \mathbf u(t), and there is no reverse effect (heating from deformation) in the product analysis. The name is inherited from the formulation "heat and mechanics on one mesh in one step". For steel, the dimensionless strength of the back-coupling

is about 1 %. For fast processes (impact, ultrasound) this correction may be significant; see Gough–Joule.

Formulation

The thermal load operator \mathbf G is assembled once by the same integral as in statics, F_i=\int c_0\,\nabla N_i\,dV, and at a step reduces to a matrix–vector multiplication. The order within a step is "heat, then mechanics" (the isothermal split of Armero–Simo, 1992): first the \theta-step of heat gives T_{n+1}, then the thermal load \mathbf f_{th}(T_{n+1}) enters the right-hand side of the mechanics step. For HHT-α the integrator takes (1+\alpha)\,\mathbf f_{th}(T_{n+1})-\alpha\,\mathbf f_{th}(T_n). With one-way coupling the scheme is unconditionally stable, and the order equals the lower of the orders of the heat and mechanics schemes.

Two modes:

• Quasi-statics without inertia: at each thermal step \mathbf K_s\,\mathbf u=\mathbf f_{\text{ext}}+\mathbf f_{th} is solved; the linear problem is factorised once. A J2 elastoplastic material is allowed (desktop only).
• Dynamics: Newmark or HHT-α integration, with a thermal step equal to several mechanical ones.

In the thermal part, radiation and \kappa(T) from nonlinear heat conduction can be used.

What the user specifies and what is obtained

The step, the end time and the initial temperature are mandatory; the stress-free temperature, the heat scheme \theta, the mechanics scheme, Rayleigh damping, the k(t) tables of thermal and mechanical loads, and probes are optional. The output is the temperature and displacement fields per frame, the series of T, \mathbf u, \dot{\mathbf u} at probes, the stress tensor with the thermal strain subtracted, equivalent stresses, the envelope, energies, and the heat balance.

solve thermoelastic transient mode=dynamic dt=1e-6 tend=2e-4 t0=20 tref=20 scheme=hht alpha=-0.2

10.7Gough–Joule effectexperimental

The reverse channel of thermoelastic coupling, "deformation → temperature": under rapid compression a body heats up, under tension it cools down. The effect is included as a volumetric source in the heat equation.

What is implemented

For adiabatic loading this gives \Delta T=-\dfrac{T_{\text{ref}}\,\beta}{\rho c_p}\,\operatorname{tr}\boldsymbol\varepsilon (Kelvin, 1853; Biot, 1956). The scheme lags by one step: after a mechanical step the nodal power is assembled and mixed into the right-hand side of the next thermal step. It works only in the quasi-static mode; dynamics with the effect is refused.

A terminological note. In the literature, the "Gough–Joule effect" is the entropic elasticity of elastomers: stretched rubber heats up, and when heated it contracts. What is implemented here is the linear thermoelastic (piezocaloric) Kelvin effect for a compressible linearly elastic body. A model of the entropic elasticity of rubber coupled to the hyperelastic laws is not implemented.

10.8CFD: incompressible flow (SU2)stable

CFD computes the steady flow of a liquid or gas around a body or inside a channel: velocity and pressure fields, wall forces, pressure losses, flow rate. The analysis runs in the external solver SU2 (Stanford University Unstructured, an open-source code, LGPL-2.1 licence) in the desktop application; CFD is not available in the browser.

Equations solved

The incompressible Navier–Stokes equations (constant density and viscosity):

In the laminar regime \mu_t=0. In turbulent regimes the Reynolds-averaged equations (RANS) are used, with the eddy viscosity \mu_t from the Menter k–\omega SST model (1994) or Spalart–Allmaras (1992). The Euler (inviscid) and compressible equations are not used. Energy is not solved: the flow is isothermal, and there is no Archimedes buoyancy.

The method is vertex-centred finite volumes on the dual mesh of the original tetrahedra (Economon et al., AIAA J., 2016). The convective fluxes use flux difference splitting (FDS) with second order by the MUSCL scheme and the Venkatakrishnan limiter; the turbulent quantities use an upwind scheme. Convergence to steady state is by implicit Euler in pseudo-time with an adaptive Courant number and, if necessary, geometric multigrid.

What the user specifies and what is obtained

Medium: air, water, oil or kerosene (density and viscosity are configured in the materials). Boundary conditions by surface numbers: a velocity inlet (magnitude, direction, turbulence intensity), an outlet with gauge pressure, a no-slip wall (with optional force monitoring), a slip wall, a symmetry plane. Unmarked faces become no-slip walls without force monitoring; the list of such faces is output in the warnings. For turbulent analyses, prismatic mesh layers at walls are needed (y⁺ of order one): the models are integrated to the wall.

The stopping criterion is set by the user: the maximum number of iterations and the tolerance. If there is a wall with force monitoring, convergence is determined by stabilisation of the force coefficient (the Cauchy criterion); for internal flows, by the relative drop of the pressure residual.

Results in SI: velocity and pressure fields, vorticity, the Q-criterion, C_p, wall shear stress, y^+, the coefficients C_d, C_l and forces (given the reference area and length), the pressure drop between inlets and outlets, mass flow rates, and the residual history. The postprocessor offers slices, streamlines, isosurfaces and particles.

Relation to the rest of the product

CFD uses the same geometry, the same surface numbering log and the same viewport. The analysis is one-way and isolated: there is no transfer of pressure to the structural analysis and no thermal coupling (CHT); you transfer forces and pressures manually.

Interface: CFD settings at the solve step, boundary conditions and the "Flow" postprocessor.

cfd fluid air
cfd set speed=15 turbulence=sst iterations=500
cfd inlet surface number 1
cfd outlet surface number 2
cfd wall surface number 3 monitor=on
solve cfd

10.9CFD: transient problems and moving bodiesbeta

The transient mode is needed when the flow does not reach a steady state (vortex shedding behind a cylinder, pulsations) or when boundaries move: wall rotation, a rotor, piston-like body motion.

What is implemented

Physical time is integrated by second-order dual time stepping: at each physical step, inner iterations are performed until convergence. The user sets the step, the duration, the number of inner iterations and the frame saving frequency. The session is saved to disk, the analysis can be cancelled and resumed from a checkpoint; frames are replayed in the postprocessor.

Motion:

• a rotating wall (the surface velocity is set by the axis and rotation frequency);
• a rotating reference frame for the fluid region (the MRF approach) — a steady averaged solution for a rotor;
• a prescribed translation of a body along one direction (linear, harmonic, expression, table) with mesh deformation by an elastic ALE scheme, first order in time; the wall velocity is controlled relative to the mesh step;
• prescribed and free rotation of a rotor with one degree of freedom (moment, damping).

The motion results are displacement, velocity, force, power and work, normalised to a selected frame.

cfd set mode=transient dt=0.001 time=0.005 inner=40 output=1
solve cfd

10.10Cavity acoustics (Helmholtz equation)experimental

Modal analysis of an air-filled (or fluid-filled) cavity: natural frequencies and pressure mode shapes. It is used to estimate the resonances of volumes — chambers, ducts, cabins.

Formulation

The Helmholtz equation for the pressure amplitude p in a harmonic process e^{i\omega t}:

On the walls the acoustically rigid condition \partial_n p=0 applies (a natural boundary condition). Finite element discretisation leads to the generalised eigenvalue problem

For a closed rigid pipe of length L, the frequencies of the longitudinal modes are f_n=n\,c_0/(2L), including the zero one (constant pressure). The zero mode is also present in the analysis.

How it is implemented

The element is a linear tetrahedron with a consistent acoustic mass. The matrices are dense, and the problem is solved by reduction to standard form via the Cholesky decomposition of the matrix \mathbf M_a and diagonalisation of the symmetric matrix by the Jacobi method. This approach is applicable to small cavities (up to several hundred nodes); for large meshes it is too slow and memory-hungry.

The user sets the speed of sound c_0 (air 343 m/s, water 1481 m/s), the density \rho_0 and the number of modes. The cavity is the body on which the tetrahedral mesh is built. The output is frequencies in Hz, pressure mode shapes by nodes, and an animation of the standing wave p(\mathbf x)\cos\omega t in the view window.

Interface: the "Vibroacoustics (Helmholtz)" section at the results step.

solve acoustics c0=343 rho0=1.204 modes=5
probe acoustic_f1

10.11Vibroacousticsexperimental

Vibroacoustics describes the interaction of a vibrating structure and the air in a cavity: the structure radiates sound, and the cavity pressure loads the structure. In the product, this analysis as a single coupled problem is absent; there are separate building blocks.

What exists in the code

• Boundary impedance (Robin condition). The damping matrix \mathbf C_a=\dfrac{\rho_0}{Z_a}\int_{\Gamma_z}N^TN\,dA with the acoustic impedance Z_a in rayl (Pa·s/m); for air the characteristic impedance is \rho_0c_0\approx413 rayl.
• The structure–cavity coupling matrix \mathbf C_{fsi}=\int_{\Gamma}\mathbf N_u^T\,\mathbf n\,N_p\,dA on the common boundary.
• The frequency response of the cavity to a source: \left(\mathbf K_a-\omega^2\mathbf M_a+i\omega\,\mathbf C_a\right)\mathbf p=\mathbf q with a solution on a frequency grid and conversion to the sound pressure level \text{SPL}=20\log_{10}(p_{\text{rms}}/p_{\text{ref}}), p_{\text{ref}}=20 µPa for air and 1 µPa for water.

These are separate functions with unit tests. A coupled "structure — cavity" system is not assembled from them, and the frequency response and impedance matrices are not connected to the analysis from the interface or DSL. Only the modal analysis of cavity acoustics with rigid walls is available.

After discretization, any linear problem reduces to a system of algebraic equations with a sparse matrix. This section describes how those systems are solved: which methods are used, why they are valid for FEM matrices, how the product checks the result, and when it honestly refuses. It describes the methods and the limits of their applicability; tuning of internal parameters is outside its scope.

The implementation uses open-source libraries: MFEM (sparse matrices and the conjugate gradient method), Eigen (sparse factorizations and small dense problems), AMGCL (algebraic multigrid preconditioner, desktop application only). On macOS, direct factorization is also performed by Apple Accelerate. All arithmetic in the core runs in double precision (IEEE 754 binary64). The SI system applies at the solver boundary: metres, newtons, pascals, kilograms. The user sees results in their own units; the conversion is performed in a single place.

The native path of the desktop application is the reference. The browser build uses the same formulas, but it is single-threaded, limited by the memory of the tab, and has no multigrid preconditioner. The differences are noted in the text.

11.1FEM systems of equations: sparsity, definiteness, constraintsstable

After assembly, a static problem has the form

The stiffness matrix \mathbf{K} has three properties on which the choice of methods rests. It is sparse: a node is coupled only to the nodes of shared elements, a row holds tens of nonzero entries, and the fraction of nonzeros falls as the model grows. It is symmetric, because the bilinear form is symmetric. It is positive semi-definite, and after the rigid body motion is removed it is positive definite: \mathbf{v}^{\mathsf T}\mathbf{K}\mathbf{v} > 0 for all \mathbf{v} \neq \mathbf{0} in the admissible subspace. For solid bodies the condition number grows as O(h^{-2}) under mesh refinement. For thin shells and beams it grows faster: roughly with the square of the ratio of membrane to bending stiffness, that is, (L/t)^2.

A single assembler combines beams, shells and solid elements into one system with six degrees of freedom per node (three displacements, three rotations). Degrees of freedom that carry no stiffness, for example the rotations of nodes surrounded only by solid elements, are eliminated automatically.

Boundary conditions and constraints

Fixed supports, linear constraint equations, rigid regions (RBE2) and bonded interfaces with exact elimination are handled by eliminating the dependent degrees of freedom, not by Lagrange multipliers:

The matrix \mathbf{T} relates the independent degrees of freedom \mathbf{u}_r to the full vector, and \mathbf{u}_0 carries the inhomogeneous parts of the constraints. The reduced matrix \mathbf{K}_r = \mathbf{T}^{\mathsf T}\mathbf{K}\mathbf{T} remains symmetric and positive definite if \mathbf{K} was so on the admissible subspace. No saddle-point system with a zero block arises. Therefore, for linear statics, methods for symmetric positive definite matrices (Cholesky, conjugate gradients) are sufficient, and indefinite methods (MINRES, factorizations with pivoting) are not required. Penalty contact adds stiffness to \mathbf{K} and also preserves positive definiteness, but it worsens the conditioning.

Scaling

Before solving, the reduced system is symmetrically scaled by its diagonal:

This is an equivalent transformation: the answer does not change, but the spread of scales between displacements and rotations, and between constraint rows and elastic rows, disappears. Rows with a non-positive diagonal are not scaled.

11.2Matrix storage and reorderingstable

Matrices are stored in the compressed sparse row (CSR) format (arrays of row offsets, column indices and values). Repeated contributions from neighbouring elements are summed during assembly. The reduced matrix is obtained as the triple product \mathbf{T}^{\mathsf T}\mathbf{K}\mathbf{T} in the same format. For direct factorizations the matrix is converted to the Eigen format, with only the lower triangle where possible.

Local element matrices are computed in parallel but accumulated into the global matrix sequentially, in element-number order. The assembled matrix therefore does not depend on the number of threads.

The number of nonzeros in the factor \mathbf{L} depends on the ordering of the unknowns. Symmetric factorizations use the approximate minimum degree (AMD) ordering: Amestoy P., Davis T., Duff I., An approximate minimum degree ordering algorithm, SIAM J. Matrix Anal. Appl., 1996. It is used by both Apple Accelerate and Eigen SimplicialLDLT. For nonsymmetric and complex systems, sparse LU uses the COLAMD column ordering. Nested dissection (METIS) and reverse Cuthill–McKee (RCM) are not used in the solvers.

For three-dimensional meshes with a good ordering, the factor memory grows roughly as O(n^{4/3}) and the factorization operation count roughly as O(n^{2}) (George A., nested dissection). The measured exponents for Eigen SimplicialLDLT on solid models are slightly higher, about n^{2.3}. The calculation screen shows a direct-solver memory forecast based on an empirical power law in the number of degrees of freedom.

11.3Direct sparse solverstable

For a symmetric positive definite matrix, a Cholesky or \mathbf{L}\mathbf{D}\mathbf{L}^{\mathsf T} factorization with a preliminary permutation \mathbf{P} is used:

On a positive definite matrix this factorization is stable without pivoting (Golub G., Van Loan C., Matrix Computations). The permutation is therefore chosen by fill-in only, not by the magnitude of the entries. The method is exact up to round-off. Conditioning affects the accuracy of the answer, but not the fact of obtaining one, as it does with iterations. A factorization failure means a loss of positive definiteness or a lack of memory.

On macOS the factorization is performed by the Apple Accelerate sparse solver (Cholesky or \mathbf{L}\mathbf{D}\mathbf{L}^{\mathsf T}), which is multithreaded and vectorized. On Windows, Linux and in the browser, Eigen SimplicialLDLT is used. It is single-threaded and noticeably slower. The direct and iterative paths give matching displacements on product models: the relative discrepancy is of order 10^{-9} (table below). If the factorization fails, fallback paths are engaged, and only when all of them are exhausted does the run return a non-convergence flag.

The factorization is performed once and reused for several right-hand sides: in eigenvalue problems, in harmonic analysis over a frequency grid, and in implicit transient runs with a constant time step.

For matrices that may be indefinite, for example \mathbf{K} + \mathbf{K}_{\sigma} in prestressed modal analysis, an \mathbf{L}\mathbf{D}\mathbf{L}^{\mathsf T} factorization without magnitude-based pivoting does not guarantee stability and by itself does not distinguish positive definiteness from indefiniteness. After factorization, the signs of the diagonal of \mathbf{D} are therefore inspected. If the matrix has ceased to be positive definite (the load is above critical), the run returns a non-convergence flag rather than a plausible frequency.

For nonsymmetric systems (the tangent matrices of some nonlinear problems) and complex systems, sparse LU with COLAMD ordering is used, see harmonic analysis and nonlinear analyses.

Measurements

A shell with a clamped edge and pressure, native path, Apple M5, Apple Accelerate (measured 13 Oct 2026):

For solid bodies the direct method is more expensive because of the greater fill-in. On Hex20 with 89 thousand active degrees of freedom, Accelerate on the reference platform solves the system in 35 s, whereas single-threaded Eigen on a Linux x86-64 bench (a third platform, not the reference, for information only) takes 734 s. For this reason, on Windows and Linux the product selects the iterative solver for solid models above a moderate size, see solver selection.

11.4Dense solver for small systemsstable

Dense solution is used for small frame and shell models and for reference checks. The algorithm is as follows: the matrix is checked for symmetry, then factorized in place by Cholesky. For a positive definite matrix, Cholesky is stable without pivoting, so no scaling or permutations are needed. If the matrix is nonsymmetric or not positive definite (a mechanism, loss of stability), the factorization falls back to LU with partial pivoting by column. Symmetry is checked explicitly, because Cholesky reads only one triangle and would silently give a wrong answer on a nonsymmetric matrix. On a singular matrix the result is non-finite, and the calling code recognizes the mechanism.

Dense solution requires O(n^2) memory and O(n^3) time, so it is suitable only for small models. On the native path, shell and frame models above a small size are switched automatically to sparse factorization; in the browser a size ceiling set by the tab memory applies. The dense and sparse paths give the same displacements: the discrepancy on a shell model is of order 10^{-11} relative.

11.5Iterative solver: CG with an algebraic multigrid preconditionerbeta

The preconditioned conjugate gradient method with preconditioner \mathbf{M} solves \tilde{\mathbf{K}}\tilde{\mathbf{u}} = \tilde{\mathbf{f}} by a sequence of steps along conjugate directions. The method is valid if both the matrix and the preconditioner are symmetric and positive definite. The convergence rate is determined by the condition number of the preconditioned matrix: the number of iterations grows as \sqrt{\kappa(\mathbf{M}^{-1}\tilde{\mathbf{K}})} (Saad Y., Iterative Methods for Sparse Linear Systems, 2003, ch. 6).

In the desktop application the preconditioner is algebraic multigrid (AMG) with smoothed aggregation from the AMGCL library (Vaněk P., Mandel J., Brezina M., Computing 56, 1996; Demidov D., AMGCL, Lobachevskii J. Math., 2019). Groups of coupled nodes are merged into aggregates, coarse levels are built, smoothing is performed on each level, and a direct solve is performed on the coarsest. One V-cycle serves as the preconditioner. Where the system is blocked, the aggregation takes into account the three degrees of freedom of a node. The smoother is symmetric, so the V-cycle is symmetric by construction, as the conjugate gradient method requires. A property of AMG is that the number of iterations barely grows under mesh refinement. On a Linux x86-64 bench (for information only), for Hex8 with 14, 105, 202 and 346 thousand active degrees of freedom, AMG-CG converges in 32, 43, 48 and 52 iterations. Plain Jacobi or Gauss–Seidel do not have this property.

The stopping criterion is the relative residual. For the AMG path it is the Euclidean norm of the residual divided by the norm of the right-hand side: \lVert \mathbf{r}_k \rVert_2 / \lVert \tilde{\mathbf{f}} \rVert_2 \le \varepsilon. For the remaining stages (the MFEM conjugate gradient method) it is the norm of the preconditioned residual divided by the initial one. The tolerance \varepsilon and the maximum number of iterations are set by the user at the "Calculation" step; the defaults are \varepsilon = 10^{-6} and 5000 iterations. The residual history is displayed in the calculation window, and the run can be cancelled.

Safeguards against incorrect operation. The product checks the preconditioner at every iteration: the inner product (\mathbf{M}^{-1}\mathbf{r}, \mathbf{r}) must remain positive. If it does not, the conjugate gradient method stops and names the cause, rather than iterating further along garbage directions. A diverging residual also aborts the AMG stage. The Euclidean norm of the conjugate gradient residual is non-monotonic and on solid bodies may temporarily grow before converging, so the abort triggers only on clear divergence. If a stage fails, the run moves to a fallback path. If no stage converges, the result is marked as not converged.

Comparison with the direct method on the reference platform (Apple M5, 2 threads, measured 21 Oct 2026 on a loaded machine, times are approximate; the same model was solved by both paths):

The time advantage of the iterative path grows with model size and becomes decisive when the direct method runs out of memory.

In the browser build AMGCL is unavailable. There the iterative path uses the conjugate gradient method with simpler preconditioners (incomplete Cholesky on small systems, Gauss–Seidel on large ones). They converge more slowly on fine meshes, so browser iteration counts must not be carried over to the desktop.

11.6Result control: residual, convergence flag, refusal instead of substitutionbeta

Every run returns a convergence flag. It is set only when the achieved residual does not exceed the tolerance and the result is finite. When it is absent, the interface shows the warning "Did not converge — the values below may be unreliable", and a script ends with the error "solver did not converge". The result is not presented as successful.

The principles followed by the core:

• Refusal instead of silent substitution. If the chosen method is inapplicable or does not fit in memory, the run returns an error stating the reason and the number of degrees of freedom, rather than quietly solving by another method with a different accuracy. Invalid solver settings are rejected at input. For thin shells the size ceiling depends on the machine memory (the limit from "Settings → Compute" is taken into account); above it the run refuses.
• Checking at every step. The conjugate gradient method checks the positivity of the preconditioner, and AMG checks residual growth. The harmonic run checks the residual at every frequency. The nonlinear solver recognizes a singular tangent and reports a diagnosis.
• Eigenvalue problems. For frames and shells on the native path, the relative residual of every requested eigenpair is checked. For solid bodies the stopping criterion is the stabilization of the Ritz values; the eigenpair residual on this path is not reported by the product.
• Equilibrium check. Support reactions and integral quantities can be requested with postprocessing commands; for frame models frame probe equilibrium is available, giving the relative residual \lVert \sum \mathbf{R} + \sum \mathbf{F} \rVert in forces and in moments.

Accuracy and reproducibility. Arithmetic is double precision in all core paths. The stiffness matrix does not depend on the number of threads. Bitwise reproducibility is not claimed: multithreaded summations and library factorizations may differ at the round-off level. Bitwise agreement between platforms does not exist and is not expected: different libraries and operation orders are used, and fused multiply–add operations are possible natively. An early comparison of the browser and native solutions on a cantilever beam gave a relative discrepancy of order 10^{-9}.

11.7Matrix-free (operator) pathbeta

In the matrix-free approach the global matrix is not assembled: the solver knows only the action \mathbf{K}\mathbf{v} on a vector. For simple linear statics of solid bodies (surface supports and loads, without rigid regions and constraints), the product uses MFEM operator assembly (partial assembly) with first- and second-order elements and the conjugate gradient method with a diagonal (Jacobi) preconditioner. The operator is physically the same as that of the assembled matrix: an early comparison of partial and full assembly on a cantilever beam gave a discrepancy of order 10^{-11}.

The geometric multigrid method (MFEM h-multigrid with Chebyshev smoothing) was tested as a research bench: when the problem size grew by a factor of 225 (189 to 42,483 degrees of freedom), the iteration count grew only by a factor of 1.9 (28 to 52), that is, the convergence does not depend on the mesh. This method is not included in the product calculation path: it requires a nested mesh hierarchy, whereas the product builds unstructured meshes.

11.8Eigenvalue problems: modal and buckling analysisstable

Modal analysis solves the generalized problem

Buckling: (\mathbf{K}_r + \lambda\,\mathbf{K}_{\sigma})\boldsymbol{\varphi} = \mathbf{0}, where \mathbf{K}_{\sigma} is the geometric stiffness due to the prestress state, and \lambda_{\text{cr}} is the smallest positive eigenvalue (the critical load multiplier). The formulations and assumptions are described in the section Linear analysis types. Here the method of solution is described.

Block inverse iteration (subspace)

For solid and mixed models, subspace iteration with the Rayleigh–Ritz procedure is used (Bathe K.-J., Finite Element Procedures, 2014, ch. 11). The iteration step:

The block contains more vectors than the number of modes requested; the small dense p \times p problem is solved by a standard method (Eigen). The iterations stop when the relative change of the requested Ritz values between steps becomes negligible. The matrix \mathbf{K}_r is factorized once, and each subsequent \mathbf{K}_r^{-1}\mathbf{b} is a pair of triangular solves. On macOS the Accelerate factorization is used for models of moderate size, and on the other platforms Eigen SimplicialLDLT. For very large models \mathbf{K}_r^{-1} is applied approximately through AMG-CG with an inner tolerance; an approximate inverse requires caution (see the limitations). A small diagonal addition (a "floor") is added to the mass matrix so that massless degrees of freedom do not break the Cholesky factorization of the small matrix \mathbf{M}_p; the spurious high modes it generates are filtered out.

The buckling problem uses the same subspace with a different form. The matrix -\mathbf{K}_{\sigma} is sign-indefinite and of low rank, so the spectrum is shifted by a constant c, and the vectors are orthonormalized in the \mathbf{K}_r inner product. The reciprocal values \theta = 1/\lambda are recovered from the shifted values.

Frames and shells: LOBPCG

For frames and shells on the native path, the block method LOBPCG is used (Knyazev A., SIAM J. Sci. Comput., 2001) with the exact preconditioner \mathbf{K}_r^{-1} from a single Cholesky factorization. In essence it is a block Lanczos method with a locally optimal direction. It works with the consistent mass matrix and does not require a spectrum shift for the buckling problem. The statics under the reference load is solved with the same factorization, so its error does not enter the critical multiplier. The stopping criterion is the relative residual of each requested pair. Comparison with the dense path on the same models gives a discrepancy of no more than 10^{-7}. For small systems the dense path remains: the generalized problem is reduced to a standard one through the Cholesky factorization of \mathbf{K}_r and solved by the Jacobi method. It is the reference, but its time grows as n^3.

Free bodies and mechanisms

For an unconstrained model the matrix \mathbf{K}_r is singular (six rigid-body modes). In that case the spectrum is shifted: (\mathbf{K}_r + \tau\,\mathbf{M}_r)\boldsymbol{\varphi} = (\lambda + \tau)\,\mathbf{M}_r\boldsymbol{\varphi} is solved with a small \tau > 0, after which the shift is subtracted. The rigid modes remain zero frequencies and are separated from the elastic ones. For frames, the command frame solve modal free=on explicitly discards the six rigid-body motion modes. If the elastic modes are lost behind a rigid mode (fewer modes than requested and the first frequency close to zero), the product rejects the result and names the cause: the model is under-constrained.

solve modal nmodes=6
solve buckling nmodes=4

11.9Complex systems of harmonic analysisbeta

The steady-state harmonic response \mathbf{u}(t) = \operatorname{Re}\left(\hat{\mathbf{u}}\,e^{i\omega t}\right) with Rayleigh damping gives, at each frequency, a complex system

The matrix is complex-symmetric but, with damping, not Hermitian, so the Cholesky method is inapplicable. Sparse complex LU is used (Eigen SparseLU, COLAMD ordering). The symbolic analysis of the structure is performed once for the whole frequency grid, and the numerical factorization at every frequency. The constraints and the reduction are the same as in statics, and the system is scaled by its diagonal. After the solve, the relative residual is computed at each frequency; if it exceeds the tolerance, the run stops and indicates the frequency. Singularity (resonance without damping) also gives an error with the frequency rather than a false answer. In a script, the modal superposition method (modal) or the direct method (direct) is selected.

solve harmonic fmin=20 fmax=400 npoints=60 method=direct rayleigh a=1.5 b=0.000002

11.10Linear algebra in nonlinear analysesstable

Nonlinear problems (geometric and material nonlinearity, contact, the arc-length method) are solved by Newton iterations. At every iteration a linear system with the tangent matrix \mathbf{K}_T(\mathbf{u})\,\Delta\mathbf{u} = -\mathbf{r}(\mathbf{u}) is solved. The tangent is factorized by a sparse direct factorization. The symbolic analysis and ordering are reused between iterations as long as the sparsity structure has not changed; a change of structure (for example, a change of the active contact) triggers a new analysis. Newton convergence criteria, load steps and the arc-length method are described in the section Nonlinear analyses and dynamics.

11.11Parallelism and platformsbeta

The desktop application uses OpenMP: parallel assembly of element matrices and a multithreaded AMG preconditioner. The Apple Accelerate direct solver is itself multithreaded. The default thread count was chosen by measurement, and the user changes it in "Settings → Compute" ("Cores for the solver"). Scaling is limited by memory bandwidth, not by the number of cores: on a model of 215 thousand degrees of freedom (Apple M5), two cores gave a speedup of 1.2–1.5, six cores about 1.8, and ten are no faster than six. The field applies to the paths that go through AMG; some native calculations are single-threaded.

The browser build is single-threaded, uses only the memory of the tab (a 32-bit WebAssembly module does not go beyond 4 GB) and has no AMG. On large models the browser hits memory and time limits long before the limits of the native path. The solver mode and size threshold are not configurable in the browser. Eigen direct factorization works there up to a memory threshold; above it, iterative CG is used.

The reference verification numbers are native ones. On the same geometry the platforms may diverge: for example, the browser and desktop meshers build different meshes. When comparing numbers, look at the platform label.

11.12Solver selection: automatic and manualbeta

The user controls the selection on two screens.

The "Calculation" step, solver panel. The type (direct or iterative), the "Prefer direct solver" checkbox, and the memory limit for the direct solver (MB). For the iterative mode, the tolerance and the maximum number of iterations are available. The panel shows the memory forecast of the direct method and the chosen route. If the forecast does not fit within the limit, the solver switches to iterative and shows a notification. With active contact the type is fixed to iterative CG.

"Settings → Compute" (desktop only): the solver mode "auto" / "direct" / "iterative", the direct-method size threshold in degrees of freedom, the thread count, the memory limit. The "auto" mode selects the direct method within the threshold and the iterative one above it. An empty threshold means automatic selection: for solid bodies it is computed from the memory limit of the calculation screen. On Windows and Linux the automatic threshold for solid models is lower, because the direct method is slower there (see the direct solver). For shells and frames the direct method is forced on automatically for the duration of the call, because iterative methods are unreliable on them.

In scripts there is no separate command for choosing the linear solver: solve uses the document settings. The settings are available from the scripting programming interface (cae.solver, cae.compute).

11.13Prototypes and unimplemented features: GPU, cloud, matrix-free P2experimental

These modules exist in the code but are not connected to the product calculation. They are listed so that the boundary of capabilities is clear.

• Matrix-free PCG for Tet10 and Hex20 on WebGPU. The module computes \mathbf{y}_e = \mathbf{K}_e\mathbf{u}_e element by element with graph coloring, preconditions with the diagonal and a two-level p-multigrid. It is covered by tests, but no calculation calls it. The "Matrix-Free P2 GPU" checkbox on the solver panel is saved in the project and changes nothing.
• WebGPU CG solver. The module solves a sparse system in single precision on the GPU with a CPU fallback variant; a development test confirms that calculations do not call it. The product core computes in double precision.
• Hybrid cloud calculation. The routing module and the mesh transfer client decide whether to hand the calculation to a remote server (by model size, analysis type, and the "Auto / Force Local / Force Cloud" setting). The server side is not connected to the calculation, and the mode selection does not change the calculation. In the current build the calculation always runs locally.
• LOBPCG with IC(0) for P2 tetrahedra. A research module of sparse modal analysis; not included in the main path.
• Geometric multigrid (GMG). Present in test benches (see above), not used in the calculation.

After a run the solver returns displacements (and rotations for beams and shells). Everything else — stresses, strains, reactions, safety factors, mesh error estimates — is built from them by the postprocessor. This section describes exactly what is computed, in which axes and units, where the postprocessor makes simplifications, and which quantities remain estimates rather than the result of a rigorous calculation.

Conventions for the whole section:

• tensors are written in Voigt order [xx,\ yy,\ zz,\ xy,\ yz,\ zx]; tension is positive;
• for stresses the shear components are tensorial, for strains \gamma_{ij} = 2\,\varepsilon_{ij} (engineering shear);
• the postprocessor of a linear elastic run works from the displacements and the element material: \boldsymbol{\varepsilon} = \mathbf{B}\,\mathbf{u}_e, \boldsymbol{\sigma} = \mathbf{C}\,\boldsymbol{\varepsilon}, where \mathbf{B} and \mathbf{C} are the same as in the element stiffness matrix;
• in nonlinear analyses the fields (Cauchy stress, plastic strain, reactions) come from the nonlinear core and are not recomputed by linear reconstruction from displacements.

12.1Result fields and unitsstable

The result of one run is a set of nodal fields. They can be colour-mapped on the model, sampled with a probe, and written to a report and to a file. The main groups:

Stress and strain components can be shown not only in global axes but also in a user coordinate system (Cartesian or cylindrical): the tensor is transformed at each node by the formula T'_{ij} = \mathbf{v}_i^{\mathsf T}\,\mathbf{T}\,\mathbf{v}_j, where \mathbf{v}_i are the unit vectors of the system at that point. The invariants (von Mises, \sigma_1, \sigma_3) do not depend on the choice of system.

Units

At the solver boundary SI applies (m, N, Pa, kg), and the model geometry is stored in millimetres. Result fields are stored in the units in which the pipeline naturally produces them:

• stresses — in pascals;
• displacements — in model units (mm), so that the deformed view adds to the original geometry without conversion;
• reactions are computed from model coordinates, so they are converted to newtons before display;
• strains and safety factors are dimensionless.

Display uses the unit system chosen in the status bar: SI, "metric engineering" (mm, N, MPa) or imperial (inch, pound-force, psi). A single table converts the values for the legend, probes, the value filter, the manual scale range and the report, so the label on the scale and the number are always consistent. In the DSL, quantities are returned in SI unless the description of the quantity says otherwise. In a VTK file the unit is part of the field name (displacement_mm, vonMises_Pa, reactionMag_N) — ParaView knows nothing about units.

12.2Stress recovery to nodesstable

A finite element gives the stress inside the element, while colour mapping and probes work with nodal values. The "element → node" transition is a separate step with its own errors.

What is computed

For a linear elastic run the stress at a point of the element is \boldsymbol{\sigma} = \mathbf{C}\,\mathbf{B}(\mathbf{x})\,\mathbf{u}_e, where for an isotropic material

The evaluation point depends on the element type:

• TET4 — the strains are constant, one stress per element;
• TET10 — inside an element with straight edges the stress varies linearly; it is evaluated at the ten nodes of the element from the gradients of the element's own shape functions;
• HEX8, HEX20, HEX27, prism, pyramid — one value at the element centre (centroid).

A node receives the arithmetic mean over the elements that contain it. The tensor components are averaged, not the equivalent quantities: von Mises, \sigma_1 and \sigma_3 are computed afterwards from the averaged tensor, because the principal stresses depend nonlinearly on the components. The weights are by element count, not by volume. In a mixed TET4/TET10 mesh (a zone of quadratic elements) each element is evaluated in its own basis, the node is averaged over its own neighbours, and the field remains single-valued at the node.

In nonlinear analyses of solids, the Cauchy stress is taken at the Gauss points of the element and averaged over the element volume in the reference configuration; then the node is averaged over its neighbours with a weight equal to that volume.

In thermoelasticity mode the thermal part \sigma^{\text{th}} = E\,\alpha\,\Delta T/(1-2\nu) is subtracted from the normal components, so the stresses shown are mechanical.

In the results panel one can switch between the "nodal" (averaged) and the "element-wise" field. The element-wise field is stepped, but it honestly shows the jumps between elements; the difference between the two modes is a rough estimate of how well the mesh resolves the gradient (the error indicator estimates it more accurately).

What matters

The boundary of materials and bodies. The nodal field is single-valued, so at the boundary of two materials (or two joined bodies) the averaging smears the physical stress jump between materials with different E. The nodal value on the boundary itself is the mean of two different quantities, not the value "on the left" or "on the right". For places where the jump is significant, look at the element-wise mode or at separate regions. For this reason a safety margin map is not built for models with several materials: the allowable stress at a boundary node is not tied to any one material. The error indicator, by contrast, recovers the field separately by material.

Peaks near the boundary. Values at the element centre underestimate the peak on a coarse mesh: the centroid lies deep inside the element, not on the outer fibre. In pure bending of a cantilever on a TET10 mesh (four element layers through the height) the nodal evaluation gave an error of 0.34 % at the outer fibre, whereas the centroid values of the same run gave 14–18 % (measured during development).

Beams are not averaged at a node. At one node of a frame, members with different internal forces meet, and the average would be incorrect by construction.

12.3Equivalent stresses and strainsstable

From the stress tensor, scalars are computed that are convenient to compare with the material limits.

Von Mises equivalent stress (distortion energy):

Principal stresses \sigma_1 \ge \sigma_2 \ge \sigma_3 are the eigenvalues of a symmetric 3\times3 matrix. They are found in closed form by Smith's trigonometric formula: let q = \operatorname{tr}\boldsymbol{\sigma}/3, p = \sqrt{\lVert \boldsymbol{\sigma} - q\mathbf{I}\rVert_F^2/6}, \mathbf{B} = (\boldsymbol{\sigma} - q\mathbf{I})/p, \varphi = \tfrac13\arccos\!\left(\det\mathbf{B}/2\right), then

For a spherical tensor (p = 0) all three equal q. As fields on the model, \sigma_1 (maximum principal) and \sigma_3 (minimum principal) are output; \sigma_2 is available in the field calculator.

Tresca stress intensity — the difference of the extreme principal stresses:

It is used in strength maps (the Tresca criterion) and in stress linearization per ASME.

Equivalent strain. For the strain tensor \boldsymbol{\varepsilon} with deviator \mathbf{e} = \boldsymbol{\varepsilon} - \tfrac13\operatorname{tr}\boldsymbol{\varepsilon}\,\mathbf{I}, the "equivalent strain" field equals the norm of the deviator:

where \varepsilon_{xy} = \gamma_{xy}/2 and so on. This is not the "plastic" equivalent strain normalized for incompressibility (\sqrt{2/3}\,\lVert\mathbf{e}\rVert): under uniaxial tension of an incompressible material the value shown is \sqrt{3/2}\approx 1.22 times the axial strain. The accumulated plastic strain \bar\varepsilon^{p} is output as a separate field from the material model of the nonlinear analysis.

12.4Strength criteria and safety factorbeta

The strength map is a postprocess over the already computed nodal stresses; changing the criterion or the limit does not require a new run. For each node the "utilization" D and the safety factor n = 1/D are computed. If D = 0, the margin is capped at the upper display ceiling.

For Drucker–Prager, \beta = (\sigma_c - \sigma_t)/(\sigma_c + \sigma_t) and \sigma_{\text{lim}} = 2\sigma_t\sigma_c/(\sigma_c + \sigma_t); the constants are chosen so that the criterion gives exactly D=1 in uniaxial tension and compression at the limits.

The three fatigue variants of the map are a simplified estimate for a pulsating cycle from zero (the amplitude equals the mean, \sigma_{\text{vm}}/2), with endurance limit \sigma_e:

For an arbitrary load cycle use the fatigue calculation from S-N curves.

If a property required by the chosen criterion is missing, the panel names the missing property rather than building a map by default.

A custom strength condition can be set in the field calculator as a "margin" field K = [\sigma]/f(\sigma).

12.5Shell resultsstable

For shells (the MITC4 quadrilateral with a membrane part and the CST+DKT triangle) stresses are output on fibres through the thickness. At the element centre the membrane part and the bending part are computed:

where \boldsymbol{\kappa} are the curvatures, m_{\alpha\beta} are the bending moments per unit length, and \sigma^m are the membrane stresses of the plane stress state (\sigma_{zz} \approx 0). The sign "+" is the top fibre (along the element normal), "−" the bottom fibre. The equivalent stress on a fibre is the plane one:

In the results panel the fibre is selected: top, bottom, or the one larger in magnitude. The components are computed in the local system of the element. Before nodal averaging, the plane tensor is rotated to global axes (\boldsymbol{\sigma}_{\text{glob}} = \mathbf{R}^{\mathsf T}\boldsymbol{\sigma}_{\text{loc}}\mathbf{R}), otherwise averaging on a curved shell would mix components in different bases; the top and bottom fibres are averaged separately.

For laminated shells the stresses are recovered ply by ply: at three points of each ply (bottom, middle, top) from the generalized strains (\boldsymbol{\varepsilon}^0, \boldsymbol{\kappa}) and the ply stiffnesses; the transverse shear stresses \tau_{13}, \tau_{23} are separate fields. The generalized strains are transferred to the nodes by patch fitting (SPR), without smoothing across a junction of different sections or a fold.

12.6Beam resultsbeta

For each beam element the postprocessor recovers the internal forces at the ends from the nodal displacements: \mathbf{f} = \mathbf{k}_{\text{loc}}\,\mathbf{T}\,\mathbf{u}_e, where \mathbf{T} is the rotation to local axes. The output comprises the axial force N, the shears Q_y, Q_z, the torque T and the bending moments M_y, M_z. Signs: N>0 is tension, T>0 follows the right-hand rule along the axis from the first node to the second. If a distributed load acts on the element, the end forces are corrected by the contribution of the equivalent loads (\mathbf{f}_{\text{int}} = \mathbf{f}_{\text{solver}} - \mathbf{F}).

The diagrams show the values at the element ends joined by line segments; at a node where two elements meet, the mean of the two end values is taken. Therefore the jump of shear force under a concentrated load is smoothed on the diagram: at the jump points look at the values at the element ends.

Section stresses are computed in two ways.

Estimate from library properties (for the summary and for finding the worst section): \sigma = N/A \pm \lvert M_y\rvert/W_z \pm \lvert M_z\rvert/W_y — a sum of magnitudes, that is, the stress at a corner point of the section. For circular and tubular sections this is a conservative estimate: the points of greatest bending in the two planes do not coincide.

Distribution over the section (the "Section" panel) — at the nodes of a grid of section points:

where the first moment of area Q and the width b of the cut-off part are computed numerically on the section grid (the Zhuravskii formula). Torsional shear stresses: circle and tube — \tau = T r/J; rectangle — Prandtl's stress function approximation; closed box — Bredt's formula \tau = T/(2A_m t); thin-walled open profiles (I-beam, T-section, channel, angle) — a linear law through the plate thickness. The equivalent stress: \sigma_{\text{vm}} = \sqrt{\sigma_{xx}^2 + 3\tau^2} with \tau = \sqrt{\tau_{xy}^2 + \tau_{xz}^2}.

12.7Stresses in local material axesbeta

For orthotropic solids (fibrous composites, wood, rolled stock) Hooke's law is written in the material axes 1,2,3, and the stresses in the same axes are more important than the global ones: strength along and across the fibres is checked with them. The orientation of the axes is set for a body by a surface, by a Cartesian, cylindrical or spherical coordinate system, plus an in-plane rotation angle.

The stiffness matrix of an orthotropic material (nine independent constants E_i, \nu_{ij}, G_{ij}) is transformed to global axes by the Bond transformation, the strains are computed through \mathbf{B}\mathbf{u} in global axes, and the stress is transformed to the material axes: \boldsymbol{\sigma}_{\text{loc}} = \mathbf{T}_\sigma\,\boldsymbol{\sigma}_{\text{glob}}, where \mathbf{T}_\sigma is the Bond matrix for stresses. The fields \sigma_{11}, \sigma_{22}, \sigma_{12} are available for colour mapping and probes (stress_local11, stress_local22, stress_local12); the values at a node are the mean over neighbouring elements.

12.8Support reactions and force balancestable

The reaction at constrained degrees of freedom is the force with which the support acts on the body. For solids it is computed through the residual of the equilibrium equation at the constrained DOFs:

where \mathbf{K}\mathbf{u} is assembled element by element (\sum_e \int_{\Omega_e} \mathbf{B}^{\mathsf T}\boldsymbol{\sigma}\,d\Omega) for the constrained nodes, and \mathbf{F} is the right-hand side that includes surface, concentrated, body and nodal loads. The reaction applies both to fixed supports (zero displacement) and to prescribed nonzero displacements; in the latter case it shows the force needed to impose the displacement.

The output comprises: the reaction for each boundary condition entry (components and magnitude), the distribution over the support nodes, and the total reaction \sum R_x, \sum R_y, \sum R_z. The balance check is to compare the total reaction with the sum of the applied loads: for statics \sum \mathbf{R} + \sum \mathbf{F} = \mathbf{0}. For frame models (beams, shells) the reactions are obtained from nodal equilibrium using the end forces and nodal loads.

A reaction is the difference of two large numbers, so any load omitted from \mathbf{F} at a constrained node goes into the answer in full. The postprocessor keeps an account of the completeness of the right-hand side: if some kind of load could not be accounted for (for example, a force on a vertex not found in the mesh), the reaction number is output with an explicit warning rather than silently.

Example

reactions
probe reaction_z

12.9Stress linearization per ASME VIII-2 (SCL)beta

In vessels and piping the strength criterion is not the peak but the integral components of stress through the wall thickness. A stress classification line (SCL) crosses the wall from the inner surface (x=0) to the outer one (x=t); the fields along it are decomposed into membrane, bending and peak parts (ASME BPVC Section VIII Division 2, Annex 5-A; EN 13445-3; GOST 34233).

The tensor components are linearized in the local system of the line (the x axis is along the SCL, the other axes are the meridian and the hoop direction), and the equivalent quantities are then computed from the linearized tensor:

on the outer surface the bending part equals -\boldsymbol{\sigma}_b^{\text{in}}. The peak part is the remainder on the surfaces: \boldsymbol{\sigma}_F = \boldsymbol{\sigma}(x_s) - \boldsymbol{\sigma}_m - \boldsymbol{\sigma}_b(x_s). The integrals are taken by the composite Simpson rule over equally spaced points along the line. For each component (P_m, P_b, P_m + P_b on each surface, F) the von Mises and Tresca equivalent stresses (the intensity \sigma_1 - \sigma_3) are output.

Allowable stresses (if the material properties are specified):

The margin is the ratio of the allowable to the actual value. The allowable S_m can also be set directly.

The line is defined by two points (inside, outside) and a third point to orient the meridian, or by the vessel axis; the default number of points along the line is 47.

Where in the interface, DSL

The "Linearization" results panel: three clicks on the model (points inside, outside, and on the meridian), a plot of the through-thickness diagram and a table of components.

scl p1=(50,0,0) p2=(62,0,0) samples=47 criterion=tresca
probe scl_pm
probe scl_pm_pb

12.10Stress concentration hotspots and the Zienkiewicz–Zhu error indicatorbeta

The stress maximum does not answer the question of whether it is a physical concentration (a hole, a fillet) or an artefact of a coarse mesh (a constraint corner). The answer comes from an estimate of the error of the finite element stress field.

Method

The Zienkiewicz–Zhu (1987) indicator for tetrahedral meshes (TET4, TET10 and mixed).

1. The finite element field \boldsymbol{\sigma}_h is discontinuous between elements. It is taken at the four Gauss points of the element (the stress superconvergence points).
2. A continuous field \boldsymbol{\sigma}^* is recovered: a polynomial is fitted over the patch of elements around a node (SPR — superconvergent patch recovery, Zienkiewicz–Zhu 1992); the value of this polynomial is taken at the node. Recovery is performed separately by material, so the jump at a material boundary is not smeared.
3. The element error is the energy norm of the difference:

The global relative error is \eta = \lVert e\rVert / \sqrt{\lVert\boldsymbol{\sigma}^*\rVert^2 + \lVert e\rVert^2}. For TET10 the difference \boldsymbol{\sigma}^* - \boldsymbol{\sigma}_h on an element is a polynomial of second degree, and the integral is taken exactly (analytical mass matrix), without quadrature error. Next to the absolute \eta_e, the error density \eta_e/\sqrt{V_e} is output — it does not depend on the element size.

From error to hotspots

1. Dörfler marking (Dörfler 1996): the minimal set of elements carrying a fraction \theta of the total squared error, \sum_{e\in M}\eta_e^2 \ge \theta \sum_e \eta_e^2; the classical choice is \theta = 0.5. The criterion is defined through the error rather than the number of elements, so it does not depend on the mesh size.
2. Merging of the marked elements into connected groups (by a shared face).
3. Filtering of noise: single elements and very small groups are more often associated with a degenerate element than with physics; they do not vanish silently but are shown in the report as discarded. If the error over the model is small and uniform, there are no hotspots, and the panel says so.
4. Ranking by error density \sqrt{\sum\eta_e^2/\sum V_e}, not by the sum: the sum puts the coarsest-meshed region first.

For each hotspot the centre, the radius of the circumscribed sphere, the volume, the share of the total error and the current element size are output. From a hotspot a refinement zone (a sphere) is built with size h_{\text{new}} = h_{\text{old}}/k (by default, half) and a margin for the mesh transition band; it can be created with one button or in all hotspots at once.

Where in the interface and DSL

The "Hotspots" results panel: a list, and the buttons "Sphere here" and "Spheres in all hotspots". In the DSL — hotspots and the automatic loop refine adaptive, for which the user sets the target tolerances: gci (the mesh convergence threshold, %) and tol (the ZZ error threshold, %).

hotspots refine=2 minsize=0.5
probe hotspot_count
refine adaptive passes=3 gci=5 tol=2

12.11Mesh convergence: Richardson extrapolation and GCIbeta

Solver convergence (the drop of the residual within one run) is one thing; convergence with the mesh (what the stress peak tends to as the mesh is refined) is another. The module answers the second question by a series of independent runs on a progressively finer mesh.

The metric is the maximum of the field inside the hotspot sphere fixed from the previous iteration (the global maximum on different meshes may lie in different places and cannot be compared). The characteristic size of the element h is taken from the mesh itself in this sphere as h = (\overline{\Delta V_e})^{1/3} (Celik et al., ASME J. Fluids Eng., 2008), not from the requested size: the mesher delivers not what was asked for, and the method works with the ratios of h.

From the last three meshes (h_1 > h_2 > h_3, peak values f_1, f_2, f_3) the following is computed:

The order p is computed signed, without the absolute value from the Celik notation: a negative p is a diagnosis of divergence, not noise. The verdict is one of:

• "converged": the limit f_{\text{ext}} \pm \text{GCI}\cdot\lvert f_3\rvert, GCI not above the target (by default 5 %);
• "converging, but the band is still wide" — refine further;
• "not converging — probably a singularity" (p \le 0): add a fillet, further refinement will not give an answer;
• "insufficient history" (fewer than three iterations) or "input unsuitable".

If the sign of the increments has changed (the series oscillates around the limit), the verdict is "refine further": extrapolation in such a series is unreliable.

12.12Fatigue from S-N curvesbeta

Life estimation in the high-cycle fatigue region by stresses: a static linear run gives the stress distribution, it is scaled by the load history, and the life is computed from the Wöhler curve at each node.

The S-N curve is one of three:

• Basquin: \sigma_a = \sigma'_f\,(2N_f)^b, whence N_f = \tfrac12(\sigma_a/\sigma'_f)^{1/b}; an optional endurance limit \sigma_e: below it N_f = \infty.
• Two-slope with a knee at the base N_e: N_f = N_e(\sigma_a/\sigma_e)^{1/b_1} for \sigma_a \ge \sigma_e and N_f = N_e(\sigma_a/\sigma_e)^{1/b_2} below the limit (for b_2 = 0 — infinite life).
• Tabular: piecewise interpolation in \log\sigma–\log N axes.

The mean stress correction \sigma_m converts the cycle (\sigma_a, \sigma_m) to an equivalent fully reversed amplitude \sigma_{a,\text{eq}}:

When \sigma_m \ge the limit of the denominator, the correction gives infinity (static failure). By default, for a compressive mean \sigma_m \le 0 the correction is not applied (conservatively, \sigma_{a,\text{eq}} = \sigma_a); accounting for the beneficial effect of compression is enabled explicitly.

The loading is specified in one of the following ways: the stress ratio R = \sigma_{\min}/\sigma_{\max} (by default R=-1); a range of scales; a history of scale factors of the static solution; a sequence of results of several steps (non-proportional loading). The history is decomposed by the rainflow method (ASTM E1049-85, §5.4.4) into closed cycles and half-cycles \{\Delta\sigma,\ \sigma_m,\ w\} (weight 1 or 0.5); the residual half-cycles enter with weight 0.5. The damage is summed by the Palmgren–Miner rule:

Here N_{\text{blocks}} is the life in repetitions of the specified load block.

The scalar stress from which the cycles are counted is von Mises, signed von Mises (the sign from the principal stress largest in magnitude) or \sigma_1. The result is the fields of damage D, \log_{10}N_f, equivalent amplitude and fatigue safety factor n = \sigma_a^{\text{allow}}(N_{\text{target}})/\sigma_{a,\text{eq}} — the ratio of the allowable amplitude at the target life to the equivalent one (capped from above), as well as the most critical node.

fatigue sigma_f=900e6 b=-0.09 sigma_e=200e6 sigma_u=600e6 correction=goodman r_ratio=0
probe fatigue_min_log_life

12.13Fracture mechanics: KI, J estimates and Paris crack growthexperimental

The module gives engineering estimates from closed-form linear elastic fracture mechanics formulas and from the stress at the location of the presumed crack. It is not a calculation of the J integral over a contour or domain around the crack tip, and not a calculation on a mesh with a crack.

Formulas.

• K_I = Y\,\sigma\sqrt{\pi a}; for a central crack 2a in an infinite plate Y=1 and K_I = \sigma\sqrt{\pi a}.
• The energy release rate (Griffith–Irwin): G = K_I^2/E', where E' = E (plane stress) or E/(1-\nu^2) (plane strain).
• Fatigue crack growth by Paris: da/dN = C\,(\Delta K)^m, \Delta K = Y\,\Delta\sigma\sqrt{\pi a}. For constant Y the integral is taken analytically: N = \dfrac{a_f^{\,1-m/2} - a_0^{\,1-m/2}}{(1-m/2)\,C\,(Y\Delta\sigma\sqrt\pi)^m} (for m=2 — N = \ln(a_f/a_0)/[C(Y\Delta\sigma\sqrt\pi)^2]).
• Crack initiation criterion: the crack is unstable at K_I \ge K_{Ic}; the critical length is a_c = \dfrac1\pi\left(\dfrac{K_{Ic}}{Y\sigma}\right)^2, the margin is n = K_{Ic}/K_I.
• For DCB (mode I) and ENF (mode II) specimens the module has beam formulas for compliance, G and critical load; they are used in the cohesive zone verification cases.

Edge-based estimate (fracture edge number N): the user specifies an edge (the location of the presumed crack), the initial and final lengths a_0 and a_f, the Paris parameters C, m, the toughness K_{Ic} and the shape factor Y. At the edge nodes the opening stress \sigma_{\text{open}} is taken — the largest principal stress from a run of the model without a crack. Then K_I = Y\,\sigma_{\text{open}}\sqrt{\pi a_0}, J = K_I^2/E', and the number of cycles is obtained by integrating the Paris law with \Delta\sigma = \sigma_{\text{open}} (a "zero to maximum" cycle). Units: K_{Ic} in MPa·√m, C for \Delta K in MPa·√m and crack length in metres; in the results panel a_0 and a_f are entered in millimetres, in the DSL — in metres.

fracture edge number 7 a0=0.002 af=0.008 c=1e-11 m=3 ki_c=50
probe sif_k1
probe paris_cycles

12.14Field calculator: user expressionsstable

The engineer defines a custom field — an expression over the nodal result fields — and it is colour-mapped on the model like a built-in one: with a legend, a probe, a value filter.

An expression is built from variables: stress components sxx … szx, vm, s1 s2 s3, strain components exx … gxz, eeq, epeq, displacements ux uy uz umag, as well as numbers, constants and functions (abs sqrt cbrt sign exp log log10 sin cos tan atan atan2 min max hypot clamp pow, arithmetic and power). Parsing is safe (no arbitrary code is executed).

Dimensions are checked at parse time. Inside an expression everything is in SI (Pa, m). Units are named quantities: the allowable is written 250*MPa, not 250e6. The expression 250e6/vm is rejected (dimension 1/Pa), as is vm + 5 (pressure cannot be added to a number). The dimension of the result is derived from the expression: Pa → "pressure" (displayed in MPa, psi or Pa according to the unit system), metre → "length", 1 → dimensionless; other powers cannot be displayed, and the expression is rejected with the name of the dimension. The optional unit= is an assertion about the dimension: if it does not match, the expression is refused.

Safety margin by a custom criterion (a separate field type): K = [\sigma]/f(\boldsymbol\sigma), where f is a criterion in units of pressure (for example vm or max(s1, -s3*ft/fc)), and [\sigma] is a constant allowable. At f=0 the margin is capped by the ceiling; at f<0 or a non-numeric value the node is undefined.

A node where the value is undefined receives "not a number" and is counted: their number and the first cause are given in the report; a silent zero is not substituted.

field define name=sigma_ratio expr="vm / (250*MPa)" unit=none
field safety name=n_shell criterion="vm" allowable="250*MPa"
probe n_shell

12.15Probes, filters and measurementsstable

In the interface a probe is placed by clicking on the model; the panel shows the value of the field at the node. Auto-probes mark the minimum and maximum of the selected field on the model. The value filter leaves on the model only the nodes outside or inside a given range. The range of the colour scale can be fixed manually (in the display units of the selected system).

In the DSL:

probe displacement_magnitude at=0,0,100
minmax vonMises
legend range min=0 max=250
result field vonMises

probe <quantity> at=… averages the field over the nodes that fall into the tolerance window tol= around the given coordinates (without interpolation; if one coordinate is given, the result is the mean over the slice); if there are no nodes in the window, the command refuses and asks to widen the tolerance. minmax returns the minimum and maximum over the whole nodal array (the value of the expression is the maximum). The available quantities are: displacements, vonMises, principalMax, principalMin, tensor components, local material stresses, reactions (reaction_x, reaction_y, reaction_z, reaction_magnitude), results of linearization, fatigue, fracture mechanics, optimization, as well as dynamic and frequency quantities. The full list is on the page /docs/dsl.

12.16Parametric DOE and surrogate optimizationexperimental

The "Parametric DOE and optimization" section is assembled from three units; none of them is yet connected to the calculation of a finite element model.

• Design of experiments. Latin hypercube sampling (LHS) with selection from several candidates by the maximum minimal distance between points (maximin); a deterministic generator with a fixable seed. The mode called "Sobol" is implemented as a Halton sequence (radical inversion over prime bases) — low-discrepancy, but not Sobol.
• Surrogate model. Ordinary kriging (a Gaussian process with a constant trend): the correlation is Gaussian or exponential, the Cholesky factorization of the correlation matrix, the prediction of mean and variance, expected improvement (EI) for Bayesian optimization. The correlation lengths are chosen by a heuristic from the mean distances between points, not by likelihood maximization.
• Constrained optimizer. An initial set of DOE points, then a local pattern search (coordinate-wise, halving the step when there is no improvement) on the penalty function f(\mathbf p) + c\,v(\mathbf p)^2, where v is the largest constraint violation. The output is the optimum, the Pareto front over the feasible points, the history and the finite-difference sensitivities \partial f/\partial p_k at the optimum. This is not an SQP method: it has no gradients or quadratic subproblems.

Constraints are set on stress, frequency, the buckling factor and displacement; targets and entity addresses follow the same numbering rules as the rest of the product (document numbering, no bbox addresses).

doe param name=thickness min=1 max=10
doe run method=lhs samples=20 objective=mass
probe doe_feasible_count

12.17Export of resultsstable

• VTK (legacy ASCII, UNSTRUCTURED_GRID) for ParaView: the nodes of the undeformed mesh, cells of all supported types (TET4, quadratic TET10, HEX8, prisms, pyramids), displacements as the vector field VECTORS (for Warp By Vector), nodal scalars (displacement magnitude, von Mises, \sigma_1, \sigma_3, safety margin, strains, contact, reactions in N) and temperature. The units are given in the field names. In the interface — the project menu, in the DSL — export vtk path=... (writing to a file at a path works in the desktop and scripting environments).
• Mesh in the format of solver packages: export mesh file="model.inp", optionally with the deformed geometry and a displacement scale.
• HTML report: run parameters, probes and a viewport snapshot; numbers are output in the chosen unit system.

12.18Legend and deformed viewstable

The colour scale on the model takes its range from the field (minimum and maximum over the nodes), and its label from the unit table of the chosen system; the range can be fixed manually (a double click on the ends of the scale or legend range min= max=). For contact the scale is categorical (status), and for the safety margin the values above the ceiling are labelled "infinity". The deformed view draws \mathbf{x} + s\,\mathbf{u} with a scale s chosen by the user; at s = 1 the true deformation is shown, at a larger s a conventional one, so take numbers from probes and the legend, not from the picture.

The same analysis in WebCAE can be carried out in three ways: with the mouse in the interface, with a script in the built-in DSL, or by delegating it to an external AI assistant through MCP. All three ways work with the same model and the same commands, so the result does not depend on how it was obtained. This section describes the workflow, the project file, the scripting language, the connection of assistants, the ready-made examples, data exchange and localization.

Why this matters to the engineer. A script turns an analysis into a reproducible document: it can be placed next to the report, rerun on a new version and the numbers compared. The project file stores the model together with the results, not only the input data. A detailed reference of DSL commands is maintained separately: DSL scripting language.

13.1Workflow: six stepsstable

The interface is organized as a chain of six steps. Each step has a tab in the top bar and its own parameter panel.

The analysis type at step 5 is chosen from the groups: structural (statics, modal analysis, buckling, harmonic response, transient dynamics), thermal, explicit dynamics and CFD. Models of beams and shells defined by the frame commands go through the same steps, but at steps 3 and 4 they show the frame panel instead of the usual panels.

The steps do not form a rigid pipeline. One can return to any available step. The steps "Geometry", "Model" and "Loads and BCs" are always open; "Mesh" and "Calculation" wait for a geometry or a mesh, and "Results" opens when there is a result of a run. Changes to the model can be undone and redone (Undo and Redo).

Autosave writes a snapshot of the model (without results) to the local storage of the browser or the application, to survive a tab reload or a crash. This is a safeguard, not a replacement for saving the project to a file.

13.2The .wce project filestable

A project is saved as a single file with the extension .wce. The former extension .caeweb continues to open; on the next save .wce is offered. The extension says nothing about the content: the reader determines the format by the signature and compression, so a renamed file opens in the same way.

Inside is the complete model, not only the input data.

• Geometry: the bytes of the source STEP, the validated B-Rep of each body, created and imported bodies, a snapshot of the topology with stable identifiers, the log of entity numbering.
• Model: materials and their assignment, sections, groups, boundary conditions, loads, contacts and constraints, mesh and analysis settings.
• Mesh: tetrahedra, hexahedra, hybrid elements and surface meshes with element quality.
• Results: statics, thermal field, explicit, harmonic response, frame results (statics, mode shapes, buckling, dynamics), the history of load steps, the archive of adaptive rebuilds.

Arrays of numbers (nodes, elements, fields) are stored in a binary block without re-encoding to text, and the whole container is compressed with gzip. The file records the data schema version. Old files are opened by a new version of the application: missing fields take defaults, obsolete entity addresses are converted to numbering on reading. A file of a newer schema is refused by an older build rather than read with data loss.

The following are not saved: the transient state of the interface (the current tab, statuses, selection), the partial series of an interrupted run, scripts and the terminal log. A project is opened from the "Project" menu, by dragging the file into the window, by a double click on the desktop (the application registers for .wce) and by the cae_open_project command for assistants. The "Project" menu also saves results separately from the model and reads them back.

13.3The unified command layer and the action logbeta

Any feature of the product is introduced as a command of the unified command layer. A ribbon button, a context menu item, a DSL line and an assistant call all invoke the same function. A form in the interface and a script do not duplicate logic, so the same operation gives the same result on any path. This is a development rule: a new capability is considered finished when it is reachable both from the interface and from a script.

Three surfaces are available. The first is the interface. The second is the DSL language, line by line. The third is the programming interface cae.* in JavaScript, through which the terminal and the MCP bridge work. In the terminal at the bottom of the window, a line beginning with cae. is executed as JavaScript; other lines are tried first as JavaScript, then as DSL. Several lines can be pasted at once.

The action log is built automatically. The terminal shows two channels: commands of the calculation core (geometry, mesh, materials, BCs, starting a run; always shown) and interface events (tab change, selection, visibility, camera; enabled by the "UI events" checkbox). Commands are recorded in the form of cae.* calls, for example cae.fix(...) or cae.setMeshOrder(...). Such a line can be copied and executed again, that is, a macro can be assembled from the log.

13.4Reproducibility of resultsbeta

To repeat an analysis, a script, the program version and the platform are needed. On one platform, rerunning the same script gives the same mesh and the same numbers: the mesher is deterministic. Between the desktop and the browser, the meshes on the same geometry may differ, so numbers from different platforms are compared with the platform label (see platform divergences).

A result carries a "solver epoch" label. When a bug is fixed in the solver because of which the numbers of a stored result could have gone outside the 5 % tolerance, the epoch number increases. On opening a project saved before such a fix, the program shows a persistent notification "results may be outdated". The result itself is not deleted: it can be compared with a recomputation. The same verdict is given by the report, the DSL (the staleResults field in the interpreter response) and MCP.

The product version is visible in the panels, in the "About" window and in "Settings → Updates", in the cae_get_state response and in the cae.version property. Clicking the version number opens the change history.

13.5The DSL scripting languagestable

The DSL is a text command language. A script is an ordinary UTF-8 file, one command per line; comments begin with # or //, blocks and loops are enclosed in { }. Each command is a word followed by positional and named arguments of the form key=value. Vectors are written as numbers separated by commas: origin=0,0,10. Numbers allow an exponent: E=200e9.

The language contains variables (let), conditions (if/else), loops (for x in [0..3], while), functions (func/return) and built-in functions (abs, sqrt, min, max, floor, ceil, round, len, str, num, basename). Output is print; writing values into the structured response is record. To protect against infinite loops, the number of execution steps is limited to one million, and the recursion depth to 64.

The commands cover the whole chain: geometry (primitive, boolean, sketch, extrude, import), material and sections (material, section, assign), mesh (mesh, refine), boundary conditions (fix, load, contact, remote), analysis (solve static|modal|buckling|thermal|nonlinear|...), results (probe, minmax, reactions, hotspots) and output (export vtk). The length unit of a script is set by the units command; loads and material constants are given in SI.

Addressing entities

Faces, edges and vertices are addressed only by the numbers from the model tree: surface number N, edge number N, vertex number N. The kind word is mandatory. Directions by bounding box (face z-) and bare indices are rejected with a syntax error. The numbers are found with the command list faces of body1 or list edges of body1. In a script the result of list can be stored in a variable: each face has .n (the number), .area, .centroid, .normal. This allows selecting faces by a criterion and substituting the number into a command: fix surface number faces[0].n. More in the subsection Topology numbering.

A script example from geometry to result

A cantilever bar with a 20×20 mm section and a length of 100 mm: fixed at one end, tension by a force of 5000 N at the other. Second-order tetrahedra, linear statics, then a check by the formula \delta = FL/(EA).

# Cantilever bar 20x20x100 mm, 5 kN tension along the z axis
units mm
material E=200e9 nu=0.3 density=7850
primitive box dx=20 dy=20 dz=100
mesh tet size=8 order=2
fix surface number 5
load surface number 6 force=5000 component=z
solve static
reactions
minmax vonMises
minmax displacement_magnitude

In this example face 5 is the end at z=0, and face 6 is the end at z=100 mm. For a box, the faces are numbered in the order x-, x+, y-, y+, z-, z+; for other bodies, look up the numbers with the command list faces of body1 before placing the BCs. The expected result: the total reaction equals the applied force in magnitude (5000 N), and the largest displacement of the end is close to the analytical \delta = FL/(EA) = 5000\cdot 0.1/(2\cdot 10^{11}\cdot 4\cdot 10^{-4}) \approx 6.25\cdot 10^{-6} m. The clamp fully prevents displacements at the end, so a local stress raiser arises near it, and the largest von Mises stress will exceed the uniform F/A = 12.5 MPa. For comparison with the analytical value, take values away from the clamp.

A script is run in the terminal of the application or from MCP with the cae_run_script command; the output of print and the values recorded by record are returned together with the execution status.

In a script, the result of list can be stored and iterated in a loop, for example to print the numbers, areas and normals of faces.

let faces = list faces of body1
for f in faces {
  print f.n f.area f.normal
}

13.6MCP server for AI assistantsbeta

MCP (Model Context Protocol) is an open protocol by which an assistant (Claude Desktop, Claude Code, Cursor and similar) calls tools. WebCAE ships a local bridge that connects the assistant to the open application window. The assistant works with a real model: it opens a project, inspects faces and edges, assigns a material, places boundary conditions, builds a mesh, starts a run, and receives the numbers and a screenshot.

The tools are typed and divided into groups: inspection (cae_inspect_surfaces, cae_inspect_edges, cae_get_state, cae_get_model_context), construction (cad_create_from_template, cad_create_from_script, cae_import_step_artifact, cae_open_project), model definition (cae_set_material, cae_add_boundary_condition), calculation and viewing (cae_solve, cae_screenshot), and the general-purpose cae_run_script, which accepts a DSL script as a string and returns the print output, the values recorded by record and the statuses. The targets of boundary conditions are addressed by the same numbers surface number N; axis pseudo-names ("top face") are rejected by the bridge. The list of tools and parameters: MCP server.

The bridge listens only on the local address. In the release desktop build it is off by default: the connection has no authorization, and a permanently open port would give the model to any process on the machine. It is enabled by the "MCP" button in the top bar or at application launch. The active application window is one: if several clients are connected, commands go to one window and the others wait.

13.7Library of examplesbeta

The application includes ready-made projects with computed results: about ninety in release 0.2.1 and more than a hundred in the working branch. They open from the "Examples" menu of the top bar, from "Settings → Open example" and by a link of the form ?example=<id>. A project opens immediately with result fields, legend and probes; there is no need to repeat the run. The examples are grouped by industry (automotive, mechanical engineering, rail transport, mining, oil and gas, power engineering, aerospace, construction, shipbuilding) and by analysis type: statics, modal analysis, buckling, large deflections, contact, RBE2 constraints, CFD (RANS), frame and shell structures. The showcase on the website shows the same models.

A separate set consists of nonlinear and multiphysics cases: braking of a pad with friction, seismic loading of a tower with RBE2 constraints, dynamic snap-through of an arch, Joule heating of a busbar with thermoelastic shock, and, in the working branch, also a hyperelastic silent block, a turbine disc sector with creep and fatigue, delamination of a two-layer DCB beam, and an acoustic cavity of a muffler. For each one a build script and a numerical passport are stored. The cases added after release 0.2.1 are absent from it.

13.8Data exchange: projects, results, geometrybeta

The "Project" menu gathers the exchange operations.

• Geometry: import of STEP and IGES (see STEP import and IGES and STL import), export to STEP and IGES when geometry is present.
• Project: "New project", "Save project…", "Open project…" (the .wce format, see above).
• Results: "Export results…" and "Import results…" work with .wce and .caeweb files, and import also accepts .vtk. This makes it possible to pass the result fields to a colleague without passing the whole project.
• VTK: export of the mesh with fields to the .vtk format read by ParaView. For explicit and CFD runs the VTK produced by the calculation module is used.
• DSL: export vtk [path=...].

On the desktop the file dialogs are the system ones. In browsers with the File System Access API (Chrome, Edge) a project is saved directly to the chosen file; in the others (Safari, Firefox) it is downloaded as an ordinary file.

13.9Mesh exchange and export of the deformed shapeexperimental

The working tree, but not the released version 0.2.1, adds the exchange of finite element meshes and the export of results to CAD-like formats. Import and export of Abaqus (.inp), Gmsh (.msh) and Nastran (.bdf) meshes, and export of the deformed shape and of a stress isosurface to STL and STEP, are opened by "Project" menu items when a mesh or a displacement field is present. More on mesh import and its limitations is in the subsection Mesh import and export.

13.10Localization and unit systemsbeta

The interface is translated into four languages: English, Russian, Spanish and Chinese. The language is determined by the page address, the lng parameter, the saved setting and the browser language; it is switched in the settings. English is the source language: if there is no translation, the English text is shown rather than an internal key.

The completeness of translation is uneven. The Spanish and Chinese locales cover almost all English strings (about 98 %). The Russian locale as of the end of October 2026 is translated at about three quarters: several panels remain untranslated, including part of the common labels, the material and section editors, planes and coordinate systems, and the tree context menu. These places are shown in English.

Units are chosen independently of language. Three input systems are available: SI (m, N, Pa, kg), metric engineering (mm, N, MPa, °C) and imperial (inch, pound-force, psi, °F). Unit labels are taken from a reference table: in the Russian locale, Cyrillic designations per GOST 8.417; in the others, international ones. At the solver boundary it is always SI: the conversion is done by a single conversion module, and the geometry is stored in millimetres.

This section answers three practical questions. Where to run an analysis: in the desktop application or in the browser, and how the results differ. What size of problems the program solves and in what time. Why the numbers can be trusted and what exactly is checked. The numbers here were measured on specific machines, and each carries a platform label. Numbers from different platforms must not be compared without a label.

There is one principle: the reference is the native desktop path. Verification, performance measurements and the analysis of discrepancies are taken on it. The browser version computes with the same core compiled to WebAssembly and serves for acquaintance and small models. If the results of the platforms differ on one model, the native result is considered correct, and the discrepancy is described explicitly.

14.1Desktop and browser: what is computed wherestable

The desktop application is released for Windows (64-bit, NSIS installer, per-user installation) and for macOS on Apple processors (a DMG image; there is no image for Intel processors). There is no Linux build. The browser version works in modern browsers with WebAssembly support; a project is saved directly to the chosen file in Chrome and Edge and downloaded as an ordinary file in Safari and Firefox.

A practical rule follows for the engineer. The browser is suitable for viewing a ready-made example, assembling and checking a small model, and testing a script. Work models, assemblies of dozens of bodies, millions of degrees of freedom, CFD and everything that must be reproducibly repeated should be computed on the desktop.

The main desktop analyses run by default on the native path. If the native path fails (no binary, a wrong response, size exceeded), the program reports the cause and does not substitute the browser path: a silent change of platform would give different numbers without warning.

14.2Typical problem sizes and timesbeta

Measurements of the native path are collected below. Each table is labelled with the machine. These are orders of magnitude for planning, not guarantees: the time depends on the processor, the machine load, the type of problem and the mesh quality. The measurements were taken at different times on different machines and cannot be merged into one table.

Linear statics of solids, native solver, 8 threads, 28 GB RAM, measured 2 Sep 2026. Second-order tetrahedra, CG with an algebraic multigrid preconditioner (AMGCL, Gauss–Seidel smoother, rigid body motion basis in the aggregation). The time is total: building the preconditioner and the solve.

The number of iterations barely grows with size: the multigrid preconditioner makes convergence practically independent of the mesh step. On a real CAD model of an antenna (741 thousand degrees of freedom, Tet10), the same method on 8 threads solves in 10.7 s, on one thread in 19.8 s; sparse direct Cholesky (Eigen) on one thread takes 843 s and 6.7 GB; on macOS the Apple Accelerate direct solver is noticeably faster. On small models (tens of thousands of degrees of freedom) the direct method is faster and more reliable than the iterative one, so the automatic selection keeps it while the model fits within the memory limit.

MITC4 shell statics, native solver, Apple M5, 24 GB, measured 13 Oct 2026. A square plate 1000×1000×10 mm, clamped along the edge, uniform pressure; the time from pressing "Solve" to the result, including model preparation and data transfer.

The solver alone takes 6–7 GB at 2 million degrees of freedom, but the whole application together with the interface at the same size takes about 15.7 GB (at 1 million about 9.6 GB), and on a machine with 24 GB swapping begins at 2 million. The size ceiling of shell statics therefore depends on physical memory: currently of the order of 1.8 million degrees of freedom on 24 GB and of the order of 1 million on 16 GB. Above the ceiling the run refuses before starting, with the name of the cause, and does not switch to the browser path.

Mesh generation, native Netgen mesher, Apple M5 (4 performance and 6 efficiency cores), 24 GB, measurements of August and October 2026 on a loaded machine.

The mesher speedup comes from simultaneously built bodies (by the number of performance cores), not from threads within one body: ten threads on one body gave only 1.2×.

What is known about the browser. The problem size is limited by 4 GB of linear memory of each module and by one solver thread. According to measurements of July 2026, the practical ceiling for second-order statics is of the order of a million tetrahedra, and for second-order bonded contacts of the order of hundreds of thousands. The exact boundary depends on the problem type, the browser and the free memory of the tab. The reference has no comparable time tables for the browser: the native path is the reference.

14.3Platform divergencesbeta

The desktop and the browser are built from the same source code but do not match bit for bit: the native build uses processor instructions with fused multiply–add and the system math library, whereas the browser build does not. Differences in the last bits of floating-point numbers sometimes lead to noticeable consequences, especially in the mesher.

• Meshes may differ. On the assembly of 271 bodies, of the 269 bodies that both platforms build, 115 (43 %) received identical meshes, while for the rest the number of tetrahedra differs, by 11.9 % in the worst body. In total the number of tetrahedra matched to within 0.01 %, because the differences in opposite directions compensate each other. Each platform on its own is deterministic: a rerun gives the same mesh.
• Failures differ. There are bodies for which the mesh is not built in the browser but is built on the desktop (one of the bodies of the same assembly, 20,599 tetrahedra on the desktop), and vice versa (another body of the same assembly). The cause is a divergence of the meshing algorithm on spline surfaces under slightly different arithmetic, not memory corruption. The practical conclusion: a mesher failure on one platform does not prove a geometry defect; try the other platform or change the element size.
• The mesh quality of the native platform is not always higher. On a test cube of 10 mm with an element size of 4 mm, the minimum SICN (signed inverse condition number of the element, see Mesh quality metrics) was 0.146 on the desktop and 0.182 in the browser. The native path remains the reference by decision, not because it is better everywhere: it matches what the desktop user gets.
• The numbers are affected by the mesh, not by the solver. The set of DSL cases with numerical acceptance passes entirely on the native mesher. When the mesher platform was changed, 111 of the 115 numerical measurements of the set matched bit for bit; the remaining four diverged by an amount from 10^{-13} to 5\cdot 10^{-2} percent, that is, at least a hundred times smaller than the tolerance.
• Features. CFD exists only on the desktop (and not in the Windows build). Verification cards taken in the browser are marked "WASM".

14.4Verification programmestable

The public verification programme is the Verification page. The catalogue holds 213 cases; 211 of them have an executed run, and all passed. Two MacNeal–Harder cases (straight and twisted cantilever beams) are published as "In development": the setup exists, but there is no run. For each executed case the page shows the setup, the reference, the obtained value, the deviation, the mesh and the solver version. For 207 cases out of 213 a DSL script is published: it can be pasted into the terminal of the application and the calculation repeated.

Examples: NAFEMS LE1, elliptic membrane, Hertz contact: sphere, thick-walled cylinder in a rigid casing (Lamé problem), buckling of a square plate, Lee frame with a limit point, cylinder in a flow at Re = 40.

There are three types of references. Analytical solutions: formulas of elasticity theory and dynamics for which the result is known in closed form. Published benchmarks: the NAFEMS and MacNeal–Harder sets. Invariants and balances: energy, force equilibrium, symmetry, independence of the rotation of the body or of the load step, where a numerical reference cannot be given.

14.5Tolerance policy ≤5 %stable

The tolerance of comparison with the reference is uniform: a deviation of no more than 5 %. The deviation is computed relative to the reference value:

where v_{\text{calc}} is the WebCAE result and v_{\text{ref}} is the reference. For cases with several quantities, the card gives the deviations for each, and the table shows one number per case.

The rule is strict: the tolerance is not relaxed to make a case pass. A deviation above 5 % is considered a defect and is eliminated in the calculation, not in the criterion. If a case cannot be brought to the tolerance, it is not closed with a wide tolerance but is published as "In development" or deferred with a reason.

Distribution of deviations over the 211 executed cases:

Deviations of hundredths of a percent are, as a rule, found in patch tests, balances, invariants and problems with exact analytical solutions. On problems with stress concentration, shells on a coarse mesh and thin bodies, the deviation is closer to the limit: the 5 % tolerance reflects the real accuracy of the method on the chosen mesh, not an excessive margin. So read the tolerance as an upper bound for the case, not as the accuracy of any model of yours.

14.6Automated regression checksbeta

Quality is held by layers of checks that are run during development. Orders of magnitude as of the end of October 2026:

• unit tests — of the order of five thousand checks in hundreds of files: formulas, matrix assembly, file parsing, commands, interface;
• type checking in two passes, separately for the product and for the tests: it catches errors that a test run misses, because the test execution environment does not check types;
• end-to-end script suites — about eighty files that run the chain from geometry to result;
• mesh corpus — 17 problems from real STEP files, native mesher; the result, the mesh quality and the connectivity are checked;
• the DSL case set — 74 cases, 115 numerical checks with tolerances, native mesher;
• tests of the native part of the application — of the order of two hundred and fifty;
• native checks of the core physics — about ninety: complex frequency responses, prescribed nonlinear paths, stress recovery in shells;
• verification cases — the same 213, with executable scripts.

A new test is accepted only after its ability to fail has been proven: a regression is returned into the product code, the test must fail, and then the change is rolled back. This protects against "green dummies" — tests that check nothing. A separate rule concerns the number of tests: the silence of a gate without the number of executed tests does not count as a result, because the suite might have stopped building.

14.7GPU and cloud modesexperimental

The calculation always runs on the processor, locally. The graphics card and a remote server do not take part in the calculation.

The code contains prototypes: a CG solver on WebGPU, a matrix-free variant for Tet10 and Hex20, and a calculation router between the browser, a local module and a remote server. They are covered by their own tests, but the product calculation does not call them, and there is no server side in the product. The calculation panel and "Settings → Compute" have corresponding switches ("WebGPU acceleration", "Matrix-Free P2 GPU", "Hybrid Cloud RPC"): at present they are saved in the project and do not change the calculation. More: Prototypes and unimplemented features.

14.8Versions and the result passportstable

The numbering is uniform for the application, the installers, the website and the verification cards: Semantic Versioning MAJOR.MINOR.PATCH, the release tag vX.Y.Z. The current version is 0.2.1. A bug fix increases the last digit; a new capability or new physics increases the middle one; while the major digit is zero, incompatible changes also go into the middle one. Planned versions on the documentation page are marked "planned" and have no dates.

Where to see the version:

• in the program: the version badge in the top and bottom bars, the "About" window, "Settings → Updates"; a click on the number opens the change history grouped by releases;
• from scripts: the cae.version property and the cae_get_state response in MCP return a string of the form WebCAE 0.2.1;
• in verification cards: the "Solver version" field.

The run passport in a card is recorded as WebCAE <version> · <short commit>. The +dirty mark after the commit means that at the time of the run the tracked files had uncommitted changes: such a number is exactly reproducible only together with those changes. Cards recorded before the introduction of uniform numbering are not rewritten, because their numbers were obtained by other builds. They show "before numbering was introduced (before 0.1.1)" or "version not recorded". Of the 211 executed cases, a passport with a version number currently exists for 14, "before numbering" for 28, and "version not recorded" for 169. The re-recording of old cards with the version recorded is not finished.

File compatibility. A project saved by an old version opens in a new one. A project of a newer schema is refused by an old version rather than read with data loss. If after an update the numbers of a stored result could have changed, the program issues the warning "results may be outdated" (see Reproducibility of results).

15.1Glossary

15.2Bibliography

Textbooks and monographs on which the formulations of this reference rely:

1. Bathe K.-J. Finite Element Procedures. 2nd ed. Watertown: K.-J. Bathe, 2014.
2. Zienkiewicz O.C., Taylor R.L., Zhu J.Z. The Finite Element Method: Its Basis and Fundamentals. 7th ed. Butterworth-Heinemann, 2013.
3. Hughes T.J.R. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Dover, 2000.
4. Belytschko T., Liu W.K., Moran B., Elkhodary K. Nonlinear Finite Elements for Continua and Structures. 2nd ed. Wiley, 2014.
5. Simo J.C., Hughes T.J.R. Computational Inelasticity. Springer, 1998.
6. Wriggers P. Computational Contact Mechanics. 2nd ed. Springer, 2006.
7. Crisfield M.A. Non-linear Finite Element Analysis of Solids and Structures. Vol. 1–2. Wiley, 1991, 1997.
8. Reddy J.N. Mechanics of Laminated Composite Plates and Shells. 2nd ed. CRC Press, 2004.
9. Clough R.W., Penzien J. Dynamics of Structures. 3rd ed. Computers & Structures, 2003.
10. Saad Y. Iterative Methods for Sparse Linear Systems. 2nd ed. SIAM, 2003.
11. Davis T.A. Direct Methods for Sparse Linear Systems. SIAM, 2006.
12. Timoshenko S.P., Goodier J.N. Theory of Elasticity. 3rd ed. McGraw-Hill, 1970.
13. Lekhnitskii S.G. Theory of Elasticity of an Anisotropic Body. Moscow: Mir, 1981.

Papers on individual methods:

1. Dvorkin E.N., Bathe K.-J. A continuum mechanics based four-node shell element for general nonlinear analysis. Engineering Computations, 1984, 1(1).
2. Batoz J.-L., Bathe K.-J., Ho L.-W. A study of three-node triangular plate bending elements. Int. J. Numer. Meth. Eng., 1980, 15.
3. Hughes T.J.R., Brezzi F. On drilling degrees of freedom. Comput. Methods Appl. Mech. Eng., 1989, 72.
4. Battini J.-M., Pacoste C. Co-rotational beam elements with warping effects in instability problems. Comput. Methods Appl. Mech. Eng., 2002, 191.
5. Felippa C.A., Haugen B. A unified formulation of small-strain corotational finite elements. Comput. Methods Appl. Mech. Eng., 2005, 194.
6. Zienkiewicz O.C., Zhu J.Z. The superconvergent patch recovery and a posteriori error estimates. Int. J. Numer. Meth. Eng., 1992, 33.
7. Dörfler W. A convergent adaptive algorithm for Poisson's equation. SIAM J. Numer. Anal., 1996, 33.
8. Celik I.B. et al. Procedure for estimation and reporting of uncertainty due to discretization in CFD applications. J. Fluids Eng., 2008, 130.
9. Popp A., Gee M.W., Wall W.A. A finite deformation mortar contact formulation using a primal–dual active set strategy. Int. J. Numer. Meth. Eng., 2009, 79.
10. Wriggers P., Zavarise G. On contact between three-dimensional beams undergoing large deflections. Commun. Numer. Meth. Eng., 1997, 13.
11. Camanho P.P., Dávila C.G. Mixed-mode decohesion finite elements for the simulation of delamination in composite materials. NASA/TM-2002-211737, 2002.
12. Benzeggagh M.L., Kenane M. Measurement of mixed-mode delamination fracture toughness. Compos. Sci. Technol., 1996, 56.
13. Cooper M.G., Mikic B.B., Yovanovich M.M. Thermal contact conductance. Int. J. Heat Mass Transfer, 1969, 12.
14. Craig R.R., Bampton M.C.C. Coupling of substructures for dynamic analyses. AIAA J., 1968, 6(7).
15. Der Kiureghian A. A response spectrum method for random vibration analysis of MDF systems. Earthquake Eng. Struct. Dyn., 1981, 9.
16. Miles J.W. On structural fatigue under random loading. J. Aeronaut. Sci., 1954, 21.
17. Hinton E., Rock T., Zienkiewicz O.C. A note on mass lumping and related processes in the finite element method. Earthquake Eng. Struct. Dyn., 1976, 4.
18. Knyazev A.V. Toward the optimal preconditioned eigensolver: LOBPCG. SIAM J. Sci. Comput., 2001, 23.
19. Vaněk P., Mandel J., Brezina M. Algebraic multigrid by smoothed aggregation for second and fourth order elliptic problems. Computing, 1996, 56.
20. Johnen A., Remacle J.-F., Geuzaine C. Geometrical validity of curvilinear finite elements. J. Comput. Phys., 2013, 233.
21. Economon T.D. et al. SU2: An open-source suite for multiphysics simulation and design. AIAA J., 2016, 54(3).

Standards and benchmark problem sets:

1. NAFEMS. The Standard NAFEMS Benchmarks. Rev. 3, 1990; Selected Benchmarks for Natural Frequency Analysis, 1987.
2. MacNeal R.H., Harder R.L. A proposed standard set of problems to test finite element accuracy. Finite Elements in Analysis and Design, 1985, 1.
3. ASME Boiler and Pressure Vessel Code, Section VIII, Division 2, Part 5.
4. ASTM E1049-85. Standard Practices for Cycle Counting in Fatigue Analysis.
5. VDI 2230 Blatt 1. Systematische Berechnung hochbeanspruchter Schraubenverbindungen.

15.3Open-source components

The computational core and model preparation use open-source libraries:

Component licences are listed on the Open source page.