Utilization of Digital Twins and Other Numerical Relatives for Efficient Monte Carlo Simulation in Structural Analysis

2.1 The (near to) perfect twin based on multi-physics models

In the present text, the concept of a digital twin is understood in the following sense: A digital twin is a numerical model capable of reproducing the state and behavior of a unique real asset in real time (or faster), with this model also being able to represent the performance of the asset for new and artificially generated conditions (i.e., in connection with extrapolated predictions). As a primary candidate for a digital twin, a complete numerical model based on first principles in terms of multi-physics modeling seems to be most relevant. Such a model will also be able to represent non-linear features of the structural behavior of the asset.

As an example, a relatively simple structure with pronounced non-linear behavior is considered: Figure 1 shows a structure composed of two truss members. The structure is subjected to a vertical load R.

If the geometry is assumed to be non-deformed, the relationship between the vertical load R and the vertical displacement is obtained as:

which for small angles can be approximated by

where α₀ is the slope angle of both truss members.

However, by accounting for changing geometry due to the vertical load, a different relationship between the vertical load, R, and displacement, r, is obtained. By consideration of geometric compatibility, equilibrium conditions and a linear stress-strain relationship, the expression for the load-displacement curve can then be derived as:

where h is the height of the truss and l is half the horizontal span length. E is the modulus of elasticity for the relevant material and A is the cross-section area of both truss members. The model uncertainty associated with this relationship is presently considered to be negligible, such that it can be assumed to represent a “digital twin” of the structure (implying that, e.g., buckling of the members themselves is not relevant due to their non-slender characteristics).

This highly nonlinear expression in Eq. (3) can be closely approximated by a third order polynomial as follows, [1]:

By inserting α₀ = h/l (also assuming small angles), this can be written as

Both of the R-r (i.e., load and displacement) relationships according to Eqs. (3) and (5) are shown in Figure 2 for a slope angle of α₀ = π/15. It is seen that they can barely be distinguished from another. (This implies that the third order representation can also be regarded as a digital twin, although not of the one-egg kind).

Both of the curves are characterized by a very non-linear behavior, where a so-called snap-through occurs when the two truss members are displaced to a completely horizontal position. After snap-through has occurred, a second equilibrium configuration is obtained for which a further increase of the vertical load can take place. However, this second equilibrium configuration will in most cases represent a “failed condition” in the sense that the structure will survive but such that an unwanted large displacement has taken place (which would, e.g., be the case if the structure represents a load-carrying roof structure or an arch system).

Up to around one quarter of the maximum load point, the load-displacement curve is quite close to being linear. Accordingly, if only empirical load-displacement data points for this interval are available, this would typically lead to the assumption that the structural behavior is linear for any load level (unless the physical behavior of the system is taken into consideration). Having available data sets for many different structures of the same type, it is very unlikely that any of the sets contain information about the post-snap interval if all the structures are still in operation.

For structures of the present type in cases where also the stress-strain behavior of the material is nonlinear, numerical solution methods will generally be required in order to compute the load-displacement curve. This will increase the computation time significantly, which will be particularly cumbersome in connection with Monte Carlo simulation procedures where a large number of repeated calculations is typically required (e.g., of the order of millions and upwards). In any case, simplified but “adequate” models need to be introduced.

2.2 More distant numerical relatives based on different kinds of simplified physics-based models

One way of being able to reduce the computation time for even larger and more complex numerical models, is to introduce a simplified representation which is no longer a twin but some more distant numerical relative. As an example, a so-called response surface model can be applied in order to overcome excessive computational efforts. Such models are frequently also referred to as “meta-models” or “cyber-physical” models. One possible approach is to use a response surface representation based, e.g., on first- or second-order polynomials as approximating functions. The parameters of these functions and their weighting coefficients are then determined, e.g., based on minimization of the mean square error. The “control points” for the approximate model are then established based on application of the physics-based model at just these points (i.e., for given input parameter values). By a proper selection of control points, the prediction error associated with the entire range of structural displacement levels can be limited in magnitude.

As examples, we consider approximation of the exact load-displacement relationship with a quadratic and also an alternative quadratic response surface model. For the former, the control points are selected as (0, 0); (0.5, 3/16) and (1.0, 0.0), where 3/16 represents the exact maximum value of the cubic function (but with the location of the maximum point shifted to an abscissa value of 0.5). For the latter, a minimum mean square error approximation within the interval 0.0–1.0. The first of these approximations is compared to the “exact physics-based model” in Figure 3.

The error associated with the second order approximations over the range from r/h = 0 to 1 is seen to be acceptable, while for the less interesting range (within the present context) from 1 to 2 it is highly inaccurate and of little use.

