© S. Torregrosa et al., Published by EDP Sciences, 2024

1 Introduction

The analytical solution of engineering models is today compromised as they become more and more complex. Since the mid-20th century [1], more and more powerful computers have been developed and democratized leading to the emergence of the so called third paradigm of science: the “virtual twins” or physics-based numerical simulations. A “virtual twin” emulates the complex behavior of a physical system using mathematical models [2]. Hence, numerical simulation has become an essential tool for scientific investigation and analysis of complex systems in engineering, facilitating the access to high-fidelity (HF) data and, thus, drastically reducing the number of experimental tests [2–5].

However, the size and complexity of the problem studied still limits the numerical simulations’ capabilities. Indeed, the intricacy of the engineering system is proportional to the computational resources needed and its linked cost. Thus, achieving HF simulations remains nowadays expensive and time-consuming [6]. HF numerical simulations are, thus, not adapted to real-time constraints or to some industrial contexts where we expect HF data quickly and with less effort [5].

Among the most known “virtual twins” one can think about computational fluid dynamics (CFD), which is studied in this paper. Since in realistic turbulent flow cases the Navier Stokes equations cannot be solved analytically, CFD offers a numerical approach to compute such flow phenomena. Through the years, CFD has been widely studied and developed leading to a high-fidelity tool with its more accurate approach: Direct Numerical Simulations (DNS) [7]. However, as highlighted just before, such methods are extremely computationally expensive and still infeasible on numerous applications [8]. Simpler CFD approaches have been developed such as large eddy simulation (LES) [9] or Reynolds Averaged Navier Stokes equations method (RANS) [10], which is nowadays the industrial standard approach. Although all the work done to improve its accuracy and cost [11,12], CFD simulations remain today either not to be trusted completely or too computationally expensive to run [8].

Moreover, the present data-based revolution leads to the fourth paradigm: the “data-intensive science” [1]. A new framework is developed where data, simulation and theory can interact and reinforce each other. Indeed, thanks to widely developed machine learning and artificial intelligence techniques, data has massively proliferated in the majority of scientific fields and data-based models (also called “digital twins”) are receiving much attention in numerous and different applications. Hence, recent interest is posed on these new tools to either replace or enhance CFD [8]. Therefore, artificial intelligence is nowadays used in numerous forms for computational efficiency but also to increase accuracy. Among those different applications one can find:

• Near wall turbulent flows modeling [13].

• Turbulence model coefficients tuning [14,15].

• Turbulence model enhancement by building a representation of closure terms [16,17].

• Surrogate model of the numerical simulation trained with CFD data [18,19].

Likewise, borrowed from computer vision and image recognition, super-resolution aims at obtaining a high resolution image output from a low-resolution image. In this sense, it can be applied in fluid mechanics simulation leading to a machine learning model taking as input a fast and computationally inexpensive simulation and returning a high-fidelity one. As any data-driven model, such approach needs an offline training stage and data base. In the literature, one can find super-resolution models trained with DNS high-resolution data and the purposely down-scaled counterpart. However, this could be judged as unreasonable since turbulent flow structures often differ significantly from their fine counterparts when computed on coarser resolutions [5,20,21]. Physics-informed machine learning models can also be found. Such models usually do not need a HF training data base and guarantee that the super-resolved fields are faithful to the physical laws and principles. However, they need an optimization step in the online stage [22,23]. Finally, correction models of the error gap between coarse and high-fidelity simulations are also developed but only for scalar variables [24,25].

In this paper we focus on this last family of approaches: data driven models combining both the “digital” and “virtual” twins that correct coarse CFD simulation fields [2]. On the one hand we want the most correct solutions, but high-fidelity data is too expensive and time-consuming.

On the other hand, we have access to computational facilities to run “as many as needed” coarse numerical simulations, but these present a non negligible deviation with the actual solution (i.e. high-fidelity data) of the system. It seems thus coherent to use both technologies together, using the advantages of each approach. The “hybrid twin” (HT) is thus born, composed of an imperfect coarse simulation of the system (the “virtual twin”) and of a data-based model emulating the ignorance gap between coarse and high- fidelity simulation fields (the “digital twin”).

Indeed, since high-fidelity simulations are the most correct ones but also too expensive, it becomes cheaper to learn the difference between some HF and the corresponding coarse solutions, allowing to correct further coarse simulations, instead of computing the corresponding HF ones. Hence, the computational complexity of the here proposed HT methodology is related to the computation of the high-fidelity simulations. Indeed, this data generation is expensive, with the complexity of usual CFD simulations. However, the data-based enrichment model construction (performed in an offline stage) that we propose is orders of magnitude cheaper. For the solutions presented here is of the order of magnitude of the hour. On the other hand, the online application of the correction data-based model is almost instantaneous, the order of magnitude of the second, which represents the key benefit of the Hybrid Twin rationale.

The question that arises now is, how do we model the ignorance, i.e. the coarse-high-fidelity simulations gap? One could think about measuring this gap as a difference. However, since we are working with fields, this approach could lead to non-physical interpolated results in, for instance, fields where a given choice of the problem parameters leads to a localized solution in different regions (e.g. fluid mechanics). A smarter solution would be to use Optimal Transport (OT) theory [26], which provides a mathematical framework to calculate distances between general objects and can be considered more physical in many fields. Such a solution has already been explored by the authors in previous works [27,28].

