Toward to an Automatic Crack Size Estimation with iFEM for Structural Health Monitoring

1. Introduction

Structural integrity and damage assessment are crucial factors to guarantee the safety of different structures and components. Safety requirements depend on the particular application and are defined in relation to the potential consequences induced by a failure or a malfunction. Aeronautical structures are particularly critical from the safety point of view since a failure can induce catastrophic consequences. For this reason, non-destructive controls (such as visual inspection, liquid penetrant, radiography, magnetic particles, etc.) are performed at scheduled intervals to early detect the presence of damages. However, these controls frequently require service interruption with a significant economic impact on the overall aircraft management. This problem is mitigated by the adoption of Structural Health Monitoring (SHM) and Prognostic Health Monitoring (PHM) techniques, which are aimed to perform a real-time assessment of the structure. The aircraft is equipped with permanently installed sensors able to represent the working conditions, then data-driven [1,2,3] or model-based [4,5] algorithms can early detect anomalies. This can increase the service time between two consecutive maintenances, decreasing the overall cost. Data-driven algorithms frequently rely on historical data to detect anomalies, while the latter approach is based on structure models in both the healthy and damaged states. Then, algorithms frequently exploit machine learning, e.g. Artificial Neural Networks (ANN), and a database of damage scenarios to perform the final structure’s assessment [6]. However, one main limitation of most of the literature’s SHM techniques is their dependency on the external load applied to the structure [7], which often hampers their applicability in real scenarios.

In recent years, there was a further shift of paradigm in the structure’s monitoring with the adoption of the Digital Twin (DT) approach [8,9,10,11]. Digital Twin is an integrated multi-physics, multi-scale, probabilistic simulation of an as-built system that uses the best available physical models, sensor updates, fleet history, etc., to mirror the life of its corresponding flying twin [12]. For each physical aeronautical structure, there is a related DT, i.e., a digital representation of the same through a high-fidelity model. The number of DT literature examples is rapidly increasing in the last years, with applications in the aeronautical field [12,13,14,15,16], marine industry (offshore wind turbine and naval transport) [17,18,19], and manufacturing processes [20,21,22]. Focusing on aeronautical applications, nowadays DTs are mainly used to model the airframe of aircraft [13,14,16] and to predict the structural fatigue life of cracked structures [18,23]. Although there is not a unique literature definition of the components of a DT [9], the following main elements can be identified: (i) a physical structure and/or system, (ii) a digital representation (master model) of the structure/system, (iii) a connection between the physical and the virtual structure/system, which is generally represented by sensors, often referred to as digital shadow. In particular, the DT model must be fast enough to run in real-time, parallelly to the operation of the real structure, based on data from sensors or simulated data to predict future states of the system. For this reason, DTs are frequently based on ANN algorithms, analytical relations, other surrogate models, or a combination of these approaches. Another fundamental aspect is that the DT must be a perfect mirror of the real structure, thus the occurrence of damage or an anomaly must be reflected in the model. This requires a constant update of the DT model, where any modification of the structure/system (also due to maintenance) is reported. In conclusion, DT allows near real-time updating of the structural model through physical measurements acquired from sensors to perform damage diagnosis and prognosis in an SHM framework.

Among the different algorithms available in the literature, the inverse Finite Element Method (iFEM) is a model-based technique to compute the structure’s displacement field from strain measurements. It requires only a mesh discretization of the structure with a definition of boundary conditions, and its efficient formulation is suitable for real-time applications. Furthermore, it does not require any material property definition or knowledge of the loading condition, only the strain measurements acquired from sensors are needed as input. More specifically, iFEM generally requires strain sensors (e.g. stain gauges or fiber optic sensors) bonded on the external sides of the structure, although applications with embedded sensors (for composite materials) are available in the literature [24]. These characteristics make the iFEM attractive for a DT and, to the best authors’ knowledge, no application has been already reported in the literature. However, one main limitation hampering the implementation of the current iFEM approach in a DT framework is the lack of model updating capability, thus this manuscript aims to cover this literature gap in view of future DT applications. In particular, current iFEM applications are always based on an initial model definition (in the undamaged configuration), without any updating to account for damage propagation or maintenance operations. This paper newly proposes an updating framework for damage identification with the iFEM, which paves the way for DT practioners to employ iFEM models rather than direct methods such as FEM and XFEM, where damage introduction and propagation are already state-of-the-art. This is on the basis that the application of the latter methods becomes significantly challenging for applications subjected to unknown and stochastic loading conditions where only strain measurements are available, which justifies the development of the inverse FEM framework.

The iFEM was originally introduced by A. Tessler et al. [25,26] and nowadays its formulation is available for beam [27,28,29,30,31,32,33] and shell structures [34,35,36,37,38,39,40,41]. It is based on the minimization of a least-square functional, which is representative of the error between the input strain measurements acquired by sensors and a numerical formulation of the same, function of the unknown nodal degrees of freedom (dof). Focusing on shell structures, three types of inverse elements are nowadays available in the literature for mesh discretization: (i) the iMIN3, a triangular element based on bilinear anisoparametric shape functions [42,43], (ii) the iQS4, a quadrilateral flat element based on bilinear anisoparametric shape functions [37,44,45], and (iii) the iCS8, an eight-nodes curved element based on quadratic isoparametric shape functions [46]. A comparative study assessed the performances of these elements also in relation to the different application scenarios [47]. In addition to the basic elements described, additional formulations have been introduced by the different authors to increase the results’ accuracy of specific classes of problems. In particular, the Refined Zig-Zag Theory (RZT) is introduced in the iFEM formulation to model the through-the-thickness displacement field of composite laminate [48,49,50,51,52], while the iFEM isogeometric analysis is beneficial in the case of large non-linear deformations [53].

Several shape sensing applications with iFEM are available in the literature, based both on numerical and experimental case studies. However, although strain sensors can be applied to the whole structure when dealing with numerical case studies [44], their number and locations are one of the main constraints for practical applications [39,54]. Hardware limitations frequently limit the number of sensors available and their installation is subordinated to practical reasons, like access to the structure. As a consequence, sensors generally cover only a limited portion of the structure and their location must be optimized, according to the different constraints, to well describe the strain field. Furthermore, the adoption of pre-extrapolation techniques was revealed an effective tool to increase the overall results’ accuracy, providing an input strain value also where sensors are not available [24,55]. In particular, this can be based on data-driven approaches, like polynomial functions or the Smoothing Element Analysis (SEA) [56,57,58,59], on a physics-based approach [60], or Gaussian Process interpolation for a statistical input strain evaluation [61]. More specifically, data-driven approaches pre-extrapolate the strain field only based on the measurements acquired from sensors, thus the sensor network must be representative of the strain field within the structure to obtain accurate results. However, in the case of high strain gradients induced by local discontinuities (e.g. holes and notches), the design of a sensor network able to well describe the strain field is a challenging task, in particular when the number of sensors is limited. Thus, the adoption of a physics-based pre-extrapolation approach combines the strain measurements from the available sensors with the physical knowledge of the discontinuity (e.g. its size and position) to better pre-extrapolate the strain field. In particular, this approach relies on the analytical stress solution of the discontinuity itself, when available, otherwise on its numerical stress solution computed with FEM.

In addition to shape sensing, iFEM has been also extended to damage detection in an SHM framework with different approaches, such as load-independent damage indices [44,54], damage parameters based on the Von-Mises strain [36,62], and Artificial Neural Networks [63]. However, these approaches are able to identify the presence and, in case, the location of damage in the structure (both metallic and composites), but not estimate its size, which is still an unexplored area of the iFEM applications. In particular, the numerical case studies proposed by Colombo et al. [44] investigate the crack detection on a metallic plate with a load-independent damage index accounting also for different crack sizes and orientations. Furthermore, in this study, a rough estimation of the crack size could be performed from the damage index pattern since all the elements around the crack itself are covered by sensors. However, this approach is only feasible in a numerical case study and cannot be performed in a real scenario. Another limitation of the previous literature works is the requirement of a damage index threshold to detect the presence of the damage itself, which selection is subordinated to the specific case study. Only Kefal et al. recently proposed a coupling of periodynamics analysis and iFEM for the crack propagation monitoring in composite plate structures [64], however, an automatic crack size estimation through the iFEM remains unexplored. For this reason, this work is aimed to define an automatic routine to perform damage size estimation in a realistic scenario within the iFEM approach. In particular, the damage is systematically introduced in the iFEM model creating different possible damage scenarios, then a maximum likelihood estimation framework selects the model that better describes the experimental strain field measured by test sensors. So, the damage is introduced in the iFEM model better approaching the actual condition of the physical structure, which could be used for future DT frameworks. The methodology is experimentally verified on an aluminum plate subjected to fatigue crack propagation. The strain field is experimentally measured by an Optical Backscatter Reflectometry fiber and from independent observations of the Digital Image Correlation, which is used to further validate the methodology.

