Theoretical Background

Theoretical background of SimSolid and its software implementation workflow, and comparison of other methods used in traditional FEA.

Overview of Initial Research

The Ritz-Galerkin method invented at the beginning of 20th century for the approximate solution of boundary value problems assumes that functions that approximate the solution are analytical functions defined on the whole domain of interest. In practical applications these functions were either trigonometric or polynomials which were infinitely smooth, i.e. they had an infinite number of derivatives. There were two main problems with such functions. Firstly, it was difficult or impossible to construct such functions that a priori meet essential boundary conditions on boundaries of arbitrary domains (in structural analysis the conditions appear as displacement constraints). And secondly, the equation system built on such functions was ill-conditioned and numerically unstable which did not allow solving real life problems with enough accuracy.

The Finite Element Method (FEM) that appeared in 1950s was just a different implementation of the classical Ritz-Galerkin approach, but it succeeded in solving of both – constraints and numerical instability issues because it consistently used functions with local supports called finite elements. Though locally the basis-functions of finite elements were infinitely differentiable standard polynomials, global basis-functions assembled from local polynomials were not smooth at all – even their first derivatives were discontinuous. The FEM’s success proved that requirements to the continuity of the approximation functions must be met only to a certain degree - just enough to provide finite energy when they are substituted into an energy functional of a boundary value problem. The spaces of such functions were introduced and investigated by Sobolev in the 1930s.

The next step in the relaxation of continuity requirements on approximation functions was the introduction of the concept of external approximations [Equation 1]. The name “external” was used in the following context. When approximation functions belong to a Sobolev space of functions with finite energy the approximation is called “internal” which means that while the approximation is refined, and the solution is converging to the exact solution, the approximation functions are always inside the Sobolev space. Alternatively, in external approximations, the approximation functions do not belong to a Sobolev space at every refinement step (they have infinite energy), but in the limit, when number of degrees of freedom (DOF) approaches infinity, the limit function must belong to the corresponding Sobolev space, i.e. it must recover the necessary smoothness properties. The abstract theory of external approximations was developed in reference [Equation 2].

The technological foundations of SimSolid have been published in reference [3]. It develops the abstract theory of external approximations. In reference [2] it was applied to the particular case of approximations by finite elements under the assumption that the elements are of absolutely arbitrary shape. In the result the necessary and sufficient condition of external approximations by finite elements has been established and convergence theorems have been proven. It was also shown that the theorems were constructive, i.e. they not only defined hallmarks of external approximations, but also provided a mechanism to build them.

Theoretical Background

An abstract boundary value problem is formulated as to find a function

U

which fulfills the equations:(1)

A U = h inside the domain Ω

(2)

L U = g at the domain boundary Γ

where

A

and

L

are differential operators.

Some boundary value problems can be equally formulated in a variational form such as to find a function $U$ which provides a functional $F (U)$ at minimum value, where the functional $F (U)$ is usually an energy functional.

In 1908 W. Ritz proposed a method of finding an approximate solution of a boundary value problem by approximating it with a linear combination of some basis-functions(3)

U_{h} = \sum_{i}^{n} a_{i} p_{i}

where

a_{i}

are unknown factors, and

p_{i}

are basis approximation functions.

The factors

a_{i}

are found by minimizing the energy functional(4)

F (\sum_{i}^{n} a_{i} p_{i}) = \min

If the boundary value problem is linear, then the minimization problem (Equation 4) can be reduced to a linear algebraic equation system with respect to the factors

a_{i}

(5)

K d = f

here K is a symmetric matrix, d is the vector of unknown factors , and f is a right-hand side of the system.

In FEM the matrix K is called a stiffness matrix, the vector f is called a load vector, and factors $a_{i}$ are called degrees of freedom.

In 1915 Galerkin proposed another approximate method of solving the boundary value problem (Equation 1)-(Equation 2). According to the Galerkin method the unknown solution U is approximated as(6)

U_{n} = U_{0} + \sum_{i}^{n} a_{i} p_{i}

Where $U_{0}$ is some function which fulfills nonhomogeneous boundary conditions (Equation 2), $p_{i}$ are analytical approximation functions which fulfill homogeneous boundary conditions, $a_{i}$ are unknown factors.

Substitution of (Equation 6) into (Equation 1) results in the residual(7)

R = A U_{0} + \sum_{i}^{n} a_{i} A p_{i} - h

The unknown factors

a_{i}

are found from the equation system(8)

\int_{Ω} R p_{i} d Ω = 0

If a boundary value problem is linear, then system (Equation 8) is a system of linear algebraic equations.

The Galerkin method does not use a variational formulation of the boundary value problem. Therefore, its applicability is much wider.

