Using scipy fsolve and odeint for simultaneous ODEs

Question

I'm trying to solve a set of transient flow equations of a body firing from a pressurised tube

from scipy.integrate import odeint
import numpy as np

# parameters
Db = 0.5 # body diameter
Lb = 2.0 # body length
mb = 1000 # body mass
Ai = 0.5 # inlet area
Ki = 1.0 # inlet loss coefficient
Li = 1.5 # inlet position
Dt = 0.6 # tube diameter
Lt = 3.0 # tube length
epsilon = 0.0001 # roughness
rho = 1000 # density
mu = 0.0011 # dynamic viscosity
par = [Db, Lb, mb, Ai, Ki, Li, Dt, Lt, epsilon, rho, mu]

# boundary conditions
pp = 5e5 # pump pressure
patm = 0.0 # atmosphere pressure
bound = [pp, patm]

# time
ts = 0.001 # time step
tmax = 2.0 # total time
t = np.linspace(0, tmax, int(tmax/ts) + 1)

# initial conditions
v0 = 0.0 # initial velocity
s0 = 0.5 # initial displacement
Qlr0 = 0.0 # initial leakage to rear
Qlf0 = 0.0 # initial leakage to front
Qf0 = 0.0 # initial front flow
initVar = [v0, s0, Qlr0, Qlf0, Qf0]

# friction factor calculation
def ff(Q, A, L):
    D = 2.0 * (A / np.pi)**0.5 # hydraulic diameter
    Re = (rho * (Q / A) * L) / mu # reynolds number
    ff = 0.25 / (np.log10((epsilon / D)+(5.74 / Re**0.9)))**2.0 # friction factor
    return ff

# main calculation
def calc_main(var, t, val):
    
    par, bound = val
    Db, Lb, mb, Ai, Ki, Li, Dt, Lt, epsilon, rho, mu = par
    pp, patm = bound
    vb, sb, Qlr, Qlf, Qf = var
    
    Ab = np.pi * (Db / 2.0)**2.0 # body area
    At = np.pi * (Dt / 2.0)**2.0 # tube area
    Alr = At - Ab # rear leak area
    Alf = At - Ab # front leak area
    Af = At # front flow area
    
    Dlr = 2.0 * (Alr / np.pi)**0.5 # hydraulic diameter of rear leak
    Dlf = 2.0 * (Alf / np.pi)**0.5 # hydraulic diameter of front leak
    Df = Dt # hydraulic diameter of front flow
    
    Lr = sb # rear position
    Lf = Lr + Lb # front position
    Llr = max(0.0, Li - Lr) # length of rear leak
    Llf = min(Lt - Li, Lf - Li) # length of front leak
    Lf = max(0.0, Lt - Lf) # length of front flow
    
    Qi = (1 / Ki)**(0.5) * Ab * ((2 * (pp - pi)) / rho)**(0.5) # inlet flow rate
    dQlrdt = Alr / (rho * Llr) * (pi - pr) * ((ff(Qlr, Alr, Llr) * Qlr**2.0) / (2 * Dlr * Alr)) # rear leak derivative
    dQlfdt = Alf / (rho * Llf) * (pi - pf) * ((ff(Qlf, Alf, Llf) * Qlf**2.0) / (2 * Dlf * Alf)) # front leak derivative
    dQfdt = Af / (rho * Lf) * (pr - patm) * ((ff(Qf, Af, Lf) * Qf**2.0) / (2 * Df * Af)) # front flow derivative
    Qf = Qi # flow continuity
    Qf = Qlr + Qlf # flow continuity
    Qlr = Ab * vb # flow continuity
    
    Fb = (pr - pf) * (pi * (Db / 2.0) ** 2.0) # force on body
    ab = Fb / mb # body acceleration
    dvbdt = ab # body velocity derivative
    dsbdt = vb # body displacement derivative
    
    return dvbdt, dsbdt, dQlrdt, dQlfdt, dQfdt

# solve
res = odeint(calc_main, initVar, t, args=([par, bound],))

