domain: math spi

text:
Kun Li
email: likun@swufe.edu.cn

Yixin Li

Liang Li

Stéphane Lanteri
email: stephane.lanteri@inria.fr

Surrogate modeling of time-domain electromagnetic wave propagation via dynamic mode decomposition and radial basis function

Keywords: Surrogate modeling, non-intrusive model order reduction, proper orthogonal decomposition, dynamic mode decomposition, radial basis function

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Numerical simulations of electromagnetic wave propagation problems primarily rely on spatial and temporal discretizations of the system of time-domain Maxwell equations using finite difference or finite element type methods. For complex and realistic 3-D situations, such a process can be computationally prohibitive, especially in view of many-query analyses (e.g., design, control, optimization, and uncertainty quantification (UQ) [START_REF] Benner | A survey of projection-based model reduction methods for parametric dynamical systems[END_REF][START_REF] Hesthaven | Certified reduced basis methods for parametrized partial differential equations[END_REF][START_REF] Peherstorfer | Survey of multifidelity methods in uncertainty propagation, inference, and optimization[END_REF][START_REF] Hesthaven | Reduced basis methods for time-dependent problems[END_REF]). Therefore, developing cost-effective and high-precision surrogate models are of great practical significance. Among the different possible approaches for building the surrogate model of a given partial differential equation (PDE) system, the data-driven based non-intrusive model order reduction (NIMOR) approaches, which are generally implemented in a fully decoupled offline-online paradigm [START_REF] Hesthaven | Non-intrusive reduced order modeling of nonlinear problems using neural networks[END_REF][START_REF] Wang | Non-intrusive reduced order modeling of unsteady flows using artificial neural networks with application to a combustion problem[END_REF],

have recently shown new promises due to their capability of handling nonlinear or/and high dimensional problems.

During the offline stage, a set of time-and parameter-independent reduced basis (RB) functions are extracted from a collection of high-fidelity (HF) solutions (also known as snapshots), which are always generated by numerical simulations or experiments at some time/parameter locations. A reduced space, which is spaned by the RB functions, can effectively capture the main dynamics of the original full-order model (FOM). A proper orthogonal decomposition (POD) method [START_REF] Berkooz | The proper orthogonal decomposition in the analysis of turbulent flows[END_REF][START_REF] Liang | Proper orthogonal decomposition and its applications -Part I: Theory[END_REF] based on a singular value decomposition (SVD) is usually used for extracting the RB functions. However, for time-domain electromagnetic wave propagation problems, the snapshot matrix has a considerable size, leading to very a time-consuming process if the POD method is used directly. To overcome this difficulty, we use a nested POD method [START_REF] Li | Non-intrusive reduced-order modeling of parameterized electromagnetic scattering problems using cubic spline interpolation[END_REF], in which the POD method is implemented separately for time and parameter spaces. As an approximation to the full-order solutions, the reduced-order solutions are recovered as a linear combination of the RB functions, where only the combination coefficients, also named the reduced-order or projection coefficients will be calculated during the online stage. In the NIMOR method, the reduced-order coefficients are approximated over the parameter domain by only using a database of reduced-order information [START_REF] Casenave | A nonintrusive reduced basis method applied to aeroacoustic simulations[END_REF].

More recently, NIMOR methods (whose mainly advantage is that they are 'equation-free' and do not require access to the governing equations of original FOM) have been developed for solving large-scale complex physical problems. For instance, Hesthaven and Ubbiali [START_REF] Hesthaven | Non-intrusive reduced order modeling of nonlinear problems using neural networks[END_REF] proposed a non-intrusive RB method combining the POD method with feed forward neural networks (FNNs) for parameterized steady-state PDEs, which is extended to time-dependent nonlinear problems [START_REF] Wang | Non-intrusive reduced order modeling of unsteady flows using artificial neural networks with application to a combustion problem[END_REF], in particular to unsteady flows. Chen and et al. [START_REF] Chen | Physics-informed machine learning for reduced-order modeling of nonlinear problems[END_REF] proposed a physics-reinforced neural network (PRNN) for nonlinear problems, which is a combination of projection driven neural network (PDNN) with physics-informed neural network (PINN) [START_REF] Raissi | Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[END_REF]. Fu and et al. [START_REF] Fu | Physics-data combined machine learning for parametric reduced-order modelling of nonlinear dynamical systems in small-data regimes[END_REF] developed a physics-data combined machine learning (PDCML) method for parametric reduced-order modelling of nonlinear dynamical systems in small-data regimes. In addition, diffierent regression and interpolation methods, such as Gaussian process regression (GPR) [START_REF] Guo | Data-driven reduced order modeling for time-dependent problems[END_REF][START_REF] Yu | Adaptive non-intrusive reduced order modeling for compressible flows[END_REF][START_REF] Li | A non-intrusive model order reduction approach for parameterized time-domain Maxwell's equations[END_REF] , radial basis function (RBF) [START_REF] Xiao | A non-intrusive reduced-order model for compressible fluid and fractured solid coupling and its application to blasting[END_REF][START_REF] Xiao | Error estimation of the parametric non-intrusive reduced order model using machine learning[END_REF] and cubic spline interpolation (CSI) [START_REF] Li | Non-intrusive reduced-order modeling of parameterized electromagnetic scattering problems using cubic spline interpolation[END_REF][START_REF] He | Model order reduction for parameterized electromagnetic problems using matrix decomposition and deep neural networks[END_REF], have also been adopted to approximate the mapping from the time/parameter to the reduced-order coefficients.

An alternative approach to construct the non-intrusive reduced-order model (ROM) is the dynamic mode decomposition (DMD) method, which is also an 'equation-free' method. The DMD method was developed by Schmid [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF] and has gained popularity in last years as ROM for dynamical systems [START_REF] Kutz | Dynamic mode decomposition: data-driven modeling of complex systems[END_REF]. The accuracy of the DMD method has been proved numerically [START_REF] Duke | An error analysis of the dynamic mode decomposition[END_REF] and theoretically [START_REF] Lu | Prediction accuracy of dynamic mode decomposition[END_REF]. Moreover, many variants of the DMD method have also been proposed to improve their efficiency and applicability [START_REF] Le Clainche | Higher order dynamic mode decomposition[END_REF][START_REF] Drmac | Data driven modal decompositions: analysis and enhancements[END_REF][START_REF] Lu | Lagrangian dynamic mode decomposition for construction of reduced-order models of advection-dominated phenomena[END_REF]. For instance, Clainchey and et al. proposed a higher order DMD (HODMD) method [START_REF] Le Clainche | Higher order dynamic mode decomposition[END_REF] that uses time-lagged snapshots and can be seen as superimposed DMD in a sliding window, which is also applied successfully to identify and extrapolate flow patterns [START_REF] Le Clainche | Higher order dynamic mode decomposition to identify and extrapolate flow patterns[END_REF]. Despite the DMD method is wildely used for reduced-order modelling of time-dependent problems, it is significantly more challenging for parameterized problems [START_REF] Gao | A reduced order method for nonlinear parameterized partial differential equations using dynamic mode decomposition coupled with k-nearest-neighbors regression[END_REF][START_REF] Hess | A data-driven surrogate modeling approach for time-dependent incompressible navier-stokes equations with dynamic mode decomposition and manifold interpolation[END_REF][START_REF] Huhn | Parametric dynamic mode decomposition for reduced order modeling[END_REF][START_REF] Lin | A dynamic mode decomposition based reduced-order model for parameterized time-dependent partial differential equations[END_REF]. For solving this limitation, Syadi, Schmid and et al. proposed a DMD-based parameter-dependent ROM framework for bifurcation analysis [START_REF] Sayadi | Parametrized data-driven decomposition for bifurcation analysis, with application to thermo-acoustically unstable systems[END_REF], where the time series HF solutions for different parameter values are "stacked" to form an augmented snapshot matrix, and the stacked DMD modes, which are calculated by the augmented snapshot matrix, are interpolated to form a new DMD mode for a new parameter. Further, Huhn and et al. developed two computationally more efficient parametric DMD approaches for diffusiondominated nonlinear dynamical problems, including the reduced eigen-pair interpolation (rEPI) and the reduced Koopman operator interpolation (rKOI). In addition, the DMD methods are also employed to build the prediction models of snapshots or latent variables at each training parameter value, and then some interpolation techniques are adopted to recover the snapshots or latent variables for some new parameter values [START_REF] Gao | A reduced order method for nonlinear parameterized partial differential equations using dynamic mode decomposition coupled with k-nearest-neighbors regression[END_REF][START_REF] Duan | Non-intrusive data-driven reduced-order modeling for timedependent parametrized problems[END_REF].