The manuscript is organized as follows. A review of the iFEM methodology is reported in Section 2, while its extension for the crack size estimation is provided in Section 3. The experimental case study is presented in Section 4 and their results and discussion are provided in Section 5. Finally, the conclusions of the work are stated in Section 6.

2. Inverse Finite Element Method review

An overview of the iFEM [25,26] is reported in this Section, while a more detailed review specific to the iQS4 element is available in [37,44] for interested readers.

Suppose a shell structure is discretized into finite elements, as in the direct FEM. However, the mesh is composed of inverse elements, the iQS4 in this case, which compute the displacement field from input strain measurements. This is done by minimizing the least-square functional of Equation (1), which is defined as the error between the input strain field acquired by sensors $\left({\cdot }^{\epsilon }\right)$ and its numerical formulation $(\cdot \left(\mathit{u}\right))$, which is is turn function of the unknown nodal displacements $\mathit{u}$. Both the input and the numerical strain fields are decoupled into three main components: the membrane $\mathit{e}$, the bending $\mathit{k}$, and the transverse shear $\mathit{g}$ strain contributions. Thus, the formulation of the i-th inverse element can be defined as:

$${\mathsf{\Phi }}_i\left({\mathit{u}}^i\right)={‖\mathit{e}\left({\mathit{u}}^i\right)-{\mathit{e}}_i^{\epsilon }‖}_{{\mathit{W}}_m^i}^2+{‖\mathit{k}\left({\mathit{u}}^i\right)-{\mathit{k}}_i^{\epsilon }‖}_{{\mathit{W}}_b^i}^2+{‖\mathit{g}\left({\mathit{u}}^i\right)-{\mathit{g}}_i^{\epsilon }‖}_{{\mathit{W}}_s^i}^2$$

(1)

where ${‖\cdot ‖}_{\mathit{W}}^2$ is the squared weighted Euclidean norm with the weight matrix $\mathit{W}$. In particular, ${\mathit{W}}_m^i$, ${\mathit{W}}_b^i$, and ${\mathit{W}}_s^i$ are diagonal matrices of weights for the membrane, bending, and transverse shear strain contributions respectively. These coefficients control the coherence between the numerical and the experimental strain measurements, in particular, in the case of sparse sensor networks. In general, a unitary reference value is associated to the elements in which the input strain field component is acquired by physical sensors, while, in other cases, the coefficients are generally reduced to small values (e.g. ${10}^{-4}$). Notice that, each matrix $\mathit{W}$ contains three weights on the main diagonal, which are related to the strain components along the x-axis, the y-axis, and the in-plane shear with respect to the element’s local reference system (Figure 2). Then, in case an element is interested by a monoaxial strain sensor, only the weight related to its direction is assigned equal to one and the others are reduced to a small value.

After a proper assembly procedure, the unknown structure’s displacement field is computed by minimizing the error functional presented in Equation (1), which will be detailed in Section 2.3 after a proper definition of the input (Section 2.1) and numerical (Section 2.2) strain fields.

2.1. Input strain formulation

In the most general case, the input strain formulation is computed from the strain measurements acquired on the structure. Sensors are generally applied on the external surfaces of the component, where their installation and maintenance are easier, although applications with embedded sensors are also possible [24].

For example, consider a couple of strain gauge rosettes applied on the two external sides of the shell as shown in Figure 1. The membrane and the bending strain components associated with the j-th sensors’ location within the i-th inverse element can be defined as:

$${\mathit{e}}_{ij}^{\epsilon }=\frac{1}{2}{\left\{\begin{array}{c}{\epsilon }_{xx}^{+}+{\epsilon }_{xx}^{-}\\ {\epsilon }_{yy}^{+}+{\epsilon }_{yy}^{-}\\ {\gamma }_{xy}^{+}+{\gamma }_{xy}^{-}\end{array}\right\}}_j {\mathit{k}}_{ij}^{\epsilon }=\frac{1}{2h}{\left\{\begin{array}{c}{\epsilon }_{xx}^{+}-{\epsilon }_{xx}^{-}\\ {\epsilon }_{yy}^{+}-{\epsilon }_{yy}^{-}\\ {\gamma }_{xy}^{+}-{\gamma }_{xy}^{-}\end{array}\right\}}_j$$

(2)

Where $2h$ is the shell thickness at the sensors’ location. The computed strain components contain the information of a plane strain tensor and, in case monoaxial strain sensors are used, only one component will be defined and the others are posed equal to zero.

Furthermore, in practical applications only few sensors are available due to costs, space constraints, and hardware limitations, limiting the definition of the input strain field. However, for an accurate iFEM computation, the input strain field should be provided on all the structure’s inverse elements and it should be representative of the strain gradient. For this reason, the sensor network should be properly designed and the element’s size can be tuned according to the expected strain gradient. Nevertheless, several elements of the mesh can be free from any input definition and, to limit their influence on the global formulation, small weighting coefficients $w_i$ ($i=m,k,s$) can be associated to these elements. To further improve the iFEM results accuracy, the input strain field can be pre-extrapolated in the element’s locations in which physical sensors are not available, obtaining the definition of the input strain field on the whole structure. This can be done with different approaches according to the specific problem, such as polynomial fitting and the Smoothing Element Analysis [56,57,58,59]. These techniques are purely based on the acquired strain measurements from sensors, thus being defined as data-based approaches, only providing a more continuous and smooth strain field in the considered domain. Then, the recent introduction of the physics-based strain pre-extrapolation technique [60] allows a further increase of the iFEM results’ accuracy in the case of sparse sensor networks on notched structures. In this case, the physical knowledge of the discontinuity together with its analytical stress function allows an accurate definition of the input strain field. This pre-extrapolation technique will be detailed in Section 3.4.2 for the particular case under analysis, i.e. a cracked plate.

2.2. Numerical strain formulation

The numerical strain formulation required by Equation (1) is based on the element’s shape functions as in any finite element approach.

A local reference system $(x,y,z)$ is defined within each inverse element, with its origin in the centroid of the element and with the z-axis defining the out-of-plane coordinate, so that $z\in [-h;+h]$, as also illustrated in Figure 2. The local coordinates are computed from the global reference system $(X, Y, Z)$, in which the structure is defined, with a proper transformation matrix, specific for each element.

Each iQS4 inverse element is composed by four nodes, each one with six degrees of freedom (dof). In particular, each element has 24 dof with 3 translations ${\mathit{u}}_q={\left\{u_x, u_y, u_z\right\}}^q$ and 3 rotations ${\mathit{\theta }}_q={\left\{{\theta }_x, {\theta }_y,{\theta }_z\right\}}^q$ for each element’s node $q=(\mathrm{1,2},\mathrm{3,4})$, which are collected into the element’s nodal displacement vector ${\mathit{u}}^i$. Then, the local displacement field within each inverse element is defined through the shape functions ${\mathit{N}}_q$, ${\mathit{L}}_q$, and ${\mathit{M}}_q$ [37] as:

$$\begin{gathered} u\left(x,y\right)={\sum }_{q=1}^4{N}_q{u}_{xq}+{\sum }_{q=1}^4{L}_q{\theta }_{zq} \\ v\left(x,y\right)={\sum }_{q=1}^4{N}_q{u}_{yq}+{\sum }_{q=1}^4{M}_q{\theta }_{zq} \\ w\left(x,y\right)={\sum }_{q=1}^4{N}_q{u}_{zq}-{\sum }_{q=1}^4{L}_q{\theta }_{xq}-{\sum }_{q=1}^4{M}_q{\theta }_{yq} \\ {\theta }_x\left(x,y\right)={\sum }_{q=1}^4{N}_q{\theta }_{xq} \\ {\theta }_y\left(x,y\right)={\sum }_{q=1}^4{N}_q{\theta }_{yq} \end{gathered}$$

