HIGH-FIDELITY DIGITAL TWINS: DETECTING AND LOCALIZING WEAKNESSES IN STRUCTURES

Given that all materials exposed to the environment and/or undergoing loads eventually age and fail, the task of trying to detect and localize weaknesses in structures is common to many fields. To mention just a few: airplanes, drones, turbines, launch pads and airport and marine infrastructure, bridges, high-rise buildings, wind turbines, and satellites. Traditionally, manual inspection was the only way of carrying out this task, aided by ultrasound, X-ray, or vibration analysis techniques. The advent of accurate, abundant and cheap sensors, together with detailed, high-fidelity computational models in an environment of digital twins has opened the possibility of enhancing and automating the detection and localization of weaknesses in structures.
From an abstract setting, it would seem that the task of determining material properties from loads and measurements is an ill-posed problem. After all, if we think of atoms, granules or some polygonal (e.g. finite element [FEM]) subdivision of space, the amount of data given resides in a much smaller space than the data sought. If we think of a cuboid domain in $d$ dimensions with $N^{d}$ subdivisions, the maximum amount of surface information/ data given is of $O(N^{d-1})$ while the data sought is of $O(N^{d})$.
Another aspect that would seem to imply that this is an ill-posed problem is the possibility that many different spatial distributions of material properties could yield very similar or equal displacements under loads.
On the other hand, the propagation of physical properties (e.g. displacements, temperature, electrical currents, etc.) through the domain obeys physical conservation laws, i.e. some partial differential equations (PDEs). This implies that the material properties that can give rise to the data measured on the boundary are restricted by these conservation laws, i.e. are constrained. This would indicate that perhaps the problem is not as ill-posed as initially thought.
As the task of damage detection is of such importance, many analytical techniques have been developed over the last decades [7, 13, 18, 12, 23, 22, 21, 8, 20, 2, 5, 11]. Some of these were developed to identify weaknesses in structures, others (e.g. [13]) to correct or update finite element models. The damage/weakness detection from measurements falls into the more general class of inverse problems where material properties are sought based on a desired cost functional [4, 14, 6, 24]. It is known that these inverse problems are ill-defined and require regularization techniques.
The analytical methods depend on the measurement device at hand, and one can classify broadly according to them. The first class of analytical methods is based on changes observed in (steady) displacements or strains [22, 8, 2, 25, 11]. The second class considers velocities or accelerations in the time domain [12, 23, 20]. The third class is based on changes observed in the frequency domain [7, 13, 18, 21, 6].
Some of the methods based on displacements, strains, velocities or accelerations used adjoint formulations [27, 22, 8, 20, 3, 25, 17] in order to obtain the gradient of the cost function with the least amount of effort. In the present case, the procedures used are also based on measured forces and displacements/strains, use adjoint formulations and smoothing of gradients to quickly localize damaged regions [1]. Unlike previous efforts, they are intended for weakness/damage detection in the context of digital twins, i.e. we assume a set of defined loadings and sensors that accompany the structure (object, product, process) throughout its lifetime in order to monitor its state. The digital twins are assumed to contain finite element discretizations/models of high fidelity, something that nowadays is common the aerospace industry. Therefore, the proposed approach fits well into the overall workflow of high-level CAD environments and high fidelity FEM models seen in the design phase.

What follows relies on the following set of assumptions:

• -

  Monitoring the weakening of a structure is carried out by applying a set of $n$ different forces ${\bf f}_{i},i=1,n$ and measuring the resulting displacements ${\bf u}^{md}_{ij},i=1,n,j=1,m$ and/or strains ${\bf s}^{ms}_{ij},i=1,n,j=1,m$ at $m$ different locations ${\bf x}_{j},j=1,m$ (the intrinsic assumption is that the forces can be standardized and perhaps even maintained throughout the life of the structure);

• -

  The weakening of a structure may occur at any location, i.e. there are no regions that are excluded for weakening; this is the most conservative assumption, and could be relaxed under certain conditions;

• -

  The sensors for displacements and strains are limited in their ability to measure by noise/signal ratios, i.e. actual displacements and strains have to be larger than a certain threshold to be of use:

  $$|{\bf u}^{m}|\geq u_{0}~{}~{},~{}~{}|{\bf s}^{m}|\geq s_{0}~{}~{}.$$

  (2.1)

