# B2000++ Solvers

## Linear Static Solver

Linear static (stationary) or linearized problems are solved with the
linear static solver. It is activated for a specific case
with the case parameter `analysis=linear`

. If several
load cases with the same constaints are present, the subcase option
(see case) can be applied during linear analysis, reducing
computation time.

## Static Nonlinear Solver (B2000++ Pro)

Nonlinear static (stationary) problems are solved with the nonlinear
static solver. It is activated for a specific case
with the case parameter `analysis=nonlinear`

.

A nonlinear problem is solved with an incremental procedure following an increment control parameter denoted \(\lambda\).

Five increment constraint strategies can by used with this solver. The figure below illustrates these by a geometric representation for a single DOF problem (DOF parameter plotted horizontally and control parameter \(\lambda\) plotted vertically):

The correction strategy is of the Newton-iteration type, optionally augmented with an acceleration method and a line search algorithm.

Optional Solver MDL Parameters

In addition to the MDL command `analysis nonlinear`

, the following
commands may be specified in the case block.

`correction_type newton|accelerated_newton`

Specifies the type of correction in a given Newton iteration. Default is

`newton`

.The

`newton`

correction type only makes use of the residuum and the tangent matrix (updated or not) to compute the correction.The

`accelerated_newton`

correction type is meaningful with a modified Newton method. It monitors whether the sequence of Newton iterations converges (whether the residuum decreases), and if not (divergence), invokes the line search algorithm to find a better solution. With line search, it is often possible to find a solution outside the convergence region. Thus, when using the MDL commands`newton_method modified`

or`newton_method delayed_modified`

, the accelerated Newton correction can converge for larger time steps than the normal Newton correction type.See also the MDL command

`line_search`

.

`correction_termination_test flux_normalised|regularised|dof_and_residue|normalised_dof_and_residue`

Specifies the type of test which indicates whether the Newton iterations of the correction phase have sufficiently converged. The tests make use of some or all of the MDL commands

`tol_residuum`

,`tol_solution`

, and`tol_work`

. Default is`flux_normalised`

.The

`flux_normalised`

correction termination test computes the residuum normalised by the energy flux for each type of field (displacements, rotations, Lagrange multipliers, etc.) separately. The`flux_normalised`

test has additional options, see Flux-normalised Correction Termination Test.The

`regularised`

correction termination test regularizes the residuum, the solution changes, and the energy change, with the goal of making the test independent of the units. Thus, changing the units of an FE model (e.g. [kg,m,s]) to another (consistent) unit system (e.g. [t,mm,s]) should not require modification of the tolerances. The following tolerances are used:`tol_residuum`

(default 1e-3),`tol_solution`

(default 1e-3), and`tol_work`

(default 1e-7).The

`dof_and_residue`

correction termination test tests the L2-norm residuum against`tol_residuum`

and the L2-norm of the solution changes w.r.t. the previous Newton iteration against`tol_solution`

. The following tolerances are used:`tol_residuum`

(default 1e-3) and`tol_solution`

(default 1e-3).The

`normalised_dof_and_residue`

correction termination test divides the L2-norm residuum by that of the first Newton iteration. The same is done for the L2-norm of the solution changes w.r.t. the previous Newton iteration. The following tolerances are used:`tol_residuum`

(default 1e-3) and`tol_solution`

(default 1e-3).

`gradients VALUE`

Controls the computation of the gradients (strains, stresses, heat transfer, etc.) for the current case. If

`VALUE`