(3)

Then, under the hypothesis of plane stress condition and after computing the partial derivatives of the shape functions, the strain field components within each element can be defined as:

$$\begin{gathered} \mathit{e}\left({\mathit{u}}^i\right)={\mathit{B}}^m{\mathit{u}}^i \\ \mathit{k}\left({\mathit{u}}^i\right)={\mathit{B}}^b{\mathit{u}}^i \\ \mathit{g}\left({\mathit{u}}^i\right)={\mathit{B}}^s{\mathit{u}}^i \end{gathered}$$

(4)

where ${\mathit{B}}^m$, ${\mathit{B}}^b$, and ${\mathit{B}}^s$ are matrices containing the derivatives of the shape functions.

Finally, the numerical strain field can be computed with the following relations in correspondence of the required $z$ coordinate:

$$\begin{gathered} \left\{\begin{array}{c}{\epsilon }_{xx}\\ {\epsilon }_{yy}\\ {\gamma }_{xy}\end{array}\right\}=\mathit{e}\left({\mathit{u}}^i\right)+z\cdot \mathit{k}\left({\mathit{u}}^i\right) \\ \left\{\begin{array}{c}{\gamma }_{xz}\\ {\gamma }_{yz}\end{array}\right\}=\mathit{g}\left({\mathit{u}}^i\right) \end{gathered}$$

(5)

Figure 2. iQS4 element with global (X,Y,Z) and local (x,y,z) reference system and the related degrees of freedoms (dof). Numbers 1 to 4 refer to the element nodes.

Figure 2. iQS4 element with global (X,Y,Z) and local (x,y,z) reference system and the related degrees of freedoms (dof). Numbers 1 to 4 refer to the element nodes.

2.3. Matrix formulation

As previously introduced, the iFEM relies on the minimization of the global formulation of Equation (1), however, limited to a single element until now. Thus, before considering the contribution of all the inverse elements discretizing the structure, Equation (1) is elaborated to achieve an efficient numerical computation. In particular, the 3 square weighted norms are expressed as normalized Euclidean norms as:

$$\begin{gathered} {‖\mathit{e}\left({\mathit{u}}^i\right)-{\mathit{e}}_i^{\epsilon }‖}_{{\mathit{W}}_m^i}^2=\frac{1}{n}{\iint }_{A_i}^{\text{}}\sum _{j=1}^n{\left(\mathit{e}{\left({\mathit{u}}^i\right)}_j-{\mathit{e}}_{ij}^{\epsilon }\right)}^T{\mathit{W}}_m^i\left(\mathit{e}{\left({\mathit{u}}^i\right)}_j-{\mathit{e}}_{ij}^{\epsilon }\right)dxdy \\ {‖\mathit{k}\left({\mathit{u}}^i\right)-{\mathit{k}}_i^{\epsilon }‖}_{{\mathit{W}}_b^i}^2=\frac{{\left(2h\right)}^2}{n}{\iint }_{A_i}^{\text{}}\sum _{j=1}^n{\left(\mathit{k}{\left({\mathit{u}}^i\right)}_j-{\mathit{k}}_{ij}^{\epsilon }\right)}^T{\mathit{W}}_b^i\left(\mathit{k}{\left({\mathit{u}}^i\right)}_j-{\mathit{k}}_{ij}^{\epsilon }\right)dxdy \\ {‖\mathit{g}\left({\mathit{u}}^i\right)-{\mathit{g}}_i^{\epsilon }‖}_{{\mathit{W}}_s^i}^2=\frac{1}{n}{\iint }_{A_i}^{\text{}}\sum _{j=1}^n{\left(\mathit{g}{\left({\mathit{u}}^i\right)}_j-{\mathit{g}}_{ij}^{\epsilon }\right)}^T{\mathit{W}}_s^i\left(\mathit{g}{\left({\mathit{u}}^i\right)}_j-{\mathit{g}}_{ij}^{\epsilon }\right)dxdy \end{gathered}$$

(6)

Where $A_i$ is the area of the considered element and $n$ the total number of input strain sensors within the same. Then, after substituting the input and the numerical strain components previously developed in Section 2.1 and Section 2.2 respectively, the least-square functional can be expressed as:

$${\mathsf{\Phi }}_i\left({\mathit{u}}^i\right)={\mathit{u}}^{i^T}{\mathit{k}}^i{\mathit{u}}^i-2{\mathit{u}}^{i^T}{\mathit{f}}^i+{\xi }^i$$

(7)

In which:

$$\begin{gathered} {\mathit{k}}^i={\iint }_{A_i}^{\text{}}({\mathit{B}}^{{\mathit{m}}^T}{\mathit{W}}_m{\mathit{B}}^{\mathit{m}}+{\left(2h\right)}^2{\mathit{B}}^{{\mathit{b}}^T}{{\mathit{W}}_b\mathit{B}}^{\mathit{b}}+{\mathit{B}}^{{\mathit{s}}^T}{{\mathit{W}}_s\mathit{B}}^{\mathit{s}})dxdy \\ {\mathit{f}}^i=\frac{1}{n}{\iint }_{A_i}^{\text{}}\sum _{j=1}^n({\mathit{B}}^{{\mathit{m}}^T}{\mathit{W}}_m{\mathit{e}}_{ij}^{\epsilon }+{\left(2h\right)}^2{\mathit{B}}^{{\mathit{b}}^T}{\mathit{W}}_b{\mathit{k}}_{ij}^{\epsilon }+{\mathit{B}}^{{\mathit{s}}^T}{\mathit{W}}_s{\mathit{g}}_{ij}^{\epsilon })dxdy \\ {\xi }^i=\frac{1}{n}{\iint }_{A_i}^{\text{}}\sum _{j=1}^n({\mathit{e}}_{ij}^{{\epsilon }^T}{{\mathit{W}}_m\mathit{e}}_{ij}^{\epsilon }+{\left(2h\right)}^2{\mathit{k}}_{ij}^{{\epsilon }^T}{{\mathit{W}}_b\mathit{k}}_{ij}^{\epsilon }+{\mathit{g}}_{ij}^{{\epsilon }^T}{{\mathit{W}}_s\mathit{g}}_{ij}^{\epsilon })dxdy \end{gathered}$$

(8)

Integrals can be efficiently computed numerically with the Gauss quadrature, for example adopting four integration points for each inverse element for an accurate result. Then, a standard assembly procedure accounts for the contribution of all structure’s inverse elements, obtaining a global formulation of Equation (7). Finally, this is minimized with respect to the global displacement field (i.e., $\partial \mathsf{\Phi }/\partial \mathit{U}=0$) and thus the problem can be written as:

$$\mathit{K}\mathit{U}=\mathit{F}$$

(9)

Where $\mathit{K}$ is a matrix linking the nodal displacement $\mathit{U}$, in global coordinates ($U_x$, $U_y$, $U_z$, ${\mathsf{\Theta }}_x$, ${\mathsf{\Theta }}_y$, and ${\mathsf{\Theta }}_z$ for each node), with the vector $\mathit{F}$, which is a function of the input strain field. Notice that the matrix $\mathit{K}$ depends on the structure’s element discretization and on the sensor network configuration, while $\mathit{F}$ only depends on the input strain measurements. However, the matrix $\mathit{K}$ is singular and it will lead to a rigid motion of the structure if unconstrained, thus, after the definition of problem-specific boundary conditions, the unconstrained (free) nodal displacements can be computed as:

$${\mathit{K}}_{\mathit{F}\mathit{F}}{\mathit{U}}_{\mathit{F}}={\mathit{F}}_{\mathit{F}}⟹{\mathit{U}}_{\mathit{F}}={{\mathit{K}}_{\mathit{F}\mathit{F}}^{-1}\mathit{F}}_{\mathit{F}}$$

(10)

Finally, after the computation of the structure’s displacement field, the numerical strain can be computed on the whole structure through Equation (5).

3. Crack size estimation technique