• -

  The type of force used to monitor the weakening of a structure is limited by practical considerations; this implies that the number of different forces is limited, and can not assume arbitrary distributions in space.

• -

  The weakening of a structure may be described by a field $\alpha({\bf x})$, where $0\leq\alpha({\bf x})\leq 1$ and $\alpha({\bf x})=0$ corresponds to total failure (no load bearing capability) while $\alpha({\bf x})=1$ is the original state;

• -

  The displacements, strains and stresses of the structure are well described by a sufficiently fine finite element discretization (e.g. trusses, beams, plates, shells, solids) [28, 26], which results in a system of equations for each load case:

  $${\bf K}{\bf u}_{i}={\bf f}_{i}~{}~{},~{}~{}i=1,n$$

  (2.2)

• where ${\bf u}_{i}$ are the displacements and ${\bf K}$ the usual stiffness matrix, which is obtained by assembling all the element matrices:

  $${\bf K}=\sum_{e=1}^{N_{e}}{\bf K}_{e}~{}~{}.$$

  (2.3)

The determination of material properties (or weaknesses) may be formulated as an optimization problem for the strength factor $\alpha({\bf x})$ as follows: Given $n$ force loadings ${\bf f}_{i},i=1,n$ and $n\cdot m$ corresponding measurements at $m$ measuring points/locations ${\bf x}_{j},j=1,m$ of their respective displacements ${\bf u}^{md}_{ij},i=1,n,j=1,m$ or strains ${\bf s}^{ms}_{ij},i=1,n,j=1,m$, obtain the spatial distribution of the strength factor $\alpha$ that minimizes the cost function:

$$I({\bf u}_{1,..,n},{\bf s}_{1,..,n},\alpha)={1\over 2}\sum_{i=1}^{n}\sum_{j=1}^{m}w^{md}_{ij}({\bf u}^{md}_{ij}-{\bf I}^{d}_{ij}{\bf u}_{i})^{2}+{1\over 2}\sum_{i=1}^{n}\sum_{j=1}^{m}w^{ms}_{ij}({\bf s}^{ms}_{ij}-{\bf I}^{s}_{ij}{\bf s}_{i})^{2}~{}~{},$$

(3.1)

subject to the finite element description (e.g. trusses, beams, plates, shells, solids) of the structure [28, 26] under consideration (i.e. the digital twin/system [18, 9]):

$${\bf K}\cdot{\bf u}_{i}={\bf f}_{i}~{}~{},~{}~{}i=1,n$$

(3.2)

where $w^{md}_{ij},w^{ms}_{ij}$ are displacement and strain weights, ${\bf I}^{d},{\bf I}^{s}$ interpolation matrices that are used to obtain the displacements and strains from the finite element mesh at the measurement locations, and ${\bf K}$ the usual stiffness matrix, which is obtained by assembling all the element matrices:

$${\bf K}=\sum_{e=1}^{N_{e}}\alpha_{e}{\bf K}_{e}~{}~{},$$

(3.3)

where the strength factor $\alpha_{e}$ of the elements has already been incorporated. We note in passing that in order to ensure that ${\bf K}$ is invertible and non-degenerate $\alpha_{e}>\epsilon>0$. Note that the optimization problem given by Eqns.(3.1)-(3.3) does not assume any specific choice of finite element basis functions, i.e. is widely applicable.

The location of a displacement or strain gauge may not coincide with any of the nodes of the finite element mesh. Therefore, in general, the displacement ${\bf u}_{k}$ at a measurement location ${\bf x}^{m}_{k}$ needs to be obtained via the interpolation matrix ${\bf I}^{d}_{k}$ as follows:

$${\bf u}_{k}({\bf x}^{m}_{k})={\bf I}^{d}_{k}({\bf x}^{m}_{k})\cdot{\bf u}~{}~{},$$

(4.1)

where ${\bf u}$ are the values of the displacements vector at all gridpoints.
In many cases it is much simpler to install strain gauges instead of displacement gauges. In this case, the strains need to be obtained from the displacement field. This can be written formally as:

$${\bf s}={\bf D}{\bf u}~{}~{},$$

(4.2)