Indeed, the Optimal Transport theory was applied by the authors in [28] to build a “digital twin”: the ignorance gap is learnt as a transport of information. In this previous work, the HT approach was applied over the error gap between experimental and numerical (the “virtual twin”) data, leading to a data-based model (the “digital twin”) able to correct a CFD simulation based on experimental knowledge. Indeed, physics-based simulations present some significant deviations when compared to measurement data. This deviations are expected to be biased since they represent the ignorance of the modeler on the subjacent physics, related to the inaccuracy in the employed models. Hence, the artificial intelligence, acting as a black box, models the physical part that is beyond the modeler’s knowledge [2,4,29].

In this new paper, the same Optimal Transport based methodology is applied, but now over the error gap between coarse and high-fidelity numerical data. Hence, the coarse data is OT-based corrected by being optimally transported to the high-fidelity data. Here, the artificial intelligence models the ignorance of the modeler over the spatial discretization of the industrial system: the meshing step. Indeed, a high-fidelity mesh leads to accurate but costly numerical results while a coarse mesh leads to poorer but cheaper ones. It is important to note that, in contrast with the numerical-experimental ignorance gap, now it is a known ignorance since the difference between the HF and coarse mesh is part of the modeler’s knowledge.

The already published HT method is a two stages approach: first, the OT-based “digital twin” is trained offline, then, this data driven correction is applied over the “virtual twin” output in an online manner. Indeed, coarse simulations and their high-fidelity counterparts are used to train the OT-based correction model. Once trained, the “digital twin” can be used in an online manner to correct further coarse simulations from the “virtual twin”.

The steps of the offline stage are as follows. First, coarse simulations (the “virtual twin”) and their high-fidelity counterparts are computed in the parametric space of our problem. Then, all the training simulations are decomposed into the sum of identical Gaussian functions (also called particles) based on a Smoothed-Particle Hydrodynamics (SPH) decomposition [30]. Next, for all the coarse- HF couples of the training set, each particle from the coarse data is matched with a particle from its corresponding high-fidelity data. Hence, the OT-based differences between coarse and HF data are computed. Finally, the OT-based “digital twin” is built by training a Neural Network architecture over the OT-based gaps.

Once the correction data-based model of the ignorance is trained, it can correct a new coarse simulation, from which we do not have access to the high-fidelity data, in an online manner. To this purpose, the “virtual twin” output is decomposed into particles and those particles transported based on the OT-based gap, interpolated by the “digital twin”. Finally, the so expected high-fidelity corresponding data is reconstructed by summing all the particles, i.e. Gaussian functions.

It can be noted that if the CFD solution of the studied problem is not regular with respect to the problem parameters’ change, the Kolmogorov n-width increases, needing for nonlinear dimensionality reduction, as Optimal Transport here proposed performs [31–33]. Additionally, if the relation between the parameters and the corresponding solution is regular, the problem becomes simpler, and the OT-based technique proposed also performs accurately.

In this article, the principal ideas of OT theory and the main steps of the methodology are shortly presented. Indeed, we focus here in the new results more than in the methodology which has already been detailed in [28]. It is important to be noted that even if the methodology is the same, the problem solved is conceptually completely different, moving from the simulation-experiment gap to the coarse-HF simulation gap.

2 Revisiting optimal transport

In this section, the OT framework is presented and the tools on which the thereafter proposed OT-based correction model is based introduced. It should be noted that this section is a non exhaustive introduction of the principal ideas of OT theory. For further documentation on this topic, [34] and the references therein can be consulted.

The initial Optimal Transport problem was introduced by Monge [35]. It consisted in finding the most optimal path to move a given quantity of soil from an initial to a target location. The cost function in the transport was defined as the distance traveled by the soil. Note that in this article we are only interested in the discrete formulation of the problem. In order to introduce the discrete problem, let us consider M factories consuming a certain resource that needs to be transported from N mines. The cost function to minimize is the square of the total Euclidean distance traveled by the resource. This discrete Optimal Transport problem is illustrated in Figure 1.

On the one hand, each factory m ∈ ⟦M⟧ (note that the notation ⟦M⟧ corresponds to {1,…,m,…,M}) is located at y_m and consumes an amount b_m of the resource. On the other hand, each mine n ∈ ⟦N⟧ is located at x_n and produces an amount a_n of this same resource. Hence, following the notion of measure, two distributions, α and β, corresponding to the resource produced and consumed respectively can be defined:

$$\alpha ={\displaystyle \sum _{n=1}^N{a}_n}{\delta }_{x_n}\text{ and }\beta ={\displaystyle \sum _{m=1}^M{b}_m}{\delta }_{y_m}$$(1)

where ${\delta }_{x_n}\text{ and }{\delta }_{y_m}$ correspond to the Dirac at locations x_n and y_m respectively.