A basic requirement of any iFEM application is that the model well describes the physical structure, which is generally true in the healthy condition of the structure. Then, in case a damage nucleates and propagates in the structure, the iFEM displacement and strain reconstruction will be different than the acquired data (often passed as input to the iFEM) due to a geometrical non-compatibility. This is exploited by different literature researches [36,44,54,62], where the strain difference between the acquired data and the iFEM reconstruction highlights the presence of a potential damage. However, these approaches are limited to damage detection and localization, without attempting a proper diagnosis. Thus, in this Section, a new strategy will be proposed to automatically estimate the damage size, in this case, the crack length on a metallic plate. In particular, the proposed approach does not require any training process and an a priori creation of a damage scenarios database for comparison is not mandatory.

3.1. Crack size estimation overview

As already anticipated, if the iFEM model will be geometrically different from the physical structure, due to damage, also the displacement and strain field reconstructions will be different from the acquired data. Thus, the idea is to systematically introduce the damage into the iFEM model to reduce its discrepancy with respect to the damaged structure. Then, the model that presents the lower discrepancy with respect to the experimental observations is the best description of the actual structural condition (e.g. its damaged state), specifically for the current mesh discretization.

The approach will be detailed for the particular case of a cracked metallic plate, however, it remains valid also for other damage mechanisms.

3.2. Crack propagation routine

Starting from an iFEM model defined in the healthy condition, the first step regards the introduction of the damage in the model itself. For the particular case of a crack in a shell structure, this can be summarized with the following procedure (Figure 3):

• The central position of the crack is identified as one node of the undamaged iFEM mesh, as highlighted by the red node in Figure 3a. This localization can be based on several approaches, such as visual inspection of the structure or automatic localization with a specific algorithm. The present research will not investigate any particular localization technique and the crack center is supposed to be known from any visual or automatic inspection, however, the approach described remains valid for any general localization method.
• The crack is introduced in the mesh by doubling the node corresponding to the crack center and modifying the connectivity matrix to unlink the elements (Figure 3b). The crack is supposed to propagate under a pure Mode I loading, thus its orientation is perpendicular to the first principal stress. Finally, the two nodes identifying the crack tips are highlighted in blue in Figure 3b. This procedure, together with the crack center identification, is detailed in Algorithm 1.
• The crack is further propagated perpendicularly to the first principal stress. The two nodes on the crack tip’s locations are doubled and the connectivity matrix is further modified to unlink the elements, as shown in Figure 3c and detailed in Algorithm 2.
• Repeat step 3 several times, opening one node per time, until the most likely crack length is reached according to the procedure detailed in Section 3.3.

For the sake of clarity, Figure 3 shows the mesh in a generic deformed state only for a better comprehension of the same. The mesh is always defined in the undeformed configuration of the structure and thus the doubled nodes are overlapped.

Algorithm 1: Crack center identification and introduction on an undamaged structure.
1:	Identify the position ($x,y,z$) of the crack center
2:	Find the maximum principal stress direction from the sensor’s measurements
3:	Find the nearest node of the mesh to the crack center → Center node $n_c={\left[x,y,z\right]}_c$
4:	Identify the two crack tip nodes $n_{T1}={\left[x,y,z\right]}_{T1}$ and ${n_{T2}=\left[x,y,z\right]}_{T2}$ as the nearest nodes to the crack center $n_c$ and perpendicular to the maximum principal stress direction
5:	Double the Center node $n_c^{*}=n_c$ → update the node matrix
6:	Update the connectivity matrix to unlink the elements in correspondence of $n_c$
7:	Compute the crack length: $2a=$ norm($n_{T1}-n_{T2}$)
8:	Compute and invert the iFEM $\mathit{K}$ matrix

Algorithm 2: Crack propagation from an existing damage condition.
1:	Identify the two nodes $n_{T11}$ and $n_{T12}$ as the nearest nodes to the crack tip $n_{T1}$ and perpendicular to the maximum principal stress direction
2:	Identify the new crack tip 1 between $n_{T11}$ and $n_{T12}$ as the farthest node from $n_c$ → $n_{T1-new}$
3:	Identify the two nodes $n_{T21}$ and $n_{T22}$ as the nearest nodes to the crack tip $n_{T2}$ and perpendicular to the maximum principal stress direction
4:	Identify the new crack tip 2 between $n_{T21}$ and $n_{T22}$ as the farthest node from $n_c$ → $n_{T2-new}$
5:	Double the crack tip nodes: $n_{T1}^{*}=n_{T1}$ and $n_{T2}^{*}=n_{T2}$ → update the node matrix
6:	Update the connectivity matrix to unlink the elements in correspondence of $n_{T1}$ and $n_{T2}$
7:	Update the crack tip location: $n_{T1}=n_{T1-new}$ and $n_{T2}=n_{T2-new}$
8:	Compute the crack length: $2a=$ norm($n_{T1}-n_{T2}$)
9:	Compute and invert the iFEM $\mathit{K}$ matrix

3.3. Damage state diagnosis

Once multiple iFEM models have been generated, the next task is identifying which model is more representative of the actual structure’s damage condition and thus, in this application, estimating the damage size. Notice that, for the sake of simplicity, the procedure described relies on a database of iFEM models containing different damage conditions (generated according to Section 3.2), however, this a priori database definition is not mandatory, as will be better discussed in Section 5.2.3.

To approach this problem, the sensors available are subdivided into two sensor networks:

• An input sensor network (${\epsilon }_{inp}$) to feed the iFEM algorithm and compute the displacement and the strain field on the whole structure
• A test sensor network (${\epsilon }_{test}$) to compare the strain field computed by the iFEM (${\epsilon }_{iFEM}$) with the experimental measurements acquired on the structure.

In the most general case, the two sensor networks can be coincident, however, to properly highlight the presence of the damage and perform a better diagnosis, the following approach is adopted (Figure 4):

• The input sensor network is composed only of sensors sufficiently far from the damaged area to acquire the far-field strain. Thus, these strain measurements will not be affected by the presence of the damage, whatever is the actual damage condition considered.
• The test sensor network is composed of sensors located in the area of influence of the damage, where the observations can be correlated to the actual damage state.

The damage-sensitive area is defined by a square domain centered on the crack location and with dimensions $8a$, where $a$ is the expected maximum semi-crack length.

Figure 4. Portion of shell structure with input and test sensor networks definition.

Figure 4. Portion of shell structure with input and test sensor networks definition.

Each strain pattern acquired from sensors is characterized by a specific crack length, which is unknown and the target of this analysis. Thus, for each acquisition, the input strain field is fed to all the iFEM models generated, which, in turn, will compute different displacement and strain reconstructions, as shown in Figure 5. Then, the computed strain fields will be compared to the test measurements and, in order to condensate the information from multiple sensors into a single scalar for each model, the likelihood index is adopted. In particular, for each iFEM model $i$, a Gaussian Likelihood index is computed as:

$$L_i=\prod _{j=1}^{N_{test}}\left(\frac{1}{\sqrt{2\pi {\sigma }^2}}\cdot \mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{{\left({\epsilon }_{test}-{\epsilon }_{iFEM}^i\right)}_j^2}{2{\sigma }^2}\right)\right)$$

(11)

Where $N_{test}$ is the total number of test sensors and ${\sigma }^2$ is a variance parameter related to the measurement uncertainties.

However, being the Likelihood the product of elements smaller than one, its final value can reach very small values which cannot be correctly computed in a numerical environment, in particular when a large number of test sensors is used. Thus, to overcome this issue, the logarithmic Gaussian Likelihood is used:

$${lnL}_i=-\frac{N_{test}}{2}\cdot \ln \left(2\pi {\sigma }^2\right)-\frac{1}{2{\sigma }^2}\sum _{j=1}^{N_{test}}{\left({\epsilon }_{test}-{\epsilon }_{iFEM}^i\right)}_j^2$$

(12)

The Likelihood (or the logarithmic Likelihood) of each model is related to its probability of well reproducing the test measurements and thus the damage state considered. So, the problem is addressed with a maximum likelihood estimation framework, where the iFEM model with the largest logarithmic Likelihood is the one that better minimizes the error with respect to the test measurements, being the diagnosis of the actual damage state.

3.4. Crack strain field pre-extrapolation