2.3 Data-driven simplified models

A numerical representation of the load-displacement relationship based on a data-driven simplified model is next considered. First, it is assumed that 10 data points in the range from r/h = 0 to 0.15 are available, which is mainly in the weakly non-linear regime. These are, e.g., obtained during normal operation of the structure. A measurement noise with a standard deviation of 10% of the measured signal is also introduced. The extrapolated second order approximation (based on regression analysis) is shown in Figure 4a together with the data points themselves. It is seen that the maximum value of the load R is significantly underpredicted by this curve. As a second approximation, the 10 data points (including noise) are next taken to lie in the range from r/h = 0 to 0.2 (i.e., into the slightly more nonlinear regime). The corresponding second order approximation is shown in Figure 4b. The maximum load is somewhat closer to the true value, but still a significant underprediction is observed.

These results are intended to illustrate the limitations of extrapolations based on data driven models unless the measurement points are available in the region with “high nonlinearity.” For structural systems, such data points are generally scarce as they represent rare events that may even correspond to failure of the structure.

2.4 Comparison of failure probabilities calculated by application of the different numerical models

By introducing a structural failure criterion for the truss in addition to joint statistical models for the inherent random variables, the failure probability corresponding to a given reference period can be computed. Presently, the failure function is expressed in terms of the maximum allowable load (i.e., Rmax), and the only random variable is the external extreme environmental load (i.e., Rex) which follows a Gumbel distribution with a mean value of 0.9 Rmax and a coefficient of variation of 10% (i.e., a standard deviation of 0.09 Rmax). In the present section, a comparison is made between structural failure probabilities, which are obtained by application of the different structural representations that were considered above (Table 1).

Not unexpectedly, the accuracy of the physics-based representations is significantly higher than the data-driven models for the present example. While the cubic response surface almost corresponds to the twin representation, the data driven model for the low loading regime could at best be referred to as a more distant relative (e.g., a half-brother or a cousin).

3.1 System modeling and reliability formulation

The failure of a structural system, e.g., offshore jacket platform is often defined as the total collapse of the structure. The collapse event can be modeled as a series system of several parallel subsystems as follows [2]:

where n is the number of components in the system, N is the number of failure modes, gFij…is the limit state function of component i for failure mode j. The system failure probability for systems like offshore jacket platforms can be accurately estimated by considering a single failure mode and expressing the system resistance R and the system load S in terms of base shear [3, 4]. The system resistance R is the ultimate capacity base shear, which is a function of system damage state’s matrix D. The system load S is the base shear load for a given environmental variable E. The probability of system overload failure for a given system damage state D is calculated shown in Eq. (7).

The performance of structural components in the system deteriorates over time due to, e.g., fatigue damage or corrosion. The system damage state’s matrix D contains the (fatigue) damage state of each component at time t, i.e., D=HF1Yt..HFiYt..HFnYt, where HFiis an indicator function that equals one if the component i fails (i.e., gFijYt≤0) and zero if otherwise. Y is a vector of random variables that influences the fatigue damage (see Chapter 3.2). The total probability theorem is then utilized to calculate the probability of system failure due to both overload and fatigue failures as follows:

Following [3], Eq. (8) can be approximated as follows:

where PFS,ointactis the system failure probability due to overload in the intact condition, PFiis the fatigue failure probability for component i, PFS,oFiis the conditional system failure probability due to overload after fatigue failure occurs at component i, and PFi∩Fjis the probability that fatigue failures occurs at components i and j before the overload failure. Eq. (9) is often referred as annual probability of system failure in the context of structural integrity management, where PFiis defined as the probability of failure at component i given survival up until year t [5]. As a first approximation, the annual probability of system failure can be calculated by keeping only the first two terms [5]:

3.2 Response surface

The system load L is a function of environmental variable vector E. In this work, the wave height H and wave period T are considered as the environmental random variables, i.e., E = [H, T]. The system load L is expressed as the base shear for a given combination of wave height and wave period. The response surface method with quadratic polynomial function is utilized to estimate the system load as follows [6]:

where a₀….a₅ are the coefficients to be determined. Probabilistic linear regression analysis is employed to obtain the coefficients and the predictive distribution of the system load L. The linear model is written as follows:

where L is a (1 × m) vector of “responses” (i.e., which here is the load), β is a (1 × m) vector of the regression coefficients (see Eq. (11)), and ε is the 1 × m vector containing the error terms. The error is assumed Normal-distributed with zero expected value and variance σe2. X is a m × p design matrix which consists of p combinations of individual terms (see Eq. (11)) and m number of samples as follows:

The unconditional predictive distribution is given by the multivariate non-central Students t-distribution i.e., L˜∣L̂∼tm−pEL˜L̂COVL˜L̂with parameters as given as follows [7]:

X̂and L̂are matrices that contain the pre-computed load points from, e.g., finite element analysis. L˜is a vector of load predictions from regression analysis for given X˜, which is calculated, e.g., from the samples of wave height and wave period. The predictive distribution of the load, pL˜L̂, can be seen as a measure of the model uncertainty associated with the response surface.

The wave height is assumed to be Weibull distribution. A special type of a conditional Weibull distribution proposed by Forristal is utilized and written as follows [8, 9]:

where h_s is the significant wave height. Probabilistic models for the wave period are less studied compared to wave height. In the present work, the wave period is assumed to follow a Lognormal distribution, which is conditional on wave height, and the parameters are defined as follows:

where b₁; b₂; b₃ and c₁; c₂; c₃ are the coefficients to be determined. Eqs. (19) and (20) ensure that the wave period is dependent on the wave height in order to avoid drawing unrealistic samples of wave height and period (e.g., very large wave heights with very small wave periods).

3.3 Structural fatigue and reliability updating with monitoring and inspection information

Fatigue failure occurs if the crack size exceeds a critical crack size, and this can be modeled by means of a limit state function as follows:

where δ_c is the critical crack size and δ(Y,t) is the crack size at time t. Here, failure of the structure will occur if this failure function becomes negative. The crack growth is modeled using Paris’ Law as follows:

where m and C are the empirical model parameters, N_s is the number of stress cycles, and ΔK is the stress range intensity factor. For through-thickness cracks on an infinite panel, the solution to Eq. (22) can be written as follows [10]:

where δ₀ is the initial crack size, and ν is the annual cycle rate. B_SIF and B_ΔS are the model uncertainties of the stress intensity factor and for the stress range calculation, respectively (see e.g., [11]). ΔSeis the so-called equivalent stress range and calculated as follows:

Y is a vector of random variables i.e., Y=δ0CBSIFBΔSlnγλ, where γ and λis the scale and shape parameter of the Weibull distributed stress range. The probability of fatigue failure is calculated as follows:

PFitis defined as the probability of annual fatigue failure given survival up until year t. Statistical dependencies between fatigue hotspots are modeled using correlation coefficients of the random variables in the Y vector. There are 6 correlation coefficients: ρδ0;ρC;ρBSIF;ρBΔS;ργ;ρλ. The coefficient ρδ0represents the statistical dependencies due to the same fabrication process. ρCindicates the dependencies due to common material characteristics. ργand ρλdescribe the statistical dependencies due to the similar loading patterns. ρBSIFand ρBΔSdepict the dependencies due to common stress intensity factor and stress range calculation.

Probabilistic models that are able to represent inspection activities in a proper way are also required. Information regarding structural performance can be obtained by carrying out inspection or structural monitoring. There are two outcomes of an inspection: no damage indication (I₁) or damage indication (I₂). The objective of inspection modeling is to obtain the marginal probability of indication (and no indication) followed by an update of the probability of system failure. By utilizing detection theory, the probability of an indication can be derived from the noise and signal distributions (see e.g., [9, 12]). Signal and noise characteristics are typically modeled by means of a Normal distribution (see e.g., [9, 13]). The updating of component fatigue failure probability is performed by utilizing Bayes’ law. Given no indication after an inspection, the probability of fatigue failure is updated as given in [14].

Furthermore, models for statistical representation of the structural monitoring methods are required. Structural health monitoring (SHM) systems can be installed to monitor specific structural properties such as, e.g., vibration or strain in the structural system. Information from a SHM can be viewed as one of the possible realizations of the model uncertainty (see e.g., [15, 16]), which is associated with the measured property such as, e.g., stress ranges. In the present work, the SHM modeling proposed by [17] is employed, i.e., three different possible SHM outcomes are considered: The outcome Z₁ corresponds to the case where monitoring indicates lower stress ranges than expected and indicates that the monitored component has a high performance. Outcome Z₃ indicates that the monitoring component has a low performance due to higher than expected stress ranges. Outcome Z₂ indicates that the monitored component performs as expected. Calculation of the updated probability of system failure is carried out as described in [14].

3.4 Quantification of the value of SIM strategies