(int) is set to 0 (default) no gradients are computed and saved, i.e. only dof solutions will be found on the database. If`VALUE` is set to -1, the gradients for the last converged step of at stage or the linear solution step are computed and saved. If ``VALUE`

is set to a positive value, the gradients will be computed at each istep step of at stage as well as for the last step of a stage.

`gradients_only YES|NO`

Specifies whether only gradients (and reaction forces) should be computed, while the solution vector is defined solely by the the initial conditions, and the boundary conditions and the current load increment. Default is

`NO`

.

`increment_control_type load|state|hyperplane|hyperelliptic|local_hyperelliptic`

Specifies the increment constraint strategy to use. The increment constraint strategy serves two purposes: It defines the increment constraint for the predictor phase and the corrector phase, and it computes the predictor solution from the solution of the last converged step. This predictor solution is then used as initial solution to the correction phase.

The following strategies are available:

`load`

means load-controlled. To avoid singularities at limit points in post-buckling analysis,`residue_function_type artificial_damping`

can be used.A state-controlled strategy is used with

`state`

.The other strategies are of the continuation type. For the

`hyperelliptic`

and`local_hyperelliptic`

strategies, the parameters`hyperelliptic_a`

and`hyperelliptic_b`

define the size of the ellipse.The default strategy is

`load`

.For more details on the kinds of constraints, refer to Implemented increment constraint strategies.

`init_step_number VALUE`

The initial step cycle number (int). Default is 1.

`line_search yes|no`

Specifies whether line search should be used for the

`accelerated_newton`

method (see the MDL command`correction_type`

). The line search method satisfies the strong Wolfe condition. When the cost of calculating the first variation is significantly lower than the cost of calculating the second variation, it is recommended to have line search enabled, as it improves the probability that the Newton iterations converge. Default is`yes`

.

`max_divergences VALUE`

Defines the maximum number of consecutive divergences (int) during the correction phase of a load step. Default is 4. When

`VALUE`

is reached, the correction phase is terminated unsuccessfully, and the step increment is re-started with a reduced step size.

`max_newton_iterations VALUE`

Specifies the maximum number of Newton iterations (int) per increment. Default is 50.

`max_steps VALUE`

Maximum number of steps (increments) (int) per stage. Default is 99999.

`newton_method conventional|modified|delayed_modified`

Specifies what kind of Newton method should be used. Default is

`conventional`

.The

`conventional`

method performs factorization for the predictor and for each Newton iteration, and thus is a full Newton method.

`For the modified`

method, factorization of the tangent matrix is performed when the correction phase has terminated unsuccessfully. Hence, several increments may be performed without any factorization.The

`delayed_modified`

method performs a factorization for the predictor and for the first Newton iteration.Note that the

`artificial_damping`

option for the MDL command`residue_function_type`

enforces factorization each time the step size changes. For the`modified`

and`delayed_modified`

Newton methods, is is recommended to use the`correction_type accelerated_newton`

option.

`newton_periodic_update VALUE`

For the modified and delayed modified, non-accelerated Newton methods, `` VALUE`` (int) specifies the maximum number of consecutive iterations without re-factorization. When `` VALUE`` is reached, a re-factorization will be performed, and the counter is reset. Note that for high values of `` VALUE`` and depending on the

`newton_method`

setting, several increments may be performed until this number of Newton iterations is reached. Default is 100.

`residue_function_type default|artificial_damping|scaled`

Specifies the type of residuum (first variation or internal forces vector) and its derivative (second variation or tangent matrix).

The

`default`

type computes the first and second variation without any modifications.In a load-controlled post-buckling analysis involving collapse, the

`artificial_damping`

type is necessary for the Newton iterations to converge. It adds viscous forces to the residuum vector (first variation). The velocities are computed from the difference of the solution to the last converged solution, divided by the step size. To the tangent matrix (second variation), the viscosity matrix divided by the step size is added. With artificial damping, singularity of the tangent stiffness matrix at limit points can be avoided and energy can be dissipated during mode jumps. See also artificial damping.The

`scaled`

type scales the DOFs by their respective diagonal element in the tangent matrix. This may be beneficial in continuation analysis (e.g.`increment_control_type local_hyperelliptic`

). Due to the scaling, it is recommended to use the default correction termination test.

`step_size_init VALUE`

Step size (float) at the beginning of the stage. Default is 0.1.

`step_size_max VALUE`

Maximum step size (float) during the stage. Default is 1.0. In continuation analysis, this is the size along the load path.

`step_size_min VALUE`

Minimum step size (float) during the stage. Default is 1e-12. In continuation analysis, this is the size along the load path.

`step_size_lambda_max VALUE`

Maximum lambda step size (float) during the stage in a continuation method analysis. Default is 1.0.

`step_size_lambda_min VALUE`

Minimum lambda step size (float) during the stage in a continuation method analysis. Default is 1e-12.

`step_size_dof_max VALUE`

Maximum degrees of freedoms step size (float) during the stage in a continuation method analysis. Default is 1.0.

`step_size_dof_min VALUE`

Minimum degrees of freedoms step size (float) during the stage in a continuation method analysis. Default is 1e-12.

`tol_residuum VALUE`

The error tolerance value (float) for the residuum (‘equilibrium’) error; used by the Newton correction termination test. The default depends on the

`correction_termination_test`

option.

`tol_solution VALUE`

The error tolerance value for the dof (‘solution’) error; used by the Newton correction termination test. The default value on the

`correction_termination_test`

option.

`tol_work VALUE`

The error tolerance value (float) for the work error; used by some Newton correction termination test. The default depends on the

`correction_termination_test`

option.

The static solver is most commonly used for stress analysis involving geometric, material, or boundary nonlinearity, and for heat transfer analysis involving nonlinear material or boundary behavior.

Note

For load-controlled quasi-static stress analysis with artificial damping, an alternative exists in the dynamic nonlinear solver, which can also be used for quasi-static stress analysis when specifying the MDL case parameters

```
analysis dynamic_nonlinear
residue_function_type artificial_damping
```

Due to its predictor and its mechanism to control the time increment by means of an error estimator, the dynamic nonlinear solver may be - depending on the problem and depending on the solution parameters - more effective than the static nonlinear solver.

The static solver can also steer a weakly-coupled multi-disciplinary analysis, this is briefly explained for the case of steady-state aero-elasticity: The spatial coupling, CFD mesh deformation, and CFD solver are wrapped inside a specially-designed and implemented nonlinear natural boundary condition object. For each Newton iteration, the static solver evaluates the natural boundary condition for the solution field; in this case, the boundary condition returns the spatially mapped aerodynamic forces that were calculated from the deformed configuration.

Example: Load-Controlled Stress Analysis with Dissipation

In stress analysis, a problem may be considered highly nonlinear if geometric instabilities, material nonlinearity such as degradation, or contact conditions are present. A typical parameter specification for such situations is:

```
analysis nonlinear
residue_function_type artificial_damping
dissipated_energy_fraction 1e-6
step_size_init 1e-2
step_size_max 1e-1
```

Load control is enabled by default. For highly nonlinear problems, the
full Newton method (default) is often more expedient than the
accelerated modified Newton method. The default maximum number of
Newton iterations `max_newton_iterations`

should be sufficient in
many cases. Lower values may cause many small time increments or even
failure to converge, while higher values may lead to increased
computation time. In case of very slow but reliable convergence,
higher values (100-200) may be necessary.

When using a load-controlled strategy, a mechanism to dissipate energy
must be enabled when instabilities occur. This is achieved with the
residue function type set to `artificial_damping`

. The amount of
damping is controlled with the `dissipated_energy_fraction`

parameter. High values accelerate the analysis, increasing the
increment sizes, but may influence the results on the other
hand. Therefore, a compromise between analysis time and accuracy has
to be found. While the parameter value is independent of the physical
units, the necessary amount of damping depends on the problem. Thus,
the default of 1e-4 may not be appropriate. It is recommended to vary
this parameter in powers of 10 (e.g. make an analysis for 1e-7, 1e-6,
1e-5, 1e-4, 1e-3) and observe its effect.

When the problem is expected to be nonlinear from the start, an
initial step size may be specified with the `step_size_init`

option. The maximum step size should always be set, using the
`step_size_max`

option, for the reason that different load paths may
be followed and collapse cannot be reliably predicted if the step size
is too large (the solver may “jump over” collapse paths). The default
of 1e-12 for the minimum step size is sufficient for most problems.

Another important option is `max_divergences`

which may be used to
trigger collapse by using smaller values (1 or 2) than the default. On
the other hand, analyses involving material degradation may benefit
from higher values (increasing the chances that the Newton iterations
will converge).

Example: Post-Buckling Analysis

A typical static potentially highly nonlinear post-buckling analysis parameter setting is

```
analysis nonlinear
increment_control_type hyperplane
step_size_init 1e-2
step_size_lambda_max 1e-1
```

The parameter `increment_control_type hyperplane`

enables a specific
method of continuation analysis, here the `hyperplane`

method.

The maximum step size in thew loading-direction can be controlled by
means of the `step_size_lambda_max`

option. Its counterpart is the
`step_size_dof_max`

option which controls the maximum step size in
the DOF direction. The `step_size_max`

option controls the
arc-length of the increments.

Example: Mildly Nonlinear Continuation Analysis

A typical static ‘mildly’ nonlinear post-buckling analysis parameter setting is

```
analysis nonlinear
increment_control_type hyperplane
correction_type accelerated_newton
newton_method delayed_modified
step_size_init 1
```

The parameter `increment_control_type hyperplane`

enables a specific
method of continuation analysis, here the `hyperplane`

method. If
the problem is only mildly nonlinear (for instance, if no limit points
are expected), there is no need to limit the step size.

In this example, the delayed-modified Newton method is used in conjunction with acceleration via line search.

The nonlinear analysis is attempted in a single increment, i.e. a
`step_size_init`

of 1.0, which might not always work. In such cases,
the default of 0.1 is more expedient.

Example: Load-Controlled Contact Analysis

```
analysis nonlinear
max_newton_iterations 100
max_divergences 10
tol_residuum 1e-3
tol_solution 1e-3
```

Load control is enabled by default. Due to the highly nonlinear contact conditions, the full Newton method (default) is often more expedient than the accelerated modified Newton method. Because the Newton iterations may converge slowly when contact is active, a high number of Newton iterations is allowed.

The maximum number of divergences is increased, this improves the chances that the Newton iterations will converge.

If the contact conditions are enforced via quadratic penalties, the
tangent stiffness matrix may be, depending on the penalty factor,
badly conditioned. In this case, the tolerances for convergence may
need to be increased using the `tol_residuum`

and `tol_solution`

options.

Example: Mildly Nonlinear Load-Controlled Analysis

A typical static ‘mildly’ nonlinear parameter setting is

```
analysis nonlinear
correction_type accelerated_newton
newton_method delayed_modified
```

Nonlinear heat-transfer analysis is often well-behaved and therefore may be carried out with the above settings.

Load-controlled analysis is enabled by default. Because the problem is considered mildly nonlinear, there is no need to limit the step size, and the delayed-modified Newton method (or the modified Newton method) can used in conjunction with acceleration via line search.

Depending on the amount of detail required, the maximum step size may be limited. Depending on the accuracy required, the tolerances for convergence may be modified.

## Dynamic Linear Solver

Linear dynamic (non-stationary) problems are solved with the linear static solver. It is activated for a specific case with the parameter

`analysis dynamic_linear`

of the case block.

Implicit time integration is performed in the linear dynamic linear solver with the Newmark method and with a fixed time step size and with one single stage.

Note

It is not recommended to solve dynamic problems with the dynamic linear solver due to the above restrictions.

By default te time (load) step is 1. However the analysis duration is usually not 1 (time unit), and an analysis stage must be defined in the case definition, scaled by the duration. Refer to the example below below for details.

In addition to the MDL command `analysis nonlinear`

, the following
commands may be specified in the case block.

`gradients VALUE`

Controls the computation of the gradients (strains, stresses, heat transfer, etc.) for the current case. If

`VALUE`

(int) is set to 0 (default) no gradients are computed and saved, i.e. only DOF solutions will be found on the database. If`VALUE`

is set to -1, the gradients for the last converged step of at stage or the linear solution step are computed and saved. Finally, if`VALUE`

is set to a positive value, the gradients will be computed at each istep step of at stage as well as for the last step of a stage.

`rayleigh_alpha|rayleigh_beta VALUE`

Specifies the coefficients (float) for Rayleigh damping (default is 0). If

`rayleigh_alpha`

and/or`rayleigh_beta`

are set, Rayleigh damping will be performed where the global viscosity matrix C is a linear combination of the global mass and linear global stiffness matrices M and K:\[C = \alpha M + \beta K\]

`step_size VALUE`

Defines the step size (float) for the stage. No default is given.

Example

A structure is fixed by means of an mdl.ebc (essential boundary condition) while the dynamic load is applied by means of a nbc (natural boundary condition), here a pressure load. The Rayleigh damping coefficients have been calculated from a free-vibration analysis. The analysis time is 10 s, and a step load is applied during the first second.

```
(duration=10.)
case 1
analysis dynamic_linear
stage 2 sfactor (duration)
end
stage 2
ebc 1
nbc 1 sfunction 't<1?1:0'
step_size (0.001*duration)
rayleigh_alpha +6.4e-02
rayleigh_beta +9.0e-06
end
adir
case 1
end
```

The analysis time of 10 s is realized by having two cases, a case (case 1) referencing a stage and a stage (stage 2) describing the pulse. The step load is realized with the expression ‘t<1?1:0’, which means that for t (the analysis time) between 0 and 1, the load is applied with the factor 1, and for all the other time values the load factor is 0.

## Dynamic Nonlinear Solver (B2000++ Pro)

Nonlinear dynamic (or non-stationary) problems are solved with the
dynamic nonlinear solver. It is activated for a specific
case with the case parameter
`analysis dynamic_nonlinear`

.

The dynamic nonlinear solver adopts for the implicit time integration the backward differential formula (BDF) method which is well suited for stiff ordinary differential equations [Sham79] like the ones encountered in stress analysis. The second-order BDF is applied by default, giving unconditional nonlinear stability. For computational efficiency, the time step size is adapted during the analysis using the Nordsieck transformation method [Nord62]. For stability and accuracy of the time integration method, the rate of increment of the step size is limited [Calvo87] and the step size is controlled using Milne’s device error estimator [Milne70].

The solver is invoked by the MDL case parameter
`analysis dynamic_nonlinear`

.

In general, the analysis duration is not 1 (second), and an analysis stage must be defined in the case definition, scaled by the duration. Refer to the first of the examples below for details.

### Applications

In stress analysis, the dynamic solver can be applied to problems involving inertia forces or viscous forces. Like the static nonlinear solver, the problem to be solved may involve any kind of geometric, material, or boundary nonlinearities. Typical uses are:

Relaxation and creep analysis, in conjunction with visco-elastic materials

Dynamic buckling analysis with or without artificial damping

Quasi-static buckling analysis with artificial damping

In heat-transfer analysis, typical applications are nonlinear loading (e.g. cyclic loading) and predicting the non-stationary behavior (i.e. the evolution of the temperature field). Like the static nonlinear solver, the problem to be solved may involve any kind of material or boundary nonlinearity.

Optional Solver MDL Parameters

In addition to the MDL command `analysis nonlinear`

, the following
commands may be specified in the case block.

`correction_termination_test flux_normalised|dof_and_residue|normalised_dof_and_residue`

Specifies the type of test which indicates whether the newton iterations have sufficiently converged. Default is

`flux_normalised`

.The

`flux_normalised`

correction termination test computes the residuum normalised by the energy flux for each type of field (displacements, rotations, Lagrange multipliers, etc.) separately. The`flux_normalised`

test has additional options, see Flux-normalised correction termination test.The

`dof_and_residue`

correction termination test tests the L2-norm residuum against`tol_residuum`

and the L2-norm of the solution changes w.r.t. the previous Newton iteration against`tol_solution`

. The following tolerances are used:`tol_residuum`

(default 1e-3) and`tol_solution`

(default 1e-3).The

`normalised_dof_and_residue`

correction termination test divides the L2-norm residuum by that of the first Newton iteration. The same is done for the L2-norm of the solution changes w.r.t. the previous Newton iteration. The following tolerances are used:`tol_residuum`

(default 1e-3) and`tol_solution`

(default 1e-3).

`gradient VALUE`

Controls the computation of the gradients (strains, stresses, heat transfer, etc.) for the current case. If `` VALUE`` (int) is set to 0 (default) no gradients are computed and saved, i.e. only DOF solutions will be found on the database. If `` VALUE`` is set to -1, the gradients for the last converged step of at stage or the linear solution step are computed and saved. If `` VALUE`` is set to a positive value, the gradients will be computed at each istep step of at stage as well as for the last step of a stage.

`line_search YES|NO`

Specifies whether line search should be used. Default is

`NO`

. This line search method satisfies the strong Wolfe condition. When the cost of calculating the first variation is significantly lower than the cost of calculating the second variation, it is recommended to activate line search, as it improves the probability that the Newton iterations converge.

`max_divergences VALUE`

Defines the maximum number of consecutive divergences (int) before the step is reduced during a load step (optional). Divergence occurs when the energy of the current correction is larger than the energy of the preceding correction. Default is 4.

`max_newton_iterations VALUE`

Maximum number of Newton iterations (int) per step. Default is 50.

`max_steps VALUE`

Maximum number of steps or increments (int) per stage. Default is 99999.

`multistep_integration_order VALUE`

Integration order (int) of the implicit transient integration scheme. Currently supported values are 1 (Euler backward difference formula), 2, 3, 4, 5, or 6. The solver is unconditionally nonlinear-stable for values of 1 and 2. Default is 2.

`newton_method modified|delayed_modified|conventional`

Specifies what kind of Newton method should be used. Default is

`conventional`

.

`newton_periodic_update VALUE`

For the modified and delayed modified Newton methods, this option specifies the maximum number of consecutive iterations (int) without re-factorization. When this number is reached, a re-factorization will be performed, and the counter is reset. Note that for high values of

`newton_periodic_update`

and depending on the`newton_method`

setting, several increments may be performed until this number of Newton iterations is reached. Default is 100.

`residue_function_type order_n|rayleigh_damping|artificial_damping`

Specifies the residue method (string). Default is

`order_n`

, for which the MDL command`multistep_integration_order`

defines the order of the implicit time-integration.The

`artificial_damping`

method works like in the static nonlinear solver; the damping is of the Rayleigh viscous type. Inertia effects are neglected, and the problem becomes a first-order ODE. The global viscosity matrix is generated in function either of the`dissipated_energy_fraction`

option or the`rayleigh_alpha`

and`rayleigh_beta`

options. In the first case, the global viscosity matrix is a scaled instance of the mass matrix, such that at the first increment, the ratio of the dissipated energy to the total energy corresponds to the value given in`dissipated_energy_fraction`

. In the second case, the global viscosity matrix is constructed in the same way as for`rayleigh_damping`

(see below).For the

`rayleigh_damping`

(second-order) method, the global viscosity matrix C is constructed from the global mass matrix M and the global linear stiffness matrix K in the following way:\[C = \alpha M + \beta K\]The coefficient \(\alpha\) can be specified with the MDL command

`rayleigh_alpha`

, and likewise, the coefficient \(\beta\) can be specified with the MDL command`rayleigh_beta`

. Note that these coefficients can be determined in the same way as explained in Artificial damping section for the static nonlinear solver.

`step_size_ini VALUE=v`

Step size (float) at the beginning of the stage. Default is 0.1.

`step_size_max VALUE`

Maximum step size (float) during the stage. Default is 1.0.

`step_size_min VALUE`

Minimum step size (float) during the stage. Default is 1e-12.

`tol_dynamic VALUE`

The dynamic nonlinear solver controls the step size with a local error estimator, which is obtained with Milne’s method. VALUE` (float) defines the maximum permitted absolute error of a degree of freedom during one time step. Default is 1e-5. The step size is adapted accordingly to this criteria. Setting

