Unsteady Cylinder Wakes from Arbitrary Bodies with Differentiable
  Physics-Assisted Neural Network
This work delineates a hybrid predictive framework configured as a
coarse-grained surrogate for reconstructing unsteady fluid flows around
multiple cylinders of diverse configurations. The presence of cylinders of
arbitrary nature causes abrupt changes in the local flow profile while globally
exhibiting a wide spectrum of dynamical wakes fluctuating in either a periodic
or chaotic manner. Consequently, the focal point of the present study is to
establish predictive frameworks that accurately reconstruct the overall fluid
velocity flowfield such that the local boundary layer profile, as well as the
wake dynamics, are both preserved for long time horizons. The hybrid framework
is realized using a base differentiable flow solver combined with a neural
network, yielding a differentiable physics-assisted neural network (DPNN). The
framework is trained using bodies with arbitrary shapes, and then it is tested
and further assessed on out-of-distribution samples. Our results indicate that
the neural network acts as a forcing function to correct the local boundary
layer profile while also remarkably improving the dissipative nature of the
flowfields. It is found that the DPNN framework clearly outperforms the
supervised learning approach while respecting the reduced feature space
dynamics. The model predictions for arbitrary bodies indicate that the Strouhal
number distribution with respect to spacing ratio exhibits similar patterns
with existing literature. In addition, our model predictions also enable us to
discover similar wake categories for flow past arbitrary bodies. For the
chaotic wakes, the present approach predicts the chaotic switch in gap flows up
to the mid-time range.
Shuvayan Brahmachary, Nils Thuerey
# Unsteady Cylinder Wakes from Arbitrary Bodies with Differentiable Physics-Assisted Neural Network

###### Abstract

This work delineates a hybrid predictive framework configured as a coarse-grained surrogate for reconstructing unsteady fluid flows around multiple cylinders of diverse configurations. The presence of cylinders of arbitrary nature causes abrupt changes in the local flow profile while globally exhibiting a wide spectrum of dynamical wakes fluctuating in either a periodic or chaotic manner. Consequently, the focal point of the present study is to establish predictive frameworks that accurately reconstruct the overall fluid velocity flowfield such that the local boundary layer profile, as well as the wake dynamics, are both preserved for long time horizons. The hybrid framework is realized using a base differentiable flow solver combined with a neural network, yielding a differentiable physics-assisted neural network (DPNN). The framework is trained using bodies with arbitrary shapes, and then it is tested and further assessed on out-of-distribution samples. Our results indicate that the neural network acts as a forcing function to correct the local boundary layer profile while also remarkably improving the dissipative nature of the flowfields. It is found that the DPNN framework clearly outperforms the supervised learning approach while respecting the reduced feature space dynamics. The model predictions for arbitrary bodies indicate that the Strouhal number distribution with respect to spacing ratio exhibits similar patterns with existing literature. In addition, our model predictions also enable us to discover similar wake categories for flow past arbitrary bodies. For the chaotic wakes, the present approach predicts the chaotic switch in gap flows up to the mid-time range.

keywords: Differentiable physics, unsteady cylinder wakes, arbitrary flows, spatio-temporal predictions +
Footnote †: journal: Journal of Computational Physics

## 1 Introduction

Machine learning has drawn significant attention towards the field of fluid simulations (Kutz, 2017; Brunton et al., 2020), particularly in the last two decades. Owing to its successful inception in the field of turbulence modeling (Ling et al., 2016; Wu et al., 2018; Duraisamy et al., 2019; Srinivasan et al., 2019), it has spawned multiple avenues such as super-resolution of fluid flows (Xie et al., 2018; Fukami et al., 2019), detection of turbulent interface (Li et al., 2020), active flow control (Pino et al., 2023), fluid-particle interaction (Davydzenka and Tahmasebi, 2022), reduced-order modeling (ROM) for predicting the dynamical state of a fluid system (Xiao et al., 2015; Pawar
et al., 2019), aerodynamic shape optimization (Chen et al., 2021), among others. To incorporate human knowledge in the form of physics models, methods were introduced that integrate them into the training procedure, _e.g.,_ physics-informed neural network (Raissi et al., 2019; Sun et al., 2020), graph neural network (Belbute-Peres et al., 2020; Brandstetter et al., 2022), and generative adversarial network (Cheng et al., 2020) with physical loss functions (Lee and You, 2019). While custom loss functions satisfy the boundary conditions and minimize the residual of the underlying governing partial differential equations, they only allows for partial coupling of the flow solver with the neural network. Another alternative that allows full integration of the flow solver within the training loop is via differentiable models, _e.g.,_ differentiable physics for turbulence modeling (List et al., 2022). This approach has demonstrated its effectiveness as a cost-efficient strategy for turbulence modeling compared to the cost-intensive direct numerical simulations. Such approaches have also shown the potential to serve as a coarse-grained surrogate for unsteady fluid flows. For instance, the works of Um et al. (2020) and Kochkov et al. (2021) demonstrated that differentiable flow solvers coupled with a neural network can yield satisfactory reconstruction of flowfields for representative test cases. However, the full merits of the hybrid strategy as a generic and accurate approach are far from established.

A primary task when evaluating the merits of a predictive framework is to identify a fluid test problem that is of significant practical interest and complexity. One such problem that has been of significant interest to the fluid dynamics community is that of flow past body (Hasegawa et al., 2020) or multiple bodies Wan and Sapsis (2018). Typical factors and challenges that plague the analysis of such fluid flows may include computational time or the turn-around time, low to large temporal scales in the flow, unsteadiness in the flow, the presence of multiple bodies, a large spectrum of wake dynamics, etc. While coarse-grained surrogates reduce the computational burden (Stachenfeld et al., 2021), the effectiveness of such a strategy to reliably reconstruct the flowfields for multiple wake regimes, besides ensuring accuracy in local boundary profiles and global wake dynamics, is yet to be examined. For instance, authors in work Morimoto et al. (2022) examine flow past two side-by-side cylinders to qualitatively evaluate the generalizability aspect of their models without considering the local fluid behavior around the body boundary. In addition, the work by Hasegawa et al. (2020) examined flow past a single bluff-body of various shapes using a convolutional neural network-based auto-encoder (CNN-AE) model with an emphasis on global fluid dynamics. In a similar vein, the investigation by Lee and You (2019) points out that their model predictions for flow past a single cylinder exhibited the largest errors in either the body boundary or the wake region. Clearly, this presents a significant limitation that should be addressed by advancing performance measuring metrics based on both the local as well as global fluid properties.

Flow past objects are of significant interest to the engineering community from the hydrodynamics as well as structural design point of view. The potential implications of the underlying instabilities (McKinley et al., 1993), wake dynamics (Williamson, 1996), wake induced vibration (Bearman, 2011) and its control (Choi et al., 2008) have led to establishing a significant body of research. One of the focal points concerning this research is understanding the spectrum of wake flow exhibited by multiple cylinders placed relative to each other. Some of the investigations that present a comprehensive study for two cylinders are (Zravkovich, 1987; Papaioannou et al., 2006; Sumner, 2010) etc. and (Lam and Cheung, 1988; Guillaume and LaRue, 1999; Bao et al., 2010; Zheng et al., 2016) for three or more cylinders. These investigations point towards the direct dependency of the resulting wake regimes with the spacing ratio (_i.e.,_, length \(L\) to diameter \(D\) ratio) \(L/D\). More recently, the investigation by (Chen et al., 2020) showed the presence of up to nine
distinct flow regimes for cylinders arranged in an equilateral-triangle position at various spacing ratios, \(1\leq L/D\leq 6\) for multiple Reynolds numbers \(50\leq Re\leq 175\). While these investigations provide fundamental insights about wake flow regimes and flow states, their analysis is limited to the use of equi-diameter circular cylinders or rectangular cylinders with varying aspect ratios. Consequently, this begs the question if the observations related to the universal non-dimensional parameters such as Strouhal number and critical spacing ratio (as reported by (Chen et al., 2020)) still hold for arbitrarily shaped cylinders.

Analyzing large datasets of wake flow exhibited by arbitrarily shaped complex configurations offers potential for novel physical insights; however, employing traditional body-conforming continuum computational fluid dynamics (CFD) flow solvers can be time-consuming and warrants special treatment for catering to changing geometries via the underlying computational grid. While non-conformal immersed boundary (IB) approaches (Mittal and Iaccarino, 2005; Brahmachary, 2019) allow for the use of a fixed Cartesian mesh, exploring a large design space of possible spacing ratios could significantly ramp up the turn-around time. Consequently, neural network-based coarse-grained surrogate models for unsteady fluid flows provide a cost-effective alternative approach (Stachenfeld et al., 2021).

