On Fostering Predictions in Data-Driven Reduced Order Model for Eulerian–Lagrangian Simulations: Decision of Sufficient Training Data

Abstract

The development of a data-driven surrogate model (SM) is extensively studied in Eulerian–Lagrangian simulations for its advantage of high computational speed. However, in the application of granular systems with violent fluid-solid flows, how to select sufficient training data to ensure consistency between the high-fidelity model and SM remains unknown and highly challenging. The accuracy of SM can be easily deteriorated due to insufficient training data. This necessitates a trial-and-error process and hinders its industrial applications. To address this issue, this study newly reveals a finding that data density is a key to sufficient training, and we propose a novel technique for deciding the sufficient training data of SM. Specifically, a feasibility index is proposed based on posterior error analysis. It is demonstrated that when the training data is determined under the proposed feasibility index $\ge$ 2, the consistency of granular dynamics between SM and the high-fidelity model can be guaranteed. Employed in a representative SM, a reduced order model (ROM), this technique enables the successful decision of sufficient training data, resulting in the remarkable predictability in violent fluid-solid flows without trial-and-error. This technique holds great potential in solving the predicament of deciding training data for data-driven models.

Keywords:

• Data-driven ROM
• Fluid-solid flows
• Sufficient training data
• Predictions

1. Introduction

In chemical engineering, multi-phase flows including solid particles are often encountered, such as fluidizations (Mori et al. 2019; Sakai et al. 2010; Takabatake and Sakai 2020), conveying (Sun and Sakai 2015), mixing (Li et al. 2021; Mori and Sakai 2022; Tsugeno et al. 2021), milling (Liu et al. 2018; Ogi et al. 2017; Tanaka et al. 2019; Yamamoto et al. 2012) and die-filling (Widartiningsih et al. 2020; Yao et al. 2018). To study such multi-phase processes, coupled discrete element method–computational fluid dynamics (DEM–CFD) is widely employed as a reliable approach for the associated Eulerian–Lagrangian simulations. Advanced applications of the DEM–CFD method have been extensively reported including but not limited to particle separation in recycling biodegradable plastics (Tsunazawa et al. 2023), wall erosion caused by particle conveying in vortex elbows (Xiao et al. 2023), chip-like particles flow in fluidized beds of solar panel modules recycling (Wang and Shen 2022), a resolved model for gas-solid-liquid three phase flows with super-quadrics particles (Washino et al. 2023), and extended DEM-based modelling for wet granulation in a rotating drum (Nakamura et al. 2022).

However, when it comes to the industrial scale, the high computational cost associated with the total amount of particles and solving complex equations in the DEM–CFD method has hindered further applications (Jiang et al. 2021; Li et al. 2022a). It is pointed out that even with the latest devices of CPUs and GPUs, days of simulations are usually required for an industrial-scale granular system (Napolitano et al. 2022; Sakai 2016). In a sand production case, it is reported necessary to employ high-performance computers (HPCs) for computations (Kazidenov et al. 2023). Thus, there has been significant interest in simplifying such simulations.

Recently, with the advances in data-science, the data-driven surrogate model (SM) has been considered highly potential to solve the problem of heavy computational cost in Eulerian–Lagrangian simulations (Li et al. 2022a). A typical representative is a data-driven reduced order model (ROM). It is a surrogate model that simplifies the model mathematically based on the dynamics identification technique. It can achieve substantial dimensionality reductions while preserving essential dynamics to reconstruct the original solutions (Karhunen 1947; Taira et al. 2017; Trefethen and Bau 1997).

The thrusts in computational speed and the physical-equation-free nature of ROMs have opened up an innovative approach for Eulerian–Lagrangian simulations. It is of profound significance in fostering real-time simulations for digital twin in the digital transformation (DX) (Chakraborty et al. 2021; Geier et al. 2023; Hartmann et al. 2021, 2018; VanDerHorn and Mahadevan 2021).

The most common approach to obtain the dynamical features for ROM is proper orthogonal decomposition (POD) (Lumley 1967). It has been widely applied to and studied in Eulerian simulations (Allery et al. 2008; Xiao et al. 2015; Xiao 2019). However, only limited studies on Lagrangian simulations have been reported (e.g., Lu and Tartakovsky 2020) due to the accuracy concern caused by complicated motion of particles. One alternative in studies of disperse phase is to map Lagrangian information onto Eulerian meshes (Higham et al. 2020; Pang and Wei 2013), but applicable parameters are limited and extra computational cost and error can be introduced in large-scale and complicated systems. Empowered by the sharply increased computational power, new approaches start to consider Lagrangian variables in ROM directly (Li et al. 2022a; Lu and Tartakovsky 2020; Schiødt et al. 2022). It has been validated efficient using Lanczos POD in mixing phenomena quantifications (Li et al. 2022b,2022c) and ROM reproductions of Eulerian–Lagrangian simulations. The implementation of Lanczos iteration (Fahl 2001; Li et al. 2022a) improves the computational efficiency drastically and makes POD more suitable for large-scale cases.