`VALUE`

to a large value will inactivate the transient integration error control.

`tol_residuum VALUE`

Error tolerance value for the residuum (‘equilibrium’) error

`VALUE`

(float) used by the Newton correction termination test. The default depends on the`correction_termination_test`

option.

`tol_solution VALUE`

The error tolerance value for the dof (‘displacement’) error

`VALUE`

(float) used by the Newton correction termination test. The default depends on the`correction_termination_test`

option.

`tol_work VALUE`

The error tolerance value for the work error

`VALUE`

(float) used by the Newton correction termination test. The default depends on the`correction_termination_test`

option.

Example: Dynamic buckling stress analysis with Rayleigh damping

In this example, a panel is compressed in the x-direction by 3mm within 0.5s. If the analysis duration is different from the value of 1 (in time units, usually seconds), a separate stage has to be defined, stage 2 in the present example.

```
ebc 1
value 0. dof [UX UY UZ RX RY RZ] nodeset nodes_left
value 0. dof [ UY UZ RX RY RZ] nodeset nodes_right
value -1. dof [UX ] nodeset nodes_right
end
end
(max_displacement=3.)
(duration=0.5)
case 1
analysis dynamic_nonlinear
stage 2 sfactor (duration)
end
stage 2
ebc 1 sfactor (max_displacement/duration)
step_size_init (0.1*duration)
tol_dynamic 1e-2
residue_function_type rayleigh_damping
rayleigh_alpha 123.6
rayleigh_beta 0.00008074
end
adirr
case 1
shell_intersection_angle 10
end
```