In this work, we build on the ideas by Um et al. (2020) and present a hybrid Differentiable Physics-Assisted Neural Network (or DPNN) framework that acts as a coarse-grained surrogate for accurate unsteady fluid simulations for low Reynolds number flows past multiple bodies. This framework requires a base solver integrated into the neural network (NN) along with ground truth data to train the resulting hybrid framework. We employ a Cartesian grid-based immersed boundary method (_i.e.,_ FoamExtend) as the _Reference_ solver. In addition, we employ an in-house developed differentiable flow solver (_i.e.,_ PhiFlow) as the _Source_ solver. For faster training and inference, the _Source_ intentionally resorts to a first-order spatial accuracy along with masked stair-step representation of the underlying body boundary. This inherently reduces the computational burden of the base differentiable flow solver. In contrast, however, the _Reference_ solver utilizes an accurate local algebraic reconstruction for fluid flow around the body boundary with second order spatial accuracy. Given this stark contrast between the _Reference_ and _Source_ solvers, besides the highly multi-modal nature of the present problem, this naturally poses a very challenging task for NNs _i.e.,_ to act as a forcing function to learn the non-trivial corrections. We comprehensively evaluate the hybrid approach for multiple wake categories while drawing necessary comparisons with various benchmarks. Moreover, our approach of utilizing arbitrarily shaped bodies in a staggered equilateral triangle position allows for a large data diversity and sheds light on the physics of the wake dynamics resulting from the arbitrary cylindrical configurations.

Our focus is thus twofold. First, we train and test the DPNN framework on a dataset comprising unsteady flow past arbitrary bodies placed in a staggered equilateral triangle position, each at a unique spacing ratio. These tests were performed once the _Reference_ solver was throughly validated and verified to generate accurate solutions. Second, after establishing the effectiveness of our approach using multiple performance measuring metrics, we evaluate the out-of-distribution generalization capabilities of our framework. These evaluations is based on configurations that align closely with the current state of existing literature. Throughout the study, we strive to understand the physical implications of the solutions derived from the predictive framework, examining both local boundary layer profiles and global wake dynamics. We also highlight the framework's performance in both physical and reduced feature spaces to demonstrate its long-term temporal stability.
## 2 Numerical setup and validation

In this section, we illustrate the computational setup underpinning the entire hybrid framework. This is followed by validation test cases undertaken to affirm the correctness of the setup as well as verify the ability to render accurate solutions.

### Reference solver

The present study employs the open source FoamExtend V 4.0 immersed boundary (IB) approach as the _Reference_ solver. As mentioned earlier, the IB-based non-conformal approach serves a dual purpose. Firstly, it greatly simplifies the use of a fixed Cartesian mesh for all underlying configurations, irrespective of the geometric complexity. This methodology circumvents the necessity for labor-intensive re-meshing for each geometric configuration and the corresponding coordinate transformation from physical to computational space. Secondly, the use of a fixed Cartesian mesh ensures seamless integration to the deep-learning pipeline wherein one can leverage the benefits of the Eulerian (_i.e._, fixed cell in space) viewpoint. FoamExtend utilizes a discrete forcing approached treatment strategy for enforcing the boundary condition for the immersed bodies and has been built on top of the well-established OpenFoam (Jasak and Tukovic, 2015) framework.

In the present work, we consider the 2D Navier Stokes equations as mentioned below:

\[\begin{cases}\frac{\partial\mathbf{u}}{\partial t}+(\mathbf{u}.\nabla)\mathbf{ u}=-\frac{1}{\rho}\nabla p+\nu\nabla^{2}\mathbf{u}\\ \nabla.\mathbf{u}=0\end{cases} \tag{1}\]

where, \(\mathbf{u}=(u,v)\) is the velocity field, \(\rho\) is the density, \(p\) is the pressure, and \(\nu\) is the dynamic viscosity. The first equation represents the momentum conservation equations, whereas the second equation represents the continuity equation, which also serves as a constraint on the velocity fields, otherwise known as the "divergence-free" constraint.

Figure 1: Computational domain including the three bodies placed in an equilateral triangle position (not to scale)
The flow scenario describes the incompressible flow past cylinders at moderate Reynolds number \(\mathrm{Re}\approx 100\). For the present case, the Reynolds number is evaluated based on the cylinder height of the primary upstream cylinder as \(Re=\frac{UH}{\nu}\), where \(U\) and \(H\) represent the freestream velocity and height of the upstream cylinder respectively. For the present setup, \(\nu\)=0.01 \(\mathrm{Ns}/m^{2}\) and \(U\)=1m/s are used. A second-order accurate Crank Nicolson scheme is used to march the solutions in time, while a second-order limited linear (central-differencing) scheme is used for convective flux discretization. Pressure-Implicit with Splitting of Operators (PISO) (Issa, 1986) is employed to solve for the discretized momentum and continuity equations in 1 predictor and 2 corrector steps. In particular, preconditioned conjugate gradient (PCG) based solver is used for \(p\) linear equations along with a diagonal-based incomplete Cholesky (DIC) precondition for symmetric matrices. In addition, for \(U\) linear equations, BiCGSTAB is used along with diagonal-based incomplete LU (DILU) preconditioned for asymmetric matrices.

### Validation: Strouhal number and force characteristics

We execute a comprehensive grid independence study to corroborate the _Reference_ solver used to acquire the ground truth data. We initiate our numerical experiments with flow past a single cylinder at Reynolds number \(Re\)=100. A computational domain of size [0,24D] \(\times\) [0,16D] is chosen using a uniformly spaced Cartesian grid. In particular, we choose medium and fine spatial resolutions to undertake a grid-independence study, _i.e.,_\(\Delta x/D\) = \(\Delta y/D\) \(\in\) {1/32, 1/50}, where \(D\) represents the diameter of the cylinder. This leads to a blockage ratio, \(B\)=\(D/H\) = 0.0625. Further, we use a constant time-step \(\Delta t\) = 0.1s for this validation test case and perform the simulations until a total time \(T\) = 200s.

The time average drag coefficient \(\bar{C_{d}}\) as well as the root mean square (r.m.s) lift coefficient \(C_{1}^{{}^{\prime}}\) obtained from the _Reference_ flow solver is compared with existing results. Specifically, Table 1 shows the \(\bar{C_{d}}\), \(C_{1}^{{}^{\prime}}\), and Strouhal number \(St\) values for flow past a single cylinder at various grid resolutions. It is found that both \(\bar{C_{d}}\) and \(C_{l}^{{}^{\prime}}\) obtained from the medium as well as the fine spatial

\begin{table}
\begin{tabular}{l c c c c} \(Re\)=100 (single cylinder) & \(\bar{C_{d}}\) & \(C_{1}^{{}^{\prime}}\) & \(St\) \\ Constant et al. (2017) (\(\Delta x=\Delta y=0.02D\)) & 1.38 & - & 0.165 \\ Constant et al. (2017) (\(\Delta x=\Delta y=0.010D\)) & 1.37 & - & 0.165 \\ Present (\(\Delta x=\Delta y=0.0312D\)) & 1.412 & 0.276 & 0.156 \\ Present (\(\Delta x=\Delta y=0.02D\)) & 1.403 & 0.285 & 0.156 \\ \end{tabular}
\end{table}
Table 1: Coefficient of lift \(C_{l}\), drag \(C_{d}\), and Strouhal number \(St\) for flow past a single cylinder at \(Re\)=100

\begin{table}
\begin{tabular}{l c c c c c c} Zheng et al. (2016) & 1.23 & 0.0 & -0.002 & 1.53 & -0.087 & 0.335 \\ Present (\(\Delta x=\Delta y=0.0312D\)) & 1.249 & 0.0 & 0.0 & 1.561 & 0.166 & 0.433 \\ \end{tabular}
\end{table}
Table 2: Coefficient of lift \(C_{l}\), drag \(C_{d}\), and Strouhal number \(St\) for flow past three cylinders at various time-step sizes
resolutions compare really well with the existing solutions of Constant et al. (2017), who also employ a non-conformal immersed boundary-based numerical solver for the computations. It must be noted that a change in the underlying grid from medium to fine spatial resolution results in an increase in the total number of control volumes \(n_{c}\) from 3,93,216 to 9,60,000. Consequently, to allow for a reasonable computational time while retaining accuracy, we choose the medium grid for the remainder of the computations in this section as well as for obtaining the ground truth solution.

Table 2 presents the mean as well as the r.m.s values for the coefficient of drag and lift for the upstream and downstream cylinder (only one) for flow past three cylinders placed in an equilateral triangle position at spacing ratio \(L/D\)=2.5. It is found that the present solutions auger really well with the ones obtained by Zheng et al. (2016), who employ a body-conformal finite-volume approach using a time step size \(\Delta t\)=0.05s. The influence of the time-step size is also analyzed in Table 3. It is found that the choice of time step size results in a marginal difference in force quantities as well as the Strouhal number. Consequently, a moderate value of \(\Delta t\)=0.1 s is adopted for the present study. It can now be remarked that the grid resolution and time-step size chosen for the present study allow for the accurate computation of flow past multiple bodies. In addition, we also qualitatively verify the wake regimes produced using the flow solver, as highlighted in B.

### Differentiable flow solver

The present study uses the open source simulation toolkit Phiflow (\(\phi_{\text{flow}}\)) as the base differentiable flow solver (or _Source_). It leverages automatic differentiation to enable an end-to-end recurrent training paradigm. While it is crucial that the base solver (_i.e.,_ Source) be differentiable, the _Reference_ solver (_i.e.,_ FoamExtend in this case) need not be differentiable. We employ the projection method to decouple the momentum and continuity equation in a two predictor-corrector step. Firstly, we use the operator splitting of the diffusion, advection, and pressure terms in the predictor phase. Specifically, the MacCormack advection step for advecting the initial velocity fields. Secondly, in the corrector phase, the spatial gradients of pressure are used to compute the divergence-free velocity field Chorin (1968).