Solving the discrete Monge problem consists, thus, in finding the map T connecting each point x_n with a single target point y_m such that the transport cost is minimized. Here, we consider the square of the L2 distance between the mine n and its corresponding factory m:

$$C_{x_n,y_m}={\Vert x_n-y_m\Vert }_2^2.$$(2)

Hence, the produced resource distribution, α, is pushed toward the consumed resource distribution, β. Since the resource cannot be destroyed or produced during its transport, transport map T:{x₁,…,x_N}→{y₁,…,y_M} must also satisfy the mass conservation:

$$\forall m\in ⟦M⟧,b_m={\displaystyle \sum _{n:T\left(x_n\right)=y_m}{a}_n},$$(3)

Note that here the map T is a surjective function. Finally, we obtain the following minimization problem:

$$\underset{\text{T}}{{\displaystyle \min }}{\displaystyle \sum _{n=1}^N{C}_{x_n,T\left(x_n\right)}}.$$(4)

In order to developed our data-driven correction model, we are interested in a simplified version of the just presented discrete Monge problem. Indeed, it is first supposed that N= M, i.e. that the number of factories and mines is the same. In addition, it is also supposed that the quantity of resource produced by every mine and consumed by every factory is the same, i.e. a_n= b_m = 1/N. Therefore, the optimization problem (4) is now a deterministic matching problem and the transport map becomes a bijective function.

Under these assumptions, the just simplified discrete Monge problem can be easily solved by linear programming. Indeed, the problem is now equivalent to an optimal matching problem between two particle clouds, as it is illustrated in Figure 2 in 2D. Note that each cloud is composed by the same number of particles, every particle has the same amount of mass and the cost is defined as the square of the L2 distance between two particles. In higher dimensions than 2D, the computational cost of resolution increases but the problem does not further complexifies.

The OT-based difference between two optimally paired particle clouds can be determined by calculating the Euclidean distances between matched particles. Indeed, the OT-based difference is the total sum of the L2 distances δ_k for all the pairs of particles.

Fig. 1 Discrete OT formulation for N= 4 mines and M=3 factories. The resource produced by the mines is: a1 =3, a2 =1, a3 = 1 and a4 = 2. The resource consumed by the factories is: b1 =4, b2 = 1 and b3 = 2. The Euclidean distance traveled by the resource is the cost to minimize.

3 Hybrid twin based on optimal transport

The “hybrid twin” approach, based on optimal transport, developed by the authors and published in [28] is here shortly reviewed. We strongly recommend the interested reader to refer to the previous article for all the details of its implementation. First, the offline stage where the “digital twin” is trained is introduced. Then, its online application, in combination with the “virtual twin”, is presented.

Let suppose a parametric space $\mathcal{W}\left({\eta }_1,\dots ,{\eta }_q,\dots ,{\eta }_Q\right)$ where η_q, q ∈ ⟦Q⟧ are the parameters and where our parametric problem is defined. Next, let consider P 3D coarse-HF simulation couples in W corresponding to the 3D simulations of the parametric problem using a coarse and a high-fidelity meshing respectively. The 3D fields are monitored in a defined 2D domain of interest Ω. Hence, each coarse and HF data sample is formally represented by a distribution ψ: Ω Є ℝ² → ℝ⁺. It can be noted that the image of ψ is supposed strictly positive.

First, the OT-based ignorance model is trained offline, based on the coarse-HF simulation couples, following the next steps (colored in blue in Fig. 3):

Pre-processing: Normalization of the distributions corresponding to the coarse simulations and the high- fidelity counterparts to obtain unitary integral distributions:

$$\rho =\frac{\psi }{ℐ}\text{ where }ℐ={\displaystyle \underset{\Omega }{\int }\psi }\text{d}\Omega \text{. }$$(5)

Particles decomposition: Every coarse and HF simulation is decompounded into a sum of N identical 2D-Gaussian functions of fixed standard deviation σ and mass 1/N. Indeed, it is important to note that the number of particles N and the standard deviation σ of each particle are hyperparameters of our methodology.

$$\begin{array}{l}\overline{\rho }(x)={\displaystyle \sum _{n=1}^N{G}_{{\mu }_n,\sigma }}(x)\text{ where }{G}_{{\mu }_n,\sigma }(x)\\ =\frac{1}{N{\sigma }^{2}2\pi }\exp \frac{-{\left(x-{\mu }_n\right)}^2}{2{\sigma }^2}.\end{array}$$(6)

Hence, the only variables are, for a given distribution, the means μ_n of each Gaussian function, i.e. N vectors of 2 components: μ_n,x and μ_n,y (because we are in 2 dimensions).

Therefore, we need to solve P × 2 minimization problems (i.e. two optimizations for each coarse-HF simulation couple) in order to place the N particles minimizing the error with respect to the original data. To this purpose a Gradient Descent approach is used to solve the optimization problem (here the subscript c denotes the coarse simulation):

$$\begin{array}{l}\underset{{\mu }_c^p}{\min }\frac{1}{2}{\Vert {\rho }_c^p-{\overline{\rho }}_c^p\Vert }_2^2\\ ={\min }_{{\mu }_c^p}\frac{1}{2}\left[{\displaystyle \sum _{i=1}^D{\left({\rho }_c^p\left(x_i\right)-{\displaystyle \sum _{n=1}^N{G}_{{\mu }_{c_n}^p,\sigma }}\left(x_i\right)\right)}^2}\right],\end{array}$$(7)