The iFEM framework proposed requires a good strain reconstruction near the damage to correctly compute the damage state. However, in this framework the input strain field provides just the far-field strain, thus, although the different iFEM models contain different damage scenarios, the crack will not open and, as a consequence, it will not induce its characteristic strain field. To overcome this issue, the recently developed physics-based strain pre-extrapolation is adopted to impose the strain field specific of each crack length. The approach is described for the specific case of a crack in an infinite plate, while its general formulation is available in [60].

3.4.1. Pre-extrapolation framework overview

As previously introduced, the physics-based strain pre-extrapolation technique is adopted to correctly compute the strain field on the whole structure and, in particular, near the crack, also in the case of a sparse sensor network. This technique accounts for the knowledge of the crack length and position, thus considering its physical strain behavior on the structure. Obviously, the correct crack scenario is not a priori known since it is the objective of the present analysis, as previously described in Section 3.3. Thus, the framework proposed (which is also summarized in Figure 5) is modified by adding the physics-based strain pre-extrapolation technique at each iteration, as detailed in Figure 6. In particular, each iFEM model considered contains a different damage scenario and thus a specific pre-extrapolated input strain field will be computed for the damage-sensitive area. Then, among all the adopted iFEM models, where each one considers a different damage scenario and its relative strain pre-extrapolation, the condition that better describes the real damage is selected with the maximum likelihood estimation framework described in Section 3.3.

3.4.2. Physics-based strain pre-extrapolation for a crack in an infinite plate

For each considered iFEM model (i.e. for each specific damage condition) the strain field is pre-extrapolated on the whole structure with the following physics-based approach [60]:

• A global pre-extrapolation, based on a data-driven approach, is computed on the whole structure to obtain the far-field strain value. In the present analysis, a bi-linear polynomial fitting is adopted, however, other techniques (for example the SEA) can be selected according to the specific case study. Notice that, each strain component (${\epsilon }_{xx}$, ${\epsilon }_{yy}$, ${\gamma }_{xy}$) is independently pre-extrapolated on each side of the structure.
• A local pre-extrapolation better predicts the strain field in the damage-sensitive area with a physics-based approach. This takes as input the far-field strain computed by the data-driven pre-extrapolation to compute the specific strain field induced by the crack.

In particular, the local physics-based approach is applied to a square damage-affected area with dimensions $8a$ (Figure 4), where $a$ is the semi-crack length considered by the current iFEM model. The pre-extrapolation is based on the Westergaard stress solution [65], which considers a crack in an infinite plate subjected to a biaxial state of stress. However, these equations are expressed in terms of stress, while only the strain is available as input data in the iFEM formulation. For this reason, the far-field stress $({\sigma }_x, {\sigma }_y, {\tau }_{xy})$ is computed from the far-field strain (${\epsilon }_{xx}$, ${\epsilon }_{yy}$, ${\gamma }_{xy}$) with the Hook’s equation under the hypothesis of plane stress condition:

$$\left\{\begin{array}{c}{\epsilon }_{xx}\\ {\epsilon }_{yy}\\ {\gamma }_{xy}\end{array}\right\}=\left[\begin{array}{ccc}1& -\nu & 0\\ -\nu & 1& 0\\ 0& 0& 2(1+\nu )\end{array}\right]\cdot \left\{\begin{array}{c}{\sigma }_x\\ {\sigma }_y\\ {\tau }_{xy}\end{array}\right\}$$

(13)

Notice that the elastic modulus of the material is not included in the equation since material properties are generally not required by the iFEM approach. Thus, the computed stress field is not exactly a stress, but it differs by a constant equal to the unknown elastic modulus. However, for proper computation of the stress field, the Poisson ratio $\nu$ is required and literature values can be adopted in case specific material properties are not available.

Then, once the far-field stress has been computed, the Westergaard stress solution can be applied. However, its original solution accounts for biaxial far-field stress (with ${\sigma }_x={\sigma }_y$), while only a monoaxial state of stress is considered under the hypothesis the crack is loaded in a pure Mode I. Different manipulation of the Westergaard solution are available in the literature to address specific problems, in particular, the solution proposed by Elfitis and Liebowitz [66] is reported in Equation (14). This set of equations accounts for a generic bi-axial state of stress, where ${\sigma }_x=k{\sigma }_y$, as illustrated in Figure 7. In particular, the monoaxial state of stress that opens the crack in a pure Mode I can be defined by imposing $S={\sigma }_y$ and $k=0$. Finally, the stress field in the damage-sensitive area can be computed as:

$$\begin{gathered} {\sigma }_x=\frac{S r}{\sqrt{r_1{r}_2}}\cos \left(\theta -\frac{{\theta }_1+{\theta }_2}{2}\right)-\frac{S a^2}{{\left(r_1{r}_2\right)}^{\frac{3}{2}}} r\sin \left(\theta \right)sin\frac{3}{2}\left({\theta }_1+{\theta }_2\right)-(1-k)S \\ {\sigma }_y=\frac{S r}{\sqrt{r_1{r}_2}}\cos \left(\theta -\frac{{\theta }_1+{\theta }_2}{2}\right)+\frac{S a^2}{{\left(r_1{r}_2\right)}^{\frac{3}{2}}} r\sin \left(\theta \right)sin\frac{3}{2}\left({\theta }_1+{\theta }_2\right) \\ {\tau }_{xy}=\frac{S a^2}{{\left(r_1{r}_2\right)}^{\frac{3}{2}}} r\sin \left(\theta \right)cos\frac{3}{2}\left({\theta }_1+{\theta }_2\right) \end{gathered}$$

(14)

where $a$ is the semi-crack length considered and $r$, $r_1$, $r_2$, $\theta$, ${\theta }_1$, and ${\theta }_2$ are linear and angular variables defining a specific point $P$ of the domain in which the stress field is computed, as reported in Figure 7. Notice that the input stress $S$ is obtained from Equation (13), which is in turn function of the strain field pre-extrapolated with a data-driven approach.

Finally, once the stress field induced by the crack has been computed, it can be transformed back into strains through Equation (13), defining the pre-extrapolated strain field to feed the iFEM algorithm.

4. Case study

The crack size estimation methodology proposed is experimentally applied and verified on an aluminum plate subjected to a fatigue crack propagation to test its effectiveness and performance in a real scenario.

This Section introduces the overall case study, in particular, the specimen with its sensor network is described in Section 4.1, the test rig with the acquisition systems is described in Section 4.2, and the iFEM models adopted are presented in Section 4.3.

4.1. Specimen and sensor network

The specimen is an aluminum plate ($E=\mathrm{70,000} MPa$ and $\nu =0.3$) with overall dimensions $650x150 mm$ and a thickness of $2 mm$ (Figure 8). At the two extremities of the specimen, aluminum tabs are connected to ensure a better stress distribution into the plate once mounted in the testing machine. Then, an artificial notch with a total length of $12 mm$ is introduced in the central point of the plate to ensure proper crack nucleation and propagation during the test.

Strains are experimentally measured through a $3 m$ long Optical Backscatter Reflectometry (OBR) fiber bounded on one side of the plate with the 3M™ DP490 epoxy adhesive. The fiber pattern has been designed to cover at the best the region of interest (ROI) of the specimen, which corresponds to the central region of the plate with dimensions $200x150 mm$. In particular, the dimension of the ROI has been selected to avoid any boundary effect induced by the tabs and the gripping system of the testing machine, obtaining a constant far-field stress distribution. This case study design mimics a portion of an aeronautical structure, where a crack can nucleate from a rivet hole. The fiber provides measurement points along both the $x$ and the $y$ axis of the plate according to the reference system in Figure 8. Then, the straight segments of fiber inside the ROI are subdivided into the input and test sensors, in particular, input sensors are sufficiently far from the crack to acquire only the far-field strain, while test sensors are affected by the crack’s strain field.

On the other side of the specimen, a fine speckle pattern was painted to create a uniformly random texture and acquire the strain field on the whole plate’s area by the Digital Image Correlation (DIC) technique. Finally, two reference surfaces are bonded at the two extremities of the ROI to measure its effective displacement with lasers during the test.

4.2. Test rig

The specimen is mounted on a hydraulic MTS monoaxial testing machine with a $100 kN$ load cell, as illustrated in Figure 9. The specimen is subjected to a tensile-tensile fatigue load with $F_{max}=20 kN$ and $R=0.1$ to nucleate and then propagate a crack from the artificial notch. A total number of $\mathrm{32,000}$ cycles have been applied subdivided into different blocks of cycles. Then, between two consecutive blocks of cycles, data have been statically acquired at the maximum load.

