Systematic Assessment of the Two-Step, One-Way Coupled Method for Computational Fluid Dynamics

Computational fluid dynamics (CFD) is a powerful tool for complex engineering design problems, which has grown rapidly in importance during the last decades. However, limited computational capacity is often a problem that forces the CFD engineer to use various approximations to generate a solution within the allowed time limits of the project in question. Hence, there is a need to find accurate and time-saving approximations for practical CFD design, particularly in engineering design optimization, where it is necessary to run many consecutive CFD solutions with different values for the design variables to find the optimum design [1]. A common simplification is to reduce the computational domain to a bare minimum and to assume a simplified inlet condition, e.g., a plug profile, to represent the rest of the computational “universe” in the upstream direction. However, a complex upstream geometry, e.g., swirler guide vanes or strong curvature, often makes it difficult to specify realistic boundary conditions. For the outlet condition, the existence of downstream geometrical complexity is less problematic since many flows have a parabolic character, which means that the equations are effectively “parabolized” so that information only propagates in the streamwise direction [2]. Hence, for the downstream boundary condition, it is usually enough to use a boundary condition that is numerically well behaved and satisfies the underlying conservation laws, e.g., a zero normal derivative condition.

For situations where the upstream flow is too complicated to motivate a plug flow profile, it is common to divide the complete physical domain into two parts and to use the result from the upstream part to define the boundary conditions for the downstream part. This method, which can be described as a Two-Step, One-Way (TSOW) coupled method, gives a strong coupling from the upstream solution to the downstream solution but no coupling and vice versa. The corresponding time-saving in engineering design optimization can be very significant compared with a direct solution of the complete problem for two reasons: first, the effort to solve the two subproblems can be smaller than that for solving the complete problem, and second the same upstream solution can be used in all design iterations, meaning that only part of the problem needs to be solved after the first iteration. By using this method, one is implicitly assuming that the flow equations are effectively parabolized. To the knowledge of the authors, this assumption has never been assessed in a systematic way.

The existing literature associated with methods to decrease computational complexity mainly deals with the simplification of features of a particular geometry [3–5] or with a 2D reduction from 3D [6,7]. Janssen et al. [4] also explore subdivision of the domain, but for the simplified version of their domain of interest. In a paper by Toffolo et al. [8], the authors exclude one part of the geometry and replace it with an equivalent acoustic excitation effect, but still most of the simplifications of geometry are based on symmetry. A very interesting paper by Schlüter et al. [9] demonstrates the reduction of computational costs by splitting the domain into three parts that use different solvers, large eddy simulation (LES) and Reynolds averaged Navier–Stokes (RANS), thus allowing the simulation of the full flow path in a gas turbine for a simple test case. The difference between the method discussed in Ref. [9] and the TSOW method is that in the former the interfaces are fully coupled, which requires more work than one-way coupling. In another paper by Hettel et al. [10], focusing on the flow in a catalytic converter, the authors exclude the monolith (the middle part of their domain that experiences laminar flow) from their simulations. The two domains are then solved in succession with the outlet velocity of the upstream domain used to calculate the inlet velocity of the downstream domain. Then, a calculation of the static pressure takes place in the downstream domain which feedbacks to the outlet pressure of the upstream domain. Similar to Ref. [9], there is a two-way coupling of the domains in Ref. [10].