where D is the number of points of the mesh where the distribution ${\rho }_v^p$ is calculated.

Once the decomposition is computed, one can introduce the matrix ${\mu }_c^p\in {ℝ}^{N\times 2}$, composed by the coordinates x and y of all the particles of the coarse simulation of the pth couple (and likewise for the high-fidelity simulation with the subscript hf):

$${\mu }_c^p=\left[\begin{array}{c}{\mu }_{c_1}^p\\ ⋮\\ {\mu }_{c_n}^p\\ ⋮\\ {\mu }_{c_N}^p\end{array}\right]=\left[\begin{array}{c}\left[{\mu }_{c_{1,x}}^p,{\mu }_{c_{1,y}}^p\right]\\ ⋮\\ \left[{\mu }_{c_{n,x}}^p,{\mu }_{c_{n,y}}^p\right]\\ ⋮\\ \left[{\mu }_{c_{N,x}}^p,{\mu }_{c_{N,y}}^p\right]\end{array}\right]\in {ℝ}^{N\times 2}.$$(8)

It is important to note that the order of the particles in this matrix ${\mu }_c^p$ is not arbitrary but will be used to represent the matching between point clouds of a given couple: the nth particle of one cloud being matched with the nth particle of the other cloud from the same couple.

P 2-dimensional matchings: The Optimal Transport behavior is calculated: for each coarse-HF simulation couple, each particle from the coarse simulation is matched with one particle from the high-fidelity counterpart. This is understood as the optimal matching problem between two N-particles clouds, where each particle is a 2D-Gaussian function represented by its μ_x and μ_y coordinates. This linear assignment problem, where the cost function for the p^th couple, $C_{c,hf}^p$, is the sum of the squared L2 distances between matched particles, can be resolved applying several algorithms. Here, the algorithm matchpairs from MATLAB is used to solve the problem [36].

$$C_{c,hf}^p\left({\varphi }_p\right)={\displaystyle \sum _{n=1}^N{\Vert {\mu }_{c_{{\varphi }_p(n)}}^p-{\mu }_{h{f}_n}^p\Vert }_2^2}$$(9)

where ϕ_p is a bijective function in the set of permutations of N elements. To each particle n of the distribution ${\rho }_c^p,{\varphi }_p:\mathbb{N}\to \mathbb{N}$ associates its new position in the sense of order in ${\mu }_c^p$.

OT-based gap: The OT-based gaps δ^p are calculated for all the P coarse-HF simulation couples by computing the difference between the coordinates of the matched N particles. Hence, the OT-based gap of the nth particle of the pth couple is noted:

$$\begin{array}{l}{\delta }_n^p=\left[{\delta }_{n,x}^p,{\delta }_{n,y}^p\right]\\ =\left[{\mu }_{c_{n,x}}^p-{\mu }_{h{f}_{n,x}}^p,{\mu }_{c_{n,y}^p}^p-{\mu }_{h{f}_{n,y}}^p\right],\end{array}$$(10)

The OT-based difference for the pth couple δ^p ∈ ℝ^N×2 is build as:

$${\delta }^p=\left[\begin{array}{c}{\delta }_1^p\\ ⋮\\ {\delta }_N^p\end{array}\right]\in {ℝ}^{N\times 2}.$$(11)

“Digital twin” training: The “digital twin” follows a Neural Network architecture. Indeed, two NN are trained, one for the x coordinates of the OT-based gap and another for the y coordinates. It can be noted that the parametric input space of the NN is built as follows. To each set of parameters p of the training set, we add the coordinates of the N particles of the pth “virtual twin” decomposition: ${\mu }_{c_x}^p\in {ℝ}^N\text{ and }{\mu }_{c_y}^p\in {ℝ}^N$ . The new parametric space is noted ${\mathcal{W}}^{\prime }\left({\eta }_1,\dots ,{\eta }_q,\dots ,{\eta }_Q,{\mu }_{c_x},{\mu }_{c_u}\right)$ . This yields a (P × N) × (Q +2) matrix X of explanatory variables. Moreover, the (P × N) × 2 matrix Y of response variables is the concatenation of the P OT-based differences δ^p, p∈ ⟦P⟧.

Then, the “digital twin” can be used in a partially online manner to correct a “virtual twin” coarse simulation, from which we do not have the high-fidelity counterpart. The online stage follows the next steps (colored in red in Fig. 3):

First, a new coarse simulation, computed in the parametric space W, is normalized and decomposed into N particles.

$$\begin{array}{l}{\mu }_c=\left[{\mu }_{c_x},{\mu }_{c_y}\right]\in {ℝ}^{N\times 2}\text{ such that }{\rho }_c\\ ={\displaystyle \sum _{n=1}^N{G}_{{\mu }_{c_n},\sigma }}(x).\end{array}$$(12)

