We investigated the two-dimensional flow around a circular cylinder as well as a square cylinder. Figure 1 illustrates the cylinder geometries. We considered a fluid with a density ρ of 1 kg/m³ and an inlet velocity of 1 m/s. Dimensions and inlet velocity were chosen to yield a scaled problem with network inputs and outputs ranging between −1 and +1 as closely as possible. The resulting flow fields featured important phenomena, such as stagnation points, high gradient boundary layers, flow separation, and wakes. The two geometries exemplified curved and cornered obstacles that impose different flow separation mechanisms triggered at the rounded or cornered walls. We employed three different versions of the cylinder geometries as presented in Table 1. Geometries G1 and G2 were used to evaluate the effect of network architecture, network shape, and PINN methods for different cylinder shapes at an elevated Reynolds number. The large scale of the cylinders allowed more accurate and spatially detailed PINN predictions. This set up represented a channel flow and allowed comparison to a PINN study of Rao et al. [26]. Geometry G3 was selected to facilitate a comparison of the PINN method with the experimental results of Lyn et al. [17] as well as a PINN study of Ang et al. [1]. Due to the small scale of the square cylinder, blockage ratio effects of geometry G3 were negligibly small; however, the spatial accuracy of the PINN method was reduced.

The associated flow field is governed by the continuity Eq 1 and the two-dimensional incompressible Reynolds-averaged Navier-Stokes Eq. 2 ∂ u i ̄ ∂ x i = 0 (1) u j ̄ ∂ u i ̄ ∂ x j + u j ′ ∂ u i ′ ∂ x j ̄ = − 1 ρ ∂ p ̄ ∂ x i + μ ρ ∂ 2 u i ̄ ∂ x j ∂ x j (2) with Reynolds-averaged velocities and pressure u i ̄ , u j ̄ , and p ̄ , constant density ρ, and constant dynamic viscosity μ. In the RANS equations, the second order moments of the turbulent velocity fluctuations, u i ′ , and u j ′ , represent the Reynolds stresses τ i j ′ ∂ τ i j ′ ∂ x j = u j ′ ∂ u i ′ ∂ x j ̄ (3) The Reynolds stresses are modeled using the Boussinesq hypothesis τ i j ′ = − u i ′ u j ′ ̄ = μ t ρ ∂ u i ̄ ∂ x j + ∂ u j ̄ ∂ x i − 2 3 k δ i j (4) where μ _t is a turbulent viscosity representing the effects of the turbulent eddies on the mean flow, k is the turbulent kinetic energy, and δ _ij is the Kronecker delta. The approach of modeling turbulence via an eddy viscosity is based on the gradient diffusion hypothesis and, hence, assumes an alignment of the turbulent transport with the negative gradient of the mean flow. The application of an isotropic turbulent viscosity, μ _t , facilitates modeling of turbulent stresses in the same way as viscous stresses. A turbulence model is required to calculate the turbulent viscosity. Here, Prandtl’s mixing-length model was selected as applied in [7] due to its stability and robustness. For this model, the turbulent viscosity is estimated using a mixing length, l _m , that represents the size of the characteristic eddies as follows: μ t = ρ l m 2 S (5) where S represents the modulus of the mean rate of strain tensor. For the two-dimensional flow considered here, S reads S = 2 ∂ u ̄ ∂ x 2 + 2 ∂ v ̄ ∂ y 2 + ∂ u ̄ ∂ y + ∂ v ̄ ∂ x 2 (6) The mixing length l _m is calculated as follows: l m = m i n κ d , 0.09 d max (7) where d is the distance from the wall, κ ≈ 0.4 is the von Kármán constant, and d _max is a characteristic maximal length scale, here taken as the maximum wall distance of 1.0 m. Using the eddy viscosity approach, a resulting viscosity, μ _res , is calculated and deployed in the RANS equations μ res = μ + μ t (8)

Two physics-informed DL methods were compared. The corresponding models were trained by minimizing a composed loss function L ( θ ) : L θ = L b θ + L e θ (9) where L b is the loss on the boundary conditions, and L e is the loss for the RANS equations. The individual loss terms were calculated using the mean squared error: L b = 1 N b ∑ n = 1 N b | U b n − U ̃ b n | 2 (10) L e = 1 N e ∑ k = 1 3 ∑ n = 1 N e | ε k n | 2 (11) where N _b and N _e represent the number of points for which the boundary conditions and RANS equations were trained. U b n = [ u b n , v b n ] are the given velocity and pressure conditions for points n on the boundaries. U ̃ b n and U ̃ d n represent the output of the PINN at the corresponding training points and ε k n is the residual of the kth equation at point n.