Nonetheless, there are also studies that favor the increase of the computational domain with the addition of an upstream geometry. In Ref. [11], the authors focus on the aerodynamic study of vehicles, where the simulation domains used are box-shaped to encompass open road conditions. However, the wind tunnels used for the experimental studies experience different conditions (e.g., blockage ratio) to open roads. The authors of Ref. [11] are in favor of adding the upstream part of the wind tunnel in order to have a closer match to the experimental studies without the need for a correction factor. A study that considers the effect that an upstream domain has on the validation of a CFD model against experiments was conducted by Longest and Vinchurkar [12]. In Ref. [12], they use two domains, one only with the double bifurcation geometry and one that extends further upstream to include a T-junction of the experimental particle delivery system. The transitional and turbulent flow effects of that upstream geometry are found to play a significant role in particle deposition profiles on that branching respiratory model. Application of the TSOW method is such a test case as Ref. [12] could prove beneficial in approximating the upstream domain flow effect without including it explicitly in the simulation of the domain where most of the interest lies. Nevertheless, there is a gap in methodically assessing this method and therefore a lack of information about its validity. Hence, it is of interest to investigate the assumptions in the TSOW method and the corresponding impact on the accuracy, compared to a direct solution over the full domain.

A swirling flow, including the swirl generator, is an excellent test case for the assessment of the TSOW method, since it is common in technical applications, it is strongly affected by the detailed geometry of the swirler, and there is no symmetry in the flow that can be exploited to reduce the geometry. The domain used in this study is the CeCOST swirl burner [13], which is an Advanced Environmental burner [14], relevant to gas turbines (see Fig. 1). The name stems from the design aiming to produce short flames in lean conditions and minimize pollutant emissions. The swirler is a crucial part of the whole setup but does not need to be modified in a typical engineering design situation where the detailed geometry of the downstream combustion chamber and the operational parameters are to be optimized. A 3D model of the swirler and a representative surface mesh is shown in Fig. 2. The swirler has a very complicated geometry and thus is quite problematic in terms of meshing/modeling. In particular, the tip of the swirler and the joining area between the guide vanes and the main duct are difficult to mesh properly and require many cells to eliminate problems with extreme angles and skewness. In transient simulations, non-uniform meshing translates into big jumps in the Courant number, which can lead to unstable solutions or divergence. If computational resources were unlimited, one could refine the whole domain to a similar size as that at the tip of the swirler, which would eliminate the problems connected to the Courant number and sudden jumps in cell size. However, for the average CFD engineer, this is not an option. The resolution of the boundary layers with so-called prism layers in such a domain is another complicating issue. As an example, the interface between swirler blades and body, as well as between the blades and the Mixing Tube (MT), in the chosen geometry needs careful attention. Achieving the right y + value on all the walls in such geometries can result in large jumps in cell size and low-quality cells in the vicinity of the prism layer.

This paper consists of an initial section that describes the numerical setup and the chosen geometry. The next section describes the methodology for the investigation of the validity of the TSOW method. This is followed by a systematic description of the results, both for the reference case with a fully coupled solution and for two versions of the Two-Step method. Finally, the results are discussed in the context of the validity of the TSOW method, followed by the conclusions that can be drawn from the study.

For the sake of clarity and to eliminate the influence from the combustion model, only isothermal non-reacting flow is studied. The commercial software STAR CCM+, version 15.06.008-R8 [16], based on a finite volume discretization was used for all the simulations. A segregated flow model, which solves pressure and velocity sequentially and iteratively, was chosen for all the simulations. A tight pressure-velocity coupling is achieved using the Rhie-Chow pressure-weighted interpolation for the flow coefficients at cell faces [17]. Second-order approximations were used for the spatial discretization of all equations. Time derivatives that occur in unsteady Reynolds averaged Navier–Stokes (URANS) simulations were discretized with a first-order time scheme, which was deemed sufficient since only the final steady-state result was of interest. Due to the complexity of the geometry, an unstructured polyhedral grid was used. One advantage of the polyhedral grid compared to a tetrahedral grid is that it obtains the same accuracy with fewer cells [16,17].

The full geometry of the CeCOST burner consists of a premixing section with a diameter of 54 mm, followed by a flow straightener with a honeycomb-porous material section (see Fig. 1). After the flow straightener follows a swirl register comprising four guide vanes with complex shape, followed by a mixing tube, 48 mm in diameter, that ends with a sudden expansion into the combustion chamber which has a height of 400 mm and is 100 mm wide. For the TSOW solution, the geometry was split in two ways, as indicated in Fig. 3. For the upstream flow, the geometry was split at the intersection between the MT and the combustion chamber (the “dump plane”). For the downstream flow, the geometry was split right after the end of the swirler blades.