Finally, the choice of Cartesian grids in both _Reference_ and _Source_ solver allows for a one-to-one mapping between each control cell. This is especially important given that the use of masked stair-step representation of the body boundary by the _Source_ results in a very approximate reconstruction of near-body fluid properties. This allows for a training setup wherein the network learns to faithfully reconstruct the local boundary flow profile based on the accurate _Reference_ solver.

### Cylinder layout

We carefully choose a setup for cylinder cluster arrangement that allows us to mimic the traditional approach adopted in open literature. However, as we seek to generate and exploit a

\begin{table}
\begin{tabular}{l c c c c c c} \(Re\)=100 (three cylinders, \(L/D=2.5\)) & \(\bar{C}_{d,1}\) & \(\bar{C}_{l,1}\) & \(C^{\prime}_{l,1}\) & \(\bar{C}_{d,2}\) & \(\bar{C}_{l,2}\) & \(C^{\prime}_{l,2}\) & \(St\) \\ Present (\(\Delta t\) = 0.25) & 1.254 & 0.0 & 0.0 & 1.584 & 0.172 & 0.49 & 0.148 \\ Present (\(\Delta t\) = 0.1) & 1.249 & 0.0 & 0.0 & 1.561 & 0.166 & 0.433 & 0.156 \\ Present (\(\Delta t\) = 0.05) & 1.244 & 0.0 & 0.0 & 1.545 & 0.166 & 0.359 & 0.156 \\ \end{tabular}
\end{table}
Table 3: Coefficient of lift \(C_{l}\), drag \(C_{d}\), and Strouhal number \(St\) for flow past three cylinders
multitude of complexities in the training data to harness the full capabilities of the present hybrid learning strategy, we introduce some randomness to the setup. Some specifications related to the cylinder arrangements are as follows:

1. **Cylinder arrangement**: We place three cylinders in an equilateral triangle position, _i.e.,_ one upstream cylinder followed by two downstream cylinders, placed apart from each other via the spacing ratio \(L/D\). This setup mimics the "regular triangular" cluster setup used in Zravkovich (1987) and has also been investigated recently by Chen et al. (2020). Additionally, this arrangement allows us to investigate _proximity_ interference (P), _wake_ interference (W), a combination of (P+W) interference as well as no-interference.
2. **Upstream composite cylinder**: The upstream cylinder (cylinder \(C_{1}\)) is framed as a composite body consisting of a base rectangle (at its center) and secondary rectangle or semi-circle on its sides (see Fig. 1). The selection between a secondary rectangle or a secondary semi-circle as one of the sides of the base rectangle is random. We define a function \(f_{1}:R_{\text{u,i}}\to S\), such that it maps a random integer \(R_{\text{u,i}}\in\{\text{0,1}\}\) to geometric shapes, \(\text{S}\in\{\text{semi-circle, rectangle}\}\), where \(i\in\{\text{1,2,3,4}\}\). The outcome can be represented as follows: \[R_{\text{u,i}}=\begin{cases}0&f_{\text{i}}=\text{semicircle}\\ 1&f_{\text{i}}=\text{rectangle}\end{cases}\] The length of the base rectangle, as well as the secondary rectangles, are randomly chosen between fixed upper and lower limits (see Table 4). The diameter of the semi-circle is equal to the edge length of the base rectangle to which it is attached.
3. **Downstream cylinders**: The choice of the downstream cylinders is either a rectangular cylinder of a certain length and height or a circular cylinder of a certain diameter. This selection is partially random. We define a function \(f_{\text{2,3}}:R_{\text{d}}\to S\), such that it maps a random integer \(R_{\text{d}}\in\{\text{0,1}\}\) to geometric shapes, \(\text{S}\in\{\text{circle, rectangle}\}\). Once the upper downstream cylinder (cylinder \(C_{2}\)) is chosen (say a rectangular cylinder), the lower downstream cylinder (cylinder \(C_{3}\)) gets fixed (a circular cylinder). This outcome can be represented as shown below.

Figure 2: Boundary reconstruction approach used by (a) FoamExtend (_Reference_ solver) (b) Phiflow (_Source_ solver)
\[R_{\rm d}=\begin{cases}0&f_{2}=\text{circle};\ \ f_{3}=\text{rectangle}\\ 1&f_{2}=\text{rectangle};\ \ f_{3}=\text{circle}\end{cases}\]

This forces the training dataset to always contain a rectangle and a circular cylinder as the two downstream cylinders. In addition, once the height of the downstream rectangular cylinder, \(H_{\rm C2}\) (or the diameter of the downstream circular cylinder, \(D_{\rm C3}\)) is chosen, the dimension (_i.e.,_ diameter or cylinder height) of the other downstream cylinder is fixed to satisfy the following criterion.

\[H_{\rm C2}+D_{\rm C3}=2H_{\rm C1} \tag{2}\]

This restricts the cylinders from overlapping. The dimensions of the rectangle and the circular cylinder are randomly chosen between a fixed upper and lower limit (see Table 4).
4. **Spacing ratio**: The distance between centre-to-centre distance between the downstream cylinders, as well as the downstream and upstream cylinders, is controlled by the spacing ratio \(L/D\). The \(L/D\) values are chosen such that it covers a wide spectrum of possible wake flow regimes (_i.e.,_\(1.2\leq L/D\leq 5.5\)). Care is taken to avoid overlapping bodies for low \(L/D\) values. Finally, to ensure an unbiased selection of the spacing ratio, the Latin Hypercube Sampling (LHS) is chosen as a design of the experiment method.

Points 2 and 3 above intentionally introduce arbitrariness to the training data to serve a dual purpose. Firstly, we analyze the role of arbitrariness in the wake flow regime to offer new insight into the present body of knowledge around flow past equi-diameter cylinders. Secondly, this introduces data diversity while also allowing us to evaluate the generalizability of the hybrid learning framework to more representative problems. Further, point 4 allows us to uncover various wake categories that have been introduced recently by Chen et al. (2020) for similar spacing ratios.

### Wake categories

This section introduces the various wake categories obtained from the _Reference_ solver at multiple spacing ratios \(L/D\) for equi-diameter cylinders placed in an equilateral triangle position. These representative spacing ratios for equi-diameter cylinders are for illustration purposes and are not part of the training data. We purposely choose equi-diameter circular cylinders so that fair comparisons can be drawn with Chen et al. (2020) who also evaluate it on identical conditions.

**Single bluff-body wake:** At extremely small spacing ratio _i.e.,_\(L/D=1.2\), the three cylinders are close to each other such that they exhibit vortices as if it were a "single bluff-body". This is further demonstrated pictorially in Fig. 30 in the Appendix B, along with other representative cases. This kind of wake is periodic in nature.

**Deflected gap wake:** At slightly larger spacing ratio _i.e.,_\(L/D=1.5\), the fluid in between the two downstream cylinders (or gap flow) switches direction to either the upper cylinder (_i.e.,_ cylinder 2) or the lower cylinder (_i.e.,_ cylinder 3), hence the name "deflected gap". This switch is constant with time. The kind of wake is periodic, with some modulation over time.

**Flip-flopping wake:** For marginally larger spacing ratio _i.e.,_\(L/D=2.25\), the gap flow in between the two downstream cylinders switches direction erratically between the upper cylinder (_i.e.,_ cylinder 2) or the lower cylinder (_i.e.,_ cylinder 3), hence the name "flip-flopping". This chaotic switch is a challenging case for any predictive framework, and the kind of wake is irregular in nature.
**Anti-phase wake:** For moderately larger spacing ratio _i.e.,_\(L/D\) = 3.5, the vortices shed by the two downstream cylinders rotate in opposing directions, hence the name "anti-phase". This wake category is periodic in nature and exhibits a symmetry about the centerline streamwise direction. For arbitrarily shaped bodies, this symmetry is lost and is consequently referred to in this work as "quasi anti-phase" wake.

**Fully developed in-phase wake:** For larger spacing ratio _i.e.,_\(L/D\) = 5, the vortices shed by the two downstream cylinders are in same direction; hence the name "in-phase". Further, the upstream cylinder also sheds vortices which merges with the vortices shed by the downstream cylinders, hence the term "fully-developed". This kind of wake is periodic in nature. For arbitrarily shaped bodies, this is referred to in this work as "quasi in-phase" wake.

\begin{table}
\begin{tabular}{l l c c c c c c} \multicolumn{2}{c}{Component} & \multicolumn{2}{c}{Length (\(L\))} & \multicolumn{2}{c}{Height (\(H\))} & \multicolumn{2}{c}{Radius (\(R\))} \\  & & & lower & upper & lower & upper & lower & upper \\  & & & limit & limit & limit & limit & limit & limit \\  & (primary) & & 0.4 & 0.5 & 0.4 & 0.5 & - & - \\ Upstream & \multicolumn{2}{c}{Base rectangle} & & & & & & \\ Cylinder & \multicolumn{2}{c}{(secondary)} & & & & & & \\  & \multicolumn{2}{c}{Rectangle} & & & & & & \\  & \multicolumn{2}{c}{Rectangle} & & & & & & \\ Downstream & \multicolumn{2}{c}{Rectangle} & & & & & & \\ Cylinder & \multicolumn{2}{c}{Circular cylinder} & & & & & & \\ \end{tabular}
\end{table}
Table 4: Upper and lower limits for the dimensions of the cylinders