Ritz and Galerkin methods proved to be effective means of solving problems in engineering and science. At the same time mathematical justification of the methods faced significant difficulties which were solved with the introduction of functional analysis as a mathematical discipline.

Modern theory of the Ritz-Galerkin method is based on the concept of the weak formulation of the boundary value problem. The weak formulation of a boundary value problem consists of finding a function from

u \in V

a corresponding Sobolev space which fulfills an abstract variational equation(9)

a (u, v) = f (v) for any function v \in V

here

V

is some subspace of a Sobolev space,

a (u, v)

is generally an unsymmetrical bilinear form which is continuous on the space product

V \times V, f (v)

is some linear form in

V

In structural analysis, the Sobolev space is a space of functions with finite strain energy.

In the Ritz-Galerkin method the space $V$ is approximated with some finite-dimensional space $X_{h}$ , and the approximate solution is found in form (Equation 3) where the functions $p_{i}$ belong to the space $X_{h}$ . Therefore, the discretized formulation of a boundary value problem is:

Find a function

U_{h} \in X_{h}

which fulfills the equation(10)

a (U_{h}, V_{h}) = f (V_{h}) for any function V_{h} \in X_{h}

Substitution of (Equation 3) into (Equation 10) results in a linear algebraic equation system from which unknown factors $a_{i}$ are found.

In the classic Ritz-Galerkin method $X_{h}$ is a space of analytical functions defined on the whole domain $Ω$ , the factors $a_{i}$ have no physical meaning. In conventional Finite Element Method $X_{h}$ is a space of piecewise polynomials and factors $a_{i}$ are values of the function $U_{h}$ in the nodes of finite elements. In structural analysis they are displacements of the nodes.

Many modifications of Ritz-Galerkin methods have been invented. They differ by the variational equations (Equation 9) and by the classes of basis-functions (Equation 3) used to approximate the solution. The same boundary value problem can have several equivalent formulations (Equation 9) which differ by the spaces $V$ .

External Approximations by Finite Elements

As already mentioned, internal finite element approximations are built on functions that belong to a corresponding Sobolev space. These functions must meet certain continuity conditions on inter-element boundaries. For instance, when 2D or 3D theory of elasticity problems are under consideration, the functions need to be continuous between finite elements. For plate bending problems not only the functions, but their first derivatives must be continuous as well.

The continuity conditions are quite restrictive. They can be met only for very simple shapes of finite elements using standard interpolation polynomials as basis-functions of finite elements. The polynomials are associated with element nodes. To provide inter-element compatibility the same interpolation functions are used to represent finite element shape. In case of curved boundaries mapping onto a canonical element is used to provide the compatibility. The geometry of finite elements and their approximation functions are tightly coupled.

To improve the approximation quality of finite elements, researchers invented incompatible finite elements. In incompatible element the interpolation basis-functions of the elements of standard shape are enriched by some other polynomials. The additional functions create discontinuity across inter-element boundaries, but incompatible finite elements often provided much better accuracy than the compatible ones. This resulted in difficulties of mathematical proof of convergence and in inconsistencies of results.

A comprehensive theory of external approximations by finite elements was developed in reference [Equation 3]. In the theory the word “finite element” was used to designate an arbitrary shaped sub-domain of the domain $Ω$ , so the definition of finite elements was not restricted anymore to canonic shapes or other shapes obtained from a canonic by mapping. The whole domain $Ω$ could be considered as one finite element, and therefore, for assemblies a part of an assembly could be one “finite element” in FEM terminology. Another assumption was that approximation functions inside the finite element could be arbitrary - not necessarily polynomials. The only requirement was that the functions belong to the corresponding Sobolev space, so they need to be sufficiently smooth inside the element.

The task was to find conditions under which the approximations built according to the assumptions above would be external approximations, i.e. they would converge to the exact solution of a boundary value problem from “outside” of a Sobolev space. The necessary and sufficient condition which provides the external approximations was found. The condition happens to be constructive – its formulation also implied the way of building finite elements that meet the condition. Convergence theorems and error estimates have been proven.

It was shown that the necessary and sufficient condition for a finite element approximation to be external is:(11)

〈 δ, γ U 〉 = 0

Here <,> is the duality pairing in certain functional spaces defined on inter-element boundaries, $δ$ and $γ$ are some operators, and $U$ are approximation functions defined inside the element.

As one can see, condition (Equation 11) does not relate to a Boundary Value Problem (BVP) formulation, or to a solution method for the BVP (Galerkin, Ritz, Trefftz, etc.). It imposes constraints on basis-functions of finite elements which guarantee that the limit approximation function belongs to the corresponding Sobolev space, so it has the necessary smoothness properties.