However, like many other data-driven approaches, ROMs also suffer the lack of training data problems, which is commonly regarded as a fatal issue in predictions (Brunton and Kutz 2019; ; Jiang et al. 2023). Especially, when the solid-fluid flows in Eulerian–Lagrangian simulations become highly fierce and complicated (e.g., in a bead mill), the assumption of consistent dynamics between ROMs and high-fidelity models (Kherad et al. 2020; Lucia et al. 2004) may not be guaranteed due to the lack of training data in dynamics identifications (e.g., POD), and the accuracy of ROMs can be significantly deteriorated even with the previously suggested local POD calculated on each parameter (Hajisharifi et al. 2023; Li et al. 2022a), particularly in the particle phase as shown in Figure 1. This is caused by the highly random and rapid motion of particles. Moreover, the random and non-linear behaviors of all particles make it rather challenging to select the training data for particle motions with significant different characteristics.

Figure 1. (Left) Comparisons of Dynamics and Typical snapshot of High-fidelity Model and ROM with insufficient training data; (Right) Snapshots of Particle & Fluid Phase at T = 2.0s in High-fidelity Model and ROM with insufficient training data.

To address the inaccurate ROMs problem and maximize the thrust in the reduction of computational cost, this research first reveals that the density of training data is crucial in terms of the temporal interpolations in ROM. Although more computational resource will be required, it can significantly improve the consistency of particle dynamics between ROM and high-fidelity model, thereby enhance the accuracy of ROM. Based on this finding, the aim is to propose a novel technique for the decision of sufficient training data in data-driven ROM for violent fluid–solid flows simulations. First, two bead mill cases are simulated by an in-house validated high-fidelity DEM–CFD code, FELMI (Mori and Sakai 2021), with different operational conditions for training data preparations. Second, a more violent case is selected as the representative of a violent granular system to conduct ROM under different selections (amount and density) of training data. Subsequently, a new feasibility index for the decision of sufficient training data is proposed and proven adequate by comparing ROM results and POD coefficients. Finally, the validation on the other case demonstrates successful reproduction and prediction when sufficient training data is determined with the proposed feasibility index $\ge$ 2.

As a result, the technique for the decision of sufficient training data in data-driven ROMs for Eulerian–Lagrangian simulations is established and validated in the considered bead mill cases. This method is designed to determine the density of training data, it is applicable to any kind of Eulerian–Lagrangian simulations including pseudo-periodic and non-periodic flows. In addition, it holds remarkable value in demonstrating the sophisticated balance between efficiency and accuracy for ROMs. Moreover, this technique guarantees consistent dynamics fundamentally with a substantially reduced computational cost, thereby high-speed & accurate ROM predictions become achievable. Notably, not only for ROMs, this method can be universally extended for evaluating the predictive accuracy of any data-driven models based on dynamics identifications.

2. Methodology

In this study, the employed high-fidelity model for the DEM–CFD method is an in-house code FELMI, and the employed data-driven ROM is LPOD-based ROM. This section provides the fundamentals of FELMI, LPOD-based ROM and the proposed technique for the decision of sufficient training data.

2.1. FELMI: DEM–CFD

In this study, the DEM–CFD method is conducted using FELMI (Mori and Sakai 2021). In detail, this method uses the Signed Distance Function (SDF), the Immersed Boundary Method (IBM), and the implicit algorithm for the drag force to couple the DEM and CFD simulations to consider all phases simultaneously. The specific numerical modeling is as follows.

2.1.1. Solid (granular) phase

The translational and rotational motion of particles (the solid phase) are described as the two following governing equations: (1) $$m\frac{d\mathbf{v}}{\mathit{\text{dt}}}=\sum {\mathbf{F}}_c+{\mathbf{F}}_d-V_p\nabla p+m\mathbf{g},$$(1) and (2) $$I\frac{d\mathit{\omega }}{\mathit{\text{dt}}}=\sum \mathbf{T},$$(2) where m, v, t, ${\mathbf{F}}_c,{\mathbf{F}}_d,$ V_p, p, g, I, $\mathit{\omega },$ and T represent the mass of the particle, the velocity of the particle, time, contact force, fluid drag force, volume of the particle, pressure, gravitational acceleration, momentum inertia, the rotational velocity of the particle, and torque, respectively.

The contact force acting on the particle is described in the normal and tangential directions: (3) $${\mathbf{F}}_c={\mathbf{F}}_{c_n}+{\mathbf{F}}_{c_t}.$$(3)