The solutions were obtained using an unsteady RANS (URANS) approach, motivated by poor convergence with a steady RANS approach for the full solution. In addition, the downstream cases, which were better conditioned due to the absence of the swirler, were also solved with a RANS approach. The turbulence was modeled using the same k–ω shear stress transport (SST) model for both RANS and URANS simulations [18,19]. The walls had a no-slip condition via blended wall functions, which require y + values for the first prism layer in the range 1 < y + < 30. For the unsteady simulations, an implicit unsteady model was used with a timestep size of 10⁻⁴ s, and a total simulation time corresponding to 7 flow-through times. For the full geometry simulation, the mean Courant number on a few cells at the swirler tip was around 20, which is also the highest value of the Courant in the simulation. In the rest of the MT, the highest value of the Courant number was less than 5, and for the combustion chamber (area of interest), it was less than 1. The outlet boundary condition was pure outflow for all the simulated cases besides the partial upstream domain that only includes the swirler and the MT (more about boundary conditions in the section Methodology). The number of cells ranged from 2 million for the partial domains, to 4 million for the full domain.

The full TSOW solution consists of two parts, one upstream and one downstream solution, where the upstream geometry overlaps the downstream geometry by 2.5 MT diameters. The TSOW method is based on two implicit assumptions: the first is that the solution in the upstream part of the domain will produce the same result as a full domain solution, except close to the outlet boundary where a simplified boundary condition must be used since the exact conditions are unknown a priori; the second implicit assumption is that the solution in the downstream part of the domain will produce the same result as a full domain solution, if the inlet boundary conditions are perfect. Moreover, since there will always be small discretization and roundoff errors, it is assumed that the effect of these errors in the boundary conditions will be negligible if the estimated error in the solution is low enough. The methodology described below is designed to assess the validity of these assumptions.

Figure 3 illustrates the geometric configurations used in this study. The reference full geometry solution is indicated as FS, and the two assessment cases are indicated as TC1 and TC2. The figure also shows the positions for the inlet boundary conditions and boundary condition transfer in the assessment cases TC1 and TC2. Case TC1 was designed to assess the validity of the second implicit assumption mentioned earlier. For this case, the upstream boundary conditions were extracted from the full solution, which presumably should yield better boundary conditions than a partial solution if the grid resolution is the same in both cases. For case TC2, which was used to assess the validity of the first assumption and to demonstrate the full TSOW method, the boundary conditions were extracted from the solution in a partial upstream domain (see Fig. 3). The overlap between the upstream and TC2 geometry limits the influence from the outlet boundary conditions of the upstream solution in the position where the inlet boundary conditions for TC2 will be extracted. Moreover, for TC2, the boundary will be sufficiently moved in the upstream direction to allow some smoothing of the flow before it enters the main region of interest. In this case, the discretization error in the boundary conditions would be different than those in the FS case and will thus differ slightly from the TC1 case.

The validity of the first assumption behind the TSOW method was assessed by a comparison of the FS to TC2 upstream. The question in this case was how sensitive the upstream solution is to the lack of feedback from the downstream flow and if the partial domain solution would yield the same result as the full domain solution, except very close to the outlet of the partial domain, where the approximate boundary conditions would introduce some errors.

The validity of the second assumption behind the TSOW method was assessed by a comparison of the FS case to the TC1 case. Moreover, the sensitivity of the solution in the downstream domain to small errors in the inlet conditions was assessed by comparison to the TC2 case. Two questions had to be answered: the first if the case with almost perfect boundary conditions (TC1) would reproduce the results from the FS; the second if the errors in the boundary conditions for case TC2 would lead to a significantly different solution than in TC1 or they would be damped, leading to similar results.