If I take out the flow equations and have pr and pf equal to pp and patm respectively, then the motion equations work fine (results below), but with the flow equations as above it complains that the pressures are undefined (understandably, I have nothing saying pf=[] or pr=[])

What I've attempted is to use a combination of fsolve and odeint to solve the simultaneous flow equations, but also solve the ODEs, to little success

from scipy.optimize import fsolve

# main calculation
def calc_main(var, t, val):
    
    par, bound = val
    Db, Lb, mb, Ai, Ki, Li, Dt, Lt, epsilon, rho, mu = par
    pp, patm = bound
    vb, sb, Qlr, Qlf, Qf = var
    
    Ab = np.pi * (Db / 2.0)**2.0 # body area
    At = np.pi * (Dt / 2.0)**2.0 # tube area
    Alr = At - Ab # rear leak area
    Alf = At - Ab # front leak area
    Af = At # front flow area
    
    Dlr = 2.0 * (Alr / np.pi)**0.5 # hydraulic diameter of rear leak
    Dlf = 2.0 * (Alf / np.pi)**0.5 # hydraulic diameter of front leak
    Df = Dt # hydraulic diameter of front flow
    
    Lr = sb # rear position
    Lf = Lr + Lb # front position
    Llr = max(0.0, Li - Lr) # length of rear leak
    Llf = min(Lt - Li, Lf - Li) # length of front leak
    Lf = max(0.0, Lt - Lf) # length of front flow
    
    def calc_fluid (var):
        
        Qi, dQlrdt, dQlfdt, dQfdt, pi, pr, pf = var
        
        eq1 = (1 / Ki)**(0.5) * Ab * ((2 * (pp - pi)) / rho)**(0.5) - Qi # inlet flow rate
        eq2 = Alr / (rho * Llr) * (pi - pr) * ((ff(Qlr, Alr, Llr) * Qlr**2.0) / (2 * Dlr * Alr)) - dQlrdt # rear leak derivative
        eq3 = Alf / (rho * Llf) * (pi - pf) * ((ff(Qlf, Alf, Llf) * Qlf**2.0) / (2 * Dlf * Alf)) - dQlfdt # front leak derivative
        eq4 = Af / (rho * Lf) * (pr - patm) * ((ff(Qf, Af, Lf) * Qf**2.0) / (2 * Df * Af)) - dQfdt # front flow derivative
        eq5 = Qi - Qf # flow continuity
        eq6 = Qlr + Qlf - Qf # flow continuity
        eq7 = Ab * vb - Qlr # flow continuity
        
        return eq1, eq2, eq3, eq4, eq5, eq6, eq7
    
    Qi, dQlrdt, dQlfdt, dQfdt, pi, pr, pf = fsolve(calc_fluid, [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1])
    
    Fb = (pr - pf) * (pi * (Db / 2.0) ** 2.0) # force on body
    ab = Fb / mb # body acceleration
    dvbdt = ab # body velocity derivative
    dsbdt = vb # body displacement derivative
    
    return dvbdt, dsbdt, dQlrdt, dQlfdt, dQfdt

I get the error "RuntimeWarning: invalid value encountered in double_scalars" for a few different lines while it runs, but I don't get told specifically which time step that happens at

I'm not convinced that the results I get are right, as I seem to now not move the body at all

Can anyone see what I've done wrong, or know of a better way to solve the simultaneous equations alongside the derivatives?

I've been shooting in the dark a little having not done anything quite like this before so it's possible I've gone completely off the rails with how these functions are meant to be used, any explanation is much appreciated!

Answer 1

Your example does has some undefined variables, so I didn't try it here.

What you have to do is to write the constraints dependent variables in a separate function and then use solve_bvp to integrate the differential equation with the boundary conditions.

The bc function is evaluated with the ya and yb being the value of your dependent variables at the beginning and the end of your interval. It will return an array of values that are 0 for the expected solution, and non-zero for other solutions.

Note

Also, you may like to know that

has a closed form solution

I am what varies in your equations, but the coefficients of your differential equations for Q then you can use the closed form, and solving the boundary conditions for the free parameters.