The quantification of the value of SIM strategies builds upon the Bayesian pre-posterior decision analysis framework as formulated by Benjamin and Cornell, [18]. A SIM strategy decision problem can be modeled by a decision tree in pre-posterior form as shown in Figure 5. The information space S consists of available information acquirement strategies i (e.g., inspection and monitoring). The outcome space O comprises the possible outcomes of a given information acquirement strategy i. The action space A consists of the possible actions that can be taken such as e.g. repair. The state space θ contains possible states such as, e.g., failure or survival.

The value of SIM strategies is quantified by utilizing the value of information and action (VOIA) analysis (see [19]). A VOIA analysis consists of a base and an enhancement scenario. The base scenario is defined as the scenario without any SHM/inspection and risk-mitigating action such as e.g., repair. There are two states considered in this system state analysis: the (collapse/no collapse) and the component state (failure/no failure). Therefore, the expected cost C₀ in the base scenario is the sum of the expected system E[C_FS] and component E[C_F,i] failure costs over the service life T_SL:

where n is the number of structural components. Procedures for calculation of E[C_FS] and E[C_F,i] are described in [14]. The failure probabilities, which are required in order to calculate these costs, are computed by means of Monte Carlo simulation.

Two different SIM strategies are analyzed as enhancement scenarios. The first strategy is to perform inspections and repair if required during the service life. The second strategy is to install the SHM system for 1 year at one component and to perform inspections and repair if required. For both strategies, a repair action is performed if the inspection indicates a damage. The expected total cost of the first strategy is calculated as the sum of expected costs of inspection E[C_i], repair E[C_R,i], system failure E[C_FS], and fatigue failure E[C_F,i] costs over the service life as follows:

where N_IC is the number of inspected components. E[C_FS] and E[C_F,i] in Eq. (27) are calculated by considering the updated probability of system and fatigue failure, respectively. E[C_R,i] is the expected repair costs over the service with N_IC repaired components.

The expected total costs for the second strategy is calculated as follows:

where ECSHMis the expected SHM costs. Further details are given in [14]. The VOIA is then calculated as follows:

3.5 Case study

A typical deepwater offshore jacket platform with 25 years of service life and located at 190 m waterdepth is utilized in the present work (see Figure 6a). The jacket platform has 200 components and each component is subjected to fatigue deterioration. In this study, each component is assumed to have exactly one hotspot for which a trough-thickness crack will result in fatigue failure. The incoming wave direction is taken as 135° and the 100-year significant wave height H_S;100y equals to 24.3 m.

System failure is defined as the collapse of the jacket platform due to overload and fatigue deterioration. In this example, it is assumed that only single component failure is possible before the overload failure and the probability of system failure is calculated by Eq. (10). The resistance R is the ultimate base shear for given damage state matrix D and is calculated by performing pushover analysis with software USFOS [20]. The system load L is approximated by utilizing response surface analysis outlined in Section 3.2 with 248 pre-computed load points. The samples of L are drawn from the predictive distribution for a given H and T. The coefficients in Eq. (19) and (20) are assumed to have the following values (based on a typical North Sea environment): b = [1.322; 0.8; 0.242] and c = [0.005; 0.09374; 0.32]. The conditional probability of system collapse (see Eq. (7)) is calculated by utilizing Monte Carlo simulation with 10⁶ samples. The response surface of the system load L is shown in Figure 6b with the coefficient of determination R² = 0.9905. Figure 6c shows the predictive distribution of L for T = 20.8 s.

The system failure probability in intact condition is 5E−6. The conditional system failure probabilities for different damaged components and its associated ultimate base shears are given in Table 2. The components in Table 2 are located on the jacket’s legs.

3.5.2 Structural integrity management (SIM)

Two SIM strategies are considered: inspection and repair and inspection with SHM and repair. For the first strategy, inspections are performed at 3 critical components 1 year before the annual system failure probability P(F_S) is estimated to exceed the threshold P_th(F_S) (i.e., constant threshold approach). The minimum system failure probability threshold during operation is set equal to 10−4which corresponds to target system failure probability recommended by JCSS [15] for structures with large consequences of failure and large relative cost of safety measure.

In order to simplify the decision analysis, a repair action is performed only if any damage is detected by inspections. This decision rule is practical and VoI-optimal, see [22]. Repaired components are assumed to behave as components with no damage indication. The probability of indication is derived from the noise and signal distributions. The noise SRis assumed to follow a Normal distribution with zero mean and a standard deviation of 0.5. The signal threshold t_hs is calibrated to the probability of false indication (PFI) of 0.01. The signal S is also Normal distributed with the following parameters:

μSt=0.8+0.1⋅δt,σSt=0.3−0.01⋅δtE30

where δ(t) is the crack size at year t.

