Python
This is the documentation for all Python solver wrappers. Currently, only solvers exist for the one-dimensional (1D) calculation of a straight flexible tube.
Tube
There are three tube solvers for a straight tube with circular cross-section.
The axial direction is along the z-axis.
All of them are 1D solvers.
This means the tube is divided in m intervals in the axial z-direction
and only variations in that direction are taken into account.
In other words, all variables are uniform over a cross-section.
They are calculated in the center of each of the m cells.
There is one flow solver SolverWrapperTubeFlow and two structural solvers, one with inertia SolverWrapperTubeStructure and one without SolverWrapperTubeRingmodel.
These solvers are very simple and provide insight in the physics, especially in the stability of fluid-structure interaction simulation.
Nonetheless, they are not meant to provide an accurate representation of reality.
Especially the pressure stabilization term in the flow solver smooths the pressure distribution.
Settings
The following parameters, listed in alphabetical order, need to be provided in the main JSON parameter file as settings of the solver wrapper.
| parameter | type | description |
|---|---|---|
debug |
bool | (optional) Default: false. For every iteration, text files are saved with the input and output data of the solver. |
delta_t |
double | Fixed time step size. This parameter is usually specified in a higher component. |
input_file |
str | (optional) Name of the input file, which may be present in the folder given in working_directory. The file contains parameters required for the solver, in JSON-format. The parameters specified in the main parameter JSON file have priority over the parameters defined in this file. |
interface_input |
list | List of dictionaries; each dictionary requires two keys: model_part and variables. The former contains the name of the ModelPart as a string. The value of the latter is a list of variables. Even if there is only one variable, a list is required. For the Python solver wrappers these variables are fixed: ['displacement'] for a flow solver and ['pressure','traction'] for a structural solver. |
interface_output |
list | Analogous to interface_input. However, the variables are different: ['pressure','traction'] for a flow solver and ['displacement'] for a structural solver. |
print_coupling_convergence |
bool | (optional) Default false. If true and if the solver coupling convergence is checked, a statement is printed when the solver converges in the first solver iteration, see solver coupling convergence. |
unsteady |
bool | (optional) Default: true. Indicates if case is steady or unsteady. |
save_restart |
int | (optional) Default: 0. Determines the time step interval upon which a pickle file case_timestep<time step>.pickle is written, to be used for restart purposes. A minus sign indicates only the file from the last interval is retained. |
time_step_start |
int | (optional) Default: 0. Time step number to (re)start a transient FSI calculation. If 0 is given, the simulation starts from scratch. Otherwise, the code looks for the pickle file case_timestep<timestep_start>.pickle to start from the corresponding time step. For a steady simulation, the value should be 0. |
working_directory |
str | Absolute path to the working directory or relative path with respect to the current directory. |
delta_t is a necessary parameter, while timestep_start and save_restart are optional, but all are usually defined in a higher component. However, they can also be given directly as parameter of the solver wrapper (e.g. for standalone testing). If they are defined both in higher component and in the solver wrapper, then the former value is used and a warning is printed.
Because these solvers are simple, it is possible to omit the use of interpolation between two solver wrappers.
In that case, the names of the ModelParts of both flow and structural solver need to be the same.
Tube flow solver
This flow solver calculates the flow inside a 1D straight and flexible tube.
The type for this solver wrapper is solver_wrappers.python.tube_flow_solver.
The required input is the radial displacement of the tube wall.
The other components of displacement are ignored.
The resulting output is the pressure on the tube wall.
Although a required variable in interface_output, traction is always returned as being zero.
From the radial displacements the area of each cross-section is calculated. The flow is governed by the continuity and momentum equation $$ \frac{\partial a}{\partial t}+\frac{\partial av}{\partial z}=0 $$ $$ \frac{\partial av}{\partial t}+\frac{\partial av^2}{\partial z}+\frac{1}{\rho_f}\left(\frac{\partial ap}{\partial z}-p\frac{\partial a}{\partial z}\right)=0 $$ with a=\pi r^2 the cross-sectional area of the tube and r the inner radius. Furthermore, t is the time, v the velocity along the axis of the tube, p the pressure and \rho_f the density of the fluid. At the start and at the end of the tube the boundary conditions have to be applied. Discretizing these equations results in a system of linear algebraic equations which can be solved for the pressure and velocity by performing Newton-Raphson iterations. Inside the solver the kinematic pressure is used, which is the pressure divided by the density of the fluid.
The boundary conditions are implemented by four additional equations: at the in- and outlet for both pressure and velocity.
At the inlet, the pressure, velocity or total pressure can be specified.
At the outlet, the pressure can be set to a fixed value or a non-reflecting boundary conditions.
These settings are specified in the input_file in the working directory.
For more information about the implementation of this solver refer to [1].
Finally, this solver also provides a Jacobian of the change in radius with pressure, which can be used as surrogate Jacobian [4].
Solver parameters
The following parameters, listed in alphabetical order, need to be specified in the main JSON file or in a file with name input_file, located in the working_directory.
Care should be taken that the values of d, e, h, l and rhof match the corresponding values of the structural solver.
| parameter | type | description |
|---|---|---|
axial_offset |
double | (optional) Default: 0. Distance over which tube is displaced axially in the coordinate system. |
d |
double | Nominal diameter of the tube. |
e |
double | Modulus of elasticity of the tube wall. |
h |
double | Thickness of the tube wall. |
inlet_boundary |
dict | Dictionary containing all information with respect to the inlet boundary condition. |
l |
double | Length of the tube. |
m |
int | Number of cells for discretization. The values are calculated in the cell centers. |
newtonmax |
double | Maximum number of Newton-Raphson iterations for the flow solver calculation. |
newtontol |
double | Relative tolerance for the Newton-Raphson iterations for the flow solver calculation. |
outlet_boundary |
dict | Dictionary containing all information with respect to the outlet boundary condition. |
preference |
double | (optional) Default: 0. Reference pressure and initial pressure. |
u0 |
double | (optional) Default: ureference. Initial velocity throughout the tube. |
ureference |
double | Reference velocity used for determination of pressure stabilization term. |
rhof |
double | Density of the fluid. |
Inlet boundary
This section describes all parameters that need to be specified in the dictionary inlet_boundary, listed in alphabetical order.
| parameter | type | description |
|---|---|---|
amplitude |
double | Amplitude of the inlet boundary condition. |
period |
double | Period of the inlet boundary condition. Period of oscillation for a periodic boundary condition, duration for a non-periodic boundary condition. Not used for a fixed value boundary condition (type 4). |
reference |
double | (optional) Reference value of inlet boundary condition. If not provided, the value of this parameter is the corresponding reference value provided above, i.e. ureference, preference or preference + rhof * ureference^2 / 2. |
type |
int | Type of inlet boundary condition. If 1, a sine wave is used with amplitude, reference and period as specified. If 2, a pulse is used with amplitude as specified and a duration equal to the parameter period. After the pulse the variable is equal to the reference value. If 3, a quadratic sine wave is used with amplitude, reference and period as specified. If 4, a fixed value equal to the sum of the reference value and the amplitude. Used for steady cases. If other, a steady increase of the value at the inlet with slope of amplitude divided by period is used. |
variable |
str | Variable upon which the inlet boundary condition is defined, either 'pressure', 'velocity' or 'total pressure'. |
Outlet boundary
This section describes all parameters that need to be specified in the dictionary outlet_boundary, listed in alphabetical order.
| parameter | type | description |
|---|---|---|
type |
int | Type of outlet boundary condition. If 1, a non-reflecting boundary condition is applied. This type cannot be used for a steady calculation. If other, a fixed value equal to the reference pressure is applied. |
Tube ringmodel solver
This structural solver calculates the deformation of the wall of a straight and flexible tube.
The type for this solver wrapper is solver_wrappers.python.ring_model_solver.
The tube is regarded as made up of independent rings and no inertia is taken into account.
Therefore, there is no dependence on previous time steps and the parameters delta_t and timestep_start are not used.
The required input is the pressure on the tube wall.
Traction is not taken into account, even though it is a required variable in interface_input.
The resulting output is the radial displacement.
For the other components of displacement, zero is returned.
This solver is not suited to calculate the propagation of a pressure pulse.
The behaviour of the elastic tube wall is described by a Hookean relation, which results in the following equation $$ a=a_0\left(\frac{p_0-2c^2_{MK}}{p_0-2c^2_{MK}}\right)^2 $$ with a=\pi r^2 the cross-sectional area of the tube and r the inner radius. Furthermore, p the pressure, p_0 the reference pressure, \rho_f the density of the fluid and c^2_{MK} the Moens-Korteweg wave speed given by $$ c^2_{MK}=\sqrt{\frac{Eh}{2\rho_f r_0}} $$ Here, E is the modulus of elasticity, h the thickness of the tube wall and r_0 the reference radius. No boundary conditions are required. Inside the solver the kinematic pressure is used, which is the pressure divided by the density of the fluid.
More information about the implementation of this solver can be found in [2].
Solver parameters
The following parameters, listed in alphabetical order, need to be specified in the main JSON file or in a file with name input_file, located in the working_directory.
Care should be taken that the values of d, e, h, l and rhof match the corresponding values of the flow solver.
| parameter | type | description |
|---|---|---|
axial_offset |
double | (optional) Default: 0. Distance over which tube is displaced axially in the coordinate system. |
d |
double | Nominal diameter of the tube. |
e |
double | Modulus of elasticity of the tube wall. |
h |
double | Thickness of tube wall. |
l |
double | Length of the tube. |
m |
int | Number of cells for discretization. The values are calculated in the cell centers. |
preference |
double | (optional) Default: 0. Reference pressure. |
residual_atol |
double | Absolute residual tolerance, only used and required for the solver coupling convergence criterion. |
rhof |
double | Density of the fluid. |
Tube structural solver
This structural solver calculates the deformation of the wall of a straight and flexible tube.
The type for this solver wrapper is solver_wrappers.python.tube_structure_solver.
In this model inertia is taken into account, but still only radial displacement is considered.
The required input is the pressure on the tube wall.
Traction is not taken into account, even though it is a required variable in interface_input.
The resulting output is the radial displacement.
For the other components of displacement, zero is returned.
The deformation of the tube in the radial direction is determined by $$ \rho_s h\frac{\partial^2 r}{\partial t^2}+b_1\frac{\partial^4 r}{\partial z^4}-b_2\frac{\partial^2 r}{\partial z^2}+b_3(r-r_o)=p-p_o $$ with \rho_s the solid density and h the thickness of the tube wall. Further, r is the inner radius, p pressure and t time. The variables r_0 and p_0 are a reference radius and pressure, respectively. The parameters b_1 and b_2 (b_1, b_2 \ge 0) account for the inner action of the bending and the axial tension in the wall, while the parameter b_3 accounts for the circumferential stress in the tube wall. For a thin-walled tube that is clamped in the axial direction, they are defined as $$ b_1=\frac{hE}{1-\nu^2}\frac{h^2}{12} \textrm{, } b_2=\frac{hE}{1-\nu^2}\frac{h^2}{12}\frac{2\nu}{r_o^2} \textrm{ and } b_3=\frac{hE}{1-\nu^2}\frac{1}{r_o^2}+\frac{hE}{1-\nu^2}\frac{h^2}{12}\frac{1}{r_o^4} $$ with E the Young's modulus and \nu Poisson's coefficient. The second term of b_3 is considered small compared to the first one because h\ll r_o and is thus neglected. Discretizing these equations results in a system of linear algebraic equations with a Jacobian that does not depend on the pressure, nor the radius. This system is solved for the radius without requiring iterations.
There are two solving options available, specified with the parameter solver. The fastest and most memory efficient is the "solve_banded" option, which writes the matrix in a diagonal ordered form.
In this way, memory use and the number of operations is reduced.
However, when the number of discretization intervals increases, the condition number of the coefficient matrix quickly increases.
For high condition number, the "direct" method is more accurate and allows deeper convergence.
This option calculates the inverse of the coefficient matrix directly, at the start of the calculation.
The inverse is only calculated once and thereafter stored. When the number of intervals is high, it is clear that this option will use a high amount of memory.
The tube is considered clamped at both ends. This boundary condition is imposed by adding four equations: two at the inlet and two at the outlet.
For more information about the implementation of this solver refer to [1]. In this work the Newmark-beta time discretization is used. Here, however, backward Euler time discretization is used as well.
For the Newmark-beta time discretization, two Newmark parameters \beta and \gamma are required, which result in an unconditionally stable integration scheme if $$ \gamma\geq\frac{1}{2} \textrm{ and } \beta\geq\frac{1}{4}\left(\frac{1}{2}+\gamma\right)^2. $$ Typical values are \gamma equal to 1/2 and \beta equal to 1/4.
A different time discretization for flow and structure can lead to difficulties for strongly coupled problems, especially looking at the resulting pressure distributions.
As most flow solvers are discretized using the backward Euler method, it is advised to choose the same method for the structural solver (the time discretization of SolverWrapperTubeFlowSolver is also backward Euler).
This avoids the occurrence of spurious oscillations of the pressure in time [3].
Finally, this solver also provides a Jacobian of the change in radius with pressure, which can be used as surrogate Jacobian [4].
Solver parameters
The following parameters, listed in alphabetical order, need to be specified in the main JSON file or in a file with name input_file, located in the working_directory.
Care should be taken that the values of d, e, h, l and rhof match the corresponding values of the flow solver.
| parameter | type | description |
|---|---|---|
axial_offset |
double | (optional) Default: 0. Distance over which tube is displaced axially in the coordinate system. |
beta |
double | (optional) Newmark parameter \beta. Only required when the Newmark-beta time discretization is used. |
d |
double | Nominal diameter of the tube. |
e |
double | Modulus of elasticity of the tube wall. |
h |
double | Thickness of the tube wall. |
gamma |
double | (optional) Newmark parameter \gamma. Only required when the Newmark-beta time discretization is used. |
l |
double | Length of the tube. |
m |
int | Number of cells for discretization. The values are calculated in the cell centers. |
nu |
double | Poisson's ratio. |
preference |
double | (optional) Default: 0. Reference pressure. |
residual_atol |
double | Absolute residual tolerance, only used and required for the solver coupling convergence criterion. |
rhof |
double | Density of the fluid. |
rhos |
double | Density of the tube wall. |
solver |
str | (optional) Default: solve_banded. Either solve_banded or direct. Specifies the solution method for the linear system of equations. |
time_disretization |
str | (optional) Default: backward Euler. Specifies the time discretization: either Newmark or backward Euler. Not case sensitive. |
Solver coupling convergence
The convergence criterion solver coupling convergence has been implemented for the Tube Python solvers. The Tube Structure solver and the Tube Ringmodel solver are linear, meaning that no solver iterations are performed. As a consequence no convergence tolerance is prescribed for these solvers.
However, when coupled with another solver, the initial residual (calculated with an updated right-hand side, but before the solution of their linear system) can still have a non-converged value.
This is exactly the residual value that is used by the convergence criterion to determine if the coupling has converged.
Therefore, a solver tolerance residual_atol has to be defined for these solvers in case the solver coupling convergence criterion is selected.
References
[1] Degroote J., Annerel S. and Vierendeels J., "Stability analysis of Gauss-Seidel iterations in a partitioned simulation of fluid-structure interaction", Computers & Structures, vol. 88, no. 5-6, pp. 263, 2010.
[2] Degroote J., Bruggeman P., Haelterman R. and Vierendeels J., "Stability of a coupling technique for partitioned solvers in FSI applications", Computers & Structures, vol. 86, no. 23–24, pp. 2224–2234, 2008.
[3] Vierendeels J., Dumont K., Dick E. and Verdonck P., "Analysis and stabilization of fluid-structure interaction algorithm for rigid-body motion", American Institute of Aeronautics and Astronautics Journal", vol. 43, no. 12, pp. 2549–2557, 2005.
[4] Delaissé N., Demeester T., Fauconnier D. and Degroote J., "Surrogate-based acceleration of quasi-Newton techniques for fluid-structure interaction simulations", Computers & Structures, vol. 260, pp. 106720, 2022.