For the first method, a traditional PINN was trained to predict the Reynolds-averaged velocity and pressure fields. The PINN reads as follows: N N θ x y = u v p (12) where θ are the trainable parameters (weights and biases) of the NN. Hence, for the traditional PINN approach, the model maps the input coordinates x and y to the output solution field components u, v, and p. The Navier-Stokes Eqs 1, 2 were used as the residual functions as follows: ε k = 1 = ∂ u i ̄ ∂ x i (13) ε k = 2,3 = u j ̄ ∂ u i ̄ ∂ x j + 1 ρ ∂ p ̄ ∂ x i − μ res ρ ∂ 2 u i ̄ ∂ x j ∂ x j (14)

For the second method, a mixed-variable approach as presented by Rao et al. [26] provided a system of equation that can be learned more easily by the PINN. The correspondingly simplified RANS equations read as follows: u j ̄ ∂ u i ̄ ∂ x j = 1 ρ ∂ σ i j ∂ x j (15) Using this method, the stress tensor σ _ij is utilized σ i j = − p δ i j + μ res ∂ u i ̄ ∂ x j + ∂ u j ̄ ∂ x i (16) The pressure p can be calculated as the trace of σ p = − 1 2 σ i i (17) As suggested by Rao et al. [26], the stream variable ψ was used to represent the velocity field. This automatically satisfies the continuity Eq 1 and amplified the learning success of the NN. The velocity components can be retrieved using u v 0 = ∇ × 0 0 ψ (18) where ∇ is the Nabla operator. A feedforward NN was trained to predict ψ, p, and the stress components σ _xx , σ _xy = σ _yx , and σ _yy . The network takes the coordinates x and y as inputs and passes the velocity components, the pressure, and the stress components back. The model reads as follows: N N θ x y = ψ p σ x x σ x y σ y y (19) In agreement with Eqs 15–17, the residuals ɛ _k are computed as ε k = 1,2 = u j ̄ ∂ u i ̄ ∂ x j − ∂ σ i j ∂ x j (20) ε k = 3,4,5 = − p + μ res ∂ u i ̄ ∂ x j + ∂ u j ̄ ∂ x i − σ i j (21) ε k = 6 = p + 1 2 σ i i (22)

For the traditional approach, the two input neurons were fed into a single fully connected feed-forward network with three output neurons as determined by Eq 12. For the mixed-variable approach, the two input neurons were fed into a single fully connected feed-forward network with five output neurons as defined by Eq 19. Figure 2 shows the two methods that were compared here.

For all networks, a tanh activation function was deployed. In a first step, the PINNs were trained using the Adam optimizer for 100,000 iterations with a decaying learning rate. The initial learning rate was set to 0.001 and was reduced by a factor of 0.9 every 2,000 iterations. In a second step, the models were trained using the L-BFGS optimizer under the predefined default settings. No loss weigthing factors were defined. The training was conducted using the tensorflow-based library DeepXDE [16] on a 16 GB NVIDIA Quadro RTX 5000.

The capability of the PINN to represent the solution of the RANS equations under consideration of the mixing length model was evaluated and compared with CFD simulations for geometries G1 and G2. For the CFD calculations, the element-based finite volume flow solver MAYA included in the Simcenter software was used. The MAYA solver provided the mixing length model with Van Driest damping. Unstructured meshes featuring prism layers on the cylinders were defined. The second order upwind discretization scheme was deployed. The incompressible calculations were carried out until reaching the minimal achievable residuals. A grid sensitivity analysis was carried out to verify the models. For the grid study, three grids were generated using a refinement factor of 1.3. Excellent agreement between all three grids was found and, consequently, the reference solution was concluded to be grid independent.

The PINN predictions for geometry G3 were compared with experimental results of Lyn et al. [17] who conducted two-component laser-Doppler measurements for a square cylinder inside a closed water channel at a Reynolds number of 2.14 × 10⁴. The time averaged data were obtained from the ERCOFTAC validation database [4].