In this study, we attempt to investigate the predictive ability of a higher-order extension of the DMD method [START_REF] Le Clainche | Higher order dynamic mode decomposition[END_REF], i.e. HODMD, for time-domain electromagnetic wave propagation from both interpolation and extrapolation perspectives. In order to extend the DMD method to parametrized problems, we propose the RBF method to interpolate the reduced-order coefficient vectors for some new parameters. The proposed NIMOR method is purely data-driven method, where the offline and online stages are fully decoupled. During the offline stage, the HF solutions as snapshots are collected, the RB bases are extracted, and the DMD models of reduced-order coefficient vectors for all training parameter values are established. In order to predict a RB solution at a new time and parameter value on the online stage, the reduced-order coefficient vector at the new time is first predicted via the surrogate models formed by the DMD models, then the new reduced-order coefficient vector in the parameter space is interpolated via the RBF method, the RB solution is finally recovered in reduced space formed by the RB bases. The main innovation of this work is the use of a multi-step procedure consisting of the nested POD, DMD and RBF methods to accurately recover field solutions from a small number of large-scale electromagnetic simulations.

The remainder of this paper is organized as follows. We briefly introduce the time-domain Maxwell's equations and the training data set in Section 2. Section 3 describes the NIMOR method including the nested POD, DMD, and RBF methods. Section 4 presents and discusses some numerical results, including the scattering of a plane wave by a dielectric disk, by a multi-layer disk, and by a 3-D dielectric sphere. Section 5 concludes this study.

Full-order solutions

Time-domain Maxwell's equations

We consider the electromagnetic scattering problems governed by the normalized timedomain Maxwell's equations with the Silver-Müller absorbing boundary condition (ABC) and the initial conditions

                   µ r ∂H(x, t) ∂t + ∇ × E(x, t) = 0, ∀(x, t) ∈ Ω × T , ε r ∂E(x, t) ∂t -∇ × H(x, t) = 0, ∀(x, t) ∈ Ω × T , L(E(x, t), H(x, t)) = L(E inc (x, t), H inc (x, t)), ∀(x, t) ∈ ∂Ω × T , E(x, 0) = E 0 (x), H(x, 0) = H 0 (x), ∀x ∈ Ω, (1)
where Ω ⊂ R 2 or R 3 is the spatial domain, T = [0, T f ] is the time interval, and ∂Ω is the boundary of Ω; E(x, t) and H(x, t) respectively denote the electric and magnetic fields; ε r and µ r are respectively the relative permittivity and permeability; the Silver-Müller ABC operator is

L(E(x, t), H(x, t)) = n × E(x, t) + Zn × (n × H(x, t))
, n denotes the unit normal vector pointing outward to ∂Ω, E inc (x, t) and H inc (x, t) stand for the incident fields, and Z = µ r /ε r ; E 0 and H 0 stand for some given functions. In this study, we consider θ

= (ε r,1 , ε r,2 , • • • , ε r,p ) ∈ P ⊂ R p
as the problem's parameters with ε r,i (i = 1, 2, • • • , p) being the relative permittivity in the i-th domain of Ω, P being the parameter domain, and p being the number of parameters.

The data set of the high-fidelity solutions

Let P tr ℏ = {θ (1) , θ (2) , • • • , θ (Np) } ⊂ P be the training parameter set with uniform or uneven sampling. For each parameter θ (i) ∈ P tr ℏ , the corresponding high-fidelity (HF) solutions or snapshots for (1) can be obtained by using classical numerical methods or physical experiments in a time partition T ℏ = {t (0) , t (1) , • • • , t (m) }. Then, we can formulate the following time snapshot matrices

S i u (i = 1, 2, • • • , N p ) in a time sampling T tr ℏ = {t (n 1 ) , t (n 2 ) , • • • , t (n N t ) } ⊆ T ℏ for each parameter sampling point S i u =        u (n 1 ) h,1 (θ (i) ) -u h,1 u (n 2 ) h,1 (θ (i) ) -u h,1 • • • u (n N t ) h,1 (θ (i) ) -u h,1 u (n 1 ) h,2 (θ (i) ) -u h,2 E (n 2 ) h,2 (θ (i) ) -u h,2 • • • u (n N t ) h,2 (θ (i) ) -u h,2 . . . . . . . . . . . . u (n 1 ) h,N h (θ (i) ) -u h,N h u (n 2 ) h,N h (θ (i) ) -u h,N h • • • u (n N t ) h,N h (θ (i) ) -u h,N h        ∈ R N h ×Nt , ( 2
)
and then the snapshot matrix for all parameter in P tr ℏ is

S u = S 1 u | S 2 u | • • • | S Np u ∈ R N h ×Ntp , u ∈ {E, H}, (3)
where u

(n j ) h,k (θ (i) ) = u h,k (t (n j ) , θ (i)
) is the HF solution, u h ∈ R N h is the mean value, N h is the number of degrees of freedom (DoF) in space, and N tp = N t • N p is the number of all snapshots.

Parameterized non-intrusive model order reduction method

In this section, we consider a possible application of the POD, DMD and RBF methods in the parameterized non-intrusive reduced-order modeling.

Nested proper orthogonal decomposition method

Based on the above data preparation, we perform a low-rank approximation by constructing a low-dimension reduced basis (RB) subspace V rb,u with rank L u ≪ min{N h , N tp } (u ∈ {E, H}), which has the lowest dimension possible, while still maintaining a good approximation of the original dataset. The proper orthogonal decomposition (POD) method with snapshots is a popular approach to form the RB subspace. However, for the parameterized electromagnetic scattering problems, the snapshot matrix S u (u ∈ {E, H}) has a considerable size because time is added as an input parameter, so a direct approach of a singular value decomposition (SVD) is very time consuming. Therefore, a nested POD method [START_REF] Li | Non-intrusive reduced-order modeling of parameterized electromagnetic scattering problems using cubic spline interpolation[END_REF] is employed to extract the timeand parameter-independent RB functions, where the POD method is implemented separately for time and parameter spaces. The detailed pseudo-code implementation of the nested POD method is shown in Algorithm 1. Figure 1 shows the framework of the classic POD and nested POD. For a study of the corresponding projection errors we refer to [START_REF] Li | Non-intrusive reduced-order modeling of parameterized electromagnetic scattering problems using cubic spline interpolation[END_REF].

Let