The OBR fiber has been acquired with a LUNA ODISI-B interrogator, providing a measurement point every $2.5 mm$ of fiber. This fiber discretization accounts for 345 measurement points for the input sensor network and 117 measurement points for the test one. Furthermore, several strain samples are statically acquired for each crack length to average the values and reduce the noise. It has to be noticed that, since the plate is loaded in tension, the bending strain components are zero (or negligibly small) and only the membrane components contribute to the displacement field of the structure. For this reason, although sensors are applied only on one side of the structure, the membrane strain components can be correctly computed.

The crack length is measured with two Dino-Lite microscopes ($1280x1024$ pixel resolution) on one side of the specimen (one at each crack tip), while, on the other side of the specimen the dual-camera ARAMIS system provided by GOM acquires the images for the DIC technique. The ARAMIS system includes 12 Megapixel cameras and a dual-led lights device. The system was calibrated to acquire the images on a sub-portion of the ROI with dimensions $140x110 mm$, corresponding to the maximum area available for this system. Finally, two MEL M5L/10 lasers acquire the effective displacement of the ROI.

4.3. iFEM models

The structure’s iFEM model is limited to the ROI, thus to a shell structure with in-plane dimensions $150x200 mm$ and a thickness of $2 mm$, as reported in Figure 10. Boundary conditions are representative of the portion of plate in tension, in particular, the lower side of the plate is constrained to avoid any free motion in the 3D space. The horizontal displacement $U_x$ is constrained just in one node to allow the transversal shrinking of the plate itself.

Then, the structure is discretized into inverse elements (the iQS4) with two different structured meshes to investigate the sensitivity of the methodology to the element’s size:

• a mesh with $3.03x3 mm$ element’s size (namely $3 mm$ model), thus the structure is discretized by a total number of 3,300 inverse elements
• a mesh with $1x1 mm$ element’s size (namely $1 mm$ model), thus the structure is discretized by a total number of 30,000 inverse elements

Notice that, although only the membrane response of the element is activated, the selection of the iQS4 element is beneficial thanks to its drilling dof ${\theta }_z$ which can avoid shear locking issues near the crack.

For both cases, the structure is initially discretized in the undamaged configuration, so without the presence of any crack, and the input and the test sensor networks are defined according to Section 4.1. Notice that, since the fiber optic spatial discretization provides a measurement point every $2.5 mm$, the $3 mm$ model contains at least one measurement in every element interested by the fiber. On the contrary, for the $1 mm$ model only one element every two/three elements is covered by an experimental input strain measurement, as also illustrated in Figure 11.

The central point of the crack is then identified on the structure by the node with coordinates $\left(x, y\right)=(75 mm, 100 mm)$, which, in the present work, is assumed known a priori. Then, multiple iFEM models are generated considering different damage conditions, as previously described in Section 3. The crack is artificially increased symmetrically along the $X$ axis with respect to the crack center since the maximum principal stress acts along the $Y$ direction. Notice that, the step increment between two consecutive crack lengths is equal to two times the element’s size considered by each model. For example, considering the $3 mm$ model, the following crack lengths are considered: 6, 12, 18, 24, 30, 36, 42, and 48 mm. While, for the $1 mm$ model the step increment is reduced to $2 mm$ and thus the following crack lengths are generated: 2, 4, 6, 8, …, 46, 48 mm. Notice that, this a priori generation of the damage conditions is not strictly necessary for the algorithm’s routine since only the required models can be generated on-demand when needed, as will be better discussed in Section 5.2.3.

5. Results and Discussion

The data experimentally acquired in the case study under analysis (Section 5.1) are systematically analyzed to evaluate the effectiveness and the performances of the proposed crack size estimation technique. In particular, the results obtained with the $3 mm$ iFEM model are analyzed in detail in Section 5.2, investigating all the aspects of the technique proposed, while a sensitivity analysis on the mesh size is performed in Section 5.3 including results for the $1 mm$ mesh model.

5.1. Fatigue propagation experimental results

The specimen was tested with a total number of 32,000 fatigue cycles subdivided into 10 different blocks. At the end of each block, both the strains and the crack length have been experimentally measured at the maximum tensile load. In particular, the crack length is reported in Table 1 as a function of the number of cycles. As already introduced, the target of the proposed methodology is estimating the crack length, which is unknown in a general application, thus the experimental measurements reported in Table 1 are exploited for validation of the methodology.

5.2. iFEM results with the 3 mm model

This Section systematically investigates all the results obtained with the $3 mm$ model. In particular, a comprehensive physical description of the Likelihood trend is reported in Section 5.2.1, assessing the basic framework of the methodology. Then, the strain field computed by the iFEM is compared with the independent measurements of DIC and the displacement field with laser measurements in Section 5.2.2. Finally, the crack size estimation from multiple measurements at different crack lengths is reported in Section 5.2.3, assessing the capability of the methodology in a real application scenario.

5.2.1. Likelihood trends

Considering the $3 mm$ mesh, 9 different models are a priori generated: the undamaged model and 8 models with different crack lengths. Crack lengths span from $6 mm$ to $48 mm$ with a step increment of $6 mm$, which corresponds to the double of the element’s size, as detailed in Section 4.3.

Let’s now consider the strain measurements experimentally acquired for a real crack length of $24.3 mm$, which was obtained after $\mathrm{20,000}$ fatigue cycles. The input strain measurements are fed to the 9 iFEM models generated, each one considering its own physics-based strain pre-extrapolation, as summarized in Figure 5 and Figure 6. Each iFEM model will compute a different strain field, which is in turn compared to the test measurements with the logarithmic Likelihood. The Likelihood is computed for a parameter $\sigma =10 \mu \epsilon$, representing a reasonable standard deviation of the strain’s uncertainties for the present application. Notice that, a different value of $\sigma$ does not affect the results of the present methodology and the maximum likelihood estimation framework would provide the same results since only a scale factor would be applied to its trend. The resulting Likelihood trend (Figure 12), which is proportional to the mean square error, resembles the capability of each model of well describing the real damage condition, in particular, the model with the highest logarithmic Likelihood considers a $24 mm$ long crack, which is almost equal to the real one. Finally, models with a different crack length have lower values of Likelihood, resembling a probability distribution.

Then, let’s analyze more in detail the trend of the Likelihood curve obtained. The two iFEM models adjacent to the maximum Likelihood point, i.e. the models with $18$ and $30 mm$ cracks, have almost the same Likelihood value which is slightly lower than the maximum one. This indicates that also these models have a good probability1 of well describing the real damage condition, although this is not the maximum one. Notice that the logarithmic Likelihood is used, thus also a small variation of the coefficient magnitude implies a considerable variation agreement between the test and the iFEM strain fields. The same concept applies also to the models with $12$ and $36 mm$ crack lengths, although they are associated to a significantly lower Likelihood. Finally, the Likelihood trend exhibits a different behavior for smaller and larger crack lengths respectively. For large values of crack length (like $42$ and $48 mm$) the Likelihood exhibits a sudden drop since the test measurements are significantly different from the strains computed by the iFEM. This behavior is induced by the relative distance between the crack tip and the test sensors, which is decreased for larger cracks. On the contrary, a small crack length may induce an almost negligible strain effect in correspondence of the test sensors. For example, the models associated to the undamaged configuration and the $6 mm$ crack length provide almost the same value of Likelihood, since their strain difference is not appreciable in correspondence of the test sensors.

Finally, the strain field computed by iFEM with different models is reported in Figure 12. Notice that the crack’s strain field and the displacement opening of the crack are obtained thanks to the physics-based strain pre-extrapolation, which imposes the strain field of each particular crack length on the input formulation, as described in Figure 6. The input measurements acquired by sensors are representative only of the far-field strain, thus, without the physics-based strain pre-extrapolation, they would not be able to compute the correct displacement and strain field given by the crack. In particular, also if a crack is present in the model, this would not open since the input strain field is not compatible with the model adopted. The pre-extrapolation technique can be avoided only in case all the elements are experimentally covered by input sensors, which is not feasible in a real scenario, in particular in correspondence of a damaged area. In conclusion, the physics-based strain pre-extrapolation is fundamental to obtain a displacement and strain field coherent with the crack length imposed in the model. However, all the results obtained by the models considered are potentially wrong since associated with a damage condition a priori set by the routine. Thus, among all the results computed, the right model is selected with the maximum likelihood estimation framework, obtaining the most accurate description of the damage state for the given mesh discretization.