Then, the “digital twin” returns the OT-based gap that needs to be added to these particles to obtain the particles of the high-fidelity counterpart. Indeed, the two NN applied over these new N inputs in W′ return the x and y components of the OT-based correction for all the particles respectively: δ ∈ ℝ^N×2. Adding this difference to the original particles leads to the the corrected positions of the particles (subscript cor):

$${\mu }_{cor}={\mu }_c+\delta .$$(13)

Next, the expected high-fidelity data ${\widehat{\rho }}_{hf}$ is reconstructed by summing these N corrected Gaussian functions:

$${\widehat{\rho }}_{hf}={\displaystyle \sum _{n=1}^N{G}_{{\mu }_{{\text{cor }}_n},\sigma }}(x).$$(14)

Finally, in order to recover ${\widehat{\psi }}_{hf}$, the correction gap of the total masses of the coarse-HF simulation couples, ${\mathcal{J}}_c-{\mathcal{J}}_{hf}$, is also interpolated by a NN architecture and, thus, the normalization step undone.

Fig. 2 2D Monge problem with N=M and an= bm = 1/N. Mines and factories are represented by red and blue circles respectively. Black arrows illustrate the optimal matching.

Fig. 3 Summary diagram of the methodology: the offline training of the “digital twin” is colored in blue and the online OT-based “hybrid twin” approach in red.

4 Results

4.1 Error evaluation

First of all, the error evaluation methodology is defined. To this purpose, let us introduce a testing data set in the parametric space W. The P_test reference solutions, i.e. the high-fidelity data, ${\psi }_{hf}^p,p\in ⟦P_{\text{test }}⟧$ of this set are compared with the corrected coarse numerical solutions ${\widehat{\psi }}_{hf}^p,p\in ⟦P_{\text{test }}⟧$ and with the “virtual twin” solutions ${\psi }_c^p,p\in ⟦P_{\text{test }}⟧$. Three error metrics are here introduced: a maximum value amplitude error, a maximum value position error and a L² – Wasserstein error. In this section, and without loss of generality, the error metrics are presented for the OT-based corrected data case.

First, the maximum value amplitude error is computed as the relative difference between the maximum value amplitude of the high-fidelity solution and the maximum value amplitude of the OT-based corrected data. Hence, the maximum value amplitude error for the pth training point writes:

$${\epsilon }_{\max }^p=100\frac{\left|\max \left({\psi }_{hf}^p\right)-\max \left({\widehat{\psi }}_{hf}^p\right)\right|}{\max \left({\psi }_{hf}^p\right)}.$$(15)

Next, the maximum value position error is calculated as the L2 norm between the positions in Ω of the maximum value of the high-fidelity simulation and of the maximum value of the “hybrid twin” solution. Note that this error is normalized with respect to l_Ω, the length of one side of Ω:

$${\epsilon }_{pos}^p=\frac{{\Vert {\mathrm{argmax}}_x\left({\psi }_{hf}^p(x)\right)-{\mathrm{argmax}}_x\left({\widehat{\psi }}_{hf}^p()\right)\Vert }_2}{l_{\Omega }}.$$(16)

Finally, the L²-Wasserstein metric $W{\left({\psi }_{hf}^p,{\widehat{\psi }}_{hf}^p\right)}_2^2$ is calculated between the reference ${\psi }_{hf}^p$ and the corrected ${\widehat{\psi }}_{hf}^p$ solutions. In order to calculate $W\left({\psi }_{hf}^p,{\widehat{\psi }}_{hf}^p\right)$, a Linear Programming methodology is followed.

The Kantorovich OT problem [37] corresponds to an infinite dimensional Linear Program. Indeed, given the distributions ${\psi }_{hf}^p$ and ${\widehat{\psi }}_{hf}^p$,defined on X and Y, the problem reads

$$W{\left({\psi }_{hf}^p,{\widehat{\psi }}_{hf}^p\right)}_2^2=\underset{{}_{\pi \in \Pi \left({\psi }_{hf}^p{\widehat{\psi }}_{hf}^p\right)}}{{\displaystyle \min }}{\displaystyle {\int }_{X\times Y}c}(x,y)\text{d}\pi (x,y),$$(17)

where c (x, y): X × Y→ ℝ is the cost function and П the set of transfer plans. The discretized measures ${\psi }_{hf}^p$ and ${\widehat{\psi }}_{hf}^p$ are defined as weighted sums of Dirac functions. The weights represent the value of the continuous measures evaluated at the corresponding nodes x_i and y_i of the mesh. Hence,

$${\psi }_{hf}^p={\displaystyle \sum _{i=1}^D{\left({\psi }_{hf}^p\right)}_i}{\delta }_{x_i}\text{ and }{\psi }_{hf}^p={\displaystyle \sum _{i=1}^D{\left({\psi }_{hf}^p\right)}_i}{\delta }_{y_i}.$$(18)

The discrete cost function is defined as

$$C_{i,j}=c\left(x_i,y_j\right)={\Vert x_i-y_j\Vert }_2^2.$$(19)

Therefore, the discrete formulation of the Kantorovich Optimal Transport problem reads