Ψ u = [ψ u,1 , ψ u,2 , • • • , ψ u,Lu ] ∈ R N h ×Lu (u ∈ {E,
H}) be the RB matrix. According to the projection theory, the RB solution u r h (t, θ), an approximation to the full-order solution u h (t, θ), can be recovered as

u h (t, θ) ≈ u r h (t, θ) = Ψ u α u (t, θ) + u h = Lu i=1 α u,i (t, θ)ϕ u,i + u h , u ∈ {E, H}, (4)
where Lu denotes the reduced-order or projection coefficient vector. In particular, the following reduced-order coefficient vectors can be calculated approximately in the training time and parameter sets

α u (t, θ) = [α u,1 (t, θ), α u,2 (t, θ), • • • , α u,Lu (t, θ)] T ∈ R
α u (t (n i ) , θ (j) ) = Ψ T u (u h (t (n i ) , θ (j) ) -u h ), (t (n i ) , θ (j) ) ∈ T tr ℏ × P tr ℏ , u ∈ {E, H}, (5)
which can be used as the training data.

Ψ i u = POD(S i u , ρ t ) with the size l i u ≪ min{N h , N t } 3: end for 4: S u = [Ψ 1 u , Ψ 2 u , • • • , Ψ Np u ] 5: Ψ u = POD( S u , ρ θ ) 6: function Ψ = POD(S, ρ) do 7: [U, Σ, V] = svd(S) 8: k = argmin{E(k) : E(k) ≥ 1 -ρ} with E(k) = k i=1 σ 2 i / r i=1 σ 2
i being the relative information content, where r is the rank of S and σ i is the singular value of S Remark 1. In Algorithm 1, the effects of (ρ t , ρ θ ) on the average projection and NIMOR errors have been investigated to provide guidance on the choice of truncation tolerances in our previous work [START_REF] Li | Non-intrusive reduced-order modeling of parameterized electromagnetic scattering problems using cubic spline interpolation[END_REF]. So we directly select the truncation tolerances and no longer test the numerical convergences in this study.

Dynamic mode decomposition method

The dynamic mode decomposition (DMD) method with data-driven framework is proposed for model reduction of time-dependent problems. Based on the database of the reduced-order coefficient vectors over a time domain for each parameter, the DMD is defined as an approximate eigendecomposition of the operator that relates one time step to the next. For each training parameter θ (j) ∈ P tr ℏ , we seek an operator (or a Koopman matrix)

K (j) u ∈ R Lu×Lu such that α u (t (n i+1 ) , θ (j) ) = K (j) u α u (t (n i ) , θ (j) ), i = 1, 2, • • • , N t -1, u ∈ {E, H}, (6)
where the operator

K (j)
u can be considered as a finite approximation of the Koopman operator [START_REF] Koopman | Hamiltonian systems and transformation in Hilbert space[END_REF], which is infinite-dimensional. In order to establish the operator K (j) u satisfying ( 6) to an acceptable approximation, we consider the time sampling T tr ℏ with equal interval ∆t. In this study, the snapshots are sampled at the first 70% of a period with equal interval ∆t for each training parameter point. In order to compute the operator K (j) u , we arrange two following sequences (whose columns are the reduced-order coefficient vectors)

X (j) αu = [α u (t (n 1 ) , θ (j) ), α u (t (n 2 ) , θ (j) ), • • • , α u (t (n N t -1 ) , θ (j) )] ∈ R Lu×(Nt-1) , (7)
and

Y (j) αu = [α u (t (n 2 ) , θ (j) ), α u (t (n 3 ) , θ (j) ), • • • , α u (t (n N t ) , θ (j) )] ∈ R Lu×(Nt-1) . (8)
One can assume that X

Y (j) αu ≈ K (j) u X (j) αu , j = 1, 2, • • • , N p , u ∈ {E, H}. (9)
Let † be the Moore-Penrose pseudoinverse, one can find the matrix K

(j) u = Y (j) αu (X (j)
αu ) † is a best-fit linear operator relating X αu to Y αu [START_REF] Penrose | On best approximate solutions of linear matrix equations[END_REF]. However, this calculation is expensive and often difficult to handle [START_REF] Andreuzzi | A dynamic mode decomposition extension for the forecasting of parametric dynamical systems[END_REF]. Therefore we adopt another method, in which the DMD can be implemented by the Algorithm 2 for the approximation of the eigenvalues and eigenvectors of K [START_REF] Chen | Physics-informed machine learning for reduced-order modeling of nonlinear problems[END_REF] In particular, if L u ≈ N t , the quality of the results obtained by the DMD method may decay [START_REF] Le Clainche | Higher order dynamic mode decomposition[END_REF]; the situation could be worse if L u < N t [START_REF] Andreuzzi | A dynamic mode decomposition extension for the forecasting of parametric dynamical systems[END_REF]. So, to avoid this situation, we consider a higher order DMD (termed as HODMD or DMD-s) method [START_REF] Le Clainche | Higher order dynamic mode decomposition[END_REF], which is able to extract more significative modes from low-dimensional data. In the HODMD method, ( 6) is modified to the following form

α u (t * , θ (j) ) ≈ α DMD u (t * , θ (j) ) = Φ (j) u (Λ (j) u ) (t * -t (n 1 ) )/∆t (Φ (j) u ) † α u (t (n 1 ) , θ (j) ), j = 1, 2, • • • , N p , u ∈ {E, H}.
α u (t (n i+s ) , θ (j) ) = K (j) u,1 α u (t (n i ) , θ (j) ) + K (j) u,2 α u (t (n i+1 ) , θ (j) ) + • • • + K (j) u,s α u (t (n i+s-1 ) , θ (j) ), (12) Algorithm 2 DMD method Input: Two matrixes X (j) αu , and Y (j) αu , j = 1, 2, • • • , N p , u ∈ {E, H} Output: DMD modes Φ (j) u and eigenvalue matrices Λ (j) u , j = 1, 2, • • • , N p , u ∈ {E, H} 1: Perform the SVD of X (j) αu ≈ U (j) u,r j Σ (j) u,r j (V (j)
u,r j ) T with r j being the rank of X (j) αu 2: Calculate the low-rank projection K

K (j) u = (U (j) u,r j ) T K (j) u U (j) u,r j = (U (j) u,r j ) T Y (j) αu (X (j) αu ) † U (j) u,r j = (U (j) u,r j ) T Y (j) αu V (j) u,r j (Σ (j) u,r j ) -1 (U (j) u,r j ) T U (j) u,r j = (U (j) u,r j ) T Y (j) αu V (j) u,r j (Σ (j) u,r j ) -1 (10) 3: Perform the eigendecomposition of K (j) u W (j) u = W (j) u Λ (j) u with Λ (j)
u being the eigenvalues matrix and W (j) u being the eigenvectors matrix 4: Obtain the DMD modes Φ

(j) u = Y (j) αu V (j) u,r j (Σ (j) u,r j ) -1 W (j) u corresponding to the particular DMD eigenvalues in Λ (j) u
where s ≥ 1 is tunable. Moreover, ( 12) can be rewritten into the following compact form [START_REF] Fu | Physics-data combined machine learning for parametric reduced-order modelling of nonlinear dynamical systems in small-data regimes[END_REF] where the modified vectors α u (t (n i ) , θ (j) ) and the modified Koopman matrix K u (j) are respectively defined as