Therefore, even before the solution method is chosen (Galerkin, Ritz, etc.), one may construct finite element spaces that have important properties. These properties can be just “good to have”, as, for instance, when solving elasticity problems, it is not required to use functions that fulfill the equilibrium in volume, but it might be useful because the use of such functions increases accuracy and reduces the number of DOF. Other properties can be crucial, for instance, only divergence-free functions can provide unconditionally stable solutions for incompressible materials.

Then condition (Equation 11) can be extended by continuity from the duality pairing into the inner product in other spaces of functions:(12)

(g, γ U) = 0

Here

g

are functions defined on inter-element boundaries, they are called boundary functions. Boundary functions are functions of surface parameters and they generate boundary DOF that are integrals of products of boundary functions onto finite element basis-functions over the finite element boundary:(13)

\int_{Γ} g_{k} γ U d Γ = 0

Here $Γ$ is the boundary of the finite element, $g_{k}$ are functions defined on the boundary of the finite element, and $U$ is a function to be approximated on the element (for example displacements in structural analysis).

For comparison, the degrees of freedom in the FEM are the value of the function

U

in the node

i

of a finite element:(14)

U (x_{i}, y_{i}, z_{i})

The functions $g_{k}$ in expression (Equation 13) are basis-functions from finite dimensional space $G_{h}$ of functions defined on element boundaries. They can be arbitrary, the only requirement is that the spaces $G_{h}$ must be dense in the space of boundary functions, i.e. they must be able to converge in the space of boundary functions. The latter is easily fulfilled in case $g_{k}$ are polynomials or piecewise polynomials defined on the element boundaries.

The functionals (Equation 13) are called boundary degrees of freedom. They do not have physical meaning and they represent approximation functions from space of finite elements compatible when the number of boundary DOF converges to infinity. The boundary DOF are responsible for meeting inter-element continuity conditions and essential boundary conditions. In adaptive solution the number of the boundary DOF is managed automatically to meet the convergence criteria.

The boundary DOF (Equation 13) are not the only DOF produced when external approximations are built. Other DOF are called internal DOF because they are associated with the element volume. Internal DOF are defined automatically when the approximation of the solution within a finite element is being built. Finally the approximation of a function on the element looks as follows:(15)

U_{h} = \sum_{i}^{n} a_{i} (U) p_{i} + \sum_{k}^{N} (\int_{Γ} g_{k} γ U d Γ) p_{k}

Here $a_{i}$ are internal DOF of the element (some factors), $p_{i}$ are basis-functions of the internal DOF, $\int_{Γ} g_{k} γ U d Γ$ are the boundary DOF, $p_{k}$ are basis-functions of the boundary DOF.

The basis-functions $p_{i}$ and $p_{k}$ constitute a finite-dimensional space $P$ of approximation functions of a finite element. It was proven that for convergence the space $P$ must be complete, for instance, in case of polynomial space it should contain all polynomials up to a certain degree assigned to an adaptive iteration.

Basis-functions of a finite element are not pre-defined because the element has an arbitrary shape. They are built on-the-fly during a solution run. What is pre-defined at an adaptive pass is the whole space

P

of approximation functions of the element. The algorithm of building basis-functions of an element at an adaptive pass works as follows:

A set of boundary functions $g_{k}$ is defined;
A complete space $P$ of approximation functions of the element is defined by choosing a complete set of generic basis-functions. In the case of polynomial spaces, a complete space of polynomials of a certain degree is specified. For instance, generic second-degree polynomials for 3D problems are:
- ${1, x, y, z, x^{2}, x y, y^{2}, x z, z^{2}, y z}$
Generic basis-functions are generated automatically on-the-fly for every sub-domain during the solution when the stiffness matrix of a sub-domain is evaluated;
The basis-functions $p_{i}$ and $p_{k}$ are found automatically by solving a certain system of linear algebraic equations.

After basis-functions of an element have been found, the element stiffness matrix and load vector are evaluated the same way as it is done in conventional FEM by integrating energy over the element volume and loads over the element boundary.

Geometry-Function Decoupling

Geometry-function decoupling is the core feature of the SimSolid technology. As one can see from the above, the basis-functions of an arbitrary element are built from generic basis functions on-the-fly during the solution. Neither element geometry representation is used in building the generic functions, nor the functions dictate the shape of the element. The only requirement to the space $P$ of approximation functions of an element is that $P$ must be a subspace of a corresponding Sobolev space associated with the formulation of boundary value problems. Therefore, any combination of generic basis-functions is allowed provided they are linearly independent.