Figure 3: DPNN-based hybrid predictive framework
## 3 Hybrid predictive framework

This section introduces the methodology behind the DPNN-based hybrid framework used in the present work. As alluded earlier, such a strategy embeds the base differentiable flow solver Phiflow with a neural network (NN) architecture, mimicking the _solver-in-the-loop_ strategy previously employed by Um et al. (2020). Such a strategy reaps the benefit of long-term training feedback while enabling an end-to-end training paradigm. We build on the above by enabling a robust and efficient coarse-grained surrogate model for unsteady flows for investigating one of the core fluid flow problems, _i.e.,_ wake flows caused by arbitrarily shaped objects. This serves a dual purpose, _i.e.,_ firstly, the arbitrarily shaped bodies allow for a rich and diverse set of datasets. Secondly, while it is well known that the shape of the body is responsible for the resulting wake characteristics, such an issue has not prompted an investigation into a large spectrum of spacing ratios to the best of our knowledge.

### Neural network architecture

In this section, we highlight the neural network (NN) architectures that were explored in the present study. We explore a range of popular architectures while keeping the total number of parameters in each network architecture within reasonable limits to facilitate a fair comparison.

The network takes as input of channel size = 3, _i.e.,_ velocity fields, \(\mathbf{u}=u,v\) along with the shape mask or marker field \(m_{f}\in R^{n_{x}\times n_{y}}\), where \(n_{x}\) and \(n_{y}\) represents the number of grid points along the longitudinal and lateral directions, respectively. The shape mask \(m_{f}\) uniquely identifies the shape of the enclosed body and is fixed for a given test configuration. While the ground truth data is obtained on a high-resolution grid of size \(n_{x},n_{y}\) = (\(768\times 512\)), the network only receives a downsampled (3x) velocity and shape mask as inputs (_i.e.,_\(96\times 64\)), for faster training. For a training sample, there are a total of \(t_{f}\) = 3000 frames or snapshots (of which the first 1500 are discarded to only consider statically stable snapshots), each representing velocity fields generated sequentially using a constant time-step of \(\Delta t\) = 0.1s (or, total time \(t\) = 300s). The entire dataset consists of 100 different experiments each at a unique spacing ratio, resulting in 100 unique shape masks \(m_{f}\). As mentioned in the earlier Section 2.4 and Fig. 1, to allow for a fair comparison, we employ the same cylinder cluster arrangement as also employed by Chen et al. (2020), with the equi-diameter cylinders employed by Chen et al. (2020) being replaced by arbitrarily shaped cylinders in the present work. Among the entire dataset, 50 are used for training, and the remaining 50 are used for testing. The training is executed using a batch size of 50 for 50 epochs using Adam (Kingma and Ba, 2014) as the optimizer. A variable learning rate \(\alpha\) is employed, which is sequentially reduced with epochs. The output from the network is constrained to be of channel size = 2, _i.e.,_ the velocity field at the next time step. In this work, we predominantly employ residual neural network (ResNet)-based architecture besides also using other popular designs. Appendix C provides the details related to each architecture.

### Present hybrid framework

The DPNN-based hybrid predictive framework is a _recursive_ learning strategy wherein the base flow solver Phiflow (\(\phi_{\text{flow}}\)) is coupled with the neural network as shown in Fig. 3. In the present study, the _Reference_ solver (_i.e.,_ FoamExtend) provides the pre-computed ground truth data. The hybrid framework receives an initial guess \(\mathbf{r}^{\prime}\) from the _Reference_ solver and feeds it to the base solver (_i.e.,_ Phiflow). The base solver then updates the velocity field to the next time-step \(\mathbf{u}^{*}\). At this
stage, no learning is involved. The output from PhiFlow is then fed to the NN architecture as input, _i.e.,_ velocity fields \(u,v\), and shape marker field \(m_{f}\). The output from the NN, \(\mathbf{f}(\mathbf{u}^{*};\theta_{W})\) acts as the forcing function which is used to correct the output from the base solver PhiFlow (\(\theta_{w}\) represents the network weights). This step is made recursive for \(m\) solver unrollment steps signifying as many forward steps by base solver in time by feeding the cumulative output back to the base solver, _i.e.,_\(\mathbf{u}^{t+1}=\mathbf{u}^{*}+\mathbf{f}(\mathbf{u}^{*};\theta_{W})\).

We remark that both the _Reference_ and the base solvers are subject to identical values for the dynamic viscosity \(\nu\), time step size \(\Delta t\), and freestream conditions, besides making sure to incorporate the same object under consideration. The critical difference here lies in their respective solution strategies and the techniques to handle body boundary conditions, _i.e.,_ while the former is spatially second-order accurate and algebraically reconstructs the fluid properties along the body boundary to locally satisfy the no-slip boundary conditions, the latter is spatially first-order accurate and resorts to a masked stair-step representation for computational efficiency. This fact serves as the key motivation to learn the global wake dynamics as well as the local boundary representation via the loss formulation, as will be discussed in the next Section.

### Loss formulation

The loss formulation provides the central goal of the learning process. In this study, this means recovering the ground truth velocity flowfields as accurately as possible. The trainable parameters of the model are iteratively updated by minimizing the loss function. This allows the optimizer to navigate the non-convex energy landscape with the goal to arrive at a global/local minima. Typical choices of a loss function \(\mathcal{L}\) use some error norm to compare network prediction \(\mathbf{u}\) with the ground truth \(\mathbf{r}\). For instance, the \(\mathcal{L}_{2}\) norm for a purely supervised training is given as:

\[\mathcal{L}=\mathcal{L}_{2}(\mathbf{r}-\mathbf{u}) \tag{3}\]

Such a formulation minimizes the Euclidean distance between \(\mathbf{u}\) and \(\mathbf{r}\). Although popular, such a loss is devoid of additional information from the physics point of view. For instance, Um et al. (2020) showed that one can formulate physics-based loss formulations by taking advantage of the temporal unrolling during the training. This was followed by an independent investigation (List et al., 2022) wherein they built a custom loss function to accurately predict turbulent flows. Consequently, a greater emphasis is placed on rendering a loss yielding a physically consistent solution. Such a physics-based loss could take the form:

\[\mathcal{L}=\frac{1}{m}\sum_{i=1}^{m}\mathcal{L}_{2}(\mathbf{r}^{t+i}- \mathbf{u}^{i}) \tag{4}\]

where \(m\) represents the number of unrollment steps used during training and \(t\) represents the instant of time. Interested readers are referred to the webbook Thuerey et al. (2021) for further details.

### Benchmarks

In addition to our present DPNN framework, we also compare it to existing benchmarks. As part of this comparative evaluation, we employ the data-driven supervised learning (SL) method. In addition, we evaluate results produced by only running the base solver _i.e.,_ _Sourcea_. For the SL approach, the network directly receives the input \(\mathbf{r}^{t}\) from the _Reference_ solver and produces the
output as the flowfield for the subsequent time step _i.e.,_\(\mathbf{f}(\mathbf{r}^{t};\theta_{W})\). As this SL framework does not contain an integrated _Source_ solver, the output generated by the network serves as the final prediction for the next time step (_i.e.,_\(\mathbf{u}^{t+1}=\mathbf{f}(\mathbf{r}^{t};\theta_{W})\)). This makes the SL approach faster. Notably, to ensure a fair comparison, the outputs from both DPNN and SL frameworks utilize identical underlying networks and receive the same input during training/testing. For the output produced from the _Source_ solver, the next time step predictions can be written as \(\mathbf{u}^{t+1}=\mathbf{u}^{*}\). Finally, For reasons of completeness and validation, we also juxtapose our predictions with those from the _Reference_ solver. Comparisons with the _Source_ solver highlight the starting point on which the DPNN solver learns to apply its corrections. Further details related to the network architecture may be found in Appendix C.

## 4 Results and Discussions

In this section, we present a detailed analysis of the results obtained from the present hybrid learning framework as well as other benchmarks. As mentioned earlier, we employ multiple neural network architectures within the hybrid learning approach. Henceforth, we begin by undertaking an analysis that compares various network architectures for a common problem in the following section.

### Role of Network Parameters in Model Performance

We start by presenting the statistical evaluations of models builds using different network architectures for all the testing samples. An effective representation of the distribution of the model performance across a large spectrum of testing samples is done via kernel density estimation (or KDE). Besides highlighting the model confidence, it further allows for the potential identification of outliers. The KDE is evaluated based on the mean absolute error, \(\mu\) of the model predictions compared to the ground truth data, for the 100th testing frame. The mean absolute error, \(\mu\) is evaluated as follows:

\[\mu=\frac{1}{n_{c}}\sum_{i=1}^{n_{c}}(|u_{\text{ref,i}}-u_{\text{ baseline,i}}|)+(|v_{\text{ref,i}}-v_{\text{baseline,i}}|) \tag{5}\]