The inlet boundary condition for the FS and the upstream domain for TC2 was a plug flow profile with an axial velocity of 6.5 m/s, which is a reasonable assumption due to the presence of the flow straightener on the upstream side of the domains (see Fig. 1). The turbulence boundary conditions for both cases were estimated in the same way according to the ERCOFTAC Best Practice Guidelines [20]. To minimize the influence on the results from the numerical mesh, the meshing for the different cases was kept as identical as possible with the built-in mesher of Star CCM+ [16]. The downstream domains for TC1 and TC2 were identical and could therefore use the same mesh. However, the meshes for the upstream domain in TC2 and FS did not overlap perfectly due to the way the automatic mesher works. This will lead to different interpolation errors for the two cases when the results are interpolated to the inlet surface mesh for the downstream domains of TC1 and TC2. The interpolation of results from the two upstream solutions to the inlet mesh of the two downstream solutions was done with the same numerical accuracy (second order) as the CFD numerical scheme. Hence, the combined error from the discretization and interpolation can be expected to be second order. The boundary values that must be extracted from the upstream solution are the velocity components, the turbulent kinetic energy, and the specific dissipation rate over the whole inlet surface.

All the simulations (FS, TC1, and TC2) were started from the same constant value initial guess. All the solutions that were used for the validation of the assumptions in the TSOW method were obtained with a URANS model. Moreover, the downstream flows in the TC1 and TC2 cases were solved with a RANS model, to demonstrate the potential for speedup. Before the decision to use the unsteady RANS approach for the assessment, a considerable effort was made, without success, to improve convergence by changing control parameters for the steady RANS solution. One possible explanation for the poor convergence could be periodic vortex shedding from the guide vanes, which due to its deterministic nature is not represented in a Reynolds averaging turbulence model. This explanation is supported by the fact that it was possible to obtain converged solutions for the two downstream flows where the swirler was excluded from the geometry.

The FS, TC1, and TC2 cases were run for 2.5 s physical time, corresponding to 7 flow-through times (the ratio of the total volume of the solution domain to the volumetric flowrate through the inlet boundary). The TC2 upstream case was run for 0.25 s, corresponding to 7 flow-through times. In all URANS cases, the collection of averages started after 4 flow-through times. To assess the convergence of the unsteady cases, monitoring points were defined at areas of interest. For the FS, this was at the start of the combustion chamber, and for the upstream solution of the TC2 case, they were distributed along the length of the MT. Convergence is assumed to be obtained when both the cumulative mean and the cumulative root mean square (RMS) values reach a steady-state for variables of interest, e.g., velocity components. Figure 4 shows the axial velocity convergence history for the FS case. The y-axis shows the normalized axial velocity U/U_max, where U is the cumulative average of the axial velocity and U_max is the maximum value of the cumulative average. The lower row shows the normalized RMS values at the same positions. Notice that the lower limit of the y-axis is different for the mean and RMS graphs. The x-axis is the normalized time t/T_F, where t is physical time and T_F is the flow-through time.

The apparent order of the scheme based on the data in Table 1 was 1.3, which is a reasonable result since the formal second order of the numerical scheme can be downgraded by grid non-uniformities and effects of the wall functions, which require that the first grid point at the wall must be at a certain distance from the wall for the wall functions to work properly. Based on these estimates, the 4M cell grid was used for the FS case. The error estimates compared to an exact solution are tabulated in Table 1. The meshes of all other cases are based on this mesh so the errors associated with them are assumed to be in the same order of magnitude.

