Real-Time Data Assimilation in Welding Operations Using Thermal Imaging and Accelerated High-Fidelity Digital Twinning

Abstract

Welding operations may be subjected to different types of defects when the process is not properly controlled and most defect detection is done a posteriori. The mechanical variables that are at the origin of these imperfections are often not observable in situ. We propose an offline/online data assimilation approach that allows for joint parameter and state estimations based on local probabilistic surrogate models and thermal imaging in real-time. Offline, the surrogate models are built from a high-fidelity thermomechanical Finite Element parametric study of the weld. The online estimations are obtained by conditioning the local models by the observed temperature and known operational parameters, thus fusing high-fidelity simulation data and experimental measurements.

1. Introduction

Welding is used extensively in the nuclear industry, for assembly and as a repair technique. It is often used in maintenance operations of different kind that involve various geometries and welding parameter settings. Very high temperatures applied in a localized zone cause expansion and nonuniform thermal contractions, resulting in plastic deformations in the welding and its surrounding areas. Thus, residual stresses and permanent deformations are produced in the welded structure. These could induce a variety of defects such as hot tearing/cracking if the process is not properly controlled. Other defects may appear such as porosity, lack of fusion or lack of penetration that need to be identified after an operation.

The detection of defects in weld beads is usually performed after the welding operation is fully performed [1]. When a defect is identified, the entire operation needs to be done again. Therefore, it would be desirable to obtain real-time estimations of the current mechanical state of the assembly using in situ measurements. With such estimations at hand, welding operations could be stopped and/or controlled whenever predicted or forecasted mechanical states are outside acceptable tolerance regions.

To this end, we propose to develop a digital twinning approach [2] whereby high-fidelity model predictions will be continuously adjusted with respect to thermal images acquired online during the welding operation. The simulation will rely on a state-of-the-art mesoscale transient thermoelastoplastic finite element model. The fusion between FEA and in situ measurements will be done by accounting for well-chosen parametric sources of uncertainties in the simulation model, leaving freedom for the digital twin to react and adapt to the sequence of thermal images. We aim for the data assimilation to be done following a statistical and, if possible, Bayesian framework [3], to enable incorporating engineering knowledge about the uncertain parameters of the model, and allow deriving credible regions for the predictions and forecasts of mechanical states.

Unfortunately, data assimilation problems, i.e., sequential inverse or sequential calibration problems, are notoriously expensive to solve when the numerical models are systems of partial differential equations [4,5,6]. In the context discussed above, the parameters of spatially detailed nonlinear FEA models would need to be optimized in real-time. This is intractable. We will therefore develop an appropriate piecewise linear meta-modeling technique to achieve real-time efficiency.

Over the last decades, surrogate modeling via model order reduction has been successfully developed for a variety of applications relying on high-fidelity modeling. This is especially true for data-driven approaches based on projecting the high-fidelity model in low-dimensional subspaces [7,8,9,10,11,12,13], among which POD-generated subspaces [14,15,16] can be seen as an extension of the principal component analysis to continuous variables. Projection-based model order reduction methods are known to reduce the computational complexity of high-fidelity models by precomputed candidate solutions corresponding to various points in the parameter domain (snapshot), and reusing the generated information to fasten online solution procedures. However, even by using efficient hyper-reduction schemes [17,18,19,20], reduced simulations of welding operations are still computationally too demanding for process monitoring. This is in particular due to (i) the lack of reducibility of moving heat source problems in general [17,21,22] and (ii) the relatively high online cost associated with hyper-reduced models (hyper-reduction generalizes well in large parameter domains owing to the fact that in the online phase, the nonlinear equations of the original simulation model are solved on a reduced integration domain [23,24]).

To circumvent these difficulties, we propose to develop an offline/online meta-modeling technique based on a mixture of probabilistic principal component analysis models (PPCA). For any given position of the heat source, the thermomechanical state, augmented by the vector of uncertain parameters, will be postulated to follow a Gaussian distribution with low-rank covariance structure. This model will be identified using the method of snapshots. Offline, we will run the high-fidelity mechanical simulations corresponding to a fine sampling of the parameter space. In a second step, the Gaussian model will be identified using the maximum likelihood approach probabilistic PCA described by Bishop [25]. Online, parameter estimation will reduce to Gaussian conditioning (i.e., the Kalman method), which can be made efficient when the covariance matrix exhibits a low-rank structure. We will show numerically that this strategy allows us to successfully, and for the first time, set up a data assimilation framework for welding operations, blending high-fidelity thermomechanical simulations and thermal imaging in real-time to predict and forecast mechanical states.

In a second stage of developments, we will treat the case where the position of the heat source is to be estimated from thermal imaging. This is of practical interest when the position of the welding torch is not accurately known or tracked during joining operations. To achieve this goal, we build, numerically and offline, a correlation model statistically linking the position of the hottest spots detected on the thermal image to the actual position of the welding torch. Then, the data assimilation method based on the PPCA can be easily adapted, using a probabilistic mixture of PPCA instead of a single Gaussian model. Of course, conditioning the mixture of PPCA remains analytically tractable, and computationally efficient as, as it will be shown, few PPCAs are associated with non-vanishing coefficients at any given time (i.e., the estimation of the position of the torch from thermal images is accurate). Even in this setting, the low-rank structure of the mixture of PPCAs ensures that data assimilation is performed in real-time, without sacrificing the accuracy provided by the high fidelity thermomechanical model.

In all cases, future states may be predicted by conditioning future statistical models to available observations. This is technically done by Gaussian (mixture) conditioning of all future mechanical states to current posterior distributions of unknown model parameters.