The duration of the analysis is given with the MDL command

```
stage 2 sfactor f
```

The amplitude of any specified boundary conditions in the stage must be
inversely-scaled by the analysis duration, this is done with the
`sfactor`

option. On the other hand, the sizes for the minimum,
maximum, and initial step are given in absolute values. Hence,
`step_size_init`

should always be defined.

The dynamic tolerance (`tol_dynamic`

) option has a significant
influence on both accuracy and numerical effort. A large dynamic
tolerance leads to relatively large time step sizes, requiring many
Newton iterations per increment. Using a dynamic tolerance that yields
smaller time step sizes which require only one or two Newton iterations
is not only more accurate but often more effective. In some cases, a
modified Newton method may further reduce the computational effort.
Conversely, when the dynamic tolerance is set to an overly conservative
value, the resulting time step size may become extremely small in the
presence of high-frequency vibrations and fast movements.

When geometric instabilities are present, Rayleigh damping should be
selected (by setting the `residue_function_type`

option to
`rayleigh_damping`

and by setting the values for `rayleigh_alpha`

and `rayleigh_beta`

) to ensure an effective solution process.

For highly nonlinear problems, the full Newton method (default) is often more expedient than the modified Newton method. For mildly nonlinear problems and for problems where the time step size is predominantly limited by the time integration error, a modified Newton method may be more effective. The default maximum number of Newton iterations (50) should be sufficient in many cases. Lower values may cause many small time increments or even failure to converge, while higher values may lead to increased computation time. In case of very slow but reliable convergence, higher values (100-200) may be necessary.