Close inspection of the transient results for the axial velocity field shows that there are signs of periodic vortex shedding from the guide vanes (see areas with negative axial velocity near the guide vanes in Fig. 5). An assumption that the asymmetry in the instantaneous flow field comes from the periodic and symmetric vortex shedding is supported by the almost perfectly symmetrical shape for the time average of the axial velocity field. This flow behavior could be the reason for the poor convergence with the steady-state approach since the vortex shedding is a deterministic phenomenon with a characteristic frequency [22] while steady RANS models are designed assuming that the fluctuations are stochastic and therefore will not properly account for the effect of periodic vortex shedding. This could potentially increase the numerical stiffness of the steady-state approach via an underestimation of the Reynolds stresses that would make the non-linear convective terms more dominant.

Figure 6 shows field plots of the time-averaged velocity components (tangential, axial, and radial) of the FS in a vertical cross section of the domain. As can be seen from the figure, there are visible signs in the MT of the helical wakes generated by the guide vanes for all three velocity components. Moreover, notice the symmetry of the velocity components around the MT center line, which comes from the symmetry of the four-bladed swirler. As a side note, for an odd-bladed swirler, one would not expect a symmetric flow in a plane through the axis. Notice the negative axial velocity near the guide vanes, in the lower part of the MT of the middle plot representing a recirculation zone. Further optimization of the swirl register could probably eliminate the periodic vortex shedding from the guide vanes, but this is out of the scope of the present paper.

It is also clear that the flow is still developing at the inlet position of the TC1 and TC2 cases (dashed line in Fig. 6). Consequently, it would not be correct to assume vanishing normal derivatives in the upstream flow simulation in the TSOW method where that surface is the outflow boundary. However, in the TSOW implementation in this paper, this problem was avoided by letting the upstream domain overlap with the downstream domain all the way to the beginning of the combustion chamber, thereby moving the outlet boundary approximately 2.5 diameters away from the inlet of the downstream domain (see Fig. 3).

Figure 7 shows the axial and tangential velocity fields for the FS and TC2 upstream cases. Overall, the two solutions agree well, supporting the assumption that the influence from the downstream boundary condition in the TC2 upstream case is limited to a small distance upstream from the outlet boundary. This means that the results that are extracted from the surface defined by the dashed line in Fig. 7 as boundary conditions to the TC2 case should be minimally influenced by the downstream condition. Moreover, the comparison with the FS case shows that the results from the solution on a partial domain (the TC2 upstream case) are very close to the FS, thereby supporting the first assumption behind the TSOW method (see Methodology).

Before the discussion of the approximate downstream solutions in TC1 and TC2, it is necessary to discuss the inlet boundary conditions in those cases. Figure 8 shows the velocity fields for the FS case on a horizontal surface at the position of the inlet of the TC1 and TC2 cases. The figure also shows the results at the same position for the TC2 upstream solution and the interpolated results that were used as inlet boundary conditions for the TC1 and TC2 downstream solutions.

It is clear from Fig. 8 that the velocity components in both the FS and the TC2 upstream cases exhibit a very complex pattern due to the distortion of the flow by the guide vanes. Notice the axial velocity maximum at the axis and the symmetric pattern for the maxima and minima of the radial and tangential velocity. For the FS and TC2 upstream cases, the results in the figure were interpolated with second-order accuracy from the corresponding solution and can be seen to agree well with each other. The data were also tabulated so they could be used for the definition of the boundary conditions in the CFD code, which included a second interpolation. As a result of the two steps of interpolation, small errors have been generated, which can be seen as a certain amount of noise in the values that were used as boundary conditions in the TC1 and TC2 simulations of the downstream flow (see plots marked TC1 and TC2 in Fig. 8).

