2. Singularity detection

2.1. Principle of the method

When the exact solution of the problem studied has singularities, the order of convergence of the finite elements solution is modified and therefore also that of the error estimator. For example, consider a discretized plane elasticity problem with triangular elements of degree \(p\).

If the exact solution \({U}_{\mathrm{ex}}\) is regular, we know that ([bib2], [bib3]):

\({\parallel u-{u}_{h}\parallel }_{\Omega }={\parallel e\parallel }_{\Omega }\le C{h}^{p}\) eq 2.1-1

Where \({\parallel e\parallel }_{\Omega }\le C{h}^{p}\) is the contribution to the energy error, which is:

\({\parallel e\parallel }_{\Omega }\le \frac{1}{2}\underset{\Omega }{\int }\varepsilon ({e}_{h})K\varepsilon ({e}_{h})d\Omega\) eq 2.1-2

On the other hand, if the exact solution has a singularity, for example if, locally in the vicinity of a point \({M}_{0}\), the field of displacement is of the form (with \(r\) and \(\theta\) polar coordinates in the vicinity of point \({M}_{0}\)):

\({U}_{\mathrm{ex}}={r}^{\alpha }V(\theta )+W\text{avec}0<\alpha <1\) eq 2.1-3

So, we show that [Strang & Fix, 1976]:

\({\parallel {e}_{h}\parallel }_{\Omega }\le C{h}^{\alpha }\) eq 2.1-4

As a result, the convergence rate of the global energy error becomes independent of the degree \(p\) of the finite elements used and the same is true of that of the measurement of the error (for example, if \(p=1\text{ou}2\) then \(\alpha =1/2\) for a crack).

In order to obtain a good prediction of optimized meshes, the previous findings lead us to use a convergence rate \({q}_{E}\) per element such that the error estimator \({\varepsilon }_{E}\) verifies:

\({\varepsilon }_{E}=O({h}^{{q}_{E}})\) eq 2.1-5

A simple way to define these local coefficients consists in taking:

• \({q}_{E}=\alpha\) if the element \(E\) is connected to a singularity of order \(\alpha\);

• \({q}_{E}=q\) for all other finite elements where \(q\) only depends on the type of finite elements used.

The method presented below therefore comprises three phases:

• detection of singular areas, in this case the singular nodes of the mesh;

• numerical evaluation of the coefficient \({q}_{E}\) for the elements connected to the nodes considered to be singular (for the other elements, \({q}_{E}=q\) is then fixed);

• calculation of the size modification coefficient \({r}_{E}\).

2.2. Detecting singular nodes

The idea is to use local errors. In fact, experiments show that these local errors have a peak in the vicinity of a singularity. For each node \(i\) of the mesh, we therefore compare the average error \({\stackrel{ˉ}{m}}^{i}\) of the elements connected to node \(i\) to the average error \(\stackrel{ˉ}{M}\) over the whole structure. Node \(i\) is considered singular if:

\({\stackrel{ˉ}{m}}^{i}\ge \beta \stackrel{ˉ}{M}\) eq 2.2.1-1

with

\({\stackrel{ˉ}{m}}^{i}=\sqrt{\frac{\sum _{\text{E connecté à i}}{{\varepsilon }_{E}}^{2}}{\sum _{\text{E connecté à i}}\mathrm{mes}(E)}}\) and \(\stackrel{ˉ}{M}=\sqrt{\frac{\sum _{\text{E}\in \text{structure}}{{\varepsilon }_{E}}^{2}}{\sum _{\text{E}\in \text{structure}}\mathrm{mes}(E)}}\) eq 2.2.1-2

where \(\beta\) is a coefficient greater than 1 and \(\mathrm{mes}(E)\) is the 2D area or the 3D volume of element \(E\). Numerical experiments have shown that singular nodes are well detected by setting \(\beta =2\) in dimension 2, \(\beta =3\) in dimension 3 for linear finite elements and \(\beta =2\) in dimension 3 for quadratic finite elements.

Note 1:

From a numerical point of view, the detection of singular nodes differs between 2D and 3D cases. The conditions given below for this detection are not based on a particular theory but rather on the experience acquired in this field by the Mechanical Laboratory and CAO in Saint-Quentin. In 2D: a node i is considered singular if it meets the following 3 conditions:

Where \({\stackrel{ˉ}{m}}_{1}^{i}\), \({\stackrel{ˉ}{m}}_{2}^{i}\), and \({\stackrel{ˉ}{m}}_{3}^{i}\) are the error averages for the elements belonging to layers 1, 2, and 3, respectively, with respect to the node i in question.

Layers are defined as follows:

• Layer 1: elements that have the node i to be tested,

• Layer 2: elements in contact (face, edge or node) with an element from layer 1,

• Layer 3: elements in contact (face, edge or node) with an element from layer 2.

In 3D: the node \(i\) is considered singular if it meets the condition \({\stackrel{ˉ}{m}}_{1}^{i}\ge \beta \stackrel{ˉ}{M}\) and if one of the nodes connected to the node \(i\) in question meets the condition \({\stackrel{ˉ}{m}}_{i}^{\text{Noeud connecté à i}}\ge \beta \stackrel{ˉ}{M}\). Unlike the 2D case, node \(i\) is singular only if one of its neighbors is also singular (forget the isolated singular nodes to keep only the singular edges).

Note 2:

In 2D, only vertex nodes are examined. In 3D, only the vertex nodes located on one edge of the structure are examined (only for reasons of calculation time; this condition could therefore be modified).

2.3. Evaluating the order of the singularity

For each singular node \(i\) detected, the order of the singularity, i.e. the value of \({q}_{E}\) that will be used for the elements connected to node \(i\), is determined by identifying the value of the energy density of the finite elements solution in the vicinity of node \(i\) with the theoretical value in the vicinity of a singular point.

2.3.1. Dimension 2 case

In this case, we calculate the average finite element energy, in disks \(A\) with center \(i\) and radius \(r\):

\({\stackrel{ˉ}{e}}_{h}(r)=\frac{1}{2\mathrm{mes}(A)}\underset{A}{\int }\varepsilon ({u}_{h})K\varepsilon ({u}_{h})\mathrm{dA}\) eq 2.3.1-1

By identifying, by a least squares method, this average energy with the theoretical value in the vicinity of a singularity of order \(\alpha\):

\(e(r)=k{r}^{2(\alpha -1)}+c\) eq 2.3.1-2

numerically, a value \(\stackrel{ˉ}{\alpha }\) close to \(\alpha\) is obtained.

Figure 2.3.1-a: Numerical evaluation of \(\alpha\)

In practice, numerical experiments show that it is sufficient to perform the identification in an area corresponding to 3 layers of elements around the singular point and to evaluate \({\stackrel{ˉ}{e}}^{h}(r)\) for 5 to 8 values of \(r\) regularly distributed in this zone (we took 10 values).

2.3.2. Dimension 3 case

In 3D the situation is more complex. Singular points are not, in most cases, isolated and it is therefore common to be in the presence of singular edges. Consider, for example, the case of a cube embedded on one face and subjected to tensile forces: all the points of the edges of the embedded face are singular [Figure 2.3.2-b].

Figure 2.3.2-b: Built-In Traction Cube

In this situation, the evaluation of the average energy in balls \(A\) of increasing radius and centered on a singular node does not make it possible to identify \({q}_{E}\). In fact, as the radius increases, the extent of the singular zone contained in ball \(A\) increases and a rapid decrease of \({\stackrel{ˉ}{e}}_{h}\) [Figure 2.3.2-c] is not obtained.

Figure 2.3.2-C:Energy in concentric balls

When the singular points are not isolated, the coefficient \({q}_{E}\) must be identified by calculating the energy density in coaxial cylinders built on the edges whose ends have been considered singular [Figure 2.3.2-d] and cf. the remark of [§ 2.2].

Figure 2.3.2-d: Energy in coaxial cylinders

2.4. Extension to areas of stress concentration

In practice, we have found that the previous method, developed on cases with singularities, also makes it possible to correctly take into account areas with high stress gradients even if mathematically these areas do not correspond to singularities.