\begin{table}
\begin{tabular}{c c c c} Architecture & No. of parameters & Mean, \(\mu_{N}\) & Standard deviation, \(\sigma_{N}\) \\ CNN & 546,478 & 0.03124 & 0.0082 \\ ResNet & 516,674 & **0.02626** & 0.00657 \\ Unet (with DS) & 520,342 & 0.07693 & 0.0138 \\ Unet & 520,342 & 0.06067 & 0.0152 \\ Dil-ResNet (with DS) & 548,066 & 0.07266 & 0.0211 \\ Dil-ResNet & 548,066 & 0.04303 & 0.01101 \\ ResNeXt & 524,322 & 0.03871 & 0.00787 \\ DenseNet & 566,578 & 0.04258 & **0.0047** \\ \end{tabular}
\end{table}
Table 5: Comparative performance among various NN architectures in terms of mean \(\mu_{N}\) and standard deviation \(\sigma_{N}\) (across 50 testing samples, at the end of 100 testing frame)
\begin{table}
\begin{tabular}{c c c} Random seed number & Mean, \(\mu\) & Standard deviation, \(\sigma\) \\
1 & 0.03029 & 0.00818 \\
2 & 0.02643 & 0.00747 \\
3 & **0.02483** & **0.00706** \\
4 & 0.02991 & 0.00896 \\
42 & 0.02626 & 0.0082 \\ \end{tabular}
\end{table}
Table 6: Comparative performance of ResNet architecture for different random seeds

Figure 4: Comparison of the mean absolute error \(\mu\) at the end of 100th testing frame for 50 testing samples obtained from different NN architectures
Figure 4 presents the approximate probability distribution via the KDE for the 50 testing samples from different network architectures. It is found that most of the approximate probability distribution is indicative of a Gaussian-like distribution. These include the KDE obtained from the convolution neural network (CNN) and the residual network-based models _i.e.,_ ResNet (He et al., 2016), together with the newer variants ResNeXt (Xie et al., 2017), DenseNet (Huang et al., 2017), and diluted residual network (Dil-ResNet) from Stachenfeld et al. (2021). On the contrary, however, the Unet (Ronneberger et al., 2015) exhibits multiple peaks in the distribution while also incurring greater sample mean error \(\mu_{N}\) (see Table 5). While attributing higher error to the architectures is not straightforward, strikingly, both the Unet with downsample (with DS) and the Dil-ResNet (with DS) architectures employ encoder blocks that compress the input data to a reduced dimension via the stride operation. This is not entirely surprising as the model struggles to retain local features of the input upon downsampling. We remark that identical choice of activation function, optimizer, learning rate \(\eta\), initial weights \(\theta_{\text{W}}\), loss function, etc. have been used for all the architectures. Consequently, it can be mentioned that the models without any compression or convolution downsample serve as a better alternative for the present choice of problems.

Neural network-based learning is a stochastic, non-linear optimization process that aims to iteratively improve the model predictions on training data while managing uncertainty and navigating through a complex, multidimensional solution space. In this context, using a gradient-based optimizer like Adam may result in a locally-optimal solution with some sensitivity to the initial guess. To check the robustness of the model, the influence of the random number seed must be evaluated. We investigate the performance of the ResNet-based architecture for multiple seeds. Figure 5 shows the approximate probability density based on the mean absolute error \(\mu\) obtained using ResNet-based architectures for different values of random number seed for all the testing samples. It can be observed that the Gaussian-like distribution of the error \(\mu\) is retained for different random seed values, besides showing consistency in the model performance in terms of the sample mean \(\mu_{N}\) (see Table 6). In light of these discussions, it can be remarked that the performance is

Figure 5: Influence of random number seed on the performance of residual network (ResNet) architecture
independent of the random number seed. Finally, we also determine the influence of the adaptive learning rate \(\eta\) that determines the convergence and consequently the nature of solutions obtained. This is briefly discussed in Appendix C.1. The remainder of the study is based on the ResNet architecture.

### Qualitative comparison of predictive models across wake categories

In this section, we compare various benchmark models, assessing their capability to accurately reproduce the dynamic wake characteristics for new, unseen testing samples with different spacing ratios. This analysis is fundamentally important because, intriguingly, different spacing ratios between the cylinders yield diverse wake behaviors, each of which necessitates faithful reproduction. For instance, the deflected gap wake should exhibit a narrow and a wide wake immediately downstream of the two subsequent cylinders, with the gap flow between these cylinders being deflected towards the narrow wake region. This deflection should remain constant with time. In contrast, for the flip-flopping or the chaotic wakes, the gap flow switches direction in an erratic manner. As such, while some wake patterns exhibit a periodic, uniform repetition, others veer into the realm of chaos. It is, therefore, a fascinating problem that tests the predictive capabilities of the models.

Five representative testing samples are chosen from the available 50 experiments, each exhibiting a unique spacing ratio and the corresponding wake dynamics. We present vorticity flowfields generated entirely from the network's output. This allows us to highlight the regions of interest, _i.e.,_ cylinder wake and near body boundary, besides also indicating the direction of rotation of the vortices. Figure 6 compares the vorticity flowfield obtained from the various benchmarks, including

Figure 6: Comparison of velocity flowfields obtained using various benchmarks for different time instances for the spacing ratio \(L/D=1.26\) (single bluff-body wake, 17th testing sample)
the _Reference_ solver at the spacing ratio \(L/D\) = 1.26. The vorticity flowfield obtained from the _Reference_ exhibits a single bluff-body wake with a certain shedding frequency. It is found that while all the baselines can reproduce the correct category of the wake, the _Source_ and the supervised learning (SL) approach very quickly deviate from the correct vortex shedding cycle. However, the present DPNN approach compares remarkably well with the _Reference_ data even for long time horizons.

Figure 7 presents the comparison of the vorticity flowfield obtained at the spacing ratio \(L/D\) = 1.57 for which the _Reference_ solver portrays deflected gap wake. As alluded to earlier, the deflected gap flow shifts toward the narrower wake region (_i.e.,_ lower downstream cylinder in this case), with the direction of the gap flow remaining stable over time. Consequently, the predictive models should preserve this phenomenon while accurately predicting the vortex-shedding cycle. It is found from Fig. 7 that the predictions from _Source_ quickly transitions to a single bluff-body wake, whereas the vorticity contours obtained from the SL approach transform into a hybrid quasi in-phase wake. In contrast, however, the current DPNN approach aligns impressively with the _Reference_ data across extended time horizons. This observation holds true for the deflected gap flow direction and the overall vortex shedding cycle.

Figure 8 showcases a comparison of the vorticity flowfield at a spacing ratio \(L/D\) = 2.34, where the _Reference_ solver displays a flip-flopping (or chaotic) wake, _i.e.,_ the gap flow in between the two downstream cylinders chaotically switches direction. This is evident in Fig. 8 wherein the gap flow pointing downwards from the frame at \(t=60\Delta t\) has transitioned to the top at \(t=120\Delta t\). The

Figure 7: Comparison of vorticity flowfields obtained using various benchmarks for different time instances for the spacing ratio \(L/D\) = 1.57 (deflected gap wake, 39th testing sample)
predictions from DPNN seem to capture this switch correctly up to frame \(t=180\Delta t\), whereas, for long temporal horizons, both the gap flow direction and the vortex shedding cycles have visible differences. These differences further augments when compared to the vorticity flowfields obtained from _Source_ and SL approaches.

Figure 9 exhibits a comparative analysis of vorticity flowfields at a spacing ratio of \(L/D\) = 3.37, which, as per the _Reference_ solver, corresponds to a fully developed quasi in-phase wake. The clockwise rotating vortices (highlighted in blue) emitted from the freestream side of the upper cylinder and the gap side of the lower cylinder, appear to be in quasi in-phase. The same can be mentioned for the anti-clockwise rotating vortices (highlighted in red) shed from the freestream side of the lower cylinder and the gap side of the upper cylinder. It is found that among all the baselines, both DPNN and SL approaches perform well, with the former preserving the vortex structure quite remarkably for long-time horizons. On the other hand, the _Source_ depicts a chaotic wake and fails to reproduce the actual wake dynamics for the given spacing ratio and geometric configuration.

Figure 10 highlights the vorticity flowfields obtained at the spacing ratio of \(L/D\) = 3.89 from various benchmarks. The flowfield obtained from the _Reference_ solver indicates a fully-developed quasi anti-phase wake. It is observed that the vortices emitted from the freestream side of the two downstream cylinders shed vortices in the opposite direction. This dynamic nature of the wake is accurately retained by the present DPNN approach. On the contrary, however, the SL approach produces a hybrid wake, _i.e.,_ wake transitions into a hybrid wake, whereas the _Source_ yields a chaotic wake.

Figure 8: Comparison of vorticity flowfields obtained using various benchmarks for different time instances for the spacing ratio \(L/D\) = 2.34 (chaotic wake, 41st testing sample)
The exercise underscored the crucial need for developing robust and precise models. These models should accurately identify the underlying wake category and deliver sustained temporal fidelity over extended periods. While the current section provides a qualitative comparison of individual baselines across multiple wake categories, the subsequent section will present a quantitative analysis of the key variables of interest.

### Quantitative comparisons in physical and reduced feature space

In the preceding section, we employed velocity and vorticity flowfields as a qualitative measure of performance accuracy. The network used in our study is trained purely on physical space, specifically the velocity flowfields as a collection of control volumes with staggered faces. Through the actions of the kernels, convolutions transform these physical quantities into an equivalent dimensional feature space, achieved through zero-padding. This process poses an engaging question - does the network retain its performance in Fourier space or the reduced feature space offered by the first three principal components?