As an initial assessment of the agreement between the approximate solutions and the FS case, the swirl number (Eq. (1)) has been plotted for all cases in Fig. 9. The x-coordinate is the axial distance from the dump plane (see Fig. 1). The first data point of TC1 and TC2 at −125 mm is from a plane just after the guide vanes. Notice that the swirl number is the same for all cases at this position, as it should be since the velocity data are almost identical for all four cases (see Fig. 8). After this point, in the MT, the two approximate solutions agree with each other but gradually deviate from the FS and TC2 upstream cases. At the end of the MT, at x = 0 mm, the swirl numbers for the TC1 and TC2 downstream cases differ by 1.8% and 2.5% compared to the FS case. Moreover, the TC2 upstream case agrees well with the FS case in a large part of the MT but near the dump plane at x = 0 mm, the TC2 upstream case starts to deviate from the FS case, because of the influence from its downstream boundary condition. In the downstream region, the TC1 and TC2 cases continue to mimic each other and initially agree well with the FS case. However, at x = 30 mm, the TC1 and TC2 cases start to deviate from the FS case.

A possible explanation for the differences between the TC1 and TC2 cases from the FS case in the combustion chamber is that the steady averaged inlet boundary values for these two cases could have a significant influence on the vortex breakdown process. The FS case at the position of the inlet in the TC1 and TC2 cases experiences an unsteady flow created by the geometry of the swirler that will influence the rest of the flow, including vortex breakdown. Notice that the two approximate downstream solutions (TC1 and TC2) are almost identical having a difference of 0.65% between them at x = 0 mm. This indicates that the small differences between the boundary conditions extracted from the FS and TC2 upstream cases have a negligible influence on the results, which supports the second assumption behind the TSOW method (see Methodology).

In the following discussion of the accuracy of the approximate solutions, the focus will be on the agreement with the FS benchmark case in the combustion chamber (see Fig. 10). Two main phenomena that need to be replicated by an approximate solution for it to be useful are the sizes and positions of the Outer Recirculation Zone (ORZ) and the Inner Recirculation Zone (IRZ). The reason is that in a reactive case the flame stabilization will be highly dependent on the recirculation of hot combustion products. Another distinct feature of the solution that also has an effect on the flame behavior in a reactive flow simulation is the spreading angle (sometimes called opening angle) of the flow when it enters the combustion chamber.

Regarding the positions of the IRZ and ORZ, indicated by negative axial velocities in Fig. 10, the two approximate solutions agree well with the FS case. The opening angle of the swirling flow was 70.74 deg for the FS case, which differed by 1.21% (0.86 deg) from the TC1 case and by 2.4% (1.68 deg) from the TC2 case. The main discrepancy between the approximate solutions TC1 and TC2 and the FS case is the height where the IRZ begins. This was quantified by using a line probe along the vertical axis to determine the height at which the first negative axial velocity value was encountered. For the FS case, the first negative axial velocity value was found at 6.7 mm above the dump plane. For the TC1 case, this height was at 4.1 mm and for TC2 2.8 mm above the dump plane. The difference between TC1 and TC2 in this case probably stems from the small differences between the inlet boundary conditions. Again, a possible reason for the discrepancy between both cases and the FS case could be that the turbulence corresponding to constant inlet values will likely be different compared to the turbulence created by the swirler, especially when phenomena like vortex shedding are present. However, both approximate cases predicted almost the same lower height for the beginning of the IRZ with approximately 1 mm difference between their values.

The pressure field needs a special discussion since it plays a special role in parabolized flow situations. As discussed by Hirsch [2], information only travels in the principal flow direction for parabolized flow, but the pressure field has a more elliptic nature and depends to some extent on the downstream boundary conditions. The pressure boundary condition in an incompressible flow is in most commercial CFD codes defined indirectly by mass conservation that is enforced through a pressure correction equation [17]. This makes it difficult to guarantee that the boundary pressure on the outlet in the TSOW method will be identical to the full solution pressure at the same position. On the other hand, the earlier comparison indicates that the solution is relatively insensitive to small errors in the pressure field at the boundaries since the results in the different cases agree so well. In Fig. 11, the pressure fields of the FS and TC1 cases are compared. The average pressure at the outlet has been subtracted from the pressure field in each case to make the comparison easier. This is an operation without impact on the solution since the equations in an incompressible flow case only depend on the pressure gradient, hence subtracting an arbitrary constant from the pressure field does not affect the results. There are some small differences between the cases in the MT and in the combustion chamber but overall, the two pressure fields are almost identical. This qualitative check indicates that an accurate result can be reached with the TSOW method using the type of indirect boundary conditions for the pressure that is implemented in most commercial codes.