The approach proposed in this paper is closely related to other data assimilation methods available in the literature. Filtering methods based on parametrized models (e.g., Kalman or particle filters) typically construct a Markov model for the propagation of parameter distributions in time, progressively assimilating data as they are made available. Our approach can be seen as a degenerate such filter whereby past data are ignored, only the current thermal image influencing the posterior distribution of unknown model parameters. Taking past data into account can be done, for instance, by concatenating mechanical states from current and past assimilation times over a sliding window (i.e., an autoregressive model). For the particular experimental setup considered in this work, taking past data into account brought no significant change in predictive parameter distributions, which justifies our development of a past-agnostic assimilation method. In terms of meta-modeling, we could have used nonlinear meta-modeling techniques such as polynomial chaos or neural network regression [26,27,28], but the choice of a linear model (PPCA) allows us to solve analytically the conditioning problem, thus bringing robustness to the method, which cannot be expected when using Markov Chain Monte Carlo solvers or sequential importance sampling [5]. Our approach is clearly inspired by the parametrized background data-weak approach [29], which adopts a variational point of view, while our method is Bayesian, thereby delivering credible intervals and not point estimates. Similar to the parametrized background data-weak method, we perform state estimation in large dimensions by using a background covariance matrix generated by parametric variations of a high-fidelity PDE system. The low rank structure of this covariance matrix is used to fasten online Gaussian conditioning, circumventing the usual $N^3$ complexity issue by making use of standard algebraic techniques [25,30].

Our paper is organized as follows. In Section 2, we present the experimental configuration of our test case. Section 3 aims to present the thermomechanical model and the inverse problem, detailing the known and unknown parameters. In Section 4, we will introduce the construction of the local surrogate models using Probabilistic PCA, and in Section 5, two different use cases are considered: a situation where the heat source position is known and another one where it needs to be estimated. Forecasting is also discussed. Finally, in Section 6 we present results for all the configurations introduced in Section 5 using noisy simulation data when the torch position is known and real experimental data when it is not. Appendix A, Appendix B and Appendix C give, respectively, details on the algebraic expressions used to accelerate the computations, some more forecasting results and material parameters of the specimen.

2. Experimental Configuration

To test the proposed approach, we need to find an application that can replicate similar welding conditions (in terms of materials and parameter variety) to those of a real maintenance operation that can be performed in a laboratory environment with proper instrumentation.

The target application is the stress prediction in a 316L stainless steel (The chemical composition and some of the temperature dependent material parameters are given in Appendix C) specimen submitted to Programmierter–Verformungs–Risstest (PVR) hot cracking tests. PVR tests, developed during the 1970s by the Austrian company Boehler, consist of making a fusion line using bead-on-plate TIG-welding with argon shielding while the specimen is uniaxially tensile loaded [31]. It allows the control of the tensile deformation progressively, which means that it can show cracks of different origin, most notably solidification and liquidation cracking. The large surface of the specimen allows easy positioning of thermocouples on both sides of it, as shown in Figure 1, and it is also possible to film the surface of the specimen with an infrared camera. Thus, process parameters and infrared images of the welding operation are the observational data from which 3D stresses and temperatures will be estimated.

The experiments are performed with a 6-axes Panasonic robotic arm and equipped with a ValkWelding torch and the tensile loading is applied with a Lloyd Instruments LS100 plus. Two parameters related to the heat source are fixed before the experiment: the speed (v) and the power (Q), which is the product of voltage (U) and current (I). The experiments are instrumented with two thermocouples and filmed with a SC7500 FLIR infrared camera. The infrared video will provide real-time measures on part of the specimen. While both sides of the specimen have thermocouples, only the surface that is not being welded is filmed. This way, reflections from the welding arc and the robot itself are avoided.

3. Digital Twin

3.1. Thermomechanical Model

3.1.1. Thermo-Elasto-Plasticity

Numerical simulation of welding is a very complex problem, as it needs to take into account a great number of parameters to represent multiscale and multiphysics phenomena, using temperature dependent material properties that are not always properly known. Usually, the interactions between the metallurgy, the heat, and the mechanical problems need to be simulated. In the case of 316L stainless steel, the metallurgic interactions are negligible [32], and thus the model is reduced to a weakly coupled nonlinear parametric thermo-elasto-viscoplasticity problem.

We consider a model of unsteady thermo-elasto-viscoplasticity over spatial domain $\Omega \in {\mathbb{R}}^3$, whose boundary will be denoted by $\partial \Omega$, and time domain $\mathcal{T}\in [0,T]$. For all $(x,t)\in \Omega \times \mathcal{T}$, the unsteady heat equation reads as

$$\rho c_p\frac{\partial T}{\partial t}-\nabla ·k\nabla T=q_d$$

(1)

where $T:\Omega \times \mathcal{T}\to \mathbb{R}$ is the temperature field, $\rho$ is the mass density, k is the thermal diffusion coefficient and $c_p$ is the specific heat capacity.

The above equation is complemented by boundary condition

$$T=T_0$$

(2)

on domain boundary $\partial {\Omega }_T$. Moreover,

$$-k\nabla T·n=\gamma (T^4-T_0^4)+\zeta (T-T_0)$$

(3)

on domain boundary $\partial {\Omega }_q$, where $\gamma$ is the radiation coefficient and $\zeta$ is the thermal convection coefficient. At time $t=0$, the temperature T is equal to $T_0$ everywhere in the computational domain.

For all $(x,t)\in \Omega \times \mathcal{T}$, the mechanical equilibrium reads as

$$\nabla ·\sigma =0$$

(4)

This equation is complemented by boundary conditions

$$u=u_d$$

(5)

on part of the domain $\partial {\Omega }_u$, and

$$\sigma ·n=0$$

(6)

on part of the domain $\partial {\Omega }_t$. The coupled thermomechanical problem is closed by the following constitutive relation:

$$\sigma =C:\left(ϵ\left(u\right)-{ϵ}^p-\alpha (T-T_0)I_d\right)$$

(7)

where $ϵ\left(u\right):=\frac{1}{2}\left(\nabla u+\nabla u^T\right)$ is the total strain, ${ϵ}^p$ is the plastic strain and $\alpha (T-T_0)I_d$ is the thermal strain, with $\alpha$ the coefficient of thermal expansion and $I_d$ the identity tensor. The plastic strain is fully determined by a Von Mises plasticity model with isotropic and kinematic hardening.

$$f(\sigma ,R,X):=\sqrt{\frac{3}{2}(\tilde{\sigma }-X):(\tilde{\sigma }-X)}-R\left(p\right)\le 0$$

(8)

$${\dot{\epsilon }}^p=\lambda \frac{df}{d\sigma }$$

(9)

$$\dot{p}=\lambda \frac{dR}{dp}$$

(10)

$$X=\frac{2}{3}C\left(p\right)\alpha$$

(11)

$$\dot{\alpha }={\dot{\epsilon }}^p-\gamma \left(p\right)\alpha \dot{p}$$

(12)

$$\lambda \ge 0\lambda f\left(\sigma \right)=0$$

(13)

where $R\left(p\right)$ is the limit of the elastic domain due to isotropic hardening, $\lambda$ is the plastic multiplier and $\tilde{\sigma }=\sigma -\frac{1}{3}Tr\left(\sigma \right)I$ is the deviatoric part of the stress tensor.

In the specimen of a PVR test, the traction boundary conditions are enforced by Dirichlet boundary conditions:

$$u(x,t)=U_d\frac{t}{t_w}\forall x\in \partial {\Omega }_u$$

(14)

The main mechanical quantity of interest is the first principal stress, denoted by ${\sigma }_I$:

$${\sigma }_I=\underset{n\in {\mathbb{R}}^3}{max}\frac{n·\sigma ·n}{{||n||}^2}$$

(15)

The first principal stress is the highest principal stress. It has a huge influence on hot cracking during the welding process.

3.1.2. Thermomechanical Load

The choice of a model for the moving heat source is a key point in numerical simulation of welding. The most commonly used one is Goldak’s double-ellipsoid model [33], represented in Figure 2. It is important to note that the heat distribution is different in the front and the rear of the heat source, thus the model depends on the current central position of the heat source $P\left(t\right)=(x\left(t\right),y\left(t\right),z\left(t\right))\in {\mathbb{R}}^3$,

$$P\left(t\right)=P_0+\left(\frac{t}{(P_F-P_0)V}-{ϵ}_P\left(t\right)\right)\left(P_F-P_0\right)$$

(16)

In the equation above, $P_0\in {\mathbb{R}}^3$ is the position of the heat source at the start of the welding operation, and $P_F\in {\mathbb{R}}^3$ is the position at the end of the operation. $T=(P_F-P_0)V$ corresponds to the total welding time, which depends linearly on the velocity V of the source. ${ϵ}_P\left(t\right)$ is a small (normalized) time delay that accounts for potential errors in the control of the welding robot (This is a very simple model, which could easily be modified to include fluctuations of the velocity field.). This parameter will be identified online.

Assuming that the heat source moves in the direction x, the final equation is