α u (t (n i+1 ) , θ (j) ) = K u (j) α u (t (n i ) , θ (j) ), i = 1, 2, • • • , N t -s, j = 1, 2, • • • , N p , u ∈ {E, H},
K u (j) =        0 I 0 • • • 0 0 0 0 I • • • 0 0 • • • • • • • • • • • • • • • • • • 0 0 0 • • • I 0 K (j) u,1 K (j) u,2 K (j) u,3 • • • K (j) u,s-1 K (j) u,s        ∈ R (sLu)×(sLu) , ( 14
)
and

α u (t (n i ) , θ (j) ) = [α u (t (n i ) , θ (j) ), α u (t (n i+1 ) , θ (j) ), • • • , α u (t (n i+s-1 ) , θ (j) )] T ∈ R (sLu) , ( 15
)
with 0 being the L u × L u null matrix and I being the L u × L u unit matrix, leading to [START_REF] Berkooz | The proper orthogonal decomposition in the analysis of turbulent flows[END_REF] and ( 8) are respectively modified as

X (j) αu = [ α u (t (n 1 ) , θ (j) ), α u (t (n 2 ) , θ (j) ), • • • , α u (t (n N t -s ) , θ (j) )] ∈ R (sLu)×(Nt-s) , ( 16
)
and

Y (j) αu = [ α u (t (n 2 ) , θ (j) ), α u (t (n 3 ) , θ (j) ), • • • , α u (t (n N t -s+1 ) , θ (j) )] ∈ R (sLu)×(Nt-s) . ( 17
)
So, we can perform the Algorithm 2 to find an approximation of the operator K u (j) in the HODMD method.

Remark 2. In Algorithm 2, the number r j is selected in terms of the rank of the matrix X

(j) αu , R j ≤ min{L u , N t -1}, to be such that EE(r j ) = R j i=r j +1 σ 2 j,i / R j i=1 σ 2 j,i ≤ ρ SVD , ( 18
)
where EE(r j ) is the relative root mean square (RMS) error, ρ SVD is the error threshold, and σ j,i is the singular value of X (j) αu sorted in decreasing order. Particularly, the reduced-order coefficient vectors are of low dimension so that ρ SVD is chosen as 1 × 10 -10 in this study. Besides, a possibly smaller number N DMD (≤ r j ) of DMD modes, after the four step of Algorithm 2, can be selected imposing that

AMP (N DMD +1) /AMP (1) < ρ DMD , ( 19
)
where AMP are the DMD amplitudes sorted in decreasing order, which can be computed by AMP = (Φ

(j) u ) † α u (t (n 1 )
, θ (j) ). For other alternative method of the computation of DMD mode amplitudes, one can see [START_REF] Le Clainche | Higher order dynamic mode decomposition[END_REF]. Remark 3. In Algorithm 2, for the HODMD method, the time complexities of steps 1 to 4 are

O(sL u (N t -s) 2 ), O(r j (sL u )(N t -s) + (N t -s)r 2 j + r 2 j ), O(r 3 j ), and O((sL u )(N t -s)r j + (sL u )r j + (sL u )r 2 j )
, respectively.

Radial basis function interpolation

Based on [START_REF] Chen | Physics-informed machine learning for reduced-order modeling of nonlinear problems[END_REF], one can predict the reduced-order coefficient vector in future time instants for the training parameter θ (j) ∈ P tr ℏ , and then obtain the RB solution via (4). However, it is inapplicable for parameterized problems. In order to resolve this problem, we propose the radial basis function (RBF) to interpolate the reduced-order coefficient vector for a new or untrained parameter θ * ∈ P\P tr ℏ . Let

F α DMD u (θ) = Np k=1 ϖ u,k φ(∥θ -θ (k) ∥), u ∈ {E, H}, (20)
where ϖ u,k is the weighting vector associated with the parameter θ (k) , and ∥• ∥ is the norm of a vector, usually chosen as Euclidean distance. Some well-known RBFs include Gaussian, linear spline, cubic spline, thin-plate spline, multi-quadratic, inverse multistory, and so on. In this work, we choose the Gaussian function, which has a form of

φ(d) = e -( d σ ) 2
with d being the radius or distance and σ being a shape parameter. In particular, choosing an optimal value of the shape parameter σ is very critical in the RBF method. Throughout this work, the shape parameter σ is selected by trials and tests. The weighting vector

ϖ u (j, i) = [ϖ u,1 (j, i), ϖ u,2 (j, i), • • • , ϖ u,Np (j, i)]
T for the j-th component of the reduced-order coefficient vector α DMD u at the n i -th time point is then determined by solving the following linear system

A RBF ϖ u (j, i) = B u,RBF (j, i), u ∈ {E, H}, (21)
where the coefficient matrix A RBF and the right hand side B u,RBF (j, i) are respectively defined as [START_REF] Duke | An error analysis of the dynamic mode decomposition[END_REF] and

A RBF =     φ(∥θ (1) -θ (1) ∥) φ(∥θ (1) -θ (2) ∥) • • • φ(∥θ (1) -θ (Np) ∥) φ(∥θ (2) -θ (1) ∥) φ(∥θ (2) -θ (2) ∥) • • • φ(∥θ (2) -θ (Np) ∥) • • • • • • • • • • • • φ(∥θ (Np) -θ (1) ∥) φ(∥θ (Np) -θ (2) ∥) • • • φ(∥θ (Np) -θ (Np) ∥)     ∈ R Np×Np ,
B u,RBF (j, i) = [α DMD u,j (t (n i ) , θ (1) ), α DMD u,j (t (n i ) , θ (2) ), • • • , α DMD u,j (t (n i ) , θ (Np) )] T ∈ R Np . ( 23
)
For a new parameter θ * ∈ P, the j-th component of the reduced-order coefficient vector α DMD u,j (t (n i ) , θ * ) at the n i -th time point can be calculated as

α DMD u,j (t (n i ) , θ * ) ≈ α RBF u,j (t (n i ) , θ * ) = Np k=1 ϖ u,k (j, i)φ(∥θ * -θ (k) ∥), u ∈ {E, H}. (24)
Finally, for the new time t * and parameter θ * , the RB solution u r h (t * , θ * ) can be quickly recovered as

u r h (t * , θ * ) = Ψ u α u (t * , θ * ) + u h ≈ Ψ u α DMD u (t * , θ * ) + u h ≈ Ψ u α RBF u (t * , θ * ) + u h , u ∈ {E, H}. (25)

Workflow of the proposed NIMOR method

Based on the discussions in the first three subsections, we propose the workflow for the construction of NIMOR method in Algorithm 3, where the offline and online stages of the POD, DMD, and RBF-based NIMOR are decoupled. Figure 2 shows the overview of the proposed method. During the offline stage, we collect suitable snapshots, obtain the RB bases, and build the DMD models for all training parameter values. During the online stage, in order to predict a RB solution at a new time and parameter value, we first obtain the reduced-order coefficient vector at the new time via the DMD models, then interpolate the new reduced-order coefficient vector in the parameter space, and finally recover the RB solution via [START_REF] Drmac | Data driven modal decompositions: analysis and enhancements[END_REF].

Algorithm 3 NIMOR for time-domain electromagnetic wave propagation problems

Offline stage

Generate the HF solutions u h (t (n i ) , θ (j) ), (t (n i ) , θ (j) ) ∈ T tr ℏ × P tr ℏ , u ∈ {E, H} Compute the POD basis functions Ψ u via Algorithm 1 Obtain the training data α u (t (n i ) , θ (j) ) via (5) Compute the DMD modes Φ Remark 4. DMD method is used to establish the surrogate models of reduced-order coefficient vectors for all training parameter values, and RBF method is adopted to interpolate reducedorder coefficient vector for a new parameter value, which is more effective than direct predicting and interpolating snapshots.

