fronts.solve_flowrate(D, i, Qb, radial, ob=1e-06, angle=6.283185307179586, height=None, itol=0.001, b_hint=None, b_bracket=None, method='implicit', maxiter=100, verbose=0)

Solve a radial problem with a fixed-flowrate boundary condition.

Given a positive function D, scalars \(\theta_i\), \(\theta_b\) and \(o_b\), and coordinate unit vector \(\mathbf{\hat{r}}\), finds a function \(\theta\) of r and t such that:

\[\begin{split}\begin{cases} \dfrac{\partial\theta}{\partial t} = \nabla\cdot\left[D(\theta)\dfrac{\partial\theta}{\partial r} \mathbf{\hat{r}}\right ] & r>r_b(t),t>0\\ \theta(r,0) = \theta_i & r>0 \\ Q(r_b(t),t) = Q_b & t>0 \\ r_b(t) = o_b\sqrt t \end{cases}\end{split}\]
  • D (callable or sympy.Expression or str or float) –

    Callable that evaluates \(D\) and its derivatives, obtained from the fronts.D module or defined in the same manner—i.e.:

    • D(theta) evaluates and returns \(D\) at theta

    • D(theta, 1) returns both the value of \(D\) and its first derivative at theta

    • D(theta, 2) returns the value of \(D\), its first derivative, and its second derivative at theta

    where theta is always a float in the latter two cases, but it may be either a single float or a NumPy array when D is called as D(theta).

    Alternatively, instead of a callable, the argument can be the expression of \(D\) in the form of a string or sympy.Expression with a single variable. In this case, the solver will differentiate and evaluate the expression as necessary.

  • i (float) – Initial condition, \(\theta_i\).

  • Qb (float) –

    Imposed flow rate of \(\theta\) at the boundary, \(Q_b\).

    The flow rate is considered in the direction of \(\mathbf{\hat{r}}\): a positive value means that \(\theta\) is flowing into the domain; negative values mean that \(\theta\) flows out of the domain.

  • radial ({'cylindrical', 'polar'}) –

    Choice of coordinate unit vector \(\mathbf{\hat{r}}\). Must be one of the following:

    • 'cylindrical' : \(\mathbf{\hat{r}}\) is the radial unit vector in a cylindrical coordinate system

    • 'polar' : \(\mathbf{\hat{r}}\) is the radial unit vector in a polar coordinate system

  • ob (float, optional) – Parameter \(o_b\), which determines the behavior of the boundary. It must be positive. The boundary acts as a line source or sink in the limit where ob tends to zero.

  • angle (float, optional) – Total angle covered by the domain. The default is \(2\pi\), which means that \(\theta\) may flow through the boundary in all directions. Must be positive and no greater than \(2\pi\).

  • height (None or float, optional) – Axial height of the domain if radial=='cylindrical'. Not allowed if radial=='polar'.

  • itol (float, optional) – Absolute tolerance for the initial condition.

  • b_hint (None or float, optional) – Optional hint to the solver. If given, it should be a number close to the expected value of \(\theta\) at the boundary (i.e. \(\theta_b\)) in the solution to be found.

  • b_bracket (None or sequence of two floats, optional) – Optional search interval that brackets the value of \(\theta_b\) in the solution. If given, the solver will use bisection to find a solution in which \(\theta_b\) falls inside that interval (a ValueError will be raised for an incorrect interval). This parameter cannot be passed together with a b_hint.

  • method ({'implicit', 'explicit'}, optional) –

    Selects the integration method used by the solver:

    • 'implicit' (default): uses a Radau IIA implicit method of order 5. A sensible default choice that will work for any problem

    • 'explicit': uses the DOP853 explicit method of order 8. As an explicit method, it trades off general solver robustness and accuracy for faster results in “well-behaved” cases. With this method, the second derivative of \(D\) is not needed

  • maxiter (int, optional) – Maximum number of iterations. A RuntimeError will be raised if the specified tolerance is not achieved within this number of iterations. Must be nonnegative.

  • verbose ({0, 1, 2}, optional) –

    Level of algorithm’s verbosity. Must be one of the following:

    • 0 (default): work silently

    • 1: display a termination report

    • 2: also display progress during iterations


solution – See Solution for a description of the solution object. Additional fields specific to this solver are included in the object:

  • o (numpy.ndarray, shape (n,)) – Final solver mesh, in terms of the Boltzmann variable.

  • niter (int) – Number of iterations required to find the solution.

  • b_bracket (sequence of two floats or None) – If available, an interval that contains the value of \(\theta_b\). May be used as the input b_bracket in a subsequent call with a smaller itol for the same problem in order to avoid reduntant iterations. Whether this interval is available or not depends on the strategy used internally by the solver; in particular, this field is never None if a b_bracket is passed when calling the function.

Return type:


See also



This function works by transforming the partial differential equation with the Boltzmann transformation using ode() and then solving the resulting ODE repeatedly with the chosen integration method as implemented in the scipy.integrate module and a custom shooting algorithm. The boundary condition is satisfied exactly as the starting point, and the algorithm iterates with different values of \(\theta\) at the boundary until it finds the solution that also verifies the initial condition within the specified tolerance. Trial values of \(\theta\) at the boundary are selected automatically by default (using heuristics, which can also take into account an optional hint if passed by the user), or by bisecting an optional search interval. This scheme assumes that \(\theta\) at the boundary varies continuously with \(\theta_i\).


[1] GERLERO, G. S.; BERLI, C. L. A.; KLER, P. A. Open-source high-performance software packages for direct and inverse solving of horizontal capillary flow. Capillarity, 2023, vol. 6, no. 2, pp. 31-40.