Flow Solver¶

The flow solver simulates laminar fluid dynamics ($Re \lesssim 10^4$) for an incompressible liquid. It is applied to phase $0$. All other phases are considered as immobile obstacles that may be affected by phase field calculations but are not transported with the liquid. Advection is included in species transport for concentration coupled simulations. If concentration dependent density changes are defined, buoyancy affects liquid flow by application of the Boussinesq approximation.

Activating the flow solver¶

The solver is activated with the flow or flow_coarse flag in the model section. If flow_coarse is chosen, an additional coarsening parameter $c$ will be asked for (which should be a low integer, preferably $2$). Fluid flow will be simulated on a coarsened grid where each cell consists of $c^d$ cells of the phase-field grid. The user needs to enter additional material data, boundary conditions for velocity and pressure and numerical parameters in the numerical section as will be discussed below.

Pressure and the velocity components averaged over the domain are given out at output time steps. Summary warnings are provided if the convergence criteria are not met.

Note that Restarts of fluid flow calculations will only work for the same cell dimension and coarsening factor of the simulation domain.

Composition driven convection¶

If any of the $d\rho/dc$ terms entered in the Material Data for Fluid Flow section is non-zero this feature will be activated. The density depends linearly on the composition and affects fluid flow according to the Boussinesq approximation.

To avoid numerical overhead, metallostatic pressures are neglected in the simulation. This means that buoyancy is evaluated relative to a reference composition. Usually this reference is the averaged composition of the liquid for the current time step unless reference concentrations are given as optional parameter for each component in the material data section. While the averaging scheme works well if concentration changes affect the simulation domain more or less uniformly, a fixed reference works better for localized concentration changes, e.g. a dendrite growing from the bottom of the simulation domain. In this case the initial concentrations for the liquid phase are probably a good choice for reference concentrations.

Boundary conditions¶

The boundary conditions for fluid flow consist of a combination of velocity and pressure boundary conditions. Basically three cases can be distinguished:

• Periodic, where the flow patterns of opposite boundaries are linked
• Conditions that fix the fluid flow across the boundary, i.e. fixed, symmetry, iso or parallel
• Conditions where pressure governs the flow across the boundary, i.e. orthogonal, free or gradient

The possible options for boundary conditions are listed in the Flow Boundary Conditions section of the input chapter.

Figure 1¶

Pressure offset for periodic pressure conditions

In the case of periodic velocity, the pressure condition should also be periodic. The user will then be asked for a pressure offset. As shown in Figure 1 this adds an overall pressure gradient to the local pressure field.

In cases of fixed fluid flow across the boundary the pressure boundary condition does not affect fluid flow and can be chosen freely or set to - (no condition). If all boundary flows are fixed outflow must match inflow.

If pressure governs the flow across boundaries (e.g. free in/outflow) at least one of the corresponding pressure conditions should be fixed (e.g. on an outflow). For other boundaries (e.g. inlets) recommended pressure conditions are fixed or von neumann.

If the fluid flow is driven by a pressure gradient, friction must be present (either from solid phase in the simulation domain or from no-slip/fixed velocity walls parallel to the gradient) to avoid constant acceleration. Frictionless (or gradient) walls along the pressure gradient will also lead to situations with constant acceleration.

In the 2-D regime the following velocity/pressure conditions will cover many cases:

• ppii/pp-- with some pressure offset and some solid phase in the simulation domain: simulation of permeability through a periodic array of grains/dendrites.

• ppfi/pp-- with zero velocity on the lower boundary (no-slip): poiseuille-like flow symmetric to the top boundary (parabolic velocity profile), e.g. for a dendrite growing perpendicular into a flow

• ppff/pp-- with zero periodic pressure offset, no-slip at the bottom and a fixed velocity parallel to the wall at the top: couette like flow (linear velocity profile)

• ffii/---- with the same, perpendicular velocity on the W and E walls: fixed volume flow

Iso boundary conditions can be replaced with symmetry conditions if that matches the geometry better. For pressure gradient driven flows without periodicity the pp../pp.. combination can be replaced with ee../ff.. or ee../nf.., and for the fixed volume flow ff../--.. can be replaced with fe../-f...

Numerical parameters¶