In the second SIM strategy, a SHM system for stress range monitoring is installed 1 year before the first inspection is performed, with a monitoring duration of 1 year (i.e., up to the time of the inspection itself). SHM performance is modeled as proposed in [17] by utilizing the stress range model uncertainty BΔS. Two thresholds distinguishing the outcomes Z₁ (low stress ranges), Z₂ (stress ranges as designed) and Z₃ (high stress ranges) are calibrated to target probabilities of P₁^T(F_i) = 1·10⁻⁴ and P₂^T(F_i) = 1·10⁻³, respectively. The target reliabilities are selected following [13] for structures with minor consequences of failure with normal and large relative cost of safety measure. The time dependent threshold’s calibration is illustrated by Figure 7. The measurement uncertainty U is assumed Normal distributed with the expected value of 1.0 and a standard deviation of 0.05. A summary of the probabilistic models is shown in Table 4.

Variable	Dimension	Distribution	Expected value	Std. deviation
SR	—	Normal	0	0.5
S	—	Normal	μS(t)	σS(t)
PFI	—	—	0.01	—
P1T(Fi)	—	—	5·10−4	—
P2T(Fi)	—	—	5·10−3	—
U	—	Normal	1.0	0.05

Table 4.

Summary of the random variables for inspection and SHM modeling.

The costs considered in this case study are consisting of inspection costs C_I, SHM costs C_SHM, repair costs C_R, system failure costs C_FS, and component failure costs C_Fi. SHM costs are further divided into investment costs C_Inv^SHM, installation costs C_Inst^SHM, and operational costs C_op^SHM. The cost model used in this example is shown in Table 5 based on [23, 24].

Type	Cost
CI	0.001
CInvSHM	1.33·10−4/channel
CInstSHM	1.33·10−4/channel
CopSHM	2·10−4/year
CR	0.01
CF,i	1
CFS	100

Table 5.

Cost models used in the case study.

3.5.3 Results

The annual component and system failure for t = 1....25 years are shown in Figure 8. The system failure probability for the intact condition is P(F_S,O|D = 0) = 5E−6. The annual system failure probability at the end of service life is 9.5E−5, which is less than the minimum operational threshold, i.e., no SIM implementations are required to achieve the minimum operational requirement. However, the decision-maker may wish to increase the structural reliability above the minimum requirement and presumably enhance the value of SIM. In this work, three different annual system failure probability thresholds are studied: 6E−5, 7E−5, and 8E−5.

For the inspection-only strategy (S = S₁), inspections and repairs are performed at three components 1 year before the system failure probability threshold P_th(F_S) is predicted to be reached. After each inspection, probabilities of fatigue and system failure are updated. Figure 9a shows the annual system failure probability for the inspection-only strategy as a function of time for a specific threshold value. The first inspection time is at t = 8 years. Increasing the threshold means that the inspections are performed later during the service life. Because of this, the inspection frequency during the service life is decreasing with a higher threshold in exchange of a higher annual system failure probability.

The second SIM strategy (S = S₂) is based on installation of a SHM system at one component to monitor stress ranges for 1 year before the first inspection. There are three possible outcomes for monitoring based on the component performance.Figure 9b shows the annual system failure probability with monitoring and inspections for a system failure probability threshold of 6E-5. Compared to Figure 9a, it is observed that the outcome of monitoring can influence the future inspection schedule. The SHM installation time differs depending on the thresholds (see Table 6). A higher threshold means that the SHM system is installed closer to the end of service life. With increasing annual system failure probability threshold, the probability of obtaining low performance outcome (Z₃) becomes higher.

SHM outcome	Annual system failure probability threshold
	6.00E-05 tSHM = 6 years	7.00E-05 tSHM = 8 years	8.00E-05 tSHM = 11 years
High performance Z1	0.175	0.0991	0.0497
As-expected performance Z2	0.553	0.5036	0.4084
Low performance Z3	0.272	0.3973	0.5419

Table 6.

The probability of SHM outcomes for different annual system probability thresholds.

t_SHM is SHM installation time for a given annual system failure probability threshold.

The values of information and action (VOIA) of the two SIM strategies are shown in Figure 10. It is observed that increasing the system failure probability threshold will reduce the value of information and action. With a higher threshold, inspection and monitoring are performed later in the service life, which reduces the benefits due to a higher annual system failure probability during the remaining service life compared to a lower system failure probability threshold. It is also observed that the VOIA of the SIM strategy SHM, inspection and repair (S₂) is higher than the inspection-only strategy (S₁) for all investigated system failure probability thresholds. This shows that information from SHM system can enhance the value of the SIM, i.e., reduce the expected total cost. In this example, the cost of system failure is dominating the expected total cost over the service life.