In detail, ${\mathbf{F}}_{c_n}$ and ${\mathbf{F}}_{c_t}$ are calculated by, (4) $$\begin{gathered} {\mathbf{F}}_{c_{\mathit{n}}}=-k{\mathbf{\delta }}_{\mathrm{n}}-\eta {\mathbf{v}}_{{\mathrm{r}}_{\mathrm{n}}}, \\ {\mathbf{F}}_{c_t}=\{\begin{array}{cc}-k{\mathit{\delta }}_t-\eta {\mathbf{v}}_{r_t},& \left|{\mathbf{F}}_{c_t}\right|<{\mu }_p\left|{\mathbf{F}}_{c_n}\right|\\ -{\mu }_p\left|{\mathbf{F}}_{c_n}\right|\frac{{\mathbf{v}}_{r_t}}{\left|{\mathbf{v}}_{r_t}\right|},& \left|{\mathbf{F}}_{c_t}\right|\ge {\mu }_p\left|{\mathbf{F}}_{c_n}\right|\end{array}, \end{gathered}$$(4) where ${\mathbf{v}}_{r_t},{\mathbf{v}}_{r_n},{\mathit{\delta }}_t,{\mathit{\delta }}_n,$ k, η, and μ_p are the relative velocity of the particle in tangential direction, the relative velocity of the particle in normal direction, the displacement of the particle in tangential direction, the displacement of the particle in normal direction, the spring constant, the damping coefficient, and the friction coefficient, respectively.

Moreover, η is given by calculating the energy dissipation as the following equation, (5) $$\eta =-2\hspace{0.17em}\mathrm{\text{ln}}\hspace{0.17em}(e)\sqrt{\frac{\mathit{\text{km}}}{{\hspace{0.17em}\mathrm{\text{ln}}\hspace{0.17em}}^2(e)+{\pi }^2}},$$(5) where e is the restitution coefficient. The drag force acting on the particle is given by, (6) $${\mathbf{F}}_d=\frac{\beta V_p}{1-\epsilon }({\mathbf{u}}_f-\mathbf{v}),$$(6) where β, ε, and ${\mathbf{u}}_f$ are the interphase momentum transfer coefficient, void fraction, and fluid velocity, respectively. Similar to the previous solid–fluid coupling simulations Sun et al. (2013); Sakai et al. (2010, 2012, 2014); Sun and Sakai (2015) a combination of equations Ergun and Orning (1949); Ergun (1952); Wen (1966) was employed for coefficient β as follows, (7) $$\beta =\{\begin{array}{cc}\frac{150{\mu }_f{(1-\epsilon )}^2}{\epsilon d^2}+1.75\rho \frac{(1-\epsilon )}{d}|{\mathbf{u}}_f-\mathbf{v}|,& \epsilon \le 0.8\\ \frac{3}{4}\frac{C_D|{\mathbf{u}}_f-\mathbf{v}|\epsilon (1-\epsilon )}{d}{\epsilon }^{-2.65},& \epsilon >0.8\end{array},$$(7) (8) $$C_D=\{\begin{array}{cc}\frac{24(1+0.15R{e}_p^{0.687})}{R{e}_p},& R{e}_p\le 1000\\ 0.44,& R{e}_p>1000\end{array},$$(8) and (9) $$R{e}_p=\frac{|{\mathbf{u}}_f-\mathbf{v}|\rho \epsilon d}{{\mu }_f},$$(9) where ρ, μ_f, d, C_D, and Re_p represent the fluid density, fluid viscosity, particle diameter, drag coefficient, and Reynolds number of the particle, respectively.

2.1.2. Fluid phase

The governing equations for fluid phase in FELMI are composed by 3-D Navier–Stokes equation and the continuity equation as follows with implementation of the local volume average method (Anderson and Jackson 1967), (10) $$\frac{\partial (\epsilon \rho {\mathbf{u}}_f)}{\partial t}+\nabla ·(\epsilon \rho {\mathbf{u}}_f{\mathbf{u}}_f)=-\epsilon \nabla p+f_s+\nabla ·(\epsilon \mathit{\tau })+\epsilon \rho \mathbf{g},$$(10) and (11) $$\frac{\partial \epsilon }{\partial t}+\nabla ·(\epsilon {\mathbf{u}}_f)=0,$$(11) where ${\mathbf{f}}_s,$ and $\mathit{\tau }$ indicate the drag force acting on the fluid, and viscous stress, respectively. The drag force in Eq. (10)(10) $$\frac{\partial (\epsilon \rho {\mathbf{u}}_f)}{\partial t}+\nabla ·(\epsilon \rho {\mathbf{u}}_f{\mathbf{u}}_f)=-\epsilon \nabla p+f_s+\nabla ·(\epsilon \mathit{\tau })+\epsilon \rho \mathbf{g},$$(10) is calculated using Eq. (6)(6) $${\mathbf{F}}_d=\frac{\beta V_p}{1-\epsilon }({\mathbf{u}}_f-\mathbf{v}),$$(6) as follows, (12) $${\mathbf{f}}_s=-\frac{\sum {\mathbf{f}}_d}{V_{\text{grid}}},$$(12) where $V_{\text{grid}}$ is the volume of the fluid grid.

2.1.3. Wall boundary

The wall boundary in FELMI is given by the combination of Signed Distance Function (SDF) (Shigeto and Sakai 2013) and Immersed Boundary Method (IBM) (Sun and Sakai 2017, 2016), which is proven accurate and robust in a much simpler manner (Mori et al. 2019; Tsugeno et al. 2021; Mori and Sakai 2021; Yao et al. 2018).

2.1.3.1. Solid phase

The wall boundary of the solid (granular) phase is modeled by the SDF as the following equation, (13) $$\varphi (\mathbf{x})=s(\mathbf{x})d(\mathbf{x}),$$(13) where x, $s(\mathbf{x}),$ and $d(\mathbf{x})$ are the position vector of the particle, the sign based on the position of the particle, and the distance between the particle and the surface of the wall boundary, respectively. For points inside of the calculation domain, the s(x) is positive, and vice versa. The contact force between the particle and the wall is given as, (14) $${\mathbf{F}}_{{\mathrm{c}}_n}=-k_{\mathit{\text{pw}}}|\nabla \varphi |{\mathbf{\delta }}_n^{\mathit{\text{SDF}}}-{\eta }_{\mathit{\text{pw}}}{\mathbf{v}}_{r_n},$$(14) and (15) $${\mathbf{F}}_{c_t}=\{\begin{array}{cc}-k_{\mathit{\text{pw}}}{\mathbf{\delta }}_t-{\eta }_{\mathit{\text{pw}}}{\mathbf{v}}_{r_t}\mathbf{ }\mathbf{,}& \left|{\mathbf{F}}_{c_t}\right|<{\mu }_{\mathit{\text{pw}}}\left|{\mathbf{F}}_{c_n}\right|\\ -{\mu }_{\mathit{\text{pw}}}\left|{\mathbf{F}}_{c_n}\right|\frac{{\mathbf{v}}_{r_t}}{\left|{\mathbf{v}}_{r_t}\right|},& \left|{\mathbf{F}}_{c_t}\right|\ge {\mu }_{\mathit{\text{pw}}}\left|{\mathbf{F}}_{c_n}\right|\end{array},$$(15) where the subscript pw represents the particle-wall physical properties.

While the overlap considered in SDF for the soft particle model, ${\mathit{\delta }}_n^{\mathit{\text{SDF}}}$ is defined as, (16) $${\mathit{\delta }}_n^{\mathit{\text{SDF}}}=(\varphi -\frac{D}{2}){\mathbf{n}}^{\mathit{\text{SDF}}},$$(16) and (17) $${\mathbf{n}}^{\mathit{\text{SDF}}}=\frac{\nabla \varphi }{|\nabla \varphi |}.$$(17)

Therefore, elastic energy is conserved in SDF.

2.1.3.2. Fluid phase

IBM is introduced to directly simulate the interaction between the solid wall and the fluid. In IBM, we first calculate the volume–weighted average velocity, u, by the following equation, (18) $$\mathbf{u}=(1-\alpha ){\mathbf{u}}_f+\alpha {\mathbf{u}}_{\mathit{\text{wb}}},$$(18) where α and ${\mathbf{u}}_{\mathit{\text{wb}}}$ are the local volume fraction calculated by SDF and the wall boundary velocity, respectively. In what’s after, the 3-D Navier-Stokes equation can be corrected as follows, (19) $$\frac{\partial (\epsilon \rho \mathbf{u})}{\partial t}+\nabla ·(\epsilon \rho \mathbf{\text{uu}})=-\epsilon \nabla p+{\mathbf{f}}_s+\nabla ·(\epsilon \mathit{\tau })+\epsilon \rho \mathbf{g}+\epsilon {\mathbf{f}}_{\mathit{\text{IB}}},$$(19) where ${\mathbf{f}}_{\mathit{\text{IB}}}$ is given by, (20) $${\mathbf{f}}_{\mathit{\text{IB}}}=\alpha \rho \frac{{\mathbf{u}}_{\mathit{\text{wb}}}-\widehat{\mathbf{u}}}{\Delta t},$$(20)

The above provides necessary information of this method, while the repeated conductions of SDF/IBM boundary conditions are out of the cope of this paper. For interested readers, more details of our original wall boundary conditions including SDF and IBM can be reached in the related papers mentioned above.

2.2. LPOD-based ROM

LPOD-based ROM conducts the simulations by (1) calculating the POD modes of different variables in each case by collecting a series of snapshots and solving the eigenvalue problem; (2) interpolating the corresponding coefficients temporally; (3) projecting the coefficients back for solution generation. The implementation of Lanczos iteration assists with the computational efficiency in POD calculations by allowing the calculation of POD modes using a subset of the total snapshots.

A series of snapshots from the results of high-fidelity model DEM–CFD or experiments (here, the FELMI is employed) is taken and written as ${\{\mathbf{W}(t_i)\}}_{i=1}^{N_t},$ where W can be any Eulerian or Lagrangian variable, $t_i\in {\left\{t_i\right\}}_{i=1}^{N_t}$ refers to the time of the snapshot(s), N_t is the total time number, and $\{·\}$ represents a set of elements. (21) $$\mathbf{W}(t_i)=\overline{\mathbf{W}}+\tilde{\mathbf{W}}(t_i),$$(21) (22) $$\overline{\mathbf{W}}=\frac{1}{N_t}\sum _{i=1}^{N_t}\mathbf{W}(t_i),$$(22) where $\overline{\mathbf{W}}$ is the mean part and $\tilde{\mathbf{W}}$ is the time-varying part. $\tilde{\mathbf{W}}(t_i)$ can be derived through solving the following minimization problem, in which ${\left\{{\mathit{\phi }}_i\right\}}_{i=1}^{N_t}$ is a set of orthogonal basis. (23) $$\begin{array}{c}\underset{{\left\{{\mathit{\phi }}_j\right\}}_{j=1}^d}{\mathrm{\text{min}}}\frac{1}{N_t}\sum _{i=1}^{N_t}{\Vert \tilde{\mathbf{W}}(t_i)-\sum _{j=1}^d\left(\tilde{\mathbf{W}}(t_i){\mathit{\phi }}_j\right){\mathit{\phi }}_j\Vert }^2,\\ s.t.({\mathit{\phi }}_i,{\mathit{\phi }}_j)={\mathit{\delta }}_{\mathit{\text{ij}}}\hspace{1em}(1\le i\le d,1\le j\le i),\end{array}$$(23) where $\left|\right|·\left|\right|$ represents the $L^2$ norm $\left|\right|f|{|}^2=(f,f).$ Given the orthogonality of ${\mathit{\phi }}_j$ ($j=1,2,\dots ,N$), the problem can be transformed into the following eigenvalue problem by introducing a Lagrangian multiplier, (24) $$\mathbf{A}\mathit{\phi }=\lambda \mathit{\phi },\mathbf{A}=\frac{1}{N_t}\sum _{i=1}^{N_t}\tilde{\mathbf{W}}(t_i){\left(\tilde{\mathbf{W}}(t_i)\right)}^T,$$(24) where, A is a covariance matrix, and $\mathit{\phi }$ is the POD modes, representing the orthogonal basis vectors, and λ is the eigenvalues, also called the proper orthogonal values (POVs), representing the energy level of the corresponding POD modes. Hereafter, the original problem has been converted to an eigenvalue problem, which is solved by an efficient LPOD technique.

For better computational efficiency in POD, Lanczos iteration is implemented as LPOD in this data-driven model. A reduced set of Krylov vectors can be conducted by exploiting the tridiagonalization of A as follows, (25) $${\mathbf{Q}}_k^T\mathbf{A}{\mathbf{Q}}_k={\mathbf{T}}_k,$$(25) (26) $${\mathbf{T}}_k=\left[\begin{array}{cc}\begin{array}{cc}{\alpha }_1& {\beta }_1\\ {\beta }_1& {\alpha }_2\\ ⋮& {\beta }_2\end{array}& \begin{array}{cc}\cdots & 0\\ {\beta }_2& 0\\ \ddots & ⋮\end{array}\\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}\cdots & {\beta }_{k-1}\\ {\beta }_{k-1}& {\alpha }_k\end{array}\end{array}\right],$$(26) (27) $${\mathbf{Q}}_k=[{\mathbf{q}}_1,{\mathbf{q}}_2,\dots ,{\mathbf{q}}_k],$$(27) where subscript k is the Lanczos-step number. A sufficiently large k can ensure the convergence of the Lanczos iterations (Li et al. 2022a). Therefore, by calculating the tridiagonalization of A, the expected eigenvalues/eigenvectors can be approximated mathematically using only a portion of snapshots. The convergence behavior of LPOD is shown in Figure 2.

Figure 2. LPOD convergence behavior of various variables.

Similar to a normal POD, the captured energy proportion (i.e., the Lanczos ratio) can be calculated as follows, (28) $$I_k(p)=\frac{\sum _{i=1}^p{\lambda }_i\prime }{\sum _{i=1}^N{\lambda }_i\prime }\approx \frac{1}{I(k)}·I(p),$$(28) where ${\lambda }_i^{\prime }$ is the i-th eigenvalue from LPOD. By choosing the Lanczos ratio, the employed modes can be determined according to the needs of computational cost reduction and accuracy. This employed modes number is called the least reconstructable mode, n_lr. It should be set as the lower limit of employed modes to guarantee reconstructability.

Given LPOD, a considered variable x at a specific time t_n can be expressed as the following approximation, (29) $$\mathbf{W}(\mathbf{x},t_n)=\overline{\mathbf{W}}(\mathbf{x})+\sum _{i=1}^N{c}_i(t_n){\mathit{\phi }}_i,$$(29) where N is the considered number of POD modes in ROM, $c_i(t_n),(i\in \{1,2,\dots ,N\}),$ represents the i-th coefficient of the corresponding POD mode at t_n. A non-intrusive ROM is developed to interpolate the coefficients as follows, (30) $$c_i(t_n)=f_i(\mathbf{c}(t_{n-1})),\hspace{1em}n\in \{1,2,\dots ,N_t\}.$$(30)

Here, f_i represents a set of interpolation functions that interpolate c_i temporally. In this study, radial basis function (RBF) is employed for the interpolation as developed and validated in previous studies (Hajisharifi et al. 2023; Li et al. 2022a). Meanwhile, other neural-network approaches can also be used for manipulating the coefficients in the reduced spaces in place of f_i.

Herein, essential fundamentals of the introduced LPOD-based ROM method are provided.

2.3. Technique for the decision of sufficient training data

To develop a technique for decision of sufficient training data, it is necessary to guarantee consistent dynamics by only using POD but not ROM so that sufficient training data can be decided without trial-and-error.

In this section, a technique is established quantitatively as the following two steps: (1) error analysis, consistency ratio for evaluating the consistency between captured dynamics and real system dynamics; (2) a feasibility index for the performance anticipation.

2.3.1. Error analysis and consistency ratio

Errors in POD-based ROM can be simply categorized as truncation error, coefficient error, and consistency error (Li et al. 2022a; Xiao et al. 2017).

Briefly, the truncation error is caused by employing insufficient POD modes. It can be reduced effectively by increasing employed POD modes (Luo et al. 2007; Xiao et al. 2015). While the coefficient error implies the error caused by the prediction of the coefficient using a regression technique (e.g., RBF or LSTM network). This error has been studied in different kinds of ROMs and its bounds has been derived (Homescu et al. 2005; Luo et al. 2007; Wu and Schaback 1993). By introducing a weight for ROM reconstructions, this error can be compensated properly (Li et al. 2022a; Xiao et al. 2015). Other probable ways to reduce this error can refer to using an appropriate model or increasing model complexity (e.g., using a multi-layers neural network). The source of consistency error refers to the inconsistency between captured dynamics and the real underlying dynamics in the system (; Duan et al. 2024; Kherad et al. 2020; Lucia et al. 2004), which is the main error to be analyzed in this research. It has been proved and demonstrated as the main error that causes the poor predictive ability of ROM (Duan et al. 2024). Due to the linear assumption in POD (Lucia et al. 2004), relatively sparse training data can significantly deteriorate the consistency between the captured dynamics and the real dynamics. Therefore, the sufficient training data is decided by both the amount and the density. In this study, within the same time range, $\Delta \tau$ (interval between training data snapshots) is used to control the amount and density of training data. A smaller sampling interval, $\Delta \tau ,$ indicates the training data with a larger amount and a higher density. Through this interval, an effective operation to decide the sufficient training data can be conducted.

To avoid any trial-and-error and potential coefficient error, a systematic error analysis posterior to POD but priori to ROM is preferable. After the calculation of POD, the dynamics are obtained and then projected directly to the training set and testing set to generate the projected solutions. The training set refers to the same data used for POD calculation, while the testing set refers to the much denser data covering the same time range (in this study, $\Delta \tau$ of testing set is set to 0.001 s). With POD calculations, they can be described as follows, (31) $${\mathbf{W}}_{\text{train}}(\mathbf{x},t_n)=\overline{\mathbf{W}}(\mathbf{x})+\sum _{i=1}^N{c}_{i,\text{train}}(t_n){\mathit{\phi }}_i,$$(31) (32) $${\mathbf{W}}_{\text{test}}(\mathbf{x},t_n)=\overline{\mathbf{W}}(\mathbf{x})+\sum _{i=1}^N{c}_{i,\text{test}}(t_n){\mathit{\phi }}_i,$$(32) where ${\mathbf{W}}_{\text{train}}(\mathbf{x},t_n)$ and ${\mathbf{W}}_{\text{test}}(\mathbf{x},t_n)$ are the projected solutions of the training set and testing set, respectively. The coefficients $c_{i,\text{train}}(t_n)$ and $c_{i,\text{test}}(t_n)$ represent the POD coefficients calculated directly using the same POD on training and testing data. Subsequently, the error on each snapshot is calculated using root mean square error (RMSE) between the projected solutions and the original solutions from the high-fidelity model. (33) $$E_{\text{train}}(t_n)=\mathit{\text{RMSE}}({\mathbf{W}}_{\text{train}}(\mathbf{x},t_n),{\mathbf{W}}_{\text{original}\_\text{train}}(\mathbf{x},t_n)),$$(33) (34) $$E_{\text{test}}(t_n)=\mathit{\text{RMSE}}({\mathbf{W}}_{\text{test}}(\mathbf{x},t_n),{\mathbf{W}}_{\text{original}\_\text{test}}(\mathbf{x},t_n)),$$(34)

With the consistency errors above, a consistency ratio is defined as follows to evaluate the consistency between the dynamics of the training set and testing set, (35) $$R_c(n)=\frac{\text{Mean}(E_{\text{test}})}{\text{Mean}(E_{\text{train}})}.$$(35)

By calculating R_c with different employed POD modes n, the consistent mode number that guarantees the consistency between the dynamics of the testing set and training set can be found. The increase of n can reduce the error in the testing set and training set simultaneously so that R_c is around 1 when the consistent dynamics hold. However, when it doesn’t hold, the increase of n can trigger the increase of R_c by increasing $E_{\text{test}}.$ This indicates that the consistency error gradually becomes dominant. In other words, the error between ROM predictions and high-fidelity simulations is increased. Note that the testing error is always higher than the training error. With sufficient snapshots considered in LPOD (e.g., 80% of total data), the value of $R_c=1.05$ is the threshold to guarantee consistency for successful reconstructions of physical variables or observations, where n can be called the consistent mode, n_c. The modes lower than n_c contribute equivalently to both the training set and testing set, while the modes higher than n_c contribute only to the training set but not the testing set. Consequently, n_c should be set as the upper limit of n to guarantee consistent dynamics.

Hereby, the conditions for consistent dynamics and successful reconstructions in ROM can be clarified as,

1. n is higher than n_lr and lower than n_c;

2. $n_c\ge n_{\mathit{\text{lr}}},$ consistent dynamics can be obtained with guaranteed reconstructability.

2.3.2. Feasibility index: Deciding sufficient training data

Because the training data is not a variable in POD, a feasibility index is presented to evaluate the feasibility of ROM with different amount/density of training data. The feasibility index is defined as follows, (36) $$K_f=\frac{n_c}{n_{\mathit{\text{lr}}}},$$(36)

This index can indicate the predictive accuracy of data-driven ROM under a set of input training data. Through the aforementioned error analysis, a suitable value of K_f indicating accurate ROMs can be derived, so that sufficient training data can be decided. With the derivation of K_f, the reconstructability, consistent dynamics and sufficient training data in POD can be all guaranteed simultaneously.

Hereby, the technique for the decision of sufficient training data can be summarized as follows,

1. Given a time range of snapshots, for $\Delta \tau$ = 0.01, 0.009, 0.008…0.002 s, conduct LPOD and get n_lr for reconstructability; (the range and step of $\Delta \tau$ can be roughly determined by ensuring a basic visual continuity in the considered systems, in this case 0.01 as an upper limit $\Delta {\tau }_{\mathit{\text{up}}},$ 0.001 as the lower limit $\Delta {\tau }_{\mathit{\text{low}}}$ and step $\Delta {\tau }_{\mathit{\text{step}}}$)

2. Calculate $E_{\text{train}},E_{\text{test}}$ and R_c to get n_c for consistent dynamics;

3. Calculate K_f for feasibility. When $K_f\ge 2.0,\Delta \tau$ is selected for the decision of sufficient training data and break.

The predictions on all variables between two adjacent output moments using data-driven ROM with sufficient training data can be achieved as illustrated in Figure 3.

Figure 3. Idea diagram of reproductions and predictions using LPOD-based ROM with sufficient training data.

3. Results and Discussion

In this chapter, two bead mill simulations are first conducted through FELMI. Second, LPOD-based ROMs are conducted on Case 1 with different training data. Subsequently, the criterion for the decision of sufficient training data is proposed and proven adequate by comparing ROM results and POD coefficients. Finally, the validation and predictions using LPOD-based ROM with sufficient training data in Case 2 is presented.

3.1. Bead mill simulations

The fluid & particles properties and calculation condition used in FELMI are shown in Table 1, the reliability of coefficient values have been studied thoroughly for the bead mill (Tanaka et al. 2019). In this study, two bead mill cases with 500,000 particles are considered, the dimensions are shown in Figure 4, these cases’ setup is listed in Table 2. One typical case, i.e. Case 1, is considered for determining an appropriate value of the feasibility index. After this value is found, Case 2 is used to validate the adequacy of this index. Table 3 lists the $\Delta \tau$ considered in trial-and-error searching of sufficient training data.

Figure 4. Dimensions of considered bead mill.

Table 1. Fluid & particle properties and calculation conditions.

Fluid
Viscosity	Variable
Density	1000 kg/m^3
Grid	67 × 95 × 67 mm^3
Particles
Amount	500,000
Density	6000 kg/m^3
Diameter	0.6 mm
Restitution coefficient	0.9
Friction coefficient	0.3
Computation
Time step	${10}^{-6}$ s
Output interval	${10}^{-3}$ s
Total time range	5 s

Table 2. Cases setup.

	Case 1	Case2
Viscosity	1 mPa $·$ s	20 mPa $·$ s
Rotational Speed	2500 rpm	1500 rpm

Table 3. Considered $\Delta \tau$ in trial-and-error.

$\Delta \tau$

0.01

0.008

0.004

0.002

3.2. Qualitative comparisons with different $\Delta \tau$

The qualitative comparisons of cases with different $\Delta \tau$ are provided through typical snapshots and POD coefficients as shown in Figures 5 and 6. They indicates qualitatively that the decrease of training data density leads to better consistent dynamics and reconstructions. Table 4 lists the training data and speed comparisons of cases with different $\Delta \tau .$ It is clear that the accuracy will increase when $\Delta \tau$ is small enough, however, the amount of training data and associated computational cost will increase as well.

Figure 5. Typical snapshots of particle phase in FELMI and ROM with different $\Delta \tau .$

Figure 6. 1st POD coefficients comparison (particle coordinate in X-direction).

Table 4. Training data and one-step speed comparison of FELMI and ROM with different $\Delta \tau$ (Sequantial on AMD Ryzen^TM Threadripper^TM Pro 5995WX).

$\Delta \tau$	Training Data	Model	DEM (s)	Reduction (*)
0.01	[-]	FELMI	10,738.138	[-]
0.01	500	ROM	0.892	12,038.271
0.008	[-]	FELMI	8647.321	[-]
0.008	625	ROM	0.943	9170.012
0.004	[-]	FELMI	4415.665	[-]
0.004	1250	ROM	2.105	2097.703
0.002	[-]	FELMI	2321.547	[-]
0.002	2500	ROM	8.948	259.460

3.3 Analysis on K_f and decision of sufficient training data

In the previous study, it was reported that a Lanczos ratio higher than 0.99 was deemed sufficient for both the fluid and particle phases (Li et al. 2022a). However, due to the violence in the bead mill, our preliminary research suggests that a Lanczos ratio greater than 0.999 (the corresponding mode number is considered n_lr) is required specifically for the particle phase. Nonetheless, even with a Lanczos ratio close to 1, it cannot guarantee a successful reconstruction in the particle phase when the training data is not dense enough.

In this section, the proposed technique is applied to Case 1 to analyse K_f for an accurate particle phase of ROM in the bead mill. Importantly, it must be conducted without ROM calculation so that the decision the sufficient training data can be made in advance without any trial-and error of ROMs. Notably, the results mainly show the calculations on particle coordinate in the X-direction, for its Lanczos-ratio convergence is slower than Y-direction (particles hold more complicated motions in the X-direction), and is similar to Z-direction due to geometrical symmetry.

Figure 7 shows the error calculated by direct projections with all POD modes employed. Note that without conducting ROM, the errors are dominated by consistency error. As the decrease of $\Delta \tau$ shrinks the band-shaped area formed by testing error, it is indicated that the captured dynamics are becoming similar to the real dynamics in the bead mill system. Also, the decrease of particle position error implies the decrease of unphysical contacts such as the particle overlaps or collision energy.

Figure 7. History of training and testing error with different $\Delta \tau$ (Particle coordinate in X-direction).

The consistency ratio is calculated as shown in Figure 8. With threshold = 1.05, n_c can be obtained.

Figure 8. Consistency ratio of particle coordinate in X-direction in cases with different $\Delta \tau .$

With the calculation of I_k and $R_c(n),$ n_lr and n_c can be determined. Subsequently, K_f can be calculated. Table 5 lists the values of K_f of cases with different $\Delta \tau .$ With examining the contours, it is suggested that a $K_f\ge 2.0$ is sufficient to guarantee reconstructability, consistent dynamics simultaneously, thereby the sufficient training data in POD can be decided.

Table 5. Feasibility index K_f for particle coordinate in X-direction with different $\Delta \tau .$

$\Delta \tau$	nlr	nc	Kf
0.01	366	317	0.8661
0.008	433	358	0.8267
0.004	582	764	1.3127
0.002	603	1322	2.1924

3.4. Validation of the proposed technique

Case 2 is calculated using LPOD-based ROM with sufficient training data for validations. With the guaranteed consistent dynamics under $K_f\ge 2.0,$ the reproductions and predictions can be realized simultaneously. Table 6 lists the related calculations. Notably, sufficient training data is decided prior by only using projections and the condition $K_f\ge 2.0$ without trial-and-error.

Table 6. Calculations for decision of sufficient training data.

Variable	Value
Selected $\Delta \tau$	0.004 s
Time range	0–5 s
Training Data (snapshots)	1250
Lanczos-step	1000
Lanczos-ratio	>0.999
nlr	323
nc	728
Kf	2.25
Actual Employed modes (Kf = 2)	646

Figure 9 shows the comparisons of FELMI and ROM solutions (reproductions and predictions) ($\Delta \tau$ is 0.004 s in LPOD, 0.001 s in ROM solutions). Remarkably, given the temporal interpolation in POD, the prediction can be done for any time within one $\Delta \tau .$ For instance, from 2.000 s to 2.004 s, only 2.000 and 2.004 are considered in LPOD, while the predictions on 2.001, 2.002, 2.003 are all achievable in ROM. Furthermore, as shown in Table 7, over 1500 times’ reduction of computational time is achieved using LPOD-based ROM with sufficient training data.

Figure 9. Reproductions and predictions of LPOD-based ROM using sufficient training data (Green box Rep.: Reproduction; Red box Pred.: Prediction).

Table 7. One-step speed comparison of FELMI and ROM with sufficient training data (Sequantial on AMD RyzenTM ThreadripperTM Pro 3995WX).

Model	DEM (s)	Reduction (*)
FELMI	1180.865	[-]
ROM	0.752	1570.299

* The online computational cost of trained ROM excluding the offline training time is only used in speed comparison.

3.5. Discussion

As the challenge of accuracy is commonly encountered in various data-driven models for both reproductions and predictions, the proposed technique offers a new approach to address the problem focusing on sufficing training data.

This technique for the decision of sufficient training data is developed without trial-and-error to guarantee consistency between captured dynamics and actual dynamics, it has the potential to be universally applied to other data-driven models. By the design based on error analysis on the training set and testing set, this technique can also be a reference of a framework for extended quantitative study on other data-driven models’ predictive performance.

Although for the considered bead mill case, a value of 2.0 is recommended for K_f to primarily focus on solving the inaccurate reconstructions in the particle phase, it is worth pointing out that different feasibility index can be applied to different variables to save more computational resource, as well as different consistency ratios, R_c can be used for different requirements of accuracy. Accordingly, a multi-timescale ROM is expected to be developed in future work to utilize the proposed technique for DX of systems using Eulerian–Lagrangian simulations.

4. Conclusions

In this study, it is firstly found that the increase of training data density has potential to solve the concern of accuracy in ROM on Lagrangian variables. Specifically, the sufficient training data for accurate data-driven ROMs needs to be decided from both data density and amount.

Therefore, a novel technique for decision of sufficient training data is established to solve the inaccurate reconstructions in particle phase. By using this technique, the sufficient training data can be decided without trial-and-error. Meanwhile, it can evaluate ROM in a temporal perspective to keep a sophisticated balance between efficiency and accuracy.

Notably, by employing this technique, the reproductions and predictions can both be achieved with a substantially reduced computational cost. Furthermore, since this technique is established through a mathematical error estimate, it can be potentially employed to other data-driven models.

Apart from solving the inaccurate reconstruction in the particle phase, it is suggested that the different sufficient training data (density and amount) can vary for different variables. Thus, by tailoring training data for each variable, towards the real-time simulations, a multi-timescale data-driven model for Eulerian–Lagrangian simulations is expected to be established in the future. Also, it is worth developing further advanced machine-learning techniques for ROM to extend the generalization capability.

Disclosure statement

The authors declare that there is no conflict of interests.

Additional information

Funding

The authors acknowledge the financial support from the New Energy and Industrial Technology Development Organization (Grant No. JPNP22005) and Japan Society for the Promotion of Science KAKENHI (Grant No. 21H04870). Also, the bead mill CAD data was provided by AIMEX CO.,Ltd. One of the authors, Kai-en Yang, was supported by the Quantum Science and TEchnology fellowship Program (Q-STEP) of the University of Tokyo.