First the user has to choose a solver scheme. Of the three choices combined, separated and piso, in most cases combined is the best choice for fluid flows through microstructures.

In cases of large velocities and large length scales the piso solver should be used. A good criterion to decide which solver to use is the Reynolds number $Re$:

where $L$ is the characteristic length (e.g. grainsize, DAS or domainsize), $v_\text{avg}$ the average velocity and $\nu$ the kinematic viscosity.

For example $L = 100 \mu m$, $\nu = 10^{-2} cm^2/s$ and $\nu_\text{avg} = 10^4 \mu m/s$ ($= 1cm/s$) would result in $Re = 1$. For low Reynolds numbers ($Re \ll 1$), viscous forces dominate the fluid flow and the combined solver should be used, for larger Reynolds numbers ($Re > 1$) inertial forces become more important and the piso solver should be considered.

After the solver choice the steady_start option can be entered to iterate to a steady fluid dynamics solution before any phasefield calculations. It should be combined with the pre_iter option to set the number of flow solver calls and ensure that a steady state has been reached. The _restart variant enables iterating to a steady flow after a restart. The ana_start option will use an analytical solution for poiseuille or couette flow determined by the boundary conditions (ignoring any solid regions) before applying the iterative solver.

As another optional parameter the flow solver verbosity level can be set. This enhances the output of the solver to help track convergence and aid in setting appropriate numerical parameters.

The verbosity level determines at which level of the solver output is produced. Figure 2 shows the basic scheme for the flow solver¹: for each time step the solution is iterated until the continuity criterion is fulfilled, and nested within this loop the impulse- and pressure equations are solved by an iterative linear solver until the convergence criteria for velocity and pressure are met.

Figure 2¶

Flow solver iteration loops

The different verbosity levels are cumulative and produce output:

Verbosity is switched off after $100,000$ lines of output to prevent clogging the file system. This limit can be changed with the verbosity_max_lines option. It is recommended to switch verbosity off once the numerical parameters are adjusted.

After the solver scheme the time step size has to be set. It is recommended to use a fix time step size. Alternatively the cfl option can be chosen. If cfl is chosen it should be used in combination with the max option to set a maximum step size and/or the first option to set the starting time step.

There are two CFL (Courant-Friedrichs-Lewy) numbers that may help to arrive at a good starting point for the time step size, one arising from the advection term in the Navier Stokes equations:

where $C_\text{adv}$ is the CFL number for advection, $\nu_x$,$\nu_y$,$\nu_z$ are components of the local velocity, $\Delta x$ is the grid distance of the flow solver and $\Delta t$ the time step. $C_\text{adv}$ is relevant for high Reynolds numbers and timesteps can be limited accordingly with the cfl option. In that case the piso solver is recommended and a limit $C_\text{adv} < 0.3$ is a conservative choice. Depending on the specific problem higher $C_\text{adv}$ may work.

For low Reynolds numbers the viscosity term becomes more important and the relevant CFL-number is

with the kinematic viscosity $\nu$. For a given viscosity and grid spacing the time step is then limited to:

Here $C_\text{visc} = 0.2$ is a conservative choice. In many cases, particularly when flow patterns change slowly with time and simulation cells, much larger time steps may converge and lead to better performance. Even time step sizes corresponding to $C_\text{visc} = 250$ or $C_\text{adv} = 30$ have been found to converge in specific cases, but this cannot be generalized and it is up to the user to find the optimal time stepping for a particular simulation.

After the time step size convergence criteria and maximal loop counts for the linear solvers need to be set. The convergence criteria limit the absolute rms (root mean square) deviations when solving the impulse (= velocity), pressure or continuity equations. The first criterion limits the velocity error in the impulse equation. When this criterion is failed at the end of a flow solver time step a warning is produced (after the very first warning only summarized warnings are produced at output time steps). These warnings can be switched off with the nowarn option.

Next the loop count for the linear solver of the impulse equation is queried (in the combined case this is the solver for the combined equation system), the default linear solver is bicgmr2, optionally one may chose bicgsafe or bicgstab. Often iteration in the range $300-500$ work very well here, but in some cases for large grids higher values (up to $3000$) may work better at the start of the iteration (or to converge steady flow patterns). In later stages of a simulation this value should rarely be reached, if it is consider shorter time steps and/or less restrictive convergence criteria.