We begin our investigation by adopting the 39th testing sample, which depicts a deflected gap wake, as discussed in Section 4.2. In Fig. 11, we present the temporal variation of the velocity components \(u,v\) at multiple spatial probes using different baselines. These probes are placed at the wake of the downstream cylinders along the centerline axis, _i.e._, \((n_{\mathrm{px}},n_{\mathrm{py}})\) = (48,32), (64,32), (80,32), where \(n_{\mathrm{px}},n_{\mathrm{py}}\) represents the nodal positioning in the \(x\) and \(y\) directions, respectively. It must be noted here that the present DPNN approach from hereon refers to \(m\)=4 steps of unrolling during training. Thus, an evaluation period up to \(t=300\Delta t\) is a substantial period wherein multiple

Figure 9: Comparison of vorticity \(\omega\) flowfields obtained using various benchmarks for different time instances for the spacing ratio \(L/D\) = 3.37 (fully developed quasi in-phase wake, 26th testing sample)
vortex-shedding cycles are captured. It is found that the present DPNN framework results in very good agreement with the modulated signals obtained for both \(u,v\) velocity components, with the agreement seemingly better for probes placed further away from the bodies. The minor offset with the data obtained from the _Reference_ solver tends to happen for higher time values \(t\). On the contrary, however, this offset is rather early for the predictions obtained from the _Source_ solver. Moreover, for higher values of time \(t\), the predictions significantly differ from the _Reference_.

Figure 12 displays the distributions of \(u\)-velocity across several cross-sections aligned in the longitudinal direction (along the rows) at multiple instances of time (along the columns), obtained from multiple benchmarks. Barring some minor differences, we found that the current framework aligns exceptionally well with data obtained from the _Reference_ solver across various spatial locations and at different time instances. Similar observations can also be drawn from Fig. 13, where the \(u\)-velocity distributions are plotted along multiple horizontal cross-sections at different instances of time. Naturally, these observations are also reflected in Fig. 14 where the mean absolute error, \(\mu\) (calculated based on Eq. 5) for the recursive predictions with respect to the ground truth data are shown. It is found that the present DPNN-based approach yields the lowest cumulative error over time, whereas the data from the _Source_ results in the highest growth of error over time. As such, it can be mentioned that the DPNN approach has successfully enabled a hybrid framework that can faithfully reproduce the dynamical behavior of the physical quantity (_i.e.,_ velocity flowfields).

We now compute the power spectral density (or PSD) based on the velocity flowfields obtained from the predictive frameworks. Figure 15(a)-(c) presents the PSD obtained from the velocity measurements at the \(n_{\mathrm{px}},n_{\mathrm{py}}\) = (64,32) probe location with respect to different time ranges. It is

Figure 10: Comparison of vorticity \(\omega\) flowfields obtained using various benchmarks for different time instances for the spacing ratio \(L/D\) = 3.89 (fully developed quasi anti-phase wake, 32nd testing sample)
Figure 11: Comparison of velocity downstream of the bodies and its variation with time (39th testing sample).
Figure 12: Comparison of the \(u\)-velocity at multiple vertical cross-sections (39th testing sample)
Figure 14: Variation of the mean absolute error from various frameworks with testing frames for the 39th testing sample

Figure 13: Comparison of the \(u\)-velocity at multiple horizontal cross-sections (39th testing sample)
found that for a relatively smaller time range of up to 100 \(\Delta t\), the PSD vs. frequency obtained from the DPNN aligns exceptionally well with the data obtained from the _Reference_ solver signifying that the present framework is able to faithfully reconstruct small as well as large scale structures in the velocity flowfield. However, minor differences in the high-frequency range are found for large recursive predictions (_i.e.,_ up to 300 \(\Delta t\)), indicating small-scale differences in the velocity flowfields. This bodes well with the observations in Fig. 11, where observable differences are found for larger recursive predictions.

We perform principal component analysis (or PCA) on the predicted velocity fields obtained from the baselines to evaluate the principal components (or PC) and their variation with respect to time \(t\). Figure 16 presents the first three PC's and their distribution with time. It is found that the temporal variation of PC 1 and PC 2 indicates a periodicity in the data obtained from the _Reference_ solver, which seems to be captured correctly in predictions obtained by DPNN. Moreover, the temporal variation in the PC 2, PC 3, and PC 3, PC 1 components are also in good overall agreement with the _Reference_ solver. The protracted vortex shedding cycle obtained from _Source_ observed in Fig. 11 is also reflected in Fig. 16 in terms of the open-ended arc indicating a longer vortex shedding time-period. As a result of these investigations, it can be remarked that the performance of the DPNN approach demonstrates consistency in both the physical and reduced feature space.

### Reconstructing flowfields: global and local perspectives

As we delve deeper into investigating the present coarse-grained surrogate model, it becomes increasingly clear that both local and global fluid structures are crucial for the overall representation of flow physics. Specifically, precisely reproducing local boundary layer phenomenon and global wake dynamics becomes essential. The intersection of these scales - the local and the global - is where true accuracy in understanding and predicting flow dynamics lies. To that extent, we investigate the ability of the learned models to correctly reconstruct the local mean velocity boundary layer profile (with momentum thickness, \(\theta\)) as well as the mean gap flow \(\bar{U_{\text{G}}}\), while also evaluating the global fluid variables in terms of kinetic energy \(KE\) and enstrophy \(\Omega\). This step is not trivial

Figure 15: Comparison of the power spectral density of the streamwise velocity for multiple time ranges (39th testing sample)
Figure 16: Comparison of the first three principal components and its variation with time (39th testing sample).
given that the _Source_ solver employs a masked stair-step representation for the underlying body boundary for computational efficiency.

It is also crucial to remark here that while three level coarsening will naturally influence the local resolution of the boundary layer profile it still retains a parabolic profile for the upstream body, as will be shown later on. Thus at a given level of grid resolution, the goal is to mitigate the difference between the under-resolved boundary layer of the _Source_ and the resolved boundary layer of the _Reference_ transferred to the _Source_ mesh. Hence, our goal is not a super-resolution task, i.e. transferring a low-resolution solution to a higher resolution, but instead to improve the solution given a fixed computational mesh.

The momentum thickness, \(\theta\) is evaluated as below,

\[\theta(x)=(\int_{0}^{\delta_{\max}}\frac{\bar{u}(x,y)}{\bar{u}_{e}(x)}\left(1- \frac{\bar{u}(x,y)}{\bar{u}_{e}(x)}\right)dy) \tag{6}\]

The time and space averaged streamwise velocity at the gap Chen et al. (2020), \(U_{G}/U\) is evaluated as follows,

\[\frac{U_{\text{G}}}{U_{\infty}}=\frac{1}{U_{\infty}}\int_{-0.5}^{0.5}\bar{u} \ d\left(\frac{y}{G}\right) \tag{7}\]

where, \(\bar{u}\) is the time-mean streamwise velocity at the gap. The kinetic energy \(KE\) and enstrophy \(\Omega\) are calculated as follows,

\[KE(t,U)=\frac{1}{2}\int_{L_{2(A)}}U(t)^{2}dA \tag{8}\]

\[\Omega(t,U)=\frac{1}{2}\int_{L_{2(A)}}\omega(t)^{2}dA \tag{9}\]

While KE is a popular choice for statistical evaluation of fluid flows, enstrophy represents the rotational energy of the flow and corresponds to the dissipative effects of the flow (Tong et al., 2015).

Figure 17 presents the boundary layer profile in the lateral \(y\) direction just above the body boundary along the cross-section drawn at the upstream cylinder for three distinct testing samples.

Figure 17: Comparison of the time-averaged boundary layer profile along the lateral \(y\) direction of the upstream body for (a) 28th (b) 40th (c) 46th testing samples
Despite the fact that downsampling high-resolution data obtained from the _Reference_ solver can distort the actual boundary profile, reconstructing the steep rise in velocity remains a critical criterion to be preserved. It is evident that the inadequate treatment of the body boundary results in a much thicker boundary layer profile obtained from the _Source_. This results in an enlarged effective diameter of the body, subsequently changing the equivalent effective spacing ratio. Such alteration points to a potential difference in the wake category (single bluff-body wake instead of deflected gap wake from _Source_ in Fig. 7). On the contrary, the boundary layer profile obtained from the DPNN approach agrees quite well with the _Reference_ data. This is also observed in Table 7, where the momentum thickness \(\theta/D\) has been reported for three representative testing samples.

Figure 18 shows the streamwise mean gap velocity profile distribution between the two downstream regions. The prominent peaks adjacent to the parabolic profile near the center are caused by the interaction of the separated shear layer from the upstream cylinder with the gap side shear

\begin{table}
\begin{tabular}{c c c c}  & \multicolumn{2}{c}{Mean streamwise velocity at the gap, \(\frac{U_{\text{G}}}{U_{\infty}}\)} \\ Testing samples no & _Reference_ & _Source_ & DPNN \\
28 & 0.896 & 0.530 & 0.857 \\
40 & 0.956 & 0.657 & 0.958 \\
46 & 0.958 & 0.728 & 0.946 \\ \end{tabular}
\end{table}
Table 8: Comparison of the time and space averaged streamwise velocity at the gap \(\frac{U_{\text{G}}}{U_{\infty}}\) for multiple testing samples