The procedure described can be repeated several times for different acquisitions performed in correspondence of different crack lengths. For example, the logarithmic Likelihood trend computed for the $32.5 mm$ crack length strain acquisition, which was obtained after $\mathrm{32,000}$ fatigue cycles, is reported in Figure 13. In this case, the experimental crack length does not correspond to any of the iFEM models considered and the maximum Likelihood estimation framework indicates the $36 mm$ crack length model as the best. This model considers a crack length immediately after the real one, however, also the $30 mm$ model is associated with a significant probability.

5.2.2. iFEM strain and displacement fields evaluation

The strain field computed by the iFEM is compared with the independent measurement of the DIC for a further assessment of the methodology, while the iFEM displacement is compared with the lasers’ measurement.

The acquisition performed at 20,000 cycles is considered since the detected crack length by the iFEM ($24 mm$) is almost equal to the target one ($24.30 mm$). The strain field along the load direction is compared in Figure 14, showing how the iFEM reconstruction is qualitatively in good agreement with the DIC observation. In particular, the DIC strain field provides a continuous measurement on the whole plate’s area, only the points in proximity of the crack are not computed due to a speckle’s degradation. This provides an independent strain field validation of the whole damage-sensitive area in locations different from the test sensors exploited by the maximum likelihood estimation framework, finally providing the effectiveness of the methodology.

The laser measurements are used to validate the displacement results, in particular the total deformation of the ROI along the load direction. The experimental measurements provide a total deformation between the unloaded and the loaded condition equal to $0.228 mm$, while the related iFEM result is equal to $0.207 mm$. The laser measurement overestimates of about $10 \%$ the iFEM numerical result and this can be attributed to a small misalignment of the lasers, thus measuring a displacement component not perfectly perpendicular to the target surface, and to experimental uncertainties.

After these considerations, the iFEM displacement and strain fields computed with the methodology proposed are accurate and representative of the actual damage condition of the structure. Thus, its adoption in a Digital-Twin framework allows a high-fidelity representation of the structure and can provide precious information about its operative conditions.

5.2.3. Automatic crack size estimation framework

The crack size estimation procedure is now applied to the whole acquisition history to assess the damage evolution.

The maximum likelihood estimation framework is applied to each strain acquisition, providing the related diagnosis state. The crack lengths detected by the iFEM are compared with the target measurements in Figure 15. In particular, the detected crack length underestimates or overestimates the target measurement since only discrete crack’ steps of $6 mm$ are possible with the present mesh discretization. Nevertheless, the maximum likelihood estimation framework always provides the correct detection giving precious diagnostic information.

However, one main limitation of the proposed technique in a real-time application is the requirement of multiple models with different damage conditions. The iFEM algorithm is particularly attractive for real-time applications since only the strain vector ${\mathit{F}}_{\mathit{F}}$ of Equation (10) is updated at each iteration, while, for a given mesh discretization and sensor network, the matrix ${\mathit{K}}_{\mathit{F}\mathit{F}}$ is constant and can be inverted just one time, resulting in a fast computation of the displacement field. In the present work, an automatic procedure to generate multiple iFEM models has been described (Section 3.2), which can be also exploited to a priori create a database of damage scenarios. However, in a real application, the database cannot consider all the possible damage conditions, since also the damage location is initially unknown. Thus, an attractive procedure is, once the crack center has been identified, generating only the required iFEM models on-demand with the procedure described in Section 3.2. This would allow a more accurate damage diagnosis since the model database may not include a damage scenario representative of the real condition. However, each time the mesh is modified by the damage introduction or evolution, the ${\mathit{K}}_{\mathit{F}\mathit{F}}$ matrix has to be recomputed and reinverted, requiring a significant computational burden. Thus, a compromise has to be found according to the specific application, as summarized with the following approaches:

• The model update and the damage diagnosis are carried out online, where the needed models are generated on-demand, although this may require a significant computation burden
• A database of damage scenarios can be generated offline considering the most relevant damage mechanisms of the structure (e.g. the crack’s nucleation from rivet holes). Then, the real-time algorithm will just load the required models, in which the matrix ${\mathit{K}}_{\mathit{F}\mathit{F}}$ has been already inverted, and perform an efficient damage diagnosis.
• The damage diagnosis procedure is carried out offline at specified intervals, in parallel with respect to the structure’s operation. With this approach, the a priori definition of the damage scenarios and the database creation is not necessary since the computation does not have strict requirements in terms of computational time. The offline procedure computes the new damage condition of the structure, which can be uploaded on the real-time Digital-Twin and used until the next model’s update. The offline damage diagnosis and model update can be carried out on-demand, for example, based on some real-time damage indices, or at specified intervals according to the expected damage growth rate.
• A combination of the approaches described can be adopted according to the specific case. For example, an initial database can be created, and then, once the damage location is identified, additional models are generated on-demand.

At this point, regardless of the specific approach adopted, let’s suppose to perform the real-time damage diagnosis and Digital-Twin update of the structure under analysis. Let’s also suppose that the structure was originally undamaged (thus without any crack) and after a visual inspection in correspondence of the 5,000 cycles acquisition the crack was discovered and its central position identified. At this point, the maximum likelihood estimation framework previously described is performed and the crack length detected is $18 mm$. In case all the 9 models previously presented are included in the computation, a total of 9 iFEM iterations are required. However, if the computation starts from the undamaged condition and the crack size is gradually increased, the maximum likelihood is correctly identified when the model with the $24 mm$ damage size is executed, thus at the 5^th iteration. More specifically, the damage size is systematically increased until the likelihood trend, after an initial increase, starts decreasing, thus the damage condition is identified by the previous iteration. Then, once the crack length in correspondence of the 5,000 cycles has been detected, the measurements acquired at 10,000 cycles are processed. The same procedure can be repeated to identify the new crack length, however, the routine can start from the damage size detected at the previous iteration, thus $18 mm$ in this case. This allows a further reduction of the number of iterations and thus computational time. Moreover, in case the damage scenarios are not loaded from a database, the model computation (i.e. the matrix inversion) of the $18$ and $24 mm$ crack length models is not necessary since performed in correspondence of the previous acquisition. Thus, each time a new model is computed for the first time, it can be stored and reused when necessary. This allows a reduction of the computation time during operation since in case the damage does not significantly propagate no new iFEM models need to be computed. The acquisition performed at $\mathrm{10,000}$ cycles will not provide any damage size increase since the likelihood computed by the $24 mm$ crack size model results lower than the $18 mm$ one. Thus, in this specific case, just two iFEM iterations are computed (with two different models) and no new model is generated and inverted, resulting in a very fast computation. The same procedure can be repeated for each iteration and the introduction of a new model, with the related matrix inversion, is carried out only when a crack length propagation is observed, as also summarized in Figure 16. Thus, according to the complexity of the model and the Digital-Twin requirements, the framework proposed can work in real-time without any a priori definition of the damage conditions, generating the required models on-demand. Finally, the numerical procedure described is also detailed by Algorithm 3.