In order to investigate the computational time savings with the TSOW method, steady-state RANS solutions for TC1 and TC2 are compared with the FS case in Fig. 12. Similar to the URANS results, the standard TSOW approach (TC2 RANS) produced almost identical results with TC1 RANS that presumably used “perfect” boundary conditions from the FS case. Hence, the second assumption behind the TSOW method that small errors in the boundary conditions have a negligible impact on the results is again supported.

In the comparison with the FS case, it was found that the lower edge of the IRZ in the two RANS solutions protrudes into the MT to about 8.9 mm below the dump plane, a 15.6 mm difference with the FS and about 13 mm compared to the TC1 and TC2 solutions with a URANS method (see Fig. 10). In addition, the position of the ORZ moved to a higher axial position and its intensity diminished in the RANS cases. In a swirl number comparison with FS (at the dump plane), TC1 RANS and TC2 RANS deviated by 16.6% and 17.6% respectively. A possible reason for the differences in the FS case could be the different modeling approaches but differences in the discretization errors could also be an explanation. Overall, both steady-state RANS cases agree quite well with the FS case and can be used for engineering design when computational speed is the highest priority.

In the present simulations, the FS case was run on a supercomputer cluster (the Kebnekaise cluster at the HPC2N supercomputer center at Umeå University) and had a total computational time of 85 h to reach 2.5 s of physical time using 168 cores. On the same hardware, the TSOW method (TC2 and TC2 upstream) was solved in 62 h. Solving only the TC2 downstream flow took 55 h. Hence, the complete TSOW method was 23 h faster than the FS. In addition, in an engineering design situation, only the downstream part would be needed for additional iterations, which would give a further 30 h speedup for every iteration. For design optimization with a response surface method and a full quadratic approximation, at least (N + 1)(N + 2)/2 CFD solutions will be required [1]. With three design variables (N = 3), this corresponds to 10 CFD evaluations, which means that the TSOW method has the potential to give a time-saving of about 270 h or about 11 calendar days for the present case.

In the present case, it turned out to be possible to use a RANS method for the downstream flow in the TSOW method, which gave a significant time-saving compared to the FS case. The TC2 RANS case was solved in 3 h with the same supercomputer hardware. This brings the total time of the TSOW method down to 10 h, or 12% of the time needed for the FS results. For further iterations where only the RANS downstream solution would be needed, the speedup would be 96% for each iteration. Moreover, the much smaller workload in the TC2 RANS made it possible to solve it on a typical CAE workstation (HP Z6 G4 Workstation) in 19 h, removing the need for supercomputer resources. In spite of the huge computational time savings of the TSOW method, confirming the desired design iteration with a full simulation is good practice.

The approximate TSOW method for CFD solutions of complex geometry problems was assessed against a benchmark solution of the cold flow in the full flow domain of the CeCOST experimental burner setup. The method solves the flow in the complete domain by dividing the domain into two parts and using results from the upstream solution as a boundary condition for the solution of the flow downstream.

The method is based on two implicit assumptions that were validated using the benchmark solution. The solutions were obtained with grids having around 2 million control volumes for the approximate solutions and 4 million control volumes for the full solution. With Richardson extrapolation, it was estimated that the grid resolution was in the asymptotic range and that the error level for the study was sufficiently low to make comparisons between different solution approaches possible.