\begin{table}
\begin{tabular}{c c c c}  & \multicolumn{2}{c}{Mean streamwise velocity at the gap, \(\frac{U_{\text{G}}}{U_{\infty}}\)} \\ Testing samples no & _Reference_ & _Source_ & DPNN \\
28 & 0.896 & 0.530 & 0.857 \\
40 & 0.956 & 0.657 & 0.958 \\
46 & 0.958 & 0.728 & 0.946 \\ \end{tabular}
\end{table}
Table 7: Comparison of the boundary layer momentum thickness \(\theta/D\) for multiple testing samples

Figure 18: Comparison of the time-averaged streamwise gap velocity profile in between the two downstream bodies for (a) 28th (b) 40th (c) 46th testing samples
layers of the downstream cylinders. This interaction penetrates the gap flow. There seems to be a very good agreement of the mean streamwise velocity profile at this gap obtained from the DPNN approach with that of the _Reference_, whereas the thicker boundary layer profile exhibited by the _Source_ results in a narrower gap profile. Finally, we also quantify the time and space averaged streamwise velocity at the gap \(\frac{U_{\Omega}}{U_{\infty}}\) in Table 8, where a very good agreement between DPNN and _Reference_ data has been noted.

To enable a comparison between global variables that represent interesting fluid phenomena, we focus on assessing the enstrophy, symbolized by \(\Omega\), and the kinetic energy, represented by \(KE\). Tables 9 and 10 present the time-averaged enstrophy \(\Omega\) and kinetic energy \(KE\) obtained from multiple benchmarks. It is found that the flowfields obtained from _Source_ exhibit considerably lower values of \(\Omega\) and \(KE\), indicating loss of rotational and kinetic energy over time. In other words, the flowfields obtained from _Source_ result in significant decay of kinetic and rotational energy due to stronger dissipation effects resulting from the coarse grid simulations. In contrast, however, the time-averaged enstrophy and kinetic energy obtained from DPNN are in excellent agreement with the data obtained from the _Reference_ solver. Given the fact that the same numerical schemes are employed in the base solver for both the _Source_ and the DPNN (i.e., the solver embedded in the loop), the network thus acts as a forcing function that counters the strong numerical dissipation introduced by the _Source_ solver. This exercise, thus, clearly points towards the ability of the DPNN approach to yield predictions that preserve the local boundary layer profile along with the global wake dynamics in a way that results in lower numerical dissipation, as achieved by the _Reference_ solver. So far, we have demonstrated the reliability of the model at an individual prediction level. Additional evaluations for all the 50 previously unseen test samples have also been performed and are shown in Appendix A, supporting the conclusions drawn above.

\begin{table}
\begin{tabular}{c c c c}  & \multicolumn{3}{c}{Mean kinetic energy, \(KE/KE_{0}\)} \\ Testing samples no & _Reference_ & _Source_ & DPNN \\
28 & 1.000 & 0.924 & 0.991 \\
40 & 1.000 & 0.929 & 0.995 \\
46 & 1.000 & 0.927 & 0.996 \\ \end{tabular}
\end{table}
Table 10: Comparison of the normalized time averaged kinetic energy \(KE/KE_{0}\) for multiple testing samples

\begin{table}
\begin{tabular}{c c c c}  & \multicolumn{3}{c}{Mean enstrophy, \(\Omega/\Omega_{0}\)} \\ Testing samples no & _Reference_ & _Source_ & DPNN \\
28 & 1.001 & 0.926 & 0.985 \\
40 & 0.998 & 0.944 & 0.979 \\
46 & 0.999 & 0.988 & 0.997 \\ \end{tabular}
\end{table}
Table 9: Comparison of the normalized time averaged enstrophy \(\Omega/\Omega_{0}\) for multiple testing samples
### Insights through model predictions

A key aspect of our research is the model's ability to produce not only robust but also practical and valuable predictions. As we advance, we will undertake a comprehensive evaluation to further substantiate the model's reliability in producing results comparable to established literature. Specifically, we will compare non-dimensional parameters such as the Strouhal number (\(St\)) and spacing ratio (\(L/D\)), which are typically associated with the flow around two or more cylinders. This broadens our evaluation beyond individual testing samples, as was conducted in previous sections. To our knowledge, this marks the first occasion when such an investigation has been carried out for arbitrarily shaped bodies.

Despite the fact that the present dataset corresponds to flowfields obtained for flow past arbitrary configurations, one can still compare the distributions of the Strouhal number \(St\) based on the dominant frequency of one of the two downstream cylinders with the spacing ratio \(L/D\). The diameter in the spacing ratio \(L/D\) is simply the projected height of the upstream cylinder. This results in the distribution of the non-dimensional number (_i.e.,_ Strouhal number and spacing ratio), which allows for comparison with (Chen et al., 2020), who report the data obtained for flow past equi-diameter cylinders computed using an immersed boundary-based flow solver. Figure 19 compares the Strouhal number based on the dominant frequency of the downstream cylinder (upper) as a function of the spacing ratio \(L/D\). Figure 19 offers multiple insights, viz., besides the observable oscillations in the Strouhal number, a good overall fit between the solutions from the two different methodologies is noted. This illustrates that despite the arbitrary nature of the embedded bodies, the underlying frequency of vortex shedding (non-dimensionalized by the corresponding diameter and freestream velocity) for a certain body at a given spacing ratio exhibits similar patterns. Besides, it is also noted that the transition from low-frequency single bluff-body wake to deflected gap wake and subsequently to high-frequency, chaotic wake happens at nearly the same spacing ratio (as compared to Chen et al. (2020)). This observation is especially important given that the Strouhal number accuracy is crucially dependent on long-term recursive predictions, and deviations cause significant differences in the vortex shedding frequency (Hasegawa et al., 2020) _e.g.,_ from _Source_ (not shown here). Consequently, it can be remarked that the present DPNN-based predictive framework allow for useful and reliable predictions enabling long-term temporal accuracy. While multiple factors, such as shear layer growth, instabilities, gap flow, etc., contribute to the frequency

Figure 19: Strouhal number along with the spacing ratio for the arbitrarily shaped cylinders and its comparison with Chen et al. 2020. Strouhal number is evaluated based on the dominant frequency of the downstream (upper) cylinder obtained for 300 frames
of vortex shedding, further investigations are needed to probe the cause of the fluctuations.

In prior sections, we emphasized the necessity of evaluating the performance of the predictive model using test samples representative of different wake categories. In this section, we expand on this by presenting additional metrics of interest. This exercise allows for inferring greater physical insights from the predictions obtained so far, thus offering a unique perspective on solutions obtained from the predictive framework. Figure 20 represents the time-averaged enstrophy \(\Omega\) evaluated over the spacial domain for 300 testing frames for all testing samples, color marked by the category of the wake. It is found that there exist clear strata or bands of regions belonging to certain wake categories, _e.g._, single bluff-body, deflected gap, and chaotic wake, whereas these bands overlap for the fully developed quasi in/anti-wake categories. This indicates that the mean rotational energy of the flowfields is lowest for the single bluff-body wake and highest for the fully-developed quasi in-phase wake. This is not entirely surprising as the vortices shed from the upstream cylinder in the fully developed wake category penetrate the gap region while merging with vortices shed from the downstream cylinder, thereby increasing the rotational energy of the flow.

\[L_{2,\phi}=\sqrt{\frac{\sum_{i=1}^{nc}(\phi_{\text{Ref},i}^{j}-\phi_{\text{ baseline},i}^{j})^{2}}{n_{x}\times n_{y}}} \tag{10}\]

Figure 21 illustrates the time-averaged \(L_{2}\) norm error computed using Eq. 10 for each testing

Figure 21: Comparison of the \(l_{2}\) norm error of the predictions obtained from DPNN

Figure 20: Comparison of the enstrophy \(\Omega\) of the predictions obtained from DPNN color-coded by the wake category. The color bars represent the region within the maximum and the minimum \(\Omega\) for each wake category.
sample based on the velocity field up to \(t_{f}\)=300 frames. The variables \(n_{x}\), \(n_{y}\) are the numbers of computational cells in the streamwise and lateral directions, and \(\phi\) represents the quantity of interest. The samples are color-marked according to the wake category they exhibit. Unsurprisingly, it is found that the testing samples belonging to the chaotic wakes incur the highest error. Additionally, the testing samples corresponding to the deflected gap wake result in the lowest error from both DPNN and Source, which can be attributed to the stable gap flow (see Fig. 22). On the contrary, the testing samples exhibiting fully developed quasi in-phase wake result in the highest error from the _Source_. To probe further, in Fig. 23, we compare the normalized enstrophy \(\Omega\) for the flowfields obtained from _Source_ and the corresponding \(L_{2}\) norm error, each in ascending order. It is found that there exists a clear relationship between the dissipation resulting from the higher enstrophy in the flowfields and the corresponding error. This suggests that the error in predictions obtained from _Source_ is predominantly influenced by the numerical dissipation.

### Generalizability towards out-of-distribution cases

A crucial component towards determining the performance of a neural network-based predictive model is its generalizability to new, unseen datasets. As such, an evaluation with test data precludes inherent biases of the training dataset. Evaluations on out-of-distribution (OOD) datasets allows for estimating the robustness of the model while potentially identifying the model's inherent bias towards certain kind of data. In addition, such practices hint towards the applicability of the model to real-world data, which are likely to be OOD. We have carefully designed experiments for OOD datasets which allow us to perform the following:

Figure 23: Comparison of the enstrophy with \(L_{2}\) norm of the predictions from _Source_ solver, sorted in ascending manner

Figure 22: Comparison of the \(l_{2}\) norm error of the predictions obtained from _Source_
* **Three equi-diameter cylinders**. This setup mimics the existing one as shown in Fig. 1 with the exception of the geometric symmetry induced by three equi-diameter cylinders at spacing ratios of \(L/D\) = [1.5, 3.5, 5.0] with each cylinder having unit diameter. This serves as an interesting problem considering that the model was never trained for symmetry (the equi-diameter cylinders result in symmetric distribution about the x-axis). Additionally, this allows for a direct comparison of the wake category with that obtained by Zheng et al. (2016). For the chosen spacing ratio, the flowfield should exhibit deflected gap wake, anti-phase wake, and fully-developed in-phase wake, respectively.
* **Two side-by-side cylinders**. This layout chooses two side-by-side equi-diameter cylinders which is contrary to the present training setup to assess the performance of the model for an entirely alternative problem of practical relevance. The spacing ratios chosen for this case are \(L/D\) = [1.5, 2.5, 4.0], which allow for direct comparison with Bao et al. (2013) and should yield flip-flopping (or chaotic) wake, in-phase wake, and anti-phase wake, respectively.

#### 4.6.1 Test case 1: three cylinders

Figure 24 presents the velocity flowfields obtained from the DPNN approach and its comparisons with the data obtained from the _Reference_ solver at multiple time instances, _i.e.,_ three equi-diameter cylinders arranged in an equilateral triangle layout for the spacing ratio \(L/D\) = 1.5. The _Reference_ (or ground truth data) solver depicts a deflected gap flow wake, with the deflected gap flow being deflected towards the narrow wake region (upper cylinder), which bodes well with the observations by Zheng et al. (2016). It is found that the predictions obtained by the present DPNN approach also result in the same wake dynamics. In addition, the vortex shedding cycle result in an excellent agreement with the data obtained from the _Reference_ solver, which instils further confidence in the long-term temporal predictions.

Figure 25 presents the velocity flowfields obtained from the present approach and its comparison with the data obtained from the _Reference_ solver at the spacing ratios \(L/D\) = 3.5 and 5, up to a total time of \(t=500\Delta t\). The wake exhibited at these spacing ratios clearly depicts anti-phase and fully developed in-phase patterns, respectively, as reported by Zheng et al. (2016). The comparison shows that DPNN is able to preserve the overall wake structures remarkably well for long time horizons, irrespective of either in-phase or fully developed anti-phase wake.

Figure 24: Comparison of velocity flowfield obtained from DPNN and _Reference_ solver at multiple time instances for flow past the three-cylinder case at \(L/D\) = 1.5
#### 4.6.2 Test case 2: two cylinders

Figure 26 highlights the velocity flowfield obtained from the present approach and its comparison with the corresponding ground truth data at multiple time instances, _i.e.,_ two side-by-side cylinders at a spacing ratio \(L/D=1.5\). The data obtained from the _Reference_ solver exhibits flip-flopping wake pattern wherein the gap flow switches direction in a chaotic manner and augers well with the observation by Bao et al. (2013). This can be clearly seen in Fig. 26, _e.g.,_ the gap flow at \(t=1\Delta t\) (directed towards the lower cylinder) switches direction at \(t=30\Delta t\) (directed towards the upper cylinder). Addressing such a tumultuous shift poses a significant challenge for any predictive framework (Vlachas et al. (2018); Doan et al. (2021)). For instance, while DPNN accurately predicts the first two shifts, its effectiveness diminishes when it comes to long-term forecasting, particularly struggling with the accurate reconstruction of chaotic wake dynamics. Nevertheless, considering the fact that the proposed framework yields satisfactory performance for predicting the chaotic dynamics for mid-time horizons, this acts as a promising direction for future research.

Figure 27 presents the velocity flowfield obtained from the present DPNN approach and its comparison with the _Reference_ data for spacing ratios \(L/D=2.5\) and \(4\) for multiple time instances. The wake dynamics exhibit in-phase and anti-phase wake, respectively, which bodes well with the observations made by Bao et al. (2013). It is found that the predictions from the present approach compare really well with the data obtained from the _Reference_ solver, especially for the anti-phase wake at \(L/D\)=4 where it is seen to preserve the symmetry of the flow about the centerline \(x\)-axis. For the in-phase wake obtained at \(L/D\)=2.5, the flowfields diverge at long-time predictions _i.e.,_\(t=500\Delta t\). Therefore, given these observations, it is encouraging to n

Figure 25: Comparison of velocity flowfield obtained from DPNN and _Reference_ solver at multiple time instances for flow past the three-cylinder case at \(L/D=3.5\) and \(5\)
Figure 27: Comparison of velocity flowfield obtained from DPNN and _Reference_ solver at multiple time instances for flow past the two cylinder case at \(L/D=1.5\)

Figure 26: Comparison of velocity flowfield obtained from DPNN and _Reference_ solver at multiple time instances for flow past the two cylinder case at \(L/D=1.5\)
approach is successful in faithfully reconstructing wake dynamics for periodic flows (both in-phase and anti-phase) over extended time horizons. Although the accuracy of predictions for chaotic wakes is limited to mid-time horizons, the progress made so far offers a solid foundation for further refinement and optimization.

## 5 Conclusions

This work implements a differentiable physics-assisted neural network (DPNN) approach as a coarse-grained surrogate for unsteady incompressible flow past arbitrarily shaped bodies for a low Reynolds number of around 100. The bodies under consideration are arranged in a staggered equilateral triangle manner at multiple spacing ratio \(1.2\leq L/D\leq 5.5\) to produce wakes of varying nature. The main conclusions from the work can be summarized as follows:

1. The DPNN-based predictive framework outperforms the standard data-driven supervised learning approach and the _Source_ as a robust coarse-grained surrogate for unsteady fluid flows past arbitrary bodies. This has been qualitatively and quantitatively verified using multiple performance measuring indices that indicate accurate reconstruction of local (boundary layer) and global (wake dynamics) fluid phenomenon and for both physical quantities and Fourier space or reduced feature space.
2. The framework has also been evaluated for out-of-distribution test samples for cases that assess the ability of the predictive framework to preserve symmetry (_i.e.,_ anti-phase wake), stable gap flow (_i.e.,_ deflected gap wake), long-term temporal stability, chaotic dynamics (_i.e.,_ flip-flopping wake). It is found that the framework is largely wake agnostic, _i.e.,_ it renders a good prediction irrespective of the wake category. However, for the chaotic or flip-flopping wake, the framework can only accurately maintain the chaotic switching of the gap flow up to a certain mid-time range.
3. The \(L_{2}\) errors computed on the predictions obtained from DPNN and _Source_ are lowest for the deflected gap wake owing to its stable gap flow as well as long vortex formation length, the latter of which occupies the bulk of the computational space. The DPNN-based predictions seem to suffer from the onset of chaotic dynamics, which later manifest into inaccurate wake dynamics. The flowfields from _Source_ seem to suffer from large dissipation due to the underlying coarse computational grid, which results in the continuous decay of kinetic/rotational energy.
4. The Strouhal number distributions with the spacing ratio for the arbitrary bodies results in good agreement with the literature based on flow past equi-diameter bodies. Also, the transition in the wake category from a single bluff-body wake to a deflected gap wake and eventually to a chaotic wake occurs at nearly the same spacing ratios. Moreover, we discover similar categories and patterns in the wake exhibited by the arbitrary bodies, except that in-phase and anti-phase wakes yield in quasi in/anti-phase wakes. Additionally, the arbitrary nature of the bodies promotes an early vortex shedding by the upstream cylinder and precludes quasi in-phase wakes from occurring at all.
5. The enstrophy calculated based on the flowfields obtained from DPNN exhibit clear strata of test samples belonging to unique wake categories. This is particularly noticeable for single bluff-body, deflected gap, and chaotic wakes, whereas these delineations tend to overlap for the fully developed quasi in/anti-phase wake categories.
In view of these discussions, we postulate that such a strategy that incorporates both a low-fidelity solver as well as the neural network as a hybrid predictive framework derives merit from each component. Such a strategy paves the way for the low-fidelity solver to incorporate as much of the underlying physics as possible while the network learns the difference between low-fidelity solver and high-fidelity data. This specific attribute is key that makes such a framework more promising than purely data-driven supervised learning.

The present approach also shows potential merits as a generic reconstruction strategy for near-body boundary fluid properties. This is especially true for non body-conformal sharp interface immersed boundary approaches for high Mach and Reynolds numbers flows (Brahmachary et al., 2021) that are prone to suffer from heat flux reconstruction at the body boundary. Moreover, the present learning framework would also be interesting avenue for flows over rough surfaces (Lee et al., 2022; Jouybari et al., 2021) with practical relevance. In addition, the generalizability of the approach opens up the potential directions for learning the kinematics of dispersed spherical particles (Ozaki and Aoyagi, 2022).

## Funding

This work was supported by the European Research Council Consolidator Grant _SpaTe_ (CoG-2019-863850)

## Data availability statement

The data and code supporting this study are openly available in disclosed upon publication.

## Declaration of interests

The authors report no conflict of interest.