where the ‘derivative matrix’ contains the local values of the derivatives of the shape-functions of ${\bf u}$. The strain at an arbitrary position ${\bf x}^{m}_{i}$ is obtained via the interpolation matrix ${\bf I}^{s}_{k}$ as follows:

$${\bf s}_{k}({\bf x}^{m}_{k})={\bf I}^{s}_{k}({\bf x}^{m}_{k})\cdot{\bf s}={\bf I}^{s}_{k}({\bf x}^{m}_{k})\cdot{\bf D}\cdot{\bf u}~{}~{}.$$

(4.3)

Note that in many cases the strains will only be defined in the elements, so that the interpolation matrices for displacements and strains may differ.

The cost function is given by Eqn.(3.1), repeated here for clarity:

$$I({\bf u}_{n},\alpha)={1\over 2}\sum_{i=1}^{n}\sum_{j=1}^{m}w^{md}_{ij}({\bf u}^{md}_{ij}-{\bf I}^{d}_{ij}{\bf u}_{i})^{2}+{1\over 2}\sum_{i=1}^{n}\sum_{j=1}^{m}w^{ms}_{ij}({\bf s}^{ms}_{ij}-{\bf I}^{s}_{ij}{\bf s}_{i})^{2}~{}~{}.$$

(5.1)

One can immediately see that the dimensions of displacements and strains are different. This implies that the weights should be chosen in order that all the dimensions coincide. The simplest way of achieving this is by making the cost function dimensionless. This implies that the displacement weights $w^{md}_{ij}$ should be of dimension [1/(displacement*displacement)] (the strains are already dimensionless). Furthermore, in order to make the procedures more generally applicable they should not depend on a particular choice of measurement units (metric, imperial, etc.). This implies that the weights for the displacements and strains should be of the order of the characteristic or measured magnitude. Several options are possible, listed below.

The gradients of the cost function with respect to $\alpha$ allow for oscillatory solutions. One must therefore smooth or ‘regularize’ the spatial distribution. This happens naturally when using few degrees of freedom, i.e. when $\alpha$ is defined via other spatial shape functions (e.g. larger spatial regions of piecewise constant $\alpha$ [22]). As the (possibly oscillatory) gradients obtained in the (many) finite elements are averaged over spatial regions, an intrinsic smoothing occurs. This is not the case if $\alpha$ and the gradient are defined and evaluated in each element separately, allowing for the largest degrees of freedom in a mesh and hence the most accurate representation. Three different types of smoothing or ‘regularization’ were considered. All of them start by performing a volume averaging from elements to points:

$$\alpha_{p}={{\sum_{e}\alpha_{e}V_{e}}\over{\sum_{e}V_{e}}}~{}~{},$$

(6.1)

where $\alpha_{p},\alpha_{e},V_{e}$ denote the value of $\alpha$ at point $p$, as well as the values of $\alpha$ in element $e$ and the volume of element $e$, and the sum extends over all the elements surrounding point $p$.

The optimization cycle outlined above can be implemented in a very efficient way if one has direct access to the source-code of finite element-based structural mechanics solvers, but is also amenable to black-box (e.g. commercial) solvers. A possible way to proceed is the following:

• -

  Output the original stiffness matrix ${\bf K}_{el}$ for each element;

• -

  For each optimization step/cycle:

• -

  With the current element values for $\alpha$: build the new stiffness matrix; this is usually done with a user-defined subroutine or module (all commercial codes allow for that);

• -

  For each load case $i$:

• -

  Solve the forward problem ($\rightarrow{\bf u}_{i}$);

• -

  Post-process the results of the forward problem in order to obtain the displacements and strains;

• -

  With the measured and computed displacements/ strains and weights: compute the cost function part of this load case ($\rightarrow I_{i}$);

• -

  Build the ‘force vector’ (i.e. the right-hand-side) for the adjoint problem by comparing the measured and computed displacements/ strains and weighting them appropriately;

• -

  Solve the adjoint problem ($\rightarrow{\tilde{\bf u}}_{i}$);

• -

  With ${\bf u}_{i},{\tilde{\bf u}}_{i}$ and the original stiffness matrix ${\bf K}_{el}$: get the gradient in each element;

• -

  Smooth the gradients (either via a ‘fake-heat-solver’ if Laplacian smoothing is desired, or via an external smoother);