This study, comprising three parts, compared the traditional method with the mixed-variable method, investigated effects of network shape and size for different geometries, and compared the predictions with measurements. For both PINN methods, several network architectures were compared. The network sizes were chosen to yield the same number N _θ of trainable parameters θ (weights and biases) inside the network, where N _θ is defined as follows: N θ = N i n N n + N n 2 N l − 1 + N out N n + N n N l (23) Here, N _out represents the number of output neurons and N _in the number of input neurons. The number of neurons per layer, N _n , and the number of hidden layers, N _l , were set to yield a progressing series of shape factor values, λ. The shape factor was defined as follows: λ = N n N l (24) A low value of λ represented a deep but narrow network architecture featuring comparatively few neurons per layer and many hidden layers. An elevated result for λ corresponded to a shallow but wide network consisting of a high amount of neurons per layer and few hidden layers.

In the first part of our study, the traditional method was tested using geometry G1. This part documented limitations of the traditional PINN method when applied to high Reynolds number flows. Table 2 lists the tested networks of the first part. In the second part, the best fitting network shapes and sizes for mixed-variable PINN models were evaluated and the results for geometry G1 and G2 were compared. The purpose of this part was to find appropriate network shapes, outline differences between the two geometries, and document capabilities and limitations of unsupervised mixed-variable models when applied to high Reynolds number flows. A total of five network architectures were deployed with corresponding shape factors ranging from 4.08 × 10⁻¹ to 1.34 × 10². Table 2 lists the evaluated networks. Furthermore, after selecting the network shapes giving the most accurate results, several networks featuring a constant shape factor and a progressing number of trainable parameters θ were defined and tested. Different shape factors were chosen for the circular cylinder and the square cylinder. Table 2 lists the networks evaluated in the second step. 10,000 randomly distributed training points were used on the boundaries as well as inside the flow domain for geometries G1 and G2. This represented 1,250 points per m² of the flow domain and 754 points per m of the boundaries. The number of points was evaluated in a prior sensitivity analysis conducted for the circular cylinder G1 using the mixed-variable PINN N50L30. The evaluation revealed no increase in accuracy by increasing the number of points, as shown in Figure 3. The density of points deployed here was 3.87 times greater than used in a study of Ang et al. [1] who reported no change in accuracy beyond 323 points per m² of the scaled flow domain for a laminar flow field.

In the third part, geometry G3 was considered and an unsupervised PINN was compared with the experimental results of Lyn et al. [17]. A search for the most suitable network shape was conducted in the same way as described above for part two. As the general method and behavior of the PINN was equivalent, the detailed results of this study are not presented here. Model N50L30 was identified as the most suitable architecture for geometry G3. The corresponding settings of the network are listed in Table 2. To account for the small scale of the cylindrical body contained in geometry G3, the number of training points inside the domain was increased to 80,000. The purpose of this part was to get an impression of the validity of the mixed-variable PINN method when no labeled training data are provided.

The traditonal PINNs failed to capture the flow field, independent of network architecture. Figure 4 displays comparative predictions for the velocity field around the circular cylinder. As a reference, the CFD results are shown as well. Several models are exhibited that represent high, moderate, and low shape factors. From the shown models, N31L76 and N268L2 both failed to capture the reference flow field accurately. Model N50L30 provided the best representation of the flow field but the wake length was overestimated and the velocity of the high speed region at the lateral sides of the cylinder was underestimated and the region separated from the cylinder wall. Overall, the flow field was not captured accurately by the traditional PINN method, neither quantitatively nor qualitatively.