Numerical experiments

In this section, some numerical experiments including the scattering of a plane wave by a dielectric disk, by a multi-layer disk and by a 3-D dielectric sphere are presented to validate the effectiveness and the accuracy of the POD, DMD and RBF-based NIMOR method. In order to evaluate the accuracy of numerical experiments, the relative error between the NIMOR solution and the HF solution for u ∈ {E, H} is defined

e u,NIMOR (t, θ) = ∥ u h (t, θ) -u r h (t, θ) ∥ R N h ∥ u h (t, θ) ∥ R N h = ∥ u h (t, θ) -(Ψ u α RBF u (t, θ) + u h ) ∥ R N h ∥ u h (t, θ) ∥ R N h , (26)
which will be compared with the relative projection error generated by the nested POD method

e u,Pro (t, θ) = ∥ (u h (t, θ) -u h ) -Ψ u Ψ T u (u h (t, θ) -u h ) ∥ R N h ∥ u h (t, θ) ∥ R N h , ( 27
)
whre || • || R N h denotes Euclidean norm. Moreover, the average relative errors on a testing time and parameter sampling T te ℏ × P te ℏ of size N te are calculated to measure the convergence of the NIMOR method, which are defined as

e u,NIMOR = (t,θ)∈T te ℏ ×P te ℏ e u,NIMOR (t, θ) N te , ( 28
)
and

e u,Pro = (t,θ)∈T te ℏ ×P te ℏ e u,Pro (t, θ) N te , ( 29
)
where the testing parameter set P te ℏ is generated via randomized latin-hypercube-sampling (LHS) method [START_REF] Sharon | Sampling: design and analysis[END_REF], and the testing time set T te ℏ is uniformly chosen within a period of the physical simulation.

The HF solutions to the scattering problem of a dielectric disk and a multi-layer disk are calculated by using a discontinuous Galerkin time-domain (DGTD) solver formulated on an unstructured simplicial mesh [START_REF] Hesthaven | Nodal discontinuous Galerkin methods: algorithms, analysis, and applications[END_REF], which combines a centered flux scheme with a second order leap-frog (LF2) time scheme. The HF solutions to the scattering problem of a 3-D dielectric sphere are obtained by using the Mie series method [START_REF] Chew | Waves and fields in inhomogenous media[END_REF], and the corresponding computer implementation has been developed by Zhu [START_REF] Zhu | Sphere scattering[END_REF]. The DGTD, Mie series and NIMOR methods are implemented in MATLAB and all simulations are run on a workstation equipped with an Intel Core i7-10700F CPU running at 2.90 GHz, and with 16 GB of RAM memory. All SVDs are implemented via the MATLAB function svd.

Scattering of a plane wave by a dielectric disk

We first consider the electromagnetic scattering of a plane wave by a dielectric disk, which is described by the following normalized form of 2-D time-domain transverse magnetic (TM) waves

     µ r ∂H(x, t) ∂t + ∇ × E(x, t) = 0, ∀(x, t) ∈ Ω × T , ε r ∂E(x, t) ∂t -∇ × H(x, t) = 0, ∀(x, t) ∈ Ω × T , ( 30
)
where E = E z (x, t) and H(x, t) = [H x (x, t), H y (x, t)] T . The differential operators in this 2-D setting are

∇ × E(x, t) = [ ∂E z (x, t) ∂y , - ∂E z (x, t) ∂x ] T , and ∇ × H(x, t) = ∂H y (x, t) ∂x - ∂H x (x, t) ∂y . ( 31
)
The incident plane waves take the form E inc z = cos(ωt -kx) and H inc y = -cos(ωt -kx) with ω = 2πf being the angular frequency of the incident wave frequency f = 300 MHz, and k being the wave number in vacuum. The computational domain is the square

Ω 2 = [-2.6 m, 2.6 m] × [-2.6 m, 2.6 m].
The radius of the disk is set to be 0.6 m, and the relative permeability is set to be 1 (non-magnetic material). The range of interest relative permittivity of the disk is [START_REF] Benner | A survey of projection-based model reduction methods for parametric dynamical systems[END_REF][START_REF] Hesthaven | Non-intrusive reduced order modeling of nonlinear problems using neural networks[END_REF], so, P = {θ : θ = ε r ∈ [1, 5]} ⊂ R. The medium exterior to the disk is assumed to be vacuum.

During the offline stage, the DGTD method is used to obtain the HF solutions in the training parameter set P tr ℏ with uniformly sampling, in which the total simulation time is set to be 50 periods (which corresponds to 50 m in normalized unit) and the time step ∆t is set to be 3.678 × 10 -3 m. For each selected training parameter, the training snapshot vectors are equally chosen from 49 m to 49.7 m (whose are the first 70% data in the last period), leading to T tr ℏ = {49.0024 m, 49.006 m, • • • , 49.6938 m, 49.6975 m}. A testing parameter set P te ℏ with size 40 generated by the LHS method and a testing time set T te ℏ with size 263 including all time points in the last period (i.e., T te ℏ = T tr ℏ ∪ {49.7012 m, • • • , 49.966 m}) are selected to evaluate the proposed method. The details for the training and testing datasets are listed in Table 1. The decay of the singular values in the nested POD method is shown in Figure 3, in which one can see that the singular values decrease rapidly indicating that a small dimension RB basis can capture the information of the original data. Particularly, the effects of (ρ t , ρ θ ) on the average projection and NIMOR errors have been investigated to provide guidance on the choice of truncation tolerances in our previous work [START_REF] Li | Non-intrusive reduced-order modeling of parameterized electromagnetic scattering problems using cubic spline interpolation[END_REF]. So we directly select the truncation tolerances and no longer test the numerical convergences in this study. With the truncation tolerances ρ t = 1 × 10 -1 and ρ θ = 1 × 10 -5 in the nested POD methd, a set of L Hx = 59, L Hy = 21, and L Ez = 22 POD basis functions are extracted. The numerical convergences of the NIMOR method with respect to the number s of time lagged vectors in [START_REF] Raissi | Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[END_REF] and the error threshold ρ DMD in [START_REF] He | Model order reduction for parameterized electromagnetic problems using matrix decomposition and deep neural networks[END_REF] are shown in Figure 4, where the shape parameter σ is set to be 0.073 in [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF]. ρ DMD = 1 × 10 -2 based on [START_REF] He | Model order reduction for parameterized electromagnetic problems using matrix decomposition and deep neural networks[END_REF]. As we can see, the dimension of the NIMOR method is further reduced. Combining Figures 4 and5, we choose the number s = 50 and the error threshold ρ DMD = 1 × 10 -2 . The corresponding average relative errors are listed in Table 2. Some of the exact and approximate (HODMD method) reduced-order coefficient vectors in T te ℏ for all training parameters are displayed Figure 6. It is clear that the accuracy of the HODMD method is very well within the testing time region (where the set [49 m, 49.7 m] is the training time region, and the set [49.7 m, 50 m] is the prediction time region) for all training parameters. This implies the potential of the HODMD method in the prediction of the RB solutions outside of the training time region.