Example: Dynamic buckling stress analysis with visco-elastic materials

In this example, the options for Rayleigh damping are omitted since all dissipation should come from the visco-elastic material. The dynamic tolerance is crucial for the accuracy of the time integration: If the time step sizes are too large, the numerical damping which is caused by the time integration error may be larger than the material damping.

```
(max_displacement=3.)
(duration=0.05)
case 1
analysis dynamic_nonlinear
stage 2 sfactor (duration)
end
stage 2
ebc 1 sfactor (max_displacement/duration)
step_size_init (0.1*duration)
tol_dynamic 1e-4
end
```

Example: Quasi-static buckling stress analysis with artificial damping

The dynamic solver can also be applied to load-controlled quasi-static analysis. In this case, the analysis duration is irrelevant, and a separate stage does not need to be defined. Also in this case, the dynamic tolerance controls the time step size and is important to computational efficiency and accuracy.

```
case 1
analysis dynamic_nonlinear
ebc 1
tol_dynamic 1e-3
residue_function_type artificial_damping
dissipated_energy_fraction 1e-3
end
```

Energy is dissipated with artificial damping; the
`residue_function_type artificial_damping`

option also tells the
solver to neglect inertia forces. The amount of damping is controlled
with the `dissipated_energy_fraction`

parameter. High values
accelerate the analysis, increasing the increment sizes, but may
influence the results on the other hand. Therefore, a compromise between
analysis time and accuracy has to be found. While the parameter value is
independent of the physical units, the necessary amount of damping
depends on the problem. Thus, the default of 1e-4 may not be
appropriate. It is recommended to vary this parameter in powers of 10
(e.g. make an analysis for 1e-7, 1e-6, 1e-5, 1e-4, 1e-3) and observe its
effect.

Example: Mildly nonlinear heat-transfer analysis

An initial temperature field is given (`dof_init`

), and a constant
temperature is applied at the boundary (`ebc 1 sfunction '1'`

). If
the problem is well behaved, it can be solved with a higher-order BDF
that is more accurate but not unconditionally stable, and a modified
Newton method can be used.

The analysis duration is one year (365 days), and the time step size is fixed to one day (84600 seconds). A very large dynamic tolerance is used to inactivate the solver’s error estimator.

```
(one_day=84600)
(duration=365*one_day)
case 1
analysis dynamic_nonlinear
stage 2 sfactor (duration)
end
stage 2
dof_init 1
ebc 1 sfunction '1'
multistep_integration_order 4
newton_method delayed_modified
step_size_min (one_day)
step_size_init (one_day)
step_size_max (one_day)
tol_dynamic 1e6
end
```

## Undamped Free Vibration Solver

Free vibration problems are solved with the frequency-independent free
vibration analysis solver.It is activated for a specific
case with the case parameter
`analysis=free_vibration`

.

Note

When executing B2000++, errors of the following kind

```
CRITICAL:all:16:22:8.124: Uncovered b2000 exception caught in the main function.
Exception: Error with the arpack subroutine aupd, info = -9.
```

may indicate that material densities are not specified or set to 0.

Warnings of the following kind

```
WARNING:solver.eigenvalue.arpack:12:36:35.234: Maximum number of iterations taken.
All possible eigenvalues (0) have been found.
```

may indicate that the number of eigemnodes `nmodes`

to be
computed option should be increased.

Optional Solver MDL Parameters

`analysis nonlinear`

, the following
commands may be specified in the case block.

`m_orthonormal yes|no`

By default the mode shapes are normalized with respect to the mass (“M-orthonormalized”), see [Bathe96]. If mode shapes should be normalized to the maximum component of each mode (infinity norm), set

`m_orthonormal`

to`no`

.

`nmodes VALUE`

Number of free-vibration modes

`VALUE`

(int) to be computed. Must be greater or equal to 2. It is often necessary to set this option to values of 10 or even 20 for the successful solution of the eigenvalue problem. Default is 2.

`shift VALUE`

The spectral frequency (not the circular frequency!) shift value

`VALUE`

(float). The computed free-vibration modes are the ones whose eigen-frequency is near this value. For complex problems, a complex spectral frequency shift value can be provided. Default is 0.

`which_eigenvalues=sm|la|ra`

Indicates where to search for eigenvalues with respect to 0 or to the current shift value.

`sm`

searches for the smallest eigenvalues around 0 or the current shift value (default).`la`

searches for the smallest eigenvalues to the left of 0 or the current shift value.`ra`

searches for the smallest eigenvalues to the right of 0 or the current shift value.

`matrices_dir VALUE`

If set to a non-empty string, the free vibration solver will output the system matrices as CSV formatted files in the directory provided by the string

`VALUE`

. The directory is created if it does not exist yet and populated by the respective system matrices. The provided path may be absolute or relative, relative directories are created in the working directory of the solver. By default an empty string is used, which deactivates the output. To activate the output but write the files directly into the working directory, provide a string like`"./"`

for this option.

## Frequency-Dependent Undamped Free-Vibration Solver (B2000++ Pro)