The mixed-variable approach was observed to be sensitive to the network shape. Too deep as well as too wide networks failed to capture the flow field. Figures 5, 6 display predicted axial velocities of several models together with comparative CFD data. From the shown models, N31L76 and N268L2 represent low and high shape factors, respectively. Both architectures failed to capture the reference flow field accurately. For the circular cylinder, model N50L30 obtained the most accurate predictions, capturing the stagnation region, the thin boundary layer on the cylinder, as well as the wake. However, the wake length was overestimated compared with the reference solution and the high speed region on the lateral sides of the cylinder was deformed. The length of the wake defined as the downstream distance behind the cylinder until 60% of the inlet velocity was recovered was predicted to be 1.32 m while the reference CFD simulation yielded a wake length of 0.27 m. Correspondingly, the PINN predicted a recirculation length of 0.16 m while the CFD computation featured a shorter recirculation region of 0.02 m. For the square cylinder, model N135L5 exhibited the most accurate results. The stagnation region, the cylinder wake, and the small area of flow separation and recirculation at the lateral sides of the cylinder was captured. However, the size of the low velocity stagnation region upstream of the cylinder was overpredicted when compared with the reference solution. For the square cylinder, the wake length was predicted more accurately by the mixed-variable PINN than for the circular cylinder as the PINN featured a wake length of 2.72 m while the CFD solution gave a length of 2.49 m. The recirculation length as predicted by the PINN was 0.52 m and the reference simulation featured a recirculation length of 0.38 m.

Figure 7 compares axial velocities of the different mixed-variable PINN models with the CFD calculations for the circular cylinder. The results are shown along several vertical lines. At x = 0, the high gradients as well as the overshoot of the developing boundary layer were best captured by N50L30 and N81L12. All models predicted a prolonged wake behind the cylinder. This resulted in a greater velocity deficit at x = 1 m and x = 2 m, accompanied by a higher free stream velocity at the lateral sides due to volume flow and momentum balances. Model N50L30 captured the reference flow most accurately. Figure 7 also displays the predicted velocity close to the square cylinder. Except for model N268L2, all models showed close agreement. At x = 0 m, model N135L5 accurately captured the recirculating flow on the cylinder surface featured by the reference solution. Furthermore, the velocity overshoot of the developing and separating boundary layer was captured. However, the thickness of the recirculating flow was slightly overestimated and, consequently, the velocity overshoot of the separating boundary layer was overestimated as well. All models accurately predicted the wake behind the square cylinder.

The trends observed for the axial velocity were also observed for the pressure field. Figures 8, 9 exhibit the pressure predictions and comparative CFD data. For the circular cylinder G1, model N31L76 predicted pressure values close to zero and did not capture the pressure field. For the square cylinder G2, model N31L76 predicted inaccurate results as well. Similarly, model N268L2 did not predict the pressure field accurately. For the circular cylinder, the positive pressure in the stagnation region of the reference solution was underestimated and a moderate negative pressure was predicted around the cylinder. For the square cylinder, model N268L2 predicted a dislocated region of positive pressure at the upstream side of the cylinder. For the circular cylinder, the best fitting model N50L30 captured the positive pressure region at the upwind side of the cylinder as well as the negative pressure at the lateral sides. However, the positive pressure on the downwind side was not predicted. For the square cylinder, the best fitting model N135L5 captured the positive pressure of the stagnation region as well as a negative pressure field around the square cylinder. However, size and magnitude of the positive pressure region as well as the negative pressure region at the upstream corners were underestimated.

The pressure values on the cylinder surface are shown in Figure 10. The deep but narrow model N31L76 predicted a vanishing pressure distribution around the entire cylinder. All models failed to predict the positive pressure region on the downstream side of the cylinder that was present in the CFD reference data. From all mixed-variable models, N50L30 provided the best representation of the reference solution. The minimum pressure peaks were more pronounced and the pressure on the downwind side of the cylinder was increased. Furthermore, the asymmetry associated with the other architectures was reduced. Figure 10 also exhibits a comparison between the pressure distributions on the square cylinder surface. Models N31L76 as well as N268L2 predicted inaccurate results. All other models captured the shift from positive to negative pressure. All networks underestimated the negative pressure maxima at the square edges that were featured by the CFD reference data. Models N50L30 and N81L12 better predicted the sharp shift from positive to negative pressure while model N135L5 more accurately matched the positive pressure magnitude at the upstream side of the square.

For the circular cylinder G1, the best fitting PINN N50L30 predicted a drag of 0.272 N while the reference CFD yielded 0.135 N. The overestimation of the drag force correlated with an underestimation of the pressure in the wake at the downwind side of the cylinder when compared to the pressure distribution of the CFD reference solution. For the square cylinder, a drag force of 0.242 N was predicted by the best fitting model N135L5 while the reference solution featured a drag of 0.473 N. The underestimation of the drag correlated with an underestimation of the pressure values on the upwind side of the square.

