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 phasefield 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 nonzero 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
orparallel
 Conditions where pressure governs the flow across the boundary, i.e.
orthogonal
,free
orgradient
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 noslip/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 2D 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 (noslip): poiseuillelike 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, noslip 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 theW
andE
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^{1}: 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  Output 

0 (off)  As no verbosity option 
1 (min)  Each flow solver call 
2  Each continuity iteration 
3  Additional output 
4  Each linear solver call 
5  Each linear solver iteration 
6 (max)  Additional 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 (CourantFriedrichsLewy) 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 CFLnumber 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 300500 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 300500. 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. 2030) 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 FlowSolver: # combinedseparatedpiso [verbosity <lvl>] # [ steady_[re]startana_[re]start [pre_iter <it>] ] combined verbosity 1 # Time stepping in FlowSolver: fix <tstep[s]> # or: cfl <CFLnumber> [first <tfirst>] [max <tmax[s]>] [min <tmin[s]>] fix 0.000005 # Convergence criterion for impulse eqn. (absolute) [mue_m/s] # (real) [warn(*)nowarn] 1e5 # Max. number of iterations for linear solver: # (int) [bicgmr2(*)bicgsafebicgstab] 500 # Convergence criterion for pressure eqn. (absolute) [Pa] # (real) [warn(*)nowarn] [adapt(*)noadapt] 1e8 # Divergence criterion for continuity equation [1/s] # (real) [warn(*)nowarn] 1E2 # 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 PVResidual 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 verbosity5
... PVsolver: pressure impulse iter: inner total Criteria: 1.00000E08 1.00000E05 500 PVLinsolv: 0.00000E+00 1.27322E+01 1 1 PVLinsolv: 5.39842E01 2.90440E+02 2 2 PVLinsolv: 2.06884E02 1.01393E+01 3 3 PVLinsolv: 2.10092E02 5.94715E+00 4 4 ... PVLinsolv: 6.15636E07 1.28979E04 19 497 PVLinsolv: 5.01661E07 8.82903E05 20 498 PVLinsolv: 5.01661E07 8.82903E05 1 499 PVLinsolv: 8.59733E07 2.31118E04 2 500 PVResidual: 7.07738E07 1.55810E04 28 500 Cont. Err: 7.00332E+01 7.07738E07 1.55810E04 1 500 ...
Example 3: tracking convergence for continuity and linear solver results for pressure and impulse for outer iterations¶
Linear solver (
combined
) output at verbosity2
... Flow: t_pf= 1.000000E04, t_flow= 5.000000E06, dt= 5.000000E06 Type: continuity pressure impulse iter. cont lin. solver Criteria: 1.00000E02 1.00000E08 1.00000E05 Cont. Err: 7.00332E+01 7.07738E07 1.55810E04 1 500 Cont. Err: 1.91968E+01 9.74374E06 2.78629E03 2 1000 Cont. Err: 9.75127E+00 4.81484E08 1.76150E04 3 1500 ... Cont. Err: 1.54485E02 7.54372E09 2.26806E06 26 9657 Cont. Err: 1.21605E02 8.22486E09 3.38942E06 27 9886 Final Res: 9.56403E03 8.94495E09 3.65755E06 28 10067 V_x = 4.3005E+01 V_z = 4.5388E13 mic/s, CFLnr = 0.8482E03 ...
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 restartfile. 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.