Linear static (stationary) or linearized problems are solved with the
linear static solver. It is activated for a specific
case with the case parameter
`analysis=linear`

.

The frequency-dependent undamped vibration solver allows for solving
frequency-dependent free vibration eigenvalue problems to determine
eigenfrequencies and eigenmodes, where the linear internal forces and
internal inertia are frequency-dependent. It is invoked by the MDL
command `analysis=frequency_dependent_free_vibration`

; this command
must be specified in the case block.

Note

When executing B2000++, errors of the following kind

```
CRITICAL:all:16:22:8.124: Uncovered b2000 exception caught in the main function.
Exception: Error with the arpack subroutine aupd, info = -9.
```

are often indicative of missing density specifications in the material definitions.

When executing B2000++, warnings of the following kind

```
WARNING:solver.eigenvalue.arpack:12:36:35.234: Maximum number of iterations taken.
All possible eigenvalues (0) have been found.
```

are indicative that the value for the `nmodes`

option should be
increased.

Optional Solver MDL Parameters

In addition to the MDL command ```
analysis
frequency_dependent_free_vibration
```

, the following commands may be
specified in the case block.

`nmodes VALUE`

Number of free-vibration modes

`VALUE`

(int) to be computed (default: 1). It is often necessary to set this option to values of 10 or even 20 for the successful solution of the eigenvalue problem.

`shift VALUE`

The spectral frequency shift value

`VALUE`

(float). The computed free-vibration modes are the ones whose eigenfrequencies are near this value (default=0). For complex problems, a complex spectral frequency shift value can be provided.

`freq_init VALUE`

The initial frequency

`VALUE`

(float) used in the first stage of the algorithm.

`tol VALUE`

The relative tolerance

`VALUE`

(float) of the computation of the exact eigenfrequency.

`tol_predictor VALUE`

The relative tolerance

`VALUE`

(float) for the re-computation of the approximated eigenfrequencies.

## Linearized Pre-Buckling Solver

Linear static (stationary) or linearized problems are solved with the linear static solver. It is activated for a specific case with the

`analysis linear`

of the case block.

The linearized pre-buckling solver is commonly used in stress analysis
of slender or thin-walled structures, it allows for the direct solving
of ‘buckling’ problems around an equilibrium position along the
nonlinear path. The solver is invoked by the MDL command
`analysis=linearised_prebuckling`

, which must be specified in the
case block.

Note

When executing B2000++, warnings of the following kind

```
WARNING:solver.eigenvalue.arpack:12:36:35.234: Maximum number of iterations taken.
All possible eigenvalues (0) have been found.
```

are indicative that the value for the `nmodes`

option should be
increased.

Optional Solver MDL Parameters

`nmodes VALUE`

Number of buckling modes v to be computed (default: 1). It is often necessary to set this option to values of 10 or even 20 for the successful solution of the eigenvalue problem.

`shift VALUE`

The spectral frequency shift value. The computed buckling modes are the ones whose eigenvalues (load factors) are near this value. The default of 0 selects the

`nmodes`

lowest buckling modes.

## Sparse Linear Equation Solvers

During the analysis of linear and nonlinear FE problems, systems of linear equations must be solved. In most cases, the matrices are sparse, and often, they are symmetric as well. The solution of these systems of linear equations constitutes an important part of the FE solving procedures and as such has a great influence on the robustness and effectiveness of any FE analysis.

B2000++ contains several sparse linear equation solvers. By default,
the specific *FE* solvers (linear, nonlinear, static, dynamic)
automatically select the proper sparse linear equation solver. This
selection can be inactivated by forcing a specific sparse linear
equation solver by means of the `sparse_linear_solver`

option of the

Optional Solver MDL Parameters

The following commands may be specified in the case block.

`sparse_linear_solver=dmumps|dmumps_dp|dll|dmf|dmf_ooc|dmf_ooc_nommap|dpastix|icg|igmres`

Selects the type of sparse linear solver (string):

`dmumps`

Selects the MUMPS direct multi-frontal solver with pivoting, for use with general symmetric matrices. This solver is the default solver. It is particularly robust and effective on matrices containing Lagrange multipliers (that is, whenever nonlinear constraints are present in the FE model). It has out-of-core support. The sequential version is available on all platforms, the parallel version is available on selected Linux distributions, such as RHEL. Fedora 21.

The parallel version is activated via the

`mpirun`

command. Example:mpirun -np 8 b2000++ DBNAMEThe following options are available:

`autospc=yes|no`

Enables fixing of singular or nearly-singular degrees-of-freedom. Default is

`no`

. This option may be useful for example in conjunction with Lagrange multipliers, or shell elements with free 6th degrees-of-freedom.The

`autospc`

option should be used with caution as it may significantly influence the solution in non-foreseeable ways, especially in free-vibration analysis and in linearized buckling analysis. For shell FE models, it is recommended to use artificial drill stiffness instead.

`autospc_value VALUE`

If

`autospc`

is enabled, use`VALUE`

(float) to fix singular or nearly-singular degrees-of-freedom. A value of 0.0 indicates that the autospc value should be determined automatically. Default is 0.0.

`autospc_threshold VALUE`

If

`autospc`

is enabled, use`VALUE`

(float) to determine singular or nearly-singular degrees-of-freedom. A value of 0.0 indicates that the autospc threshold will be chosen automatically. Default is 0.0.

`compute_null_space_vector=n`

If

`autospc`

is enabled, compute a maximum of`n`

(int) first null space vectors (zero-energy modes) and store them on the database as`DISP_NS`

datasets(not documented). Default is 0.Use of this option requires sometimes to adjust

`autospc_threshold`

.The following example is a single 2D quadrilateral element. No degrees-of-freedom are constrained, hence 3 zero-energy modes are present.

nodes 1 0.0 0.0 0.0 2 1.0 0.0 0.0 3 1.0 1.0 0.0 4 0.0 1.0 0.0 end elements type Q4.S.2D.TL mid 1 thickness 0.1 1 1 2 3 4 end material 1 type isotropic e 73.e3 nu 0.3 end case 1 analysis linear sparse_linear_solver dmumps autospc yes autospc_threshold 1.e-15 compute_null_space_vector 3 end adir case 1 end

`mumps_ooc=yes|no`

Enables the out-of-core solver.

`mumps_ooc_tmpdir=d`

When the out-of-core solver is active,

`d`

(string) ist the directory path name where the temporary files should be created. By default, the`/tmp`

directory is used.