The correlation of shape factor and mean squared error (MSE) is shown in Figure 11 for the circular cylinder G1. As seen, the axial velocity in the x-direction achieved the highest overall MSE. The relationship between MSE and shape factor was similar for all flow field variables. As λ tended to zero, the error increased markedly. For the networks with shape factors above 1.67 × 10⁰, the error grew with a decreasing slope. The obtained data suggests an optimal network shape for networks with shape factors close to 1.67 × 10⁰. For the best fitting model, N50L30, the losses for the two momentum equations were 9.23 × 10⁻⁶ and 6.99 × 10⁻⁶ for the x- and y-directions, respectively. The achieved losses indicate that the governing equations were satisfied. The correlation between MSE and shape factor was different for the two cylinder geometries. Figure 11 also displays the MSE of the different network shapes for the square cylinder G2. As seen, the error tended to a minimum for shape factors between 6.75 × 10⁰ and 2.70 × 10¹. For the square cylinder, the best shape factor was greater than for the circular cylinder. Furthermore, close to its minimum, the growth of the MSE was less and featured a broader region of minimal MSE. However, the models featuring the lowest NMSEs, i.e., N50L30 and N81L12 also featured artificial velocity peaks at the upstream corners of the square and, hence, were discarded. For model N135L5, featuring a shape factor of 2.70 × 10¹, the achieved losses for the momentum equations were 3.23 × 10⁻⁵ and 2.47 × 10⁻⁵. The minimal global MSEs for u, v, and p were lower for the square cylinder while the corresponding losses of the best model were lower for the circular cylinder.

Figure 12 exhibits the MSEs obtained by the different network sizes featuring the best fitting shape factor 1.67 × 10⁰ for the circular cylinder. As seen, the network N20L12, featuring 4,780 trainable parameters, yielded the greatest errors. The corresponding predictions of the velocity field were inaccurate and comparable to that of model N268L2 shown in Figure 5. The other models, featuring a greater number of trainable parameters, yielded similar errors. Model N50L30, featuring 74,350 trainable parameters, was well within the size independent range and neither a reduction nor an increase of the number of trainable parameters led to a further reduction of the global MSE, indicating the superiority of this architecture. Figure 12 also exhibits the results of the evaluation of the network size for the square cylinder. In contrast to the circular cylinder, the global MSEs decreased over the entire range of tested network sizes. For model N54L2, featuring 3,402 trainable parameters, the highest errors were recorded, and the corresponding predictions were comparable to that of model N31L76 shown in Figure 6. Increasing the number of trainable parameters from 74,520 to 133,326 led to a reduction of 20%–30% of the global NMSE for model N162L6. However, the predictions showed no qualitative improvement and featured artificial velocity peaks at the upstream corners of the square. Hence, the larger model did not lead to more favorable predictions of the flow field.

Figure 13A) exhibits the axial velocity field for geometry G3 as predicted by model N81L12. As seen, the PINN captured a stagnation region of lower velocity at the upwind edge of the square as well as a wake at the downwind side of the cylinder. Furthermore, the high velocity areas extending from the upwind corners of the square were captured. Overall, the flow field was well captured. A comparison of the axial velocity as predicted by the PINN and the values obtained from the time averaged laser-Doppler measurements is shown in Figure 13B) for several vertical lines in the wake of the square cylinder. In the PINN predictions, the velocity deficit inside the wake was less pronounced at x/c = 1.0 and, due to momentum and mass conservation, a lower velocity was predicted at the outer extent of the shear layer at y/c = 1.0. Overall, there was a favorable agreement between the PINN predictions and the measured results, considering the small scale of the cylinder, the three-dimensional nature of the reference flow, and the absence of any simulated or measured data during training of the neural net.

We compared different network architectures using mixed-variable physics-informed deep learning and traditional PINNs with CFD and measured reference data. The models were trained to solve the two-dimensional RANS equations for a turbulent flow around a circular cylinder and a square cylinder. The mixing length turbulence model was deployed. The main findings are summarized as follows:

• For the elevated Reynolds number flow considered here, the superiority of the mixed-variable approach of Rao et al. [26] was confirmed. The traditional PINN method failed to capture the flow field accurately, independent of the network architecture.