During the online stage, the reduced-order solutions are visualized for some testing parameters including θ (1) = 1.215, θ (2) = 2.215, θ (3) = 3.215, and θ (4) = 4.215. Figure 7 shows the time evolution fields H y and E z (HF and RB solutions) in T te ℏ at a given space point for above four testing parameters. It can be found that the RB solutions and the HF solutions match each other very well. Furthermore, the corresponding time evolution of the relative projection and NIMOR errors for the four testing parameters are displayed in Figure 8. In order to have a look on visual effects of electromagnetic fields, over the Fourier domain during the last oscillation period of the incident wave, we display in Figure 9 the total field of |H x | calculated by the NIMOR method for the parameter θ (3) = 3.215, which coincides with the HF solution. Finally, Table 3 shows the average computing time of the NIMOR and HF solutions. Although we spend a long time preparing the training data for the reduced-order modeling in the offine stage, it is worthwhile especially when solving the electromagnetic field value under multiple parameters, because the online computing time of the NIMOR method is greatly shortened.

Scattering of a plane wave by a multi-layer disk

Next, we consider a multivariate parameter problem, which is induced by a multi-layer disk illuminated by an incident plane wave, as shown in Figure 10. The computational domain is artificially truncated by a square with 6.4 m side length. The radius r i of each layer of the multi-layer disk and the corresponding relative permittivity range are summarized in Table 4. The medium exterior to the multi-layer disk is also assumed to be vacuum. So, the material

Layer i P (i) µ r,i r i 1 ε r,1 ∈ [5.0, 5.6] 1 0.15 Multi-layer disk 2 ε r,2 ∈ [3.25, 3.75] 1 0.3 3 ε r,3 ∈ [2.0, 2.5] 1 0.45 4 ε r,4 ∈ [1.25, 1.75] 1 0.6
parameters can be represented as a 4-D vector θ = (ε r,1 , ε r,2 , ε r,3 , ε r,4 ) ∈ P = P (1) × P (2) × P (3) × P (4) . Other physical information is the same as the case in section 4.1.

During the offline stage, the HF solutions, in the training parameter set P tr ℏ chosen by a Smolyak sparse grid method with approximation level L = 3, are generated by using the DGTD solver. The Smolyak sparse grid method is often used for the integration or interpolation of Figure 9: Comparison of the total fields of |Hx| for disk (top), multi-layer disk (middle), and sphere (bottom) for the testing parameter θ (3) , including the HF solutions (left), the RB solutions (middle), and the absolute errors (right).

Figure 10: . Geometry of the multi-layer disk.