`mumps_ooc_prefix=d`

When the out-of-core solver is active, specifies

`d`

(string) as prefix for the temporary files.

`mumps_icntl VALUE`

Sets one or several values in the ICNTL (int) array, which is used to control the behavior of the MUMPS solver. The string

`VALUE`

consists of comma-separated`key:value`

pairs, where`key`

is the index into the ICNTL array, starting from 1, and`value`

is an integer number.The following example enforces the use of the SCOTCH library for pivoting by setting “ICNTL(7)=3”:

case 1 analysis linear sparse_linear_solver dmumps mumps_icntl "7:3" endTo also use iterative refinement with a fixed number of 3 iterations use:

mumps_icntl "7:3,10:-3"Only in rare cases it is necessary to set ICNTL, the defaults that are set by B2000++ are usually satisfactory or are supported directly. While ICNTL allows to fine-tune the workings of MUMPS, not all options will work with B2000++, and the options

`5`

(matrix format) and`18`

(matrix distribution) are forbidden, as they are determined by B2000++ itself. Refer to the MUMPS user’s guide, section 5.1, for an explanation of ICNTL.

`mumps_cntl VALUE`

Sets one or several values in the CNTL (floating-point) array, which is used to control the behavior of the MUMPS solver. The string

`VALUE`

consists of comma-separated`key:value`

pairs, where`key`

is the index into the CNTL array, starting from 1, and`value`

is a floating-point number.Only in rare cases it is required to set CNTL, the defaults that are set by b2000++ are usually satisfactory or are supported directly (e.g.

`autospc_threshold`

and`autospc_value`

).One setting that may be of particular interest is

`CNTL(3)`

, which controls the threshold for identifying null pivots. Since MUMPS 5.4.0 the criterion for the default setting`CNTL(3)=0`

defines a threshold relative to the sup norm of the preprocessed matrix based on the machine precision and a factor dependent on the number of variables on the deepest elimination branch. When provided a positive value, this one is used instead of the machine precision and factor for the calculation and a negative value is used directly as an absolute threshold.If the simulation aborts due to a null pivot it may help to set the threshold to a value smaller than machine precision. As it is possible to correctly represent 18 decimal digits in the 80 bit extended precision of x86 processors, we suggest to try

`mumps_cntl "3:1.e-18"`

as a starting point if the default threshold does not work.Note: In older versions of MUMPS (before 5.4.0) the default value was a factor of 1.e-5 below machine precision. To emulate that behaviour set

`mumps_cntl "3:1.1e-21"`

Refer to the MUMPS user’s guide, section 5.2, for an explanation of CNTL.

`dmumps_dp`

Direct multi-frontal solver without pivoting, for use with symmetric positive definite matrices, making use of the MUMPS library. It has out-of-core support. Since the matrices have to be symmetric positive definite, this solver cannot be used in FE analysis involving Lagrange multipliers. This solver is moderately faster than

`dmumps`

and supports the same set of options.

`dll`

Sequential direct super-nodal left-looking sparse solver for symmetric problems. The matrices need not be positive definite, however, no zero elements on the diagonal must be present.

`dmf`

Parallel direct multi-frontal sparse solver for symmetric problems. The matrices need not be positive definite, however, no zero elements on the diagonal must be present.

`dmf_ooc`

Direct multi-frontal, out-of-core, sparse symmetric solver. This solver is similar to

`dmf`

, but can also be used for FE analyses where the problem’s memory footprint exceeds the available memory. The out-of-core data is stored in a temporary file and is mapped into memory using UNIX’smmapmechanism.

`dmf_ooc_nommap`

Similar to

`dmf_ooc`

, but less efficient. Makes use of a regular file instead of ammap’edfile. To be used on platforms wheremmapis not reliable.

`dpastix`

Parallel direct left-looking super-nodal solver with static pivoting, for use with positive-definite symmetric matrices, making use of the the PASTIX library.

In B2000++, the solver is executed with multi-threading enabled and MPI and out-of-core support disabled. This solver is available for Centos.

`icg`

or`igmres`

Iterative conjugate-gradient (

`icg`

) and GMRES (`igmres`

) solver for use with symmetric resp. non-symmetric problems, making use of the GMM++ library. Note that the successful use of iterative solvers largely depends on the type of problem and the preconditioner.The following options are available:

`sls_tol_residuum`

Specifies the tolerance (default 1.0e-8) for the residuum. When the residuum becomes lower than this tolerance, the iterative solution process is considered converged.

`sls_max_iter`

Specifies maximum number of iterations. By default, an infinite amount of iterations is permitted.

`sls_precond`

Specifies the type of preconditioner to use. For the symmetric solver (

`icg`

) , the`ildlt`

(default),`ildltt`

, and`diagonal`

preconditioners are available. For the unsymmetric solver (`igmres`

), the`ilu`

(default),`ilut`

,`ilutp`

,`diagonal`

, and`mr_approx_inverse`

preconditioners are available.`sls_restart`

Specifies after how many iterations the preconditioner should be applied again. Default is 200.

`sls_verbosity`

Specifies the amount of information that GMM++ should print to the console. A value of 0 (default) means no information, a value of 1 will print the iterations, and a value of 2 will print iterations and sub-iterations.

`unsymmetric=yes|no`

Specifies whether an unsymmetric sparse linear solver should be used. If the problem is per se symmetric, default is

`no`

. This option is normally not necessary as B2000++ automatically detects whether the problem is unsymmetric and selects the appropriate sparse linear solver.`complex yes|no`

Specifies whether a complex sparse linear solver should be used. For non-complex problems, default is

`no`

. This option is normally not necessary as B2000++ automatically detects whether the problem is complex and selects the appropriate sparse linear solver.

## Linear and Nonlinear Constraint Control

Constraints originate from prescribed conditions such as mdl.ebc, linc, field_transfer or rigid-coupling (RBE elements) or distributed-coupling (FMDE) elements.

B2000++ supports different ways for imposing constraints on a FE problem. Constraints may be linear or nonlinear, and may be imposed by boundary conditions and by Finite Elements. For many problems, the default settings (reduction for linear constraints and Augmented Lagrange method or Lagrange method for nonlinear constraints) will suffice to effectively obtain a solution. In this chapter, the MDL commands for overriding these settings are described.

Optional Solver MDL Parameters

The following commands may be specified in the case block.

`constraint_method=reduction|lagrange|penalty|augmented_lagrange`