The geometry-function decoupling proved to be the key feature which provides better performance, better accuracy, robustness, less computer resources, less modeling errors. The following substantial benefits can be realized when finding an accurate solution for a specific problem, or managing adaptive solutions:

It is possible to build special approximations that make approximate solutions of boundary value problems unconditionally stable. For instance, when parts made of incompressible materials are simulated, SimSolid uses divergence-free functions which exactly meet the incompressibility condition. Here is an example of some generic divergence-free 3D functions of 3rd degree, where $(u, v, w)$ are displacement components:
- $\begin{array}{l} u = - x z^{2}, v = y z^{2}, w = 0 \\ u = - 3 x z^{2}, v = 0, w = z^{3} \\ u = - 2 x y z, v = y^{2} z, w = 0 \\ u = - x y z, v = 0, w = y z^{2} \\ u = - x y^{2}, v = 0, w = y^{2} z \end{array}$
Neighboring parts may have approximation functions of different classes. For instance, in case an assembly contains parts made of compressible and incompressible materials (rubber insertions, or cavities with liquid) the approximation functions for the incompressible material are built as special divergence-free functions. On neighboring parts with compressible material regular functions like standard polynomials are used.
It is always possible to use basis-functions that a-priori fulfill the governing equations of boundary value problems which provide better accuracy and reduce the number of DOF. For instance, thermo-elastic problems are solved using a complete polynomial solution of the corresponding governing equations:(16)
$\begin{array}{l} (λ + μ) \frac{\partial ε}{\partial x} + μ Δ μ = \frac{α Ε}{1 - 2 v} \frac{\partial T}{\partial x} \\ (λ + μ) \frac{\partial ε}{\partial y} + μ Δ μ = \frac{α Ε}{1 - 2 v} \frac{\partial T}{\partial y} \\ (λ + μ) \frac{\partial ε}{\partial z} + μ Δ μ = \frac{α Ε}{1 - 2 v} \frac{\partial T}{\partial z} \end{array}$

here $(u, v, w)$ are displacement components,

$\begin{matrix} ε = \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} \\ λ = \frac{E v}{(1 + v) (1 - 2 v)}, μ = \frac{E}{2 (1 + v)} \\ Δ = \frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}} + \frac{\partial^{2}}{\partial z^{2}} \end{matrix}$
Further, $α$ is the thermal expansion coefficient; $E$ is Young’s modulus, $v$ is Poisson’s ratio, and $T$ is the temperature field. The equation system (Equation 16) is non-homogeneous. For instance, when:
$E = 1, v = 0.25, α = 1$
and temperature field is described by a monomial:
$T = a x^{m} y^{n} z^{p}$
the solution of the non-homogeneous problem for a=1, m=0, n=2, p=3 is:(17)
$u = 0, v = 0.1667 y z^{5}, w = 0.4167 y^{2} z^{4} - 0.02778 z^{6}$
Here is an example of a polynomial solution of the homogeneous equations (14):(18)
$u = 20 x^{4} z - 20 x^{2} z^{3}, v = 20 x^{4} y z - 20 x y z^{3}, w = 8 x^{5} - 60 x^{3} z^{2}$

When solving a thermo-elastic problem, polynomial approximations of temperature $T$ are imported from thermal analysis, functions of type (Equation 17) are generated for every element, and generic functions of type (Equation 18) are used to build basis-functions of elements.

For heat transfer problems harmonic polynomials are used as basis-functions which precisely fulfill the corresponding equation of heat transfer. Here are some generic harmonic functions of degree 3:
$\begin{array}{l} f = x^{3} - 3 x z^{2} \\ f = x^{2} y - y z^{2} \\ f = x y^{2} - x z^{2} \\ f = y^{3} - 3 y z^{2} \\ f = 3 x z^{2} - z^{3} \end{array}$
The approximations are always built in the physical coordinate space without mapping onto a canonic shape. Therefore, the properties of generic basis-functions are preserved throughout the solution which eliminates a substantial source of approximation errors.
A complete set of basis-functions is always used to approximate solutions on a sub-domain. Completeness means that no functions are missing from a space of a certain degree. For instance, if the solution is approximated with harmonic polynomials of degree 5, then all harmonic generic polynomials of degree 5 are included into the approximation space of a sub-domain. This provides high accuracy, ease of building p-adaptive solutions globally and locally, and ease of implementation of new types of problem-specific basis-functions.
Geometry-functions decoupling allows effectively handle assemblies of parts with incomparable geometries in terms of size and shape (multi-scale assemblies).
Local effects like concentrated forces, cracks, stress concentration, etc., can be easily simulated by enriching the approximation space of sub-domains with special functions that possess corresponding characteristics associated with the feature.