Algorithm 3: Real-time crack size estimation procedure without a priori creation of a database of damage scenarios.
1:	Define the undamaged iFEM model (mesh, connectivity, BC, etc.)
2:	Set the initial state variables:  $iFE{M}_{model}\left(1\right)$ = undamaged iFEM model (mesh, connectivity, BC, etc.)  ${2a}_{model}\left(1\right)=0$ (zero initial crack length)  $count=0$ (time instant counter)  $2a_{detected}\left(0\right)=0$ (zero initial detected crack length)
3:	Identify the crack center
4:	Define the input (${\epsilon }_{inp}$) and the test (${\epsilon }_{test}$) sensor networks
5:	while experimental data acquisition ON do
6:	$count=count+1$
7:	Acquire the new experimental measurements ${\epsilon }_{inp}\left(count\right)$ and ${\epsilon }_{test}\left(count\right)$
8:	$inde{x}_{model}=$ find ($2a_{model}=={2a}_{detected}(count-1))$ (index of the previous iteration iFEM model)
9:	$iter=1$ (iteration counter)
10:	($\mathit{U}\left(iter\right), {\mathit{\epsilon }}_{iFEM}\left(iter\right)$, $LogL\left(iter\right)$) = Compute the iFEM results ($iFE{M}_{model}$($inde{x}_{model})$, ${\epsilon }_{inp}\left(count\right)$, ${\epsilon }_{test}\left(count\right)$) (Algorithm 4)
11:	$condition=true$
12:	while $condition==true$ do
13:	$inde{x}_{model}=$ $inde{x}_{model}+1$
14:	$iter=iter+1$
15:	if $iFE{M}_{model}\left(inde{x}_{model}\right)$ does not exist do
16:	if $a_{model}\left(inde{x}_{model}-1\right)==0$ do
17:	$iFE{M}_{model}\left(inde{x}_{model}\right)$ = crack introduction from $iFE{M}_{model}\left(inde{x}_{model}-1\right)$ (Algorithm 1)
18:	else
19:	$iFE{M}_{model}\left(inde{x}_{model}\right)$ = crack propagation from $iFE{M}_{model}\left(inde{x}_{model}-1\right)$ (Algorithm 2)
20:	end
21:	Get the new ${2a}_{model}\left(inde{x}_{model}\right)$ from $iFE{M}_{model}\left(inde{x}_{model}\right)$
22:	end
23:	($\mathit{U}\left(iter\right), {\mathit{\epsilon }}_{iFEM}\left(iter\right)$, $LnL\left(iter\right)$) = Compute the iFEM results ($iFE{M}_{model}$($inde{x}_{model})$, ${\epsilon }_{inp}\left(count\right)$, ${\epsilon }_{test}\left(count\right)$) (Algorithm 4)
24:	if $LnL\left(iter-1\right)>LnL\left(iter\right)$ do
25:	${2a}_{detected}\left(count\right)={2a}_{model}\left(inde{x}_{model}-1\right)$ (store the detected crack length)
26:	${\mathit{U}}_{final}\left(count\right)=\mathit{U}(iter-1)$
27:	${\mathit{\epsilon }}_{iFEM}^{final}\left(count\right)={\mathit{\epsilon }}_{iFEM}(iter-1)$
28:	$condition=false$ (exit from the while loop)
29:	end
30:	end
31:	end

Algorithm 4: iFEM computation subroutine with physics-based strain pre-extrapolation.
1:	Pre-extrapolate the strain field with the physics-based approach considering the ${2a}_{model}\left(inde{x}_{model}\right)$ crack length
2:	Define the iFEM input vector $\mathit{F}$
3:	Compute the iFEM displacement field ${\mathit{U}}_{\mathit{F}}$ using Equation (10)
4:	Compute the iFEM strain field ${\mathit{\epsilon }}_{iFEM}$ using Equation (5)
5:	Compute the Logarithmic Likelihood ($LnL$) between ${\mathit{\epsilon }}_{iFEM}$ and ${\mathit{\epsilon }}_{test}$ using Equation (12).

5.3. iFEM results with the 1 mm model

The same methodology is now applied with a finer mesh with a $1 mm$ element size. The likelihood trends are briefly described in Section 5.3.1, while the final crack size estimation routine is presented in Section 5.3.2.

5.3.1. Likelihood trends

The $1 mm$ mesh allows a better discretization of the crack length and so a total number of $25$ models are generated. The first model is related to the undamaged structure (no crack), while the other $24$ consider different crack lengths up to $48 mm$ with a step size of $2 mm$, corresponding to the double of the element size.

The likelihood trends computed for the 20,000 and the 32,000 cycles acquisition are shown in Figure 17 for comparison with respect to Section 5.2.1. In particular, the maximum likelihood is obtained in correspondence of a model that well describes the target damage state, assessing the capability of the technique for a finer mesh. In particular, all the conclusions obtained for the $3 mm$ model are still valid, but with a more accurate detection of the crack length.

5.3.2. Crack size estimation results

The procedure detailed by Algorithm 3 can be now applied to all the acquisitions, resembling a real-time application with an on-demand generation of the damaged scenarios. The results are reported in Figure 18a, showing a fine discretization of the target crack length. More specifically, the root mean square error between the detected crack length and its true value decreases from $2.05 mm$, for the $3 mm$ mesh, to $0.97 mm$ for the $1 mm$ case. Then, the number of iterations required by the algorithm is shown in Figure 18b. In particular, the computation of a new iFEM model is performed in almost every acquisition since the fine mesh is able to appreciate the damage propagation that occurred. Obviously, a finer mesh requires a higher computational time for two main reasons: (i) the higher number of elements and degree of freedoms; (ii) the higher number of model updates required. Thus, in practical applications, a compromise on the mesh size must be chosen and, in case, the most relevant damage scenarios can be generated a priori and stored in a numerical database to speed up the overall procedure.

6. Conclusions

This paper newly extends the load-independent iFEM approach to estimate the structural damage size leveraging on strain measurements acquired on the structure itself.

In case damage occurs on a structure, the iFEM model is no more compatible with the physical structure, inducing a wrong displacement and strain field reconstruction. Thus, the iFEM approach proposed in this work systematically introduces the damage in the iFEM model to decrease this discrepancy. Several iFEM models are generated to account for different possible damage scenarios, thus considering different damage sizes, then a maximum likelihood estimation framework selects the model that better approaches the real damage condition. This framework minimizes the error between the iFEM strain reconstruction and the measurements from test sensors, which are representative of the real structure’s damaged condition. In addition, the methodology presented overcomes some limitations of previous literature works on the damage detection with the iFEM, in particular: (i) it can account for a lower number of strain sensors, (ii) it does not require the definition of a threshold for the damage index evaluation.

The methodology is experimentally validated on an aluminum plate subjected to fatigue crack propagation, where the objective of the analysis is the crack length size estimation. The iFEM methodology exploits also the recently developed physics-based strain pre-extrapolation to obtain the correct strain gradient induced by the crack and, at the same time, relies on a sparse sensor network constituted of an OBR fiber optic. The methodology proposed provides a crack length estimation very close to the one experimentally measured by visual inspection, although only discrete values are possible according to the mesh discretization adopted. The iFEM strain reconstruction is further validated with the independent observations of DIC, showing how the iFEM is able to well describe the strain field near the crack, even in presence of a very sparse sensor network.

Nevertheless, the approach proposed has some limitations which will be further investigated and solved in future works by the authors. First, the present application assumes the crack is detected and its location is a priori known from a visual inspection or other identification techniques. Second, the proposed methodology investigates only regular crack growths in which the crack’s direction is constant, while different case studies (e.g., with composite materials) may exhibit different behaviors. Finally, the proposed approach does not require a priori knowledge of the possible damage conditions (i.e., of the possible crack positions and lengths) since the required models are generated on-demand by the algorithm itself, which is beneficial compared to other diagnosis approaches, e.g. based on machine learning techniques, where the a priori generation of different damage conditions (i.e., a database of damage scenarios) is a fundamental step. However, each time a new damage condition is generated, the iFEM matrix ${\mathit{K}}_{\mathit{F}\mathit{F}}$ should be computed and inverted, resulting in a computationally demanding task. This point in the most significant limitation of the present technique and, in view of realizing a computationally efficient algorithm that can work in real-time applications also for complex structures, different strategies can be adopted. For example, the most relevant damage scenarios can be a priori generated and then specific models, more representative of the actual damage condition, can be computed on-demand. As an alternative, model updating can be run simultaneously to structure operation and, once the most likely model is selected, it will be provided as input to the real-time routine.

These approaches, opportunely combined, can result in a fast routine that can be exploited to develop an iFEM Digital-Twin (iFEM-DT). This is a high-fidelity model of the real structure, which is constantly updated according to the damage evolution, maintenance operations, and other structure modification. Thus, the automatic damage size estimation based on iFEM proposed in this work appears a good candidate to develop new DT frameworks.

In conclusion, this work introduces the possibility to include the damage information in the iFEM model, which can be exploited in an SHM or a DT framework. Future works of the authors will be dedicated to further developing the proposed methodology and assess its capabilities in different scenarios.