Selects the type of method used to impose linear constraints on the equations.

`reduction`

The system of linear equations is modified such that the constraints are eliminated. The constraints are satisfied ‘exactly’ (within numerical precision). This is the default method, and it is applicable to linear and nonlinear problems, but only for

linear constraints.

`augmented_lagrange`

A combination of Lagrange multipliers and penalties is used to enforce the constraints. The MDL command

`constraint_penalty`

specifies the penalty factor (default value 1, it is rarely necessary to change this value). The penalty factor is solely used to ensure the positive-definiteness of the matrix and does not impair the matrix’ condition number. The constraints are, due to the Lagrange multipliers, satisfied ‘exactly’ (within numerical precision).Linear systems containing Lagrange multipliers are best solved using a pivoting solver. It is recommended to use

`sparse_linear_solver=dmumps`

, which is the default.

`penalty`

With the penalty method, the number of degrees-of-freedom remains unchanged. The constraints are imposed approximately, where the MDL command

`constraint_penalty`

specifies the penalty factor (default value is 1e10). A high penalty factor ensures a close approximation of the constraints but makes the matrix badly conditioned; while a low penalty factor may introduce artifacts in the solution of the FE analysis, especially for eigenmode problems.

`lagrange`

To the system of linear equations, Lagrange multiplicators are added. The constraints are satisfied ‘exactly’ (within numerical precision).

Linear systems containing Lagrange multipliers are best solved using a pivoting solver. It is recommended to use

`sparse_linear_solver dmumps`

, which is the default.

`nonlinear_constraint_method=lagrange|penalty|augmented_lagrange`

Selects the type of method used to impose nonlinear (linearized) constraints on the equations.

`augmented_lagrange`

A combination of Lagrange multipliers and penalties is used to enforce the constraints. The MDL adir command

`constraint_penalty`

specifies the penalty factor (default value 1, it is rarely necessary to change this value). The penalty factor is solely used to ensure the positive-definiteness of the matrix and does not impair the matrix’ condition number. The constraints are, due to the Lagrange multipliers, satisfied exactly.Linear systems containing Lagrange multipliers are best solved using a pivoting solver. It is recommended to use

`sparse_linear_solver dmumps`

, which is the default.

`penalty`

With the penalty method, the number of degrees-of-freeom remains unchanged. The constraints are imposed approximately, where the MDL adir command

`constraint_penalty`

specifies the penalty factor (default value is 1e10). A high penalty factor ensures a close approximation of the constraints but makes the matrix badly conditioned; while a low penalty factor may introduce artefacts in the solution of the FE analysis, especially for eigenmode problems.

`lagrange`

To the system of linear equations, Lagrange multiplicators are added. The constraints are satisfied exactly. Since for the static case, the matrix is no longer positive definite, an unsymmetric solver should be used.

Linear systems containing Lagrange multipliers are best solved using a pivoting solver. It is recommended to use

`sparse_linear_solver=dmumps`

, which is the default.

## Flux-normalised Correction Termination Test

This is the default test to determine convergence of the Newton iterations. The flux can be considered a measure of the individual internal and external forces in the FE model. For each Newton iteration, the instantaneous flux is normalised by the average flux. The average flux is computed from the converged solutions of the previous steps and from the solution of the current Newton iteration.

For all fluxes, and for the Newton correction of the solution, the inf-norm (maximum) is used.

The flux normalization is performed per field, because each field may have a different unit and a different range of values. The finite elements present in the analysis model define what fields are active. The following fields are commonly used in elastic analysis and heat-transfer analysis:

Field name |
Description |
---|---|

DISPLACEMENT |
The displacement degrees-of-freedom. |

ROTATION |
The rotational degrees-of-freedom in shell and beam elements and in rigid-coupling elements. |

LAGRANGE |
The degrees-of-freedom associated to Lagrange multipliers. |

TEMPERATURE |
The degrees-of-freedom in heat-transfer elements. |

The following case parameters are active for all fields:

`max_newton_iterations VALUE`

Specifies the maximum number of Newton iterations

`VALUE`

(int) per increment. Default is 50.

`max_divergences VALUE`

Defines the maximum number of consecutive divergences

`VALUE`

(int) during the correction phase of a load step. Default is 4. When this number is reached, the correction phase is terminated unsuccessfully, and the step increment is re-started with a reduced step size.

`min_newton_iterations_check_divergence VALUE`

Checks divergence only after

`VALUE`

(int) Newton iterations. Default is 4.

The following options can be specified either for all fields or per field. In the second case, the field name must be appended to the relevant case parameter. Example:

```
case 1
ebc 1
analysis= nonlinear
tol_residuum= 1.e-3 # for all fields
tol_residuum.rotation= 1.e-5 # only for the rotation field
end
```

`tol_residuum v`

Threshold for convergence: The maximum flux divided by the average flux must be smaller than

`VALUE`

. Default is 1e-5.

`tol_solution v`

Threshold for convergence: The maximum of the Newton correction, divided by the maximum of the current step increment must be smaller than

`VALUE`

. Default is 1e-5.

`tol_residuum_linear_problem v`

Threshold for linear convergence: If the maximum flux divided by the average flux is smaller than

`VALUE`

, linear convergence is assumed, without consideration of`tol_solution`

. Default is 1e-8.

`nb_iter_nonquadratic_convergence v`

Number of Newton iterations that will be executed until the convergence history will be examined to determine whether quadratic convergence has been obtained. Default is 9 iterations.

`tol_residuum_nonquadratic_convergence v`

If convergence is not quadratic, after

`nb_iter_nonquadratic_convergence`

Newton iterations,`tol_residuum`

is replaced by`VALUE`

. Default is 2e-2.

`average_flux v`

Use

`VALUE`

to normalise the flux, instead of the averaged flux which is computed from the analysis.

`initial_average_flux v`

Use

`VALUE`

as start value to compute the averaged flux. Default is 0.01.

`zero_flux_criterion v`

If the magnitude of the instantaneous flux divided by the averaged flux is smaller then

`VALUE`

, the flux will be considered zero, and`tol_solution_zero_flux`

will be used instead of`tol_solution`

. Default is 1e-5.

`tol_solution_zero_flux v`

Replaces

`tol_solution`

for fields having zero flux. Default is 1e-5.

`zero_flux_relative_criterion v`

For the computation of the normalised flux, discard regions in the model where the flux is below

`VALUE`

. Default 1e-5.