multivariate functions, and can be used to deal with the so-called curse of dimensionality problem in high dimensional space, which is briefly described in the Appendix. The total simulation time is also set to be 50 m, and the time step ∆t is set to be 3.800 × 10 -3 m. Similarly, one can equally choose the training snapshots from 49 m to 49.7 m for each training parameter, and the corresponding training time set is T tr ℏ = {49.000 m, 49.004 m, • • • , 49.6945 m, 49.6984 m}. A testing parameter set P te ℏ with size 81 generated by the LHS method and a testing time set T te ℏ with size 253 including all time points in the last period are chosen to evaluate the NIMOR method. Table 1 summarizes the training and testing sets. With the truncation errors ρ t = 1 × 10 -1 and ρ θ = 1 × -5 for the nested POD method, L Hx = 16 POD basis functions are extracted for H x , L Hy = 15 for H y , and L Ez = 15 for E z . As with the previous problem, Figure 11 shows the numerical histories of the NIMOR method with respect to the number s and the error threshold ρ DMD , in which we set the shape parameter σ = 1. It can be seen that the errors e u,NIMOR are almost unchanged for varying s and ρ DMD for this case. The main reason is that the first few terms of amplitude play a key role. Figure 12 ( very close to the projection errors e u,Pro . In this test, s and ρ DMD are respectively set to be 5 and 1 × 10 -1 , and the corresponding sizes of the DMD modes for all training parameters, which are shown in Figure 12 (b), are all reduced to N DMD = 2. The corresponding average relative errors are presented in Table 2, and some of the approximation and exact reduced-order coefficient vectors, in the testing time set T te ℏ (where the set [49 m, 49.7 m] is the training time region, and the set [49.7 m, 50 m] is the prediction time region) for all some parameters, are shown in Figure 13. We can draw conclusions similar to Figure 6 from Figure 13.

During the online stage, the RB solutions are calculated by using the NIMOR method for four testing parameters θ (1) = {(5.125, 3.375, 2.125, 1.375)}, θ (2) = {(5.425, 3.625, 2.425, 1.625)}, θ (3) = {(5.125, 3.625, 2.125, 1.625)}, and θ (4) = {(5.425, 3.375, 2.425, 1.375)}. An evidence of the validity of the NIMOR method is shown in Figure 14, which reports the comparison between the time evolution of RB solutions at a given space point and their HF solutions. The time evolution of the relative projection and NIMOR errors are displayed in Figure 15. The total field of |H x | calculated by the NIMOR method for the parameter θ (3) = {(5.125, 3.625, 2.125, 1.625)} in the Fourier domain during the last period of wave oscillation is displayed in Figure 9. The time performance of the porposed NIMOR is presented in Table 3. Based on the numerical performance of the NIMOR method proposed here, this method is effective for multivariate parameter problem.

Scattering of a plane wave by a dielectric sphere

Finally, we consider the electromagnetic scattering of a plane wave by a 3-D dielectric sphere, where the plane wave is polarized in the x direction propagating in the +z direction. The incident field in frequency domain is expressed as

E inc = A 0 e -ikz n x , and H inc = E inc x η e -ikz n y , ( 32
)
where the magnitude A 0 is 1 V/m, the wave number k = ω √ µ r ε r with ω = 2πf being the angular frequency of the incident wave frequency f = 300 MHz, and i is the imaginary unit; n x = (1, 0, 0) T , and n y = (0, 1, 0) T . In particular, one can obtain the following time evolution solutions once obtaining the frequency-domain solutions

u(t, θ) = R(u(θ)e iωt ), θ ∈ P, u ∈ {E, H}, (33)
where R(• ) denotes the real part of a vector. The computational domain is the square

Ω 2 = [-0.2 m, 2 m] 3 .
The radius of the sphere is set to be 1 m. The range of corresponding relative permittivity is [START_REF] Hesthaven | Certified reduced basis methods for parametrized partial differential equations[END_REF][START_REF] Wang | Non-intrusive reduced order modeling of unsteady flows using artificial neural networks with application to a combustion problem[END_REF],

resulting in P = {θ : θ = ε r ∈ [2, 6]} ⊂ R.
The medium exterior to the sphere is assumed to be vacuum. During the offline stage, the HF solutions, at the 81 uniformly distributed parameter points, are obtained by using the Mie series method with ∆t = 1.0101 × 10 -2 m in the first period, which are implemented by Zhu [START_REF] Zhu | Sphere scattering[END_REF] for the frequency-domain solutions. In this test, the number of spatial grid points in x, y, and z directions is all 45, i.e., N h = 45 3 = 91125. The training time set is evenly chosen as T tr ℏ = {0 m, 1.0101 × 10 -2 m, set of L Hx = 38, L Hy = 37, L Hz = 34, L Ex = 36, L Ey = 32, and L Ez = 31 POD basis functions are extracted. Furthermore, s, ρ DMD , and σ are respectively set to be 5, 1 × 10 -1 , and 0.039, resulting in N DMD = 2. Table 2 shows the corresponding average relative errors, and Figure 16 and 17 present some of the approximation and exact reduced-order coefficient vectors in the testing set T te ℏ × P te ℏ for all training parameters, where the set [0 m, 0.7 m] is the training time region, and the set [0.7 m, 1 m] is the prediction time region.

During the online stage, in order to show the performance of NIMOR method for 3-D problem, we compare the RB solutions with the HF solutions generated by the Mie series method for four testing parameters θ (1) = 2.145, θ (2) = 3.145, θ (3) = 4.145, and θ (4) = 5.145. Figure 18 reports the comparison between the time evolution of the RB solutions at a given space point for the four testing parameters and the corresponding HF solutions. Figure 19 shows the time evolution of the corresponding relative projection and NIMOR errors. Figure 9 displays the total field of |H x | calculated by the NIMOR method for the parameter θ (3) = 4.145 in the Fourier domain during the first period of wave oscillation. Table 3 presents the average computing time for the offline and online stages.

Conclusion

In this paper, we proposed and studied an 'equation-free' NIMOR method for time-domain electromagnetic wave propagation. From the HF database on a given training parameter set with uniform or uneven sampling, a set of RB functions were extracted by using the nested POD method, which performs a prior dimensionality reduction. The DMD method, which is defined by finding the best-fit linear model of the Koopman operator, was then used to predict the lowdimensional reduced-order cofficient vectors in future time instants on this training parameter set, where the dimension of NIMOR method is further reduced. The RBF interpolation was finally used to estimate the reduced-order cofficient vectors at a new specified parameter value from the DMD calculations. The main contribution of this work is that we exploited the POD, DMD and RBF techniques to obtain a quick approximation for new specified parameters in future time instants. The accuracy and effectiveness of the NIMOR method were illustrated numerically by considering three classical test cases: the scattering of a plane wave by a dielectric disk, by a multi-layer disk, and by a 3-D dielectric sphere. Future research direction include the usage of other interpolation methods, e.g., manifold interpolation [START_REF] Hess | A data-driven surrogate modeling approach for time-dependent incompressible navier-stokes equations with dynamic mode decomposition and manifold interpolation[END_REF], to obtain the reduced DMD operator, or the reduction of geometric parameters.

Appendix: Smolyak sparse grid

The sampling process of the training parameter points over the parameter space P ∈ R p is briefly introduced. Let Q i k be an interpolation operator in the kth-direction of p-dimensional space, and the interpolation operator of Smolyak algorithm in whole p-dimensional space can be then expressed as

Q(L, p) = L-p+1≤|i|≤L (-1) L-|i| p -1 L -|i| (Q i 1 ⊗ Q i 2 ⊗ • • • ⊗ Q ip ), ( 34
)
where ⊗ is the tensor product operator, i k (k = 1, 2, • • • , p) denotes the approximation level of sparse grid in the kth-direction, L is the entire approximation level of sparse grid, and

|i| = i 1 + i 2 + • • • + i p .
The set of Smolyak sparse sampling points in (34) can be then derived as

S(L, p) = L-p+1≤|i|≤L (S i 1 ⊗ S i 2 ⊗ • • • ⊗ S ip ), (35)
with S i k = {α (1) , α (2) ,

• • • , α (M i k ) },
where S i k is the vector of sampling points with approximation level i k in the kth-direction, and M i k denotes the corresponding number of sampling points. Particularly, one can choose the Clenshaw-Curtis points as the sampling points on each dimension. Examples of Smolyak grids under approximation levels L = 1, L = 2, and L = 3 are illustrated in Figure 20 for 2-D case. For comparison, a tensor-product grid of 5 2 points is also shown. For a more detailed introduction of this technnique, one can see [START_REF] Judd | Smolyak method for solving dynamic economic models: Lagrange interpolation, anisotropic grid and adaptive domain[END_REF].

9 :

 9 Ψ = U(:, 1 : k) 10: end function

Figure 1 :

 1 Figure 1: Classic POD and nested POD.

  αu are connected by the following linear dynamics

u

  based on the data sets X (j) αu and Y (j) αu . Performing Algorithm 2, the eigenvalues and the corresponding eigenvectors of K (j) u can be approximated by Λ (j) u and the columns of Φ (j) u . Then, one can predict the solution for a new time t * in the future as follows

  F α DMD u (θ) denotes the interpolation function representing α DMD u (t, θ), which is a linear combination of the N p RBFs φ(• ) and has the following form

u

  (j = 1, 2, • • • , N p ) and the eigenvalue matrices Λ (j) u via Algorithm 2 with (13) (HODMD) Online stage Compute the DMD predictions α DMD u (t * , θ (j) ) for a new time t * via (11), θ (j) ∈ P tr ℏ Calculate the weighting vectors ϖ u (j, i) via (21) to (23) Calculate the RBF interpolation α RBF u (t * , θ * ) for a new parameter θ * ∈ P Evaluate the RB solution u r h (t * , θ * ) via (25)

Figure 2 :

 2 Figure 2: Overview of the proposed NIMOR method.

Figure 3 :

 3 Figure 3: Scattering of a plane wave by a dielectric disk: the decay of the singular values in the nested POD method. The first to third columns are the singular values of the 10th, 30th and 50th training parameter in the first step of the nested POD method. The fourth column is the singular values in the second step of the nested POD method.

Figure 5 Figure 4 :

 54 Figure 4: Scattering of a plane wave by a dielectric disk: the numerical convergences of (a) e E,NIMOR , and (b) e H,NIMOR on the testing set T te ℏ × P te ℏ with vary s and ρDMD, where eu,Pro is the average projection error of the HF solutions for E or H onto the corresponding RB subspace.

Figure 5 :

 5 Figure 5: Scattering of a plane wave by a dielectric disk: the size of the DMD modes with s = 20 and ρDMD = 1 × 10 -4 , s = 50 and ρDMD = 1 × 10 -2 , and the dimension of RB basis for Hx, Hy and Ez.

Figure 6 :

 6 Scattering of a plane wave by a dielectric disk: some of the exact and approximate (HODMD method) reduced-order coefficient vectors α u,l (t, θ) (l = 1, 2, 4, 6) of (a) Hx, (b) Hy and (c) Ez in T te ℏ for all training parameters. The dashed line indicates the end of the training time region, which is t = 49.7 m.

Figure 8 :

 8 Figure 8: Scattering of a plane wave by a dielectric disk: comparison of the time evolution of the relative errors of E (left) and H (right) for the four testing parameters in the testing time region T te ℏ .

Figure 11 :

 11 Figure 11: Scattering of a plane wave by a multi-layer disk: the numerical histories of (a) e E,NIMOR , and (b) e H,NIMOR on the testing set T te ℏ × P te ℏ with vary s and ρDMD, where eu,Pro is the average projection error of the HF solutions for E or H onto the corresponding RB subspace.

Figure 12 :

 12 Figure 12: Scattering of a plane wave by a multi-layer disk: (a) the DMD amplitudes for the first training parameter when s = 5, and (b) the size of the DMD modes for all training parameters when s = 5 and ρDMD = 1 × 10 -1 .

Figure 13 :

 13 Scattering of a plane wave by a multi-layer disk: some of the exact (-) and approximate (o) (HODMD method) reduced-order coefficient vectors α u,l (t, θ) (l = 1, 5, 10) of (a) Hx, (b) Hy and (c) Ez in T te ℏ for some training parameters. Different colors indicate different parameters. The dashed line indicates the end of the training time region, which is t = 49.7 m.

Figure 15 :

 15 Figure 15: Scattering of a plane wave by a multi-layer disk: comparison of the time evolution of the relative errors of E (left) and H (right) for the four testing parameters in the testing time region T te ℏ .

Figure 16 :

 16 Scattering of a plane wave by a dielectric sphere: some of the exact and approximate (HODMD method) reduced-order coefficient vectors α u,l (t, θ) (l = 1, 2, 4, 6) of (a) Ex, (b) Ey and (c) Ez in T te ℏ for all training parameters. The dashed line indicates the end of the training time region, which is t = 0.7 m.

Figure 17 :Figure 18 :

 1718 Figure 18: Scattering of a plane wave by a dielectric sphere: comparison of the time evolution of the fields Hx (top, left), Hy (middle, left), Hz (bottom, left), Ex (top, right), Ey (middle, right), Ez (bottom, right) at a given space point for the four testing parameters in the testing time region T te ℏ .

Figure 19 :

 19 Figure 19: Scattering of a plane wave by a dielectric sphere: comparison of the time evolution of the relative errors of E (left) and H (right) for the four testing parameters in the testing time region T te ℏ .

Figure 20 : 2 -

 202 Figure 20: 2-D Smolyak sparse grids versus a full tensor product grid.

Table 1 :

 1 Settings for the training and testing datasets.

	Case	Method	Data set	Training set	Testing set
			Parameter sampling 65, uniform	40, random (LHS)
	Disk (2-D)	DGTD	Time sampling	190, uniform	263, uniform
			Size	12350	10520
			Parameter sampling 137, uneven (Smolyak, 4-D, L = 3) 81, random (LHS)
	Multi-layer disk (2-D)	DGTD	Time sampling	183, uniform	253, uniform
			Size	25071	20493
			Parameter sampling 81, uniform	40, random (LHS)
	Sphere (3-D)	Mie series Time sampling	70, uniform	100, uniform
			Size	5670	4000

Table 2 :

 2 Comparison of the average relative projection and NIMOR errors on the testing set T te ℏ × P te ℏ .

	Case	e E,Pro	e E,NIMOR	e H,Pro	e H,NIMOR
	Disk	1.178 × 10 -2	1.348 × 10 -2	1.072 × 10 -2	1.388 × 10 -2
	Multi-layer disk	7.622 × 10 -3	7.690 × 10 -3	6.727 × 10 -3	6.775 × 10 -3
	Sphere	2.451 × 10 -3	5.908 × 10 -3	2.479 × 10 -3	5.269 × 10 -3

Table 3 :

 3 Comparison of the average computing time of the RB and HF solutions. The unit of time cost is second.

	Case	Offlin stage	Online stage
		(HF solutions, Nested POD, DMD training) (DMD prediction, RBF, one run for new paramter)
	Disk	2.765 × 10

+4

, 2.015 × 10 +1 , 2.920 × 10 +1 6.690 × 10 -1 , 5.089 × 10 -1 , 4.218 × 10 -2 ≪ 1 Multi-layer disk 6.184 × 10 +4 , 1.071 × 10 +2 , 5.224 × 10 +1 7.040 × 10 -1 , 5.313 × 10 -1 , 5.680 × 10 -2 ≪ 1 Sphere 2.541 × 10 +5 , 2.194 × 10 +2 , 6.900 × 10 +1 7.800 × 10 -1 , 8.066 × 10 -1 , 5.098 × 10 -2 ≪ 1

Table 4 :

 4 Distribution and range of material parameters.

  • • • , 6.9697 × 10 -1 m}. The testing parameter set P te ℏ with size 40 is chosen by LHS method, and the testing time set T te ℏ is set to be {0 m, 1.0101 × 10 -2 m, • • • , 1 m}. With the truncation errors ρ t = 1 × 10 -1 and ρ θ = 1 × 10 -5 , a

	100
	0
	-100
	49.1	49.2	49.3	49.4	49.5	49.6	49.7	49.8	49.9
				Time (m)
	10
	0
	-10
	49.1	49.2	49.3	49.4	49.5	49.6	49.7	49.8	49.9
				Time (m)
	2
	0
	-2
	49.1	49.2	49.3	49.4	49.5	49.6	49.7	49.8	49.9
				Time (m)
				(a)
	200
	0
	-200
	49.1	49.2	49.3	49.4	49.5	49.6	49.7	49.8	49.9
				Time (m)
	20
	0
	-20
	49.1	49.2	49.3	49.4	49.5	49.6	49.7	49.8	49.9
				Time (m)
	5
	0
	-5
	49.1	49.2	49.3	49.4	49.5	49.6	49.7	49.8	49.9
				Time (m)
				(b)
	200
	0
	-200
	49.1	49.2	49.3	49.4	49.5	49.6	49.7	49.8	49.9
				Time (m)
	20
	0
	-20
	49.1	49.2	49.3	49.4	49.5	49.6	49.7	49.8	49.9
				Time (m)
	49.1	49.2	49.3	49.4	49.5	49.6	49.7	49.8	49.9
				Time (m)

Acknowledgments

The authors are very grateful for the constructive advices from anonymous reviewers and the help of the editors.

Acknowledgements

This research was supported by the National Natural Science Foundation of China (Grant No. 12101511) and by the Fundamental Research Funds for the Central Universities (Grant No. JBK22YJ37).

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.


Yanshu Wang

Anne Sergent

Gilles Bernard-Michel

Patrick Le Quéré

Numerical assessment of entrainment coefficient of turbulent jet in case of hydrogen release in a two-vented cavity

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

NUMERICAL ASSESSMENT OF ENTRAINMENT COEFFICIENT OF TURBULENT JET IN CASE OF HYDROGEN RELEASE IN A TWO-VENTED CAVITY

ABSTRACT

Hydrogen is considered as one of the green and alternative energy sources due to its reduced emission of CO2. As its high risk of detonation, which is often associated with a mixture of hydrogen and air, appropriate security measures are required for its usage. In accidental situations, hydrogen escapes and rises in form of a plume-jet that entrains ambient air and potentially forms a dangerous flammable mixture. Predictive models have therefore to be developed to assess its potential use in specific configurations like highly confined environments (garage, containment building …). Macroscopic models as those proposed by Linden et al. (1990) based on Morton's turbulent jet model (Morton et al., 1956) for building applications have been used to predict hydrogen concentrations in vented cavities. The determination of entrainment coefficient 𝛼, which is supposed to be constant in the Morton et al.'s model, plays an important role in the estimation of hydrogen concentration.

We present numerical results from direct numerical simulations (DNS) of an air-hydrogen buoyant jet in a two-vented cavity, configuration dedicated to analysis of hydrogen distribution in accidental situations for highly confined environments. We compare the result of 3D numerical simulation with the simplified models proposed by Linden et al. (1990) based on Morton's turbulent jet model (Morton et al. 1956) used to predict hydrogen concentrations in vented cavities in industrial context. By integration of time-averaged fields obtained by DNS simulation, we associate the profiles of characteristic quantities of the jet with its predicted values, to study the choice of entrainment coefficient 𝛼 in Morton's model. The coefficient is determined by three methods with different values which gives indication about its value depended on different flow regime. By comparing the results, we find that 𝛼-constant models present restriction in application of simplified models used in industrial context.

Keywords: Hydrogen safety, DNS, turbulent jet, ventilation model, entrainment coefficient

  1,2* Yanshu Wang,2,3 Anne Sergent, 1 Gilles Bernard-Michel and 2 Patrick Le Quéré 1 DEN/DANS/DM2S/STMF/LIEFT, Commissariat à l'énergie atomique et aux énergies alternatives (CEA), Saclay, France 2 Laboratoire d'Informatique pour la Mécanique et les Sciences de l'Ingénieur (LIMSI), CNRS, UPR 3251, Orsay, France

	3 UFR Ingénierie, Sorbonne University, Paris, France
	*Corresponding author e-mail: yanshu.wang@cea.fr