After that a convergence criterion for the pressure equation is entered. In an adaptive scheme this value is reduced in successive continuity iterations unless the noadapt option is used. For the separated and piso solvers a loop count for the linear solver of the pressure equation needs to be entered, typical limits should be in the range $300-500$. The piso solver also demands an "inner" iteration count for the number of corrector cycles.

Next one needs to enter a criterion for the continuity equation, often $10^{-2}$ is a good choice here. The unit is $1/s$ as this limits a volume flow density. To visualize this, a value of $10^{-2}$ in one simulation cell over one second would mean that the cells outflow exceeds the inflow by $1%$ of the cell volume.

The next parameter limits the number of continuity iterations (indicated by the green arrow in Figure 2. In some cases higher values (e.g. $20-30$) here may help convergence at the beginning of a simulation or a steady flow calculation, but if many iterations are needed in later stages of the simulation, shorter flow solver time steps should be considered.

For the combined and separated solvers an additional underrelaxation parameter α is needed. For the separated solver it is possible to set this value to $1$. For the combined solver α must be less than $1$ and should typically be in the range $0.7 \lt \alpha \lt 0.97$.

Example 1¶

Numerical parameters for fluid flow

...
# Other numerical parameters
# ==========================
# (*) marks default settings for options
# Selection of the Flow-Solver:
# combined|separated|piso  [verbosity <lvl>]
#     [ steady_[re]start|ana_[re]start [pre_iter <it>] ]
combined verbosity 1
# Time stepping in Flow-Solver:  fix <tstep[s]>
# or:  cfl <CFL-number> [first <tfirst>] [max <tmax[s]>] [min <tmin[s]>]
fix 0.000005
# Convergence criterion for impulse eqn. (absolute) [mue_m/s]
#    (real) [warn(*)|nowarn]
1e-5
# Max. number of iterations for linear solver:
#   (int) [bicgmr2(*)|bicgsafe|bicgstab]
500
# Convergence criterion for pressure eqn. (absolute) [Pa]
#    (real) [warn(*)|nowarn] [adapt(*)|noadapt]
1e-8
# Divergence criterion for continuity equation [1/s]
#    (real) [warn(*)|nowarn]
1E-2
# Max. iterations of the flow solver
50
# Factor for underrelaxation of pressure 0.0 < alpha < 1.0
0.97
# Phase minimum?
...

Adjusting numerical parameters¶

For a given fluid flow simulation the errors in velocity, pressure and continuity are correlated, so a matching set of convergence criteria needs to be found. Here it helps to track the convergence by switching on flow solver verbosity. Example 2 shows the convergence of the linear combined solver at the beginning of the calculation. The criteria match the input of Example 1. The top shows the convergence of the linear solver at verbosity 5, in the combined case pressure and impulse equation are solved together, so pressure and velocity converge side by side. At the head of the pressure and impulse column the criteria are shown that the linear solver should reach. These match the parameters from the input file. The 4^th and 5^th column show the number of "inner" iterations since the last restart and the total iterations in this solver call. Here the linear solver is stopped after $500$ cycles and the solution is further converged in subsequent calls to the linear solver. The convergence reached at this stage is shown in the PV-Residual line. In the piso and separated cases the linear solvers for pressure and impulse are separate so either the impulse or the pressure column is printed. Interspersed between output from the linear solver are output lines resulting from lower verbosity levels.

The bottom shows how to track convergence of the continuity equation with the outer convergence cycles at verbosity $2$. With this verbosity the 4^th column shows the loop count for continuity and the 5^th column the accumulated linear solver iterations. As one can see when comparing the Final Resolutions line to the Criteria line the pressure criterion is sufficient to reach the desired accuracy for the continuity and the impulse criterion is chosen to match at this point. Here after $28$ cycles the continuity criterion (first column) is reached. The high number of iterations in this example is due to the high accelerations at the start of the simulation.

At the beginning of each fluid flow time step the current times for phasefield (t_pf) and fluid flow simulation (t_flow) as well as the timestepsize (dt) are printed, at the end the velocities averaged over the liquid volume and $C_\text{adv}$ are printed.

Example 2: tracking convergence for pressure and impulse per linear solver iteration¶

Linear solver (combined) output at verbosity 5

...
   PV-solver:                      pressure        impulse iter: inner       total
    Criteria:                   1.00000E-08    1.00000E-05                     500
  PV-Linsolv:                   0.00000E+00    1.27322E+01           1           1
  PV-Linsolv:                   5.39842E-01    2.90440E+02           2           2
  PV-Linsolv:                   2.06884E-02    1.01393E+01           3           3
  PV-Linsolv:                   2.10092E-02    5.94715E+00           4           4
...
  PV-Linsolv:                   6.15636E-07    1.28979E-04          19         497
  PV-Linsolv:                   5.01661E-07    8.82903E-05          20         498
  PV-Linsolv:                   5.01661E-07    8.82903E-05           1         499
  PV-Linsolv:                   8.59733E-07    2.31118E-04           2         500
 PV-Residual:                   7.07738E-07    1.55810E-04          28         500
   Cont. Err:    7.00332E+01    7.07738E-07    1.55810E-04           1         500
...

Example 3: tracking convergence for continuity and linear solver results for pressure and impulse for outer iterations¶

Linear solver (combined) output at verbosity 2

...
Flow: t_pf=  1.000000E-04, t_flow=  5.000000E-06, dt=   5.000000E-06
        Type:     continuity       pressure        impulse  iter. cont lin. solver
    Criteria:    1.00000E-02    1.00000E-08    1.00000E-05
   Cont. Err:    7.00332E+01    7.07738E-07    1.55810E-04           1         500
   Cont. Err:    1.91968E+01    9.74374E-06    2.78629E-03           2        1000
   Cont. Err:    9.75127E+00    4.81484E-08    1.76150E-04           3        1500
...
   Cont. Err:    1.54485E-02    7.54372E-09    2.26806E-06          26        9657
   Cont. Err:    1.21605E-02    8.22486E-09    3.38942E-06          27        9886
   Final Res:    9.56403E-03    8.94495E-09    3.65755E-06          28       10067
V_x = 4.3005E+01 V_z = 4.5388E-13 mic/s, CFL-nr = 0.8482E-03
...

Optimizing flow solver performance¶

Fluid flow simulations are time intensive, so some care should be taken to optimize them for performance. There are four approaches to speed up your calculation:

• Employing parallelisation, see input section Parallel Execution
• Solve fluid flow on a coarser grid (see input section Model, current topic Activating the flow solver)
• Optimizing the time step size
• Optimizing solver type and parameters

When optimizing time step size and solver parameters it should be noted, that in most simulations it takes some time to establish a steady or near steady state, e.g. a pressure driven flow that accelerates until an equilibrium between accelerating pressure and frictional drag is reached. If you start with a fully liquid domain the ana_start option (see current topic Numerical Parameters) is a shortcut to establish a steady flow for some boundary conditions. If multiple simulations are done for parameter variations consider starting from a common restart-file. Optimal time step sizes are not necessarily the longest possible time steps, better minimize (linear) solver cycles per simulation time (on average).

When tuning the convergence parameters (see current topic Numerical Parameters) less restrictive convergence conditions and/or fewer cycles benefit performance only up to a point since reduced accuracy in one time step adversely affects performance of subsequent time steps. Also if accuracy is repeatedly reduced in successive time steps this will quickly lead to divergence.

Usually each fluid flow time step is converged until the convergence conditions for continuity and impulse equation are met, but it is possible to use solver loop counts as limiting element. A good example of this is using the piso solver with two inner (correction) cycles and only one solver iteration. The convergence of continuity is then governed by the condition for the pressure equation. The error for continuity and impulse equation can then be tracked by the warning messages. When using the solver in this way choose short time steps to ensure convergence.

Besides numerical parameters also the choice of boundary conditions (see input section Flow Boundary Conditions) has a major impact on flow solver performance. Often similar conditions can be modelled with different sets of boundary conditions. Preferable choices of boundary conditions are discussed in in the current topic Numerical Parameters).

If the Boussinesq model is activated and concentration changes are localized consider setting a fixed reference concentration to avoid large scale changes in the pressure field (see the current topic Composition driven convection.