$$W{\left({\psi }_{hf}^p,{\widehat{\psi }}_{hf}^p\right)}_2^2=\underset{\pi \in \text{Π}\left({\psi }_{hf}^p{\widehat{\psi }}_{hf}^p\right)}{\min }\underset{i,j}{{{\displaystyle \sum }}^{\text{}}}{C}_{i,j}{\pi }_{i,j},$$(20)

where π_i,j represents the quantity of mass transported from π_i towards y_j. The set of transfer plans reads

$$\begin{array}{l}\Pi \left({\psi }_{hf}^p,{\psi }_{hf}^p\right)=\{\pi =\left({\pi }_{ij}\right)\mid {\displaystyle \sum _j{\pi }_{ij}}={\left({\psi }_{hf}^p\right)}_i,\\ {\displaystyle \sum _i{\pi }_{ij}}={\left({\widehat{\psi }}_{hf}^p\right)}_j\}\end{array}$$(21)

It should be noted that, in order to analyse the three different errors, the value of the error metric presented corresponds to the mean value of the p ∈ ⟦P_test⟧ points of the test set. Hence, the three error mean metrics are noted:

$${\epsilon }_{\max },{\epsilon }_{\text{pos }}\text{ and }{W}_2^2.$$(22)

4.2 Fluid dynamics problem

In this section, the methodology developed is applied to a fluid dynamics problem. A 3D steady turbulent flow into a channel facing a backward ramp is studied, as it is illustrated in Figure 4. As indicated in (23), the fluid, considered as incompressible, has a uniform velocity profile at the inlet domain Ω_Inlet with a inlet velocity of v_Inlet.

The geometry used in this paper is closed to the Ahmed’s study geometry. The slant angle of our geometry corresponds nearly to the minimum of drag found in the Ahmed’s study geometry [38]. A non slip condition is imposed on the walls Ω_Wall and a zero gradient condition is imposed on the outlet section Ω_Outlet. Therefore, the problem writes:

$$\{\begin{array}{ll}v\cdot \nabla v=-\frac{1}{\rho }\nabla p+v{\nabla }^{2}v\hfill & \text{ in }{\Omega }_{\text{Channel }},\hfill \\ \nabla \cdot v=0\hfill & \text{ in }{\Omega }_{\text{Channel }},\hfill \\ v(x=0,y,z)=v_{\text{Inlet }}x\hfill & \text{ in }{\Omega }_{\text{Inlet }},\hfill \\ v=0\hfill & \text{ on }{\Omega }_{\text{Wall }},\hfill \\ \nabla v\cdot n=0\hfill & \text{ on }{\Omega }_{\text{Outlet }},\hfill \end{array}$$(23)

where ρ is the density, v the kinematic viscosity, n is the outward normal from Ω_Outlet and x the elementary vector of the x axis. The turbulence model chosen is k − ω SST with a stepwise switch wall function:

$$\{\begin{array}{l}\omega ={\omega }_{vis}=\frac{6v_w}{{\beta }_1{y}^2}\text{\hspace{1em} if }{y}^{+}\le y_{lam}^{+}\hfill \\ \omega ={\omega }_{log}=\frac{\sqrt{k}}{C_{\mu }\kappa y}\text{\hspace{1em} if }{y}^{+}>y_{lam}^{+},\hfill \end{array}$$(24)

where ω is the specific dissipation rate, k the turbulent kinetic energy, y the wall-normal distance, C_μ and β₁ model constants, v_w the kinematic viscosity of fluid near wall, κthe von Karman constant, y⁺ the estimated wall-normal distance of the cell center in wall units and $y_{lam}^{+}$ the estimated intersection of the viscous and inertial sub-layers in wall units.

The geometry is parameterized as it can be seen in Figure 5. The numerical values chosen for each parameter are gathered in Table 1. Hence, the parameters defining the parametric space are v_Inlet, v and α: $\mathcal{W}\left(v_{\text{Inlet }},v,\alpha \right)\in {ℝ}^3$. It can be noted that by determining α, and since h₃ is fixed, h₁ and h₂ are thus also fixed by:

$$h_2=\frac{l_2}{\tan (90-\alpha )}\text{ and }{h}_1=H-h_3-h_2.$$(25)

Likewise, by fixing L, l₁ and l₂, l₃ is thus also fixed by:

$$l_3=L-l_2-l_1.$$(26)

In order to solve the problem, the channel is meshed with an hexahedral mesh. On the one hand, we mesh the geometry with a high-fidelity mesh, as it can be seen in Figure 6, for the high-fidelity data (i.e. the reference simulations). On the other hand, we mesh the geometry with a coarse mesh, as it can be seen in Figure 7, for the coarse data (i.e. the “virtual twin” simulations). It can be noted that in Figure 5, the number of the nodes per edge for the coarse and HF meshes are indicated in red and blue respectively.

The Computational Fluid Dynamics OpenFOAM code is used to solve both finite volumes problems. The SIMPLE solver is selected to solve the Navier-Stokes equations. The convergence of the simulations are assured by controlling the residual convergence. Moreover, the norm of the velocity field is monitored on a plane of interest PoI perpendicular to the channel at x = l_PoI, as it is represented in Figure 4. In order to compare the high-fidelity simulation with the coarse one, the norm of two velocity fields are presented. First, a longitudinal cut of the channel is shown in Figure 8 for both the HF mesh (top) and the coarse one (bottom). It can be noticed that the flow behavior (recirculation bubble) of the Ahmed body is observed [38]. Indeed, here, for a low angle a, the flow remains attached over the ramp. Then, a perpendicular cut of the channel at the plane of interest PoI is shown in Figure 9 for both the HF mesh (left) and the coarse one (right). It can be noted that the poor mesh leads to different physical results with respect to the high-fidelity mesh.

Here, the high-fidelity results, and their coherence with the fine mesh and the numerical model chosen, are analyzed form a physical point of view. First, the dimensionless wall distance parameter y⁺ is analyzed. One can interpret y⁺ as a local Reynolds number, such as its magnitude can be expected to determine the relative importance of viscous and turbulment processes. This parameter writes

$$y^{+}=\frac{u_{\tau }z}{v}\text{ with }{u}_{\tau }=\sqrt{\frac{{\tau }_w}{\rho }}\text{ and }{\tau }_w=\rho \nu {\left(\frac{dU}{dz}\right)}_{\text{wall }}\text{, }$$(27)

where u_τ is the friction velocity, z the absolute distance from the wall, U the velocity magnitude and τ_w the wall shear stress.

Turbulence models are usually not solved at the wall. Indeed, at some distance to the wall, boundary conditions are applied following the so called wall functions. Depending on the turbulence model used, different wall functions are applied. Hence, depending on the turbulence model and the corresponding wall function, a value of y⁺ at the wall mesh cell is required. As indicated before, here the k − ω SST model is applied using a wall function. In such a case, the OpenFOAM documentation recommends a y⁺ value between 30 and 300. If we calculate the y⁺ value across our parametric space, as it is shown for one solution of the problem in Figure 10, we get values from 50 to 150, which validates our high-fidelity mesh.

Then, the boundary layer velocity profiles are presented in Figure 11. The x coordinate of the velocity field is represented along a straight line, perpendicular to the wall, through the boundary layer for different x positions. On the one hand, upstream of the ramp, it can be observed how the boundary layer grows and merges with the opposite wall boundary layer leading to a fully developed velocity profile. On the other hand, over the ramp, it can be observed how the boundary layer thickness grows progressively til reaching a separation point where a negative x coordinate of the velocity field appears, leading to the boundary layer separation and to a flow recirculation zone.

Next, the x coordinate of the velocity field is analyzed with a contour plot, as it can be seen in Figure 12. The darkest blue contour zone represents the negative values of the field. Again, the flow recirculation over the ramp can be observed. Moreover, a big flow recirculation, downstream of the ramp, due to the previous flow detachment, can easily be seen. Further downstream, it can be observed how the flow reattaches to the wall.

Finally, the flow is studied using streamlines colored by the x coordinate of the velocity field, as it is presented in Figure 13 for 2 different solutions of the problem. On the top solution, both recirculation zones can again be observed. Moreover, the plane of interest, where the coarse-CFD solution is corrected, is represented by a violet line. It should be noted that, depending on the solution of the problem, the recirculation zone can be included or not in this plane of interest.

The resolution of the problems with both meshes is now presented in the plane of interest PoI (i.e. x= l_PoI = 100 cm) for the test data set, as it can be seen in Figure 14. It can be observed that the flow behavior depends on the mesh. Indeed, a difference in the amplitude, position and shape of the velocity field can clearly be seen.

Then, the trained OT-based “digital twin” is applied to the “virtual twin” simulations of the test set. The database is composed by 100 simulations: 90 for the training and 10 for the test. Note that the training set is sampled following a classic Latin Hypercube Sampling strategy [39,40]. The results are presented in Figure 15. It can observed that the OT-based correction leads to a solution very close to the high-fidelity data. In order to quantify the accuracy improvement, the original coarse simulations and the corrected ones are compared with the high-fidelity data. The three error metrics are presented in the Table 2. It can be observed that the original ignorance gap between the “virtual twin” and the high-fidelity simulations is considerably reduced thanks to the OT-based ignorance model.

Fig. 4 Problem geometry schema.

Fig. 5 Parameterized geometry: the number of the nodes per edge for the coarse and HF meshes are indicated in red and blue respectively.

Fig. 6 High-fidelity mesh.

Fig. 7 Coarse mesh.

Fig. 8 Norm of the velocity field in a longitudinal cut of the channel for both the HF mesh (top) and the coarse one (bottom).

Fig. 9 Norm of the velocity field in a perpendicular cut of the channel at the plane of interest PoI for both the HF mesh (left) and the coarse one (right).

Fig. 10 Contour plot of the y+ coefficient for a given solution of the problem. The high-fidelity mesh is plotted over the contour plot.

Fig. 11 Boundary layer velocity profiles at different x values for a given solution of the problem: (a) upstream of the ramp and (b) over the ramp.

Fig. 12 Contour plot of the x coordinate of the velocity field for a given solution of the problem.

Fig. 13 Streamlines of the flow, colored by the intensity of the x coordinate of the velocity field, for 2 different solutions of the problem. The plane of interest is represented by a violet straight line at x = lPoI.

Fig. 14 Test set data points: (a) OpenFOAM solutions of the problem (23) corresponding to the high-fidelity simultaion. (b) OpenFOAM solutions of the problem (23) corresponding to the coarse simulation.

Fig. 15 “Hybrid twin” results in the test set: “virtual twin” simulations corrected by the OT-based “digital twin”.

5 Conclusion

FD is today a widely used tool for scientific investigation and analysis of complex systems in engineering, allowing a numerical resolution of realistic turbulent flows which cannot be solved analytically. Although all the work done to improve its accuracy and cost, CFD simulations remain today either not to be trusted completely or too computationally expensive to run in many industrial contexts. The “hybrid twin” rationale brings a solution to this problem by allowing a correction of fast but inaccurate coarse simulations, bringing them close to precise but costly high-fidelity simulations. However, in fields such as fluid dynamics, filling the ignorance gap of the mesher by classical data-driven models leads to non physical results. Combining the “hybrid twin” rationale with the simplified Optimal Transport Monge problem, our approach leads to an OT-based “digital twin” able to correct “virtual twin” coarse simulations from an OT point of view. This OT-based “hybrid twin” methodology was already presented in a previous paper of the authors in order to correct the simulation-experiment ignorance gap. In this new paper, the same approach has been applied to a completely different scenario, correcting the coarse-HF simulation ignorance gap. Therefore, the previously proposed OT-based “hybrid twin” methodology is shown to be able to correct coarse numerical simulations giving solutions very close to the high-fidelity counterpart data and leading, thus, to a faster and cheaper access to data almost as much accurate as precise but costly numerical solutions.

Acknowledgments

We thank Angelo Pasquale of PIMM Laboratory at Arts et Métiers Institute of Technology for providing helpful assistance and advice on the CFD simulations with OpenFOAM. The research work was carried out at Stellantis as part of a CIFRE (Conventions Industrielles de Formation par la REcherche) thesis.

Funding

This research received no external funding.

Conflicts of interest

The authors have nothing to disclose.

Data availability statement

The data presented in this study are available on request from the corresponding author.

Author contribution statement

Conceptualization, S.T., V.C., A.A., V.H., and F.C.; Methodology, S.T., V.C., A.A. and F.C.; Software, S.T. V.C. and A.A.; Validation, S.T.; Formal Analysis, S.T. and V.C.; Investigation, S.T. and V.C.; Resources, V.H. and F.C.; Data Curation, S.T. and V.C.; Writing—Original draft preparation, S.T.; Writing— Review and editing, S.T., V.C., A.A., V.H., and F.C.; Visualization, S.T.; Supervision, V.H. and F.C.; Project administration, V.H. and F.C.; Funding acquisition, V.H. and F.C.; All authors have read and agreed to the published version of the manuscript.

References

All Tables

All Figures

	Fig. 1 Discrete OT formulation for N= 4 mines and M=3 factories. The resource produced by the mines is: a1 =3, a2 =1, a3 = 1 and a4 = 2. The resource consumed by the factories is: b1 =4, b2 = 1 and b3 = 2. The Euclidean distance traveled by the resource is the cost to minimize.
In the text

	Fig. 2 2D Monge problem with N=M and an= bm = 1/N. Mines and factories are represented by red and blue circles respectively. Black arrows illustrate the optimal matching.
In the text

	Fig. 3 Summary diagram of the methodology: the offline training of the “digital twin” is colored in blue and the online OT-based “hybrid twin” approach in red.
In the text

	Fig. 4 Problem geometry schema.
In the text

	Fig. 5 Parameterized geometry: the number of the nodes per edge for the coarse and HF meshes are indicated in red and blue respectively.
In the text

	Fig. 6 High-fidelity mesh.
In the text

	Fig. 7 Coarse mesh.
In the text

	Fig. 8 Norm of the velocity field in a longitudinal cut of the channel for both the HF mesh (top) and the coarse one (bottom).
In the text

	Fig. 9 Norm of the velocity field in a perpendicular cut of the channel at the plane of interest PoI for both the HF mesh (left) and the coarse one (right).
In the text

	Fig. 10 Contour plot of the y+ coefficient for a given solution of the problem. The high-fidelity mesh is plotted over the contour plot.
In the text

	Fig. 11 Boundary layer velocity profiles at different x values for a given solution of the problem: (a) upstream of the ramp and (b) over the ramp.
In the text

	Fig. 12 Contour plot of the x coordinate of the velocity field for a given solution of the problem.
In the text

	Fig. 13 Streamlines of the flow, colored by the intensity of the x coordinate of the velocity field, for 2 different solutions of the problem. The plane of interest is represented by a violet straight line at x = lPoI.
In the text

	Fig. 14 Test set data points: (a) OpenFOAM solutions of the problem (23) corresponding to the high-fidelity simultaion. (b) OpenFOAM solutions of the problem (23) corresponding to the coarse simulation.
In the text

	Fig. 15 “Hybrid twin” results in the test set: “virtual twin” simulations corrected by the OT-based “digital twin”.
In the text