• -

  Send the cost function and the smoothed gradients to the optimizer;

• -

  Update $\alpha$

The aim of choosing a minimal yet optimal set of forces to monitor the weakening of structures is to be able to determine the field $\alpha({\bf x})$ as best as possible. From structural mechanics, this implies that one should avoid regions where the strains are very small or vanish. For these regions, $\alpha({\bf x})$ can assume an arbitrary value without having any effect on the overall displacements or strains. Therefore, the forces should be chosen such that the number of regions with very small or vanishing strains should be minimized.
As stated before, for practical reasons the number and type of possible loads is limited. This ‘limitation of the search space for loads’ opens up the possibility of a simple algorithm to determine the optimal choice. For each of the $n$ possible loads ${\bf f}_{i},i=1,n$, obtain the resulting displacement and strains, and record all elements for which the strains are above a minimal sensor threshold $s_{0}$.
If only one force is to be applied, the obvious choice is to select the one that produces the largest area with strains that are above a minimal threshold. Having selected this force, the regions that have already been affected by this force (i.e. with strains that are above the minimal threshold) are excluded from further consideration. The next best force is then again the one that is able to measure the largest area with strains that are above a minimal threshold. And so on recursively.

Let us assume that a certain part $\Omega^{w}$ of the structure has weakened. This could be a region of several elements, or a single element. This will lead to a change in the stiffness matrix, and a resulting change in displacements and strains. The aim is to be able to record and identify the spatial location of this weakening with the minimum number of sensors. The change in displacements or strains due to a weakening requires the evaluation of the derivative

$$D_{\alpha}({\bf x}_{i},{\bf x}_{j})={{\partial{\bf u}({\bf x}_{i})}\over{\partial\alpha({\bf x}_{j})}}$$

(9.1)

for all possible combination of locations ${\bf x}_{i},{\bf x}_{j}$. In the most general case ${\bf x}_{i}$ is arbitrary, i.e. it could be any node or element of the mesh. However, if sensors can only be placed in certain regions of the domain, the location of ${\bf x}_{i}$ can be reduced significantly.
Two possible ways were explored to obtain $D_{\alpha}({\bf x}_{i},{\bf x}_{j})$: forward-based and adjoint-based.

All the numerical examples were carried out using two finite element codes. The first, FEELAST [16], is a finite element code based on simple linear (truss), triangular (plate) and tetrahedral (volume) elements with constant material properties per element that only solves the linear elasticity equations. The second, CALCULIX [10], is a general, open source finite element code for structural mechanical applications with many element types, material models and options. The optimization loops were steered via a simple shell-script for the adjoint-based gradient descent method. In all cases, a ‘target’ distribution of $\alpha({\bf x})$ was given, together with defined external forces ${\bf f}_{\Gamma}$. The problem was then solved, i.e. the displacements ${\bf u}({\bf x})$ and strains ${\bf s}({\bf x})$ were obtained and recorded at the ‘measurement locations’ ${\bf x}_{j},~{}j=1,m$. This then yielded the ‘measurement pair’ ${\bf f},{\bf u}_{j},~{}j=1,m$ or ${\bf f},{\bf s}_{j},~{}j=1,m$ that was used to determine the material strength distributions $\alpha({\bf x})$ in the field.

An adjoint-based procedure to determine weaknesses, or, more generally the material properties of structures has been presented. Given a series of force and displacement/strain measurements, the material properties are obtained by minimizing the weighted differences between the measured and computed values. In a subsequent step techniques to optimize the number of loadings and sensors have been proposed and tested.
Several examples show the viability, accuracy and efficiency of the proposed methodology.
We consider this a first step that demonstrates the viability of the adjoint-based methodology for system identification and its use for high fidelity digital twins [19, 9].
Many questions remain open, of which we just mention some obvious ones:

• -

  Will these techniques work for nonlinear problems ?

• -

  Which sensor resolution is required to obtain reliable results ?

• -

  Will these techniques work under uncertain measurements ? [3]

• -

  Can one detect faulty sensors in a systematic way ?

This work has been partially supported by NSF grant DMS-2110263 and the Air Force Office of Scientific Research (AFOSR) under Award NO: FA9550-22-1-0248.