$$q(x,y,z,t)=\left\{\genfrac{}{}{0pt}{}{\frac{12\sqrt{3}\eta Q}{(a_r+a_f)bc\sqrt{{\pi }^3}}exp\left(\frac{-3{(x-X)}^2}{a_f^2}+\frac{-3{(y-Y)}^2}{b^2}+\frac{-3{(z-Z)}^2}{c^2}\right),x\ge X.}{\frac{12\sqrt{3}\eta Q}{(a_r+a_f)bc\sqrt{{\pi }^3}}exp\left(\frac{-3{(x-X)}^2}{a_r^2}+\frac{-3{(y-Y)}^2}{b^2}+\frac{-3{(z-Z)}^2}{c^2}\right),x<X.}\right.$$

(17)

where $a_r$, $a_f$, b and c are unknown parameters describing the geometry of the double ellipsoid (In practice, we calibrate $a_r$, b, c, $\eta$ and K where K is the ratio between $a_r$ and $a_f$: $a_r=K{a}_f$), $Q=UI$ is the power, with U the voltage, I the current and $\eta$ the efficiency.

3.1.3. Space and Time Discretization

The thermomechanical problem is discretized in space using the standard $P1$ Lagrange finite element method. In the following, $T\left(t\right)$ and ${\sigma }^{\square }\left(t\right)$ will respectively denote the vector of finite element nodal temperature and the vector of components of the stress tensor at the quadrature points.

The thermomechanical problem is discretized in time using the standard backward Euler finite difference scheme.

3.2. “Truth” Online Inverse Problem

We now describe the problem of data assimilation and online forecasting for the welding operation.

3.2.1. Parametrized Probabilistic Setting

The power Q of the heat source and its velocity V are supposed to be controlled with a good degree of accuracy. The mechanical load $U_d$ is also well controlled during the experimental procedure. In order to build the “truth” digital twin, we further assume that the thermal capacity and diffusivity of the material are well characterized. The mechanical behaviour of the structure (thermal expansion, elasticity and plasticity) are also assumed to be well characterized, qualitatively and quantitatively. This is consistent with EDF’s decades of experience in modeling, characterizing and simulating such welding processes. Finally, the thermal and mechanical boundary conditions are reasonably well quantified in our experimental setting.

However, several sources of uncertainties negatively affect the predictive capabilities of the simulation model:

• The position P of the center of the heat source is not perfectly well known. We will consider that parameter ${ϵ}_t$ is random, i.e., is a source of epistemic uncertainty.
• The spatial distribution of the heat source is not known with precision. We surmise that that the main contribution to the overall uncertainty of the model is the spatial length-scale of the Goldak model.

Therefore, the parameter space $\mathcal{P}=\mathcal{M}\times \Theta$ space comprises two distinct blocks:

• The set of known parameters $\mu =\{Q,V\}\in \mathcal{M}$, which will vary from one welding operation to the next, but can be controlled with precision.
• The set of unknown parameters $\theta =\{a_r,b,c,\eta ,K,\}\in \Theta$, which are given a probability distribution $\theta \sim p_{\theta }$ that encodes the prior knowledge available about these parameters.

3.2.2. “Truth” Online Bayesian Conditioning

Data will be assimilated at homogeneously distributed times $\overline{\mathcal{T}}=\{t^0=0,t^1=\Delta t,t^2=2\Delta t,\dots t^{N_t}=N_t\Delta t\}$. The assimilation time step $\Delta t$ is adjusted so that the number of assimilation steps $N_t$ is the same for all simulations, independently of the velocity V of the heat source,

$$\Delta t=\frac{(P_F-P_0)}{N_t}V$$

(18)

Online at time $t^k\in \overline{\mathcal{T}}$, we assume that surface temperatures are measured noisily (i.e., with the thermal camera). We will write

$$d_T^k=H_T{T}_{|\theta ,\mu }^k+{ϵ}_s^k$$

(19)

where H is a fixed Boolean operator acting on the vector of finite element temperature nodal values $T^k$ at time step k. The additive measurement noise is supposed to be zero-mean Gaussian distributed ${ϵ}_T\sim \mathcal{N}(0,{\Sigma }_M={\sigma }_M^2{I}_d)$. The sources of the temperature measure uncertainty are varied and might include light reflections or lack of knowledge of material parameters like the emissivity at different temperatures [34]. The measure error parameter ${\sigma }_M^2$ will be calibrated empirically comparing the measure from two sensors available in EDF’s welding lab.

Assuming that the successive noise vectors ${\left\{{ϵ}^k\right\}}_{k\in [0,t]\cap \overline{\mathcal{T}}}$ are independent, the statistical inverse problem to be solved online is the following.

At time $t\in \overline{\mathcal{T}}$, the posterior distribution of unknown parameters given the past measurements is

$$p(\theta |d^k,{\mu }_k)\propto {\Pi }_{k^{\prime }\le k}{\mathcal{L}}^{k^{\prime }}(\theta ;d^{k^{\prime }})p_{\theta }\left(\theta \right)$$

(20)

where

$${\mathcal{L}}^k(\theta ;d^k)=\mathcal{N}(d^k;H_T{T}_{|\theta ,\mu }^k,{\Sigma }_T)$$

(21)

Unfortunately, evaluating the likelihood function in real-time is unfeasible, even when using standard Markov assumptions.

4. Real-Time Predictions with Gaussian Mixture of Local Surrogate Models

In the previous section, the truth inverse problem was presented, and it was concluded that it is not well suited for real-time applications. To overcome this, in this section the construction of Gaussian local surrogate models from snapshots of a Finite Elements parametric study is discussed.

4.1. Local Multiphysics PCA Model

Our proposal is to make local surrogate models for a linear joint representation of the state, known parameters and unknown parameters for every position of the heat source at the assimilation times. The positions are fixed in a grid between the initial position $P_0$ and the final position $P_F$ with a regular separation $\Delta P$. The parametric solutions are clustered by the position $P\left(t\right)$ of the heat source and associated to the position $P^k=k\Delta P$ if $|P\left(t\right)-P^k|<\frac{\Delta P}{2}$. The local model at position $P^k$ reads as

$$s^k=\left(\begin{array}{c}{T}^k\\ {\sigma }^k\\ {\mu }^k\\ {\theta }^k\end{array}\right)={\overline{s}}^k+{\Phi }^k{\alpha }^k+{ϵ}_s^k$$

(22)

with ${\alpha }^k\sim \mathcal{N}(0,I_d)$ and ${ϵ}_s^k\sim \mathcal{N}(0,{\sigma }^2{I}_d)$. Operator ${\Phi }^k$—a decoder—is fixed for each position and possesses $n_{\varphi }$ columns. It is obtained by using all parametric solutions of the welding problem corresponding to position $P^k$. This probabilistic model encodes the dependency between all the state variables and the known and unknown parameters in the form of a multivariate Gaussian:

$$s^k\sim \mathcal{N}({\overline{s}}^k,{\Phi }^k{\Phi }^{kT}+{\sigma }_s^2{I}_d)$$

(23)

The choice of using local models aims at avoiding the use of a global reduced basis, as it is well known that moving heat sources generate irreducible parametric solutions [17]. Notice as well that there is no Markovian assumption in the model. In other words, measurements at times ${\left\{t_l\right\}}_{l<k}$ do not provide any new information about the state at time k.

Using Gaussian assumptions on ${\alpha }^k$ and ${ϵ}_s$ follows a certain logic as well. With $s^k$ being Gaussian, we can deduce that $d^k$, the observed data at time k is also Gaussian distributed. This means that $p(s^k|d^k,\mu )$ will also be Gaussian and can be calculated analytically, allowing us to greatly accelerate the computation times.

4.2. Maximum Likelihood PCA

In order to calibrate the surrogate model, we assume that a snapshot of $n_s$ extended states is available, which we denote by $\mathcal{S}=\{S^1,S^2,\dots ,S^{n_s}\}$. This snapshot should cover the time and parameter domains appropriately.

In 1999, Tipping and Bishop [35] showed that a probabilistic formulation of PCA can be obtained from a Gaussian latent variable model, where the basis vectors are maximum-likelihood estimates. Given Equation (22), there are explicit expressions for the maximum likelihood estimates of the parameters:

$$\Phi =U_q\sqrt{\left({\Lambda }_q-{\sigma }^2{I}_d\right)}R$$

(24)

where $U_q$ contains the left singular vectors associated to the greatest $q<n_s$ singular values of a Singular Value Decomposition (SVD) of the snapshot matrix $\mathcal{S}$. The diagonal matrix ${\Lambda }_q$ is also obtained from the SVD, and it contains the greatest q singular values in decreasing order. R is an orthogonal rotation matrix that, in practice, will be chosen to be the identity matrix. Finally, ${\sigma }^2$ is related to the truncation error:

$${\displaystyle {\sigma }^2=\frac{1}{n_s-q}\sum _{j=q+1}^{n_s}{\lambda }_j^2.}$$

(25)

where ${\lambda }_j$ are the smaller singular values.

This means that the generative PCA model can be easily computed from a SVD of the output of a parametric study of the thermomechanical finite elements model.

The snapshots $\mathcal{S}=\{S^1,S^2,\dots ,S^{N_t}\}$ are generated using a straightforward procedure. All the parameters in $\mathcal{P}$ are assumed to be uniformly distributed over a hyper-cube and will be sampled using Latin Hypercube Sampling [36]. The minimum and maximum values for each parameter are given by experts.

The samples are generated as follows. For every $(\mu ,\theta )$ sampled with LHS the finite element simulation is ran over $\mathcal{T}\left(\mu \right)=[0,T(\mu \left)\right]$, delivering a history of $N_t$ snapshots. Additional postprocessing could be performed depending on the physics that are represented in the state vector s.

5. Online Prediction

In this section, we treat the resolution of the inverse problem. Two cases have to be considered depending on whether the torch position is known or not: The first case, when the heat source position is known at all times, is straightforward and the predictions are obtained using a single well-identified PPCA model. In the second case, the torch position needs to be estimated with some uncertainty. The predictions are now computed from a mixture of PPCAs. Finally, the forecast of future states is discussed.

5.1. Known Position

Let us assume that the heat source position is known at all times. This happens when all measurements are synchronized and knowing how much time has passed from the starting point allows perfect knowledge of the position $P^k$ via the speed or when the robotic arm is equipped with a sensor measuring the advanced distance and this signal is available. This corresponds to a situation where ${ϵ}_P\left(t\right)=0$ in Equation (16).

When the position $P^k$ is known, the closest model is also known and it is enough to estimate the current state and unknown parameters. This is done by computing the posterior distribution $p(s^k|d^k,\mu )$, which is Gaussian as seen in Section 4. It can be computed analytically and it is completely determined by its mean $m_{s^k|d^k,\mu }$ and covariance ${\Sigma }_{s^k|d^k,\mu }$ and it. Introducing the notation ${\Sigma }_{s^k}={\Phi }^k{\Phi }^{kT}+{\sigma }_T^2{I}_d$, the mean and covariance are

$$\begin{array}{cc}\hfill {\mu }_{s^k|d^k,\mu }& ={\overline{s}}^k+{\Sigma }_{s^k}{H}^{kT}{\left(H^k{\Sigma }_{s^k}{H}^{kT}+{\sigma }_m^2{I}_d\right)}^{-1}(d^k-{\overline{d}}^k)\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill {\Sigma }_{s^k|d^k,\mu }& ={\Sigma }_{s^k}-{\Sigma }_{s^k}{H}^{kT}{\left(H^k{\Sigma }_{s^k}{H}^{kT}+{\sigma }_m^2{I}_d\right)}^{-1}{H}^k{\Sigma }_{s^k}\hfill \end{array}$$

(27)

where $H^k$ is a Boolean operator acting as the observation function.

The evaluation of this expression is very slow due to the size of the matrices involved. To avoid this cost, we propose to compute the posterior distribution of the reduced coordinates ${\alpha }^k$ instead. Once again, the posterior distribution $p({\alpha }^k|d^k,\mu )$ is Gaussian, so it is determined by its mean $m_{{\alpha }^k|d^k}$ and covariance ${\Sigma }_{{\alpha }^k|d^k}$:

$$\begin{array}{cc}\hfill m_{{\alpha }^k|d^k,\mu }& ={\Phi }^{kT}{H}^{kT}{\left(H^k{\Sigma }_{s^k}{H}^{kT}+{\sigma }_m^2{I}_d\right)}^{-1}(d^k-{\overline{d}}^k)\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill {\Sigma }_{{\alpha }^k|d^k,\mu }& =I_d-{\Phi }^{kT}{H}^{kT}{\left(H^k{\Sigma }_{{\alpha }^k}{H}^{kT}+{\sigma }_m^2{I}_d\right)}^{-1}{H}^k{\Phi }^k\hfill \end{array}$$

(29)

The mean of the posterior distribution of $s^k$ can be then deduced using the reduced basis ${\Phi }^k$ as follows:

$${\mu }_{s^k|d^k,\mu }={\overline{s}}^k+{\Phi }^k{m}_{{\alpha }^k|d^k}$$

(30)

As for the covariance matrix, it can be calculated by

$${\Sigma }_{s^k|d^k,\mu }={\Phi }^k{\Sigma }_{{\alpha }^k|d^k,\mu }{\Phi }^{kT}+{ϵ}_s^k$$

(31)

The whole covariance matrix might be impossible to store in RAM if the number of degrees of freedom is large. In this case, only the diagonal of the matrix product ${\Phi }^k{\Sigma }_{{\alpha }^k|d^k,\mu }{\Phi }^{kT}$ is calculated.

For further details on the acceleration of these computations, see Appendix A.

5.2. Unknown Position

As it is the case in EDF’s lab, we may not have access to the exact heat source position for lack of synchronicity between the measure sensors. In this case, ${ϵ}_P\left(t\right)>0$ in Equation (16) and the torch position needs to be estimated from the video frame, which adds more uncertainty to the model. The filmed side of the specimen is the opposite side of the weld and there exists a delay between the highest temperature on this side and the torch position. Furthermore, this delay is dependent on material and operational parameters, such as the speed.

Using the finite elements simulations used for the prior generation of the local PPCA surrogate models, we created a Gaussian surrogate model that links the position of the highest temperature measured on the camera side of the specimen $x^k$ and the known parameters $\mu$ with the position of the heat source on the welded surface of the specimen $P^k$. The model was fitted using a Matérn kernel with $\nu =\frac{3}{2}$ as a covariance function. The output of the model is the probability distribution of the estimated heat source position ${\widehat{P}}^k$. We sample this distribution to find possible torch positions.

Not being able to determine the exact position, the data assimilation needs to be adapted using a mixture of PPCAs. The multiphysics state for this video frame is named $s^{\widehat{k}}$ to indicate an unknown position. For each sample of the heat source position we assign it the closest local PPCA model and a discrete assignment probability distribution is calculated to find the active coefficients of the mixture. Assuming a total of J local models are active, $p(s^{\widehat{k}}|d^k,\mu )$ is computed as a Gaussian mixture of the estimations of each local model weighted by the discrete assignment probability $p_j\forall j\in J$:

$$\begin{array}{cc}\hfill m_{s^{\widehat{k}}|d^k,\mu }& =\sum _{j\in J}{p}_j{m}_{s^j|d^k,\mu }\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill {\Sigma }_{s^{\widehat{k}}|d^k,\mu }& =\sum _{j\in J}{p}_j{\Sigma }_{s^j|d^k,\mu }\hfill \end{array}$$

(33)

The torch position estimation is sufficiently accurate to ensure that the number of active PPCA models is not too large. The assimilation is still performed in real-time due to the very quick evaluation of each individual PPCA.

5.3. Forecasting of Future States

Despite the lack of Markovian assumption in the model, it can be used to forecast future states without observed experimental data by integrating previously estimated unknown parameter values. Indeed, the unknown parameters are not changed online and, assuming that they were estimated for a position $P^k$, ${\theta }^k$ could be used to obtain accurate predictions of future states.

The position of the future state is defined by a displacement of $j\Delta P$ in the direction of the weld, where $\Delta P$ is the regular separator in the positions grid. This opens two possibilities whether the torch position was known or was estimated from the data frame. If the heat source position was known, the prediction is computed using the PPCA model for position $P^{k+j}=(k+j)\Delta P$, this is $p({\tilde{s}}^{k+j}|\mu ,{\theta }^k)$.

The other possibility is that the current position was estimated from the experimental data. In this case, the displacement will be added to the position samples so that the uncertainty on the torch position is maintained. Then, the prediction will be calculated using the mixture of PPCAs.

6. Results

In this section, we will show the use of the proposed models. First of all, the experimental configuration of a PVR hot cracking test is explained. Then, details about the parametric finite element study from which the snapshots are generated will be given.

For the numerical tests, two sources of data are considered to show the case where the torch position is well known and the case where and estimation of its position is needed. The first test will use noisy synthetic data obtained from a previously calibrated Finite Element solution as input. The joint multiphysics state and unknown parameters estimation are obtained using a single PPCA model. Then, we will consider a frame of an infrared video of a PVR experiment. The torch position will be deduced from the image and the estimation computed from a mixture of PPCAs.

Last, the forecasting capabilities of the model will be showcased for the prediction of a state situated 12 mm further than the last studied position.

6.1. Experimental Procedure

A PVR hot cracking test was done at EDF R&D welding lab to obtain experimental data. The 316L steel specimen is 200 mm long and 3.5 mm thick. Its width is 80 mm on both ends and 40 mm in the center. The fusion line is 130 mm long. The specimen was placed in a tensile testing machine which was configured to augment the tensile deformation speed progressively from 0 mm/s to a maximum speed of 0.333 mm/s. The welding parameters used in the experiment are shown in

Table 1. PVR experiment parameter values.

Current	Voltage	Travel Speed	Shielding Gas	Maximal Stroke Rate
81 A	8.4 V	2.00 mm/s	Argon	20 mm/min

. Thus, the known parameters are v = 2 mm/s and $Q_0=U\times I=8.4\times 81=680.4$ W. The welded specimen is shown in Figure 3.

The experiment is instrumented with two type K thermocouples and a SC7500 FLIR infrared camera. The thermocouples, one on each side of the specimen (see Figure 1), are placed 110 mm from the bottom and 4 mm to the left from the center. The camera films the surface that is not being welded in order to avoid reflections from the welding arc. Figure 4 shows the projection of a frame of the video on the FE mesh. The resolution of the camera images is 320 × 256 pixels.

Previous uses of this configuration of thermocouples and infrared camera helped with the calibration of the parameter ${\sigma }_m^2$. To estimate the value of ${\sigma }_m^2$, we compare the measures of both sensors on a single point over time. The aim is to obtain an estimation of deviations between measures. Figure 5 shows the comparison between the signals. Here, both measures seem very close but show differences in some areas. Due to the configuration of the camera, only temperatures above 400 °C should be considered. One more thing to note is the peak that appears in the camera measure during the cooling phase, which was probably caused by the reflection of light. Work is being done at the lab to improve the quality of instrumentation and avoid this kind of issue. Considering the thermocouple and camera measures as ${\Theta }_T\left(t\right)$ and ${\Theta }_C\left(t\right)$, respectively, the value of ${\sigma }_m^2$ is estimated as

$${\sigma }_m^2=Var({\Theta }_T\left(t\right)-{\Theta }_C\left(t\right))=323.880246$$

(34)

which corresponds to a standard deviation of ~18 °C.

The thermocouple measures, and not the camera data, are also used to calibrate the unknown parameter values for a high fidelity Finite Elements simulation of the experiment that will be used as reference. This is a standard procedure that is done for every simulation of an experiment at EDF and it involves the resolution of an optimization problem. The result of the calibration of $\theta$ is shown in

Table 2. Deterministic calibration of the unknown parameters.

	${\mathit{a}}_{\mathit{f}}$	b	c	$\mathit{\eta }$	K
Calibrated value	6.657 mm	3 mm	1.5 mm	0.9	1.15

.

6.2. Prior Generation

In Section 4, we briefly discussed the generation of the prior that is obtained by running finite element simulations for a Latin Hypercube Sampling of the parameter space $\mathcal{P}$, where every parameter is assumed to be uniformly distributed. The minimum and maximum values for each parameter are issued from EDF’s decades of experience in both welding and numerical analysis of welding. These values are shown in

Table 3. Minimum and maximum values for the known parameters $\mu$, as determined by the experts.

	Heat Source Speed	Voltage	Current
Min. value	1 mm/s	8 V	70 A
Max. value	3 mm/s	12 V	90 A

for the known parameters and in

Table 4. Minimum and maximum values for the known parameters $\theta$, as determined by the experts.

	${\mathit{a}}_{\mathit{f}}$	b	c	$\mathit{\eta }$	K
Min. value	3 mm	1 mm	0.5 mm	0.75	1.15
Max. value	12 mm	7 mm	2.5 mm	0.95	2.5

for the unknown parameters.

A total of 128 simulations of PVR experiments were run with code_aster [37], EDF’s open-source thermomechanical simulation software. The selected physical fields for the multiphysics state vector are the temperature (T) and the maximum principal stress ${\sigma }_I$. The maximum principal stress is postprocessed from the stress tensor. It was chosen because it is used in a hot crack criterion developed in a previous PhD [38].

6.3. Numerical Tests

6.3.1. Estimation of the Heat Source Position

Before launching the two tests, the heat source position is identified using the surrogate model presented in Section 5. Figure 6 shows the discrete probability of assignment from 1000 samples of the position distribution.

The most probable basis is 130, which corresponds to 65 mm from the starting position. This position is the one that will be used as known position in the synthetic data test.

6.3.2. Tests with Noisy Synthetic Data

In this first test, the input data come from the calibrated finite elements simulation of the experiment, with $\mu =(2.0,680.4)$. The snapshot used as input is the one corresponding to the highest probability. White noise of the same amplitude as the measure error estimated for the camera has been added to the data. The observation function $H^k$ restricts the view to a region “seen” by the camera, and it is represented in Figure 7.

The mean of the posterior distribution $m_{s_k|d^k}^k$ contains the estimation of the temperature and maximum principal stress fields, as well as the unknown parameters. We can compare these results to the calibrated finite elements simulation, from which the input data was taken.

The results are compared along a line of 130 mm at the center of both sides of the specimen, which coincides with the the weld line on the torch side. The 95% confidence interval of the estimation is also plotted around the posterior mean in Figure 8 and Figure 9. The reconstruction of the temperature and stress fields is very accurate, with a global relative error of 1.3693% for the temperature and 6.3419% for the principal stress. This relative error is calculated as

$$\begin{array}{cc}\hfill e_T^k& =\frac{||T_{|\mu ,\theta }^k-s_{|T}^k{||}_2}{||T_{|\mu ,\theta }^k{||}_2}\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill e_{{\sigma }_I}^k& =\frac{||{\sigma }_{I|\mu ,\theta }^k-s_{|{\sigma }_I}^k{||}_2}{||{\sigma }_{I|\mu ,\theta }^k{||}_2}\hfill \end{array}$$

(36)

where $T_{|\mu ,\theta }^k$ and ${\sigma }_{I|\mu ,\theta }^k$ are the nodal temperature and principal stress simulation values, respectively, and $s_{|T}^k$ and $s_{|{\sigma }_I}^k$ are the estimated nodal temperature and principal stress, respectively.

The confidence intervals are thin for the observed temperature data, as it is expected, but it is larger around the peak on the non-observed surface. For the principal stress estimations, although no data was seen, the estimated posterior mean is close to the simulation on both surfaces.

The estimation of the unknown parameters $\theta$ is very close to the results of the deterministic calibration, with a relative error inferior to $10\%$, meaning that the parameter calibration successfully identified the theta values used for the simulation. A comparison of the values obtained by the deterministic calibration and the mean posterior theta estimation is given in

Table 5. Posterior estimation of the unknown parameters $\theta$ obtained from noisy synthetic data.

	${\mathit{a}}_{\mathit{f}}$	b	c	$\mathit{\eta }$	K
Calibrated value	6.657	3	1.5	0.9	1.15
Estimated value	6.18868114	2.84188617	1.58588521	0.89049624	1.16543712

.

6.3.3. Tests with Real Experimental Data

In this test, we will use the camera frame projected onto the mesh as input data for the model. The camera is configured to capture temperatures above 400 °C, so the observation function $H^k$ is modified to only use data above 400 °C, which is the area represented in Figure 10.

The estimation is calculated as a mixture of the results given by the local PCA models in Figure 6. The amount of active local models is small and each individual conditioning is computed very fast. All the confidence intervals are the 95% confidence intervals.

On the camera side, the estimations can be compared to the measured temperature. In Figure 11, we can see that the model estimates a temperature that follows the experimentally measured one. Interestingly, the estimation deviates from the simulation, as seen in Figure 12, where the temperature is higher for the FE results. We interpret this difference in the estimation and the FE simulation as a model correction given by the partial observation of the temperature. We remind the reader that the FE model was calibrated using only the thermocouple measures and not the infrared camera. This difference between the FE simulation and the measured data in Figure 12 may indicate that the calibration of $\theta$ using the thermocouples is not good enough. This is supported by the fact that the estimated efficiency $\eta$, shown in

Table 6. Posterior estimation of the unknown parameters $\theta$ obtained from infrared experimental data.

	${\mathit{a}}_{\mathit{f}}$	b	c	$\mathit{\eta }$	K
Calibrated value	6.657	3	1.5	0.9	1.15
Estimated value	6.40815805	3.8549261	1.47004239	0.73615823	1.09763587

along with the other unknown parameters, is smaller at ~0.74 while the initial calibration estimated it at 0.9. A smaller efficiency would transfer less temperature to the specimen, explaining the lower estimated temperature.

Note that on the camera side the confidence interval is very small for the temperature estimation, as it is expected from observed data, while on the torch side the confidence interval is larger, especially around the peak position. No mechanical data is observed, and thus the maximum principal stress confidence intervals reflect the uncertainty of the posterior estimation (see Figure 13).

6.4. Forecasting

In this last example, the forecasting capabilities of the model will be shown. Let us assume that we want to estimate the temperature and maximum principal stress fields in a future position that has not been observed yet. This position is situated 12 mm further than the one studied for the previous tests. All the information available is the value of the known parameters and the posterior estimation of the unknown parameters with their associated posterior covariance.

Choosing the active coefficient of the mixture of PPCAs for this position $P^j=(k+12)\Delta P$ is the first problem to solve. When an image is available, the surrogate model for position estimation takes the highest temperature identified in the image and returns the probability distribution of the torch position. In the forecasting problem, we have no video frame to find the highest temperature, so it is assumed to be shifted by the 12 mm and then it is used to obtain the discrete assignment probability. The multiphysics state ${\tilde{s}}^j$ will be predicted by conditioning its distribution by the known parameters $\mu$ and the estimated unknown parameters ${\theta }^k$ with its associated posterior variance.

As with previous tests, the predicted estimations, $p({\tilde{s}}^j|\mu ,{\theta }^k)$, are compared to the experimental data (when it is possible) and to a Finite Elements simulation over a line on both sides of the specimen. Figure 14 shows a frame of the video where the highest temperature position is situated 12 mm to the right of the highest temperature position found in Figure 10. The measures corresponding to that frame are compared to the predictions in Figure 15. The estimation is very close to the camera data but the 95% confidence interval is very large compared to the one in Figure 11, a case where the data was observed, reducing the uncertainty.

The Finite Elements simulation was used as reference in Figure 16 and Figure 17, where we observe generally large confidence intervals on both sides of the specimen. Overall, knowledge of the known and unknown parameters is enough to obtain predictions that represent the behaviour of the temperature and maximum principal stress fields but with a high degree of uncertainty, which can be greatly reduced by observation of the thermal images. This is not the case when only the known parameters are observed. Additional figures are shown in the Appendix B to support this claim.

7. Conclusions

We have proposed a novel data assimilation methodology for automatized welding operations monitored using thermal imaging. The approach is based on digital twinning, a physically detailed thermomechanical finite element model assimilating online data and predicting unseen mechanical quantities of interest such as stress fields. Sources of variability are clearly identified and modeled using appropriate random parameters, the distribution of which are designed based on expert knowledge. Data assimilation then consists in finding the posterior distribution of the uncertain parameters given the sequence of thermal images available at current process time.

The data assimilation framework is made in real-time by deploying and offline–online meta-modeling technology. Sparse linear models are postulated for every position of the welding torch, linking observations and hidden mechanical quantities to the parameters of the welding process and to the random parameters. The coefficients of the linear models are identified using the method of snapshots, using hundreds of prior high-fidelity computations. Online, predicting unseen mechanical states operation reduces to simple Gaussian conditioning with a background covariance matrix exhibiting a low-rank structure. This is made computationally efficient using standard algebraic manipulations.

Thanks to this model/data fusion technology, we have shown that we were, for the first time, able to predict mechanical stress fields during welding in real-time, the predictions being continuously adjusted based on thermal imaging. We have shown that our digital twin may produce predictions that differ significantly from those produced by a high-fidelity model that is precalibrated using standard thermal sensing via a set of thermocouples. We interpret this finding as an online model correction, the data acquired online being richer, corresponding to thermal sensing closer to the regions of intense thermal gradients.

Finally, we were able to forecast stress predictions into the future of welding operations, by simply using current posterior distributions of unknown parameters and perform uncertainty propagation. We have shown that such predictions where reasonably sharp and could be used to stop welding operations in advance, when forecasting inadmissible levels of stresses.

The developments presented in this paper are expected to be the foundations for further theoretical and applied research. Future experiments may include additional instrumentation to obtain more data such as measures of the strain field by digital image correlation. In terms of algorithmic efficiency, the proposed meta-modeling methodology is piecewise linear, homogeneously along the trajectory of the welding torch. We are now investigating a data-driven clustering approach of the welding parameter domain, which we expect will help us minimize the memory required to store the various meta-models. In terms of application, the present developments focus on the prediction of stresses during the operations, but could easily be extended to the prediction of residual stresses. The exploration of parameter spaces describing geometrical variations of joining operations is also of high practical interest. Finally, the proposed approach may constitute the computational core of a model-based control technology aimed at adjusting the process parameters online in order to ensure that the joining operations produce assemblies that are of acceptable mechanical qualities.