The first assumption behind the method is that the solution in a partial upstream domain will be the same as the benchmark solution. The agreement between the two solutions was found to be excellent, except very close to the outflow boundary in the partial domain solution. The overlap between the upstream and the downstream domain means that the boundary condition for this flow will be extracted in a part of the domain where the agreement is excellent. This overlap is crucial for the successful application of the TSOW method as it allows the extraction of data from a surface that is not influenced by the peculiarities of boundary conditions (e.g., zero normal derivatives at the outlet). Hence, the first assumption is supported by the numerical results.

The second assumption is that the flow in the downstream solution will be the same in the common domain with the benchmark solution if the boundary conditions are equal to the results of the benchmark solution at the position of the inflow boundary. This assumption was tested first by using the results from the benchmark case as boundary conditions for the downstream solution and second by using the results from the solution with the upstream partial domain as boundary conditions for the downstream solution. In the comparison, both downstream solutions agreed excellently with each other, indicating that small errors in the boundary conditions are damped. In the comparison with the benchmark solution, the overall agreement for the swirl number through the whole domain was relatively good although the difference grew with increasing flow distance toward the outlet. It is hypothesized that this difference is an effect of the steady-state inlet conditions that were used for the two downstream solutions and not an effect of the TSOW method. Using steady-state boundary conditions is equivalent to a filtering of the benchmark solution that removes the effect from the periodic vortex shedding in the swirler.

In the development of the benchmark solution, it was necessary to use a URANS method to obtain converged results. It is believed that the reason was periodic flow separation in the swirler that is poorly represented by a steady RANS turbulence model. In the URANS simulation, these large-scale deterministic motions were resolved which improved convergence and made it possible to obtain stable averages. In the assessment of the TSOW method, URANS was used for all simulations, to avoid influences of different turbulent models in the different cases. However, after the assessment, additional runs were made with a steady RANS model for the two downstream flows and convergence was achieved. It is believed that the reason for the improved convergence was the use of steady-state boundary conditions instead of the time-periodic flow that was generated by the swirler in the benchmark case and the upstream part of the TSOW solution. The two approximate downstream RANS solutions deviated more from the benchmark solution than the URANS solutions for TC1 and TC2. However, since the boundary conditions were steady in both cases, this difference is most likely a direct consequence of the different modeling approaches (URANS vs. RANS).

The goal with the TSOW method is to save time and to make it possible to solve complex problems on less powerful computers. In the comparison of time consumption, it was found that the complete TSOW method with URANS took 22 h less than the benchmark solution. Furthermore, it was shown that the method can give massive time savings in practical engineering design optimization where only the downstream part will be modified. For this case and with a relatively modest number of design variables, savings of the order of whole calendar weeks can be obtained. Moreover, with the faster steady RANS method for the downstream flow, the computational time was less than 10% of the benchmark solution. This speedup made it possible to use a standard engineering workstation to solve a problem that required an advanced supercomputer to solve for the benchmark case.

The simulations were performed on resources (Kebnekaise) provided by the Swedish National Infrastructure for Computing (SNIC) at High Performance Computing Center North (HPC2N) partially funded by the Swedish Research Council through grant agreement no. 2016-07213. The authors thank all the staff of HPC2N for their technical assistance. The financial support from the Swedish Biomass Gasification Centre SFC, which is funded by the Swedish Energy Agency and a consortium of companies through project number 34721-3, is acknowledged.

There are no conflicts of interest. This article does not include research in which human participants were involved. The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

Informed consent is not applicable. This article does not include any research in which animal participants were involved.

t =

physical time, s

A =

area, m2

R =

maximum radius, mm

S =

swirl number

U =

mean axial velocity, m/s

W =

tangential velocity, m/s

T_F =

flow-through time, s

U_max =

maximum mean axial velocity, m/s

U_RMS =

RMS axial velocity, m/s

U_{RMS max} =

maximum RMS axial velocity, m/s

FS =

full solution

TC1 =

two-step case 1

TC2 =

two-step case 2

TSOW =

two-step one-way coupled

MT =

mixing tube