Kernel: Python 3 (system-wide)
In [1]:
import sympy as sp import numpy as np import matplotlib.pyplot as plt from scipy import integrate from scipy.integrate import solve_ivp,trapezoid import math import mpmath as mp import os from scipy.interpolate import interp1d
Physical parameters
In [2]:
## All air-mucus parameters correspond to Table 2 in ## Dietze and Ruyer-Quil (J. Fluid Mech., 762, 2015, 68). ep = 1 ## aspect ratio nG = 256 ## spatial grid points d0 = 0.86 ## average interface position L = 7.73 ## Domain length R = 3e-4 ## Tube radius gamma_am = 0.025 ## air-mucus surface tension mu_a = 1.8e-5 ## air viscosity mu_m = 0.01 ## mucus viscosity rho_a = 1.2 ## air density rho_m = 1098.3 ## mucus density Lamnda_R = L
In [3]:
## calculating non-dimensional parameters alpha_star = 2*np.pi*d0/(Lamnda_R) lambda_star = 1/alpha_star tau = (6*mu_m*R*d0/gamma_am)/(alpha_star**4*(1/alpha_star**2-1)*(1/d0-1)**2*(1/d0**2-1)) ## tau correspond to eq. 3.1 in Dietze et al. 2015 u_char = gamma_am/mu_m PI_rho = rho_a/rho_m PI_mu = mu_a/mu_m Re_m = rho_m*R**2/(mu_m*tau) Re_a = Re_m*PI_rho/PI_mu Ca = 1
Setting up the numerical solver for the WRIBL equations
In [4]:
## grid n = nG h = (L-0) / n lesz = np.arange(0,L,h)
In [5]:
## Derivative matrices ## D1: central difference D1c =np.zeros((n,n)) ## boundary nodes: periodic D1c[0,0]= 0 D1c[0,1]=1 D1c[0,-1]=-1 D1c[-1,0]= 1 D1c[-1,-2]=-1 D1c[-1,-1]=0 ## internal nodes for i in range(1, n-1): D1c[i,i-1] = -1 D1c[i,i] = 0 D1c[i,i+1]=1 D1c = 1/(2*h)*D1c ## D2: central difference D2c =np.zeros((n, n)) ## boundary nodes: periodic D2c[0, 0] = -2 D2c[0, n-1] = 1 D2c[0, 1] = 1 D2c[n-1, n-1] = -2 D2c[n-1, n-2] = 1 D2c[n-1, 0] = 1 ## internal nodes for i in range(1, n-1): D2c[i, i-1] = 1 D2c[i, i] = -2 D2c[i, i+1] = 1 D2c = 1/(h**2)*D2c ## D1: forward difference D1f =np.zeros((n,n)) ## boundary nodes: periodic D1f[0,0]= -1 D1f[0,1]=1 D1f[-1,0]= 1 D1f[-1,-1]=-1 ## internal nodes for i in range(1, n-1): D1f[i,i] = -1 D1f[i,i+1]=1 D1f = 1/(h)*D1f ## av matrix: M1av = np.zeros((n,n)) ## boundary nodes: periodic M1av[0,0]= 1 M1av[0,1]=1 M1av[-1,0]= 1 M1av[-1,-1]=1 ## internal nodes for i in range(1, n-1): M1av[i,i] = 1 M1av[i,i+1]=1 M1av = 1/(2)*M1av
In [6]:
## grid for interpolation function d_i = np.linspace(0.1,0.999,10000)
All coefficients with tilde in the paper are presented here with a subscript 'p', e.g., Saa​~​ as Saap
Coefficients associated with cilia (subscript 'c') in the paper are presented denoted with 'w' here, e.g., Facc​ as Faww
In [7]:
## defining symbolic variables and functions r,t,z = sp.symbols('r t z ',real = True) d = sp.Function('d',real=True)(z,t) Q_a = sp.Function('Q_a',real=True)(z,t) Q_m = sp.Function('Q_m',real=True)(z,t) Q_t = sp.Function('Q_t',real=True)(t) U_w= sp.Function('U_w',real=True)(z,t) kappa= sp.Function('kappa',real=True)(z,t) dpdz_d = sp.Function('dpdz_d',real=True)(z,t) Saa_s = sp.Function('Saa_s',real=True)(d) Sam_s = sp.Function('Sam_s',real=True)(d) Saap_s = sp.Function('Saap_s',real=True)(d) Samp_s = sp.Function('Samp_s',real=True)(d) Sma_s = sp.Function('Sma_s',real=True)(d) Smm_s = sp.Function('Smm_s',real=True)(d) Smap_s = sp.Function('Smap_s',real=True)(d) Smmp_s = sp.Function('Smmp_s',real=True)(d)
In [8]:
## Reading WRIBL coefficients from the text files Aa1= open("./Two_way_coefficients/Saa.txt","r") Aa1s = Aa1.read() def Saa(d): return eval(Aa1s) Saa_in = Saa(d_i) interp_func_Saa = interp1d(d_i,Saa_in,kind='cubic') def Saa_i(d): return interp_func_Saa(d) Aap1= open("./Two_way_coefficients/Saap.txt","r") Aap1s = Aap1.read() def Saa_p(d): return eval(Aap1s) Saa_p_in = Saa_p(d_i) interp_func_Saa_p = interp1d(d_i,Saa_p_in,kind='cubic') def Saa_pi(d): return interp_func_Saa_p(d) Am1= open("./Two_way_coefficients/Sma.txt","r") Am1s = Am1.read() def Sma(d): return eval(Am1s) Sma_in = Sma(d_i) interp_func_Sma = interp1d(d_i,Sma_in,kind='cubic') def Sma_i(d): return interp_func_Sma(d) Amp1= open("./Two_way_coefficients/Smap.txt","r") Amp1s = Amp1.read() def Sma_p(d): return eval(Amp1s) Sma_p_in = Sma_p(d_i) interp_func_Sma_p = interp1d(d_i,Sma_p_in,kind='cubic') def Sma_pi(d): return interp_func_Sma_p(d) Aa2= open("./Two_way_coefficients/Sam.txt","r") Aa2s = Aa2.read() def Sam(d): return eval(Aa2s) Sam_in = Sam(d_i) interp_func_Sam = interp1d(d_i,Sam_in,kind='cubic') def Sam_i(d): return interp_func_Sam(d) Aap2= open("./Two_way_coefficients/Samp.txt","r") Aap2s = Aap2.read() def Sam_p(d): return eval(Aap2s) Sam_p_in = Sam_p(d_i) interp_func_Sam_p = interp1d(d_i,Sam_p_in,kind='cubic') def Sam_pi(d): return interp_func_Sam_p(d) Am2= open("./Two_way_coefficients/Smm.txt","r") Am2s = Am2.read() def Smm(d): return eval(Am2s) Smm_in = Smm(d_i) interp_func_Smm = interp1d(d_i,Smm_in,kind='cubic') def Smm_i(d): return interp_func_Smm(d) Amp2= open("./Two_way_coefficients/Smmp.txt","r") Amp2s = Amp2.read() def Smm_p(d): return eval(Amp2s) Smm_p_in = Smm_p(d_i) interp_func_Smm_p = interp1d(d_i,Smm_p_in,kind='cubic') def Smm_pi(d): return interp_func_Smm_p(d) Aa4= open("./Two_way_coefficients/Faaa.txt","r") Aa4s = Aa4.read() def Faaa(d): return eval(Aa4s) Faaa_in = Faaa(d_i) interp_func_Faaa = interp1d(d_i,Faaa_in,kind='cubic') def Faaa_i(d): return interp_func_Faaa(d) Aap4= open("./Two_way_coefficients/Faaap.txt","r") Aap4s = Aap4.read() def Faaa_p(d): return eval(Aap4s) Faaap_in = Faaa_p(d_i) interp_func_Faaap = interp1d(d_i,Faaap_in,kind='cubic') def Faaa_pi(d): return interp_func_Faaap(d) Am4= open("./Two_way_coefficients/Fmaa.txt","r") Am4s = Am4.read() def Fmaa(d): return eval(Am4s) Fmaa_in = Fmaa(d_i) interp_func_Fmaa = interp1d(d_i,Fmaa_in,kind='cubic') def Fmaa_i(d): return interp_func_Fmaa(d) Amp4= open("./Two_way_coefficients/Fmaap.txt","r") Amp4s = Amp4.read() def Fmaa_p(d): return eval(Amp4s) Fmaap_in = Fmaa_p(d_i) interp_func_Fmaap = interp1d(d_i,Fmaap_in,kind='cubic') def Fmaa_pi(d): return interp_func_Fmaap(d) Aa5= open("./Two_way_coefficients/Faam.txt","r") Aa5s = Aa5.read() def Faam(d): return eval(Aa5s) Faam_in = Faam(d_i) interp_func_Faam = interp1d(d_i,Faam_in,kind='cubic') def Faam_i(d): return interp_func_Faam(d) Aap5= open("./Two_way_coefficients/Faamp.txt","r") Aap5s = Aap5.read() def Faam_p(d): return eval(Aap5s) Faamp_in = Faam_p(d_i) interp_func_Faamp = interp1d(d_i,Faamp_in,kind='cubic') def Faam_pi(d): return interp_func_Faamp(d) Am5= open("./Two_way_coefficients/Fmam.txt","r") Am5s = Am5.read() def Fmam(d): return eval(Am5s) Fmam_in = Fmam(d_i) interp_func_Fmam = interp1d(d_i,Fmam_in,kind='cubic') def Fmam_i(d): return interp_func_Fmam(d) Amp5= open("./Two_way_coefficients/Fmamp.txt","r") Amp5s = Amp5.read() def Fmam_p(d): return eval(Amp5s) Fmamp_in = Fmam_p(d_i) interp_func_Fmamp = interp1d(d_i,Fmamp_in,kind='cubic') def Fmam_pi(d): return interp_func_Fmamp(d) Aa6= open("./Two_way_coefficients/Fama.txt","r") Aa6s = Aa6.read() def Fama(d): return eval(Aa6s) Fama_in = Fama(d_i) interp_func_Fama = interp1d(d_i,Fama_in,kind='cubic') def Fama_i(d): return interp_func_Fama(d) Aap6= open("./Two_way_coefficients/Famap.txt","r") Aap6s = Aap6.read() def Fama_p(d): return eval(Aap6s) Famap_in = Fama_p(d_i) interp_func_Famap = interp1d(d_i,Famap_in,kind='cubic') def Fama_pi(d): return interp_func_Famap(d) Am6= open("./Two_way_coefficients/Fmma.txt","r") Am6s = Am6.read() def Fmma(d): return eval(Am6s) Fmma_in = Fmma(d_i) interp_func_Fmma = interp1d(d_i,Fmma_in,kind='cubic') def Fmma_i(d): return interp_func_Fmma(d) Amp6= open("./Two_way_coefficients/Fmmap.txt","r") Amp6s = Amp6.read() def Fmma_p(d): return eval(Amp6s) Fmmap_in = Fmma_p(d_i) interp_func_Fmmap = interp1d(d_i,Fmmap_in,kind='cubic') def Fmma_pi(d): return interp_func_Fmmap(d) Aa7= open("./Two_way_coefficients/Famm.txt","r") Aa7s = Aa7.read() def Famm(d): return eval(Aa7s) Famm_in = Famm(d_i) interp_func_Famm = interp1d(d_i,Famm_in,kind='cubic') def Famm_i(d): return interp_func_Famm(d) Aap7= open("./Two_way_coefficients/Fammp.txt","r") Aap7s = Aap7.read() def Famm_p(d): return eval(Aap7s) Fammp_in = Famm_p(d_i) interp_func_Fammp = interp1d(d_i,Fammp_in,kind='cubic') def Famm_pi(d): return interp_func_Fammp(d) Am7= open("./Two_way_coefficients/Fmmm.txt","r") Am7s = Am7.read() def Fmmm(d): return eval(Am7s) Fmmm_in = Fmmm(d_i) interp_func_Fmmm = interp1d(d_i,Fmmm_in,kind='cubic') def Fmmm_i(d): return interp_func_Fmmm(d) Amp7= open("./Two_way_coefficients/Fmmmp.txt","r") Amp7s = Amp7.read() def Fmmm_p(d): return eval(Amp7s) Fmmmp_in = Fmmm_p(d_i) interp_func_Fmmmp = interp1d(d_i,Fmmmp_in,kind='cubic') def Fmmm_pi(d): return interp_func_Fmmmp(d) Aa13= open("./Two_way_coefficients/Gaaa.txt","r") Aa13s = Aa13.read() def Gaaa(d): return eval(Aa13s) Gaaa_in = Gaaa(d_i) interp_func_Gaaa = interp1d(d_i,Gaaa_in,kind='cubic') def Gaaa_i(d): return interp_func_Gaaa(d) Aap13= open("./Two_way_coefficients/Gaaap.txt","r") Aap13s = Aap13.read() def Gaaa_p(d): return eval(Aap13s) Gaaap_in = Gaaa_p(d_i) interp_func_Gaaap = interp1d(d_i,Gaaap_in,kind='cubic') def Gaaa_pi(d): return interp_func_Gaaap(d) Am13= open("./Two_way_coefficients/Gmaa.txt","r") Am13s = Am13.read() def Gmaa(d): return eval(Am13s) Gmaa_in = Gmaa(d_i) interp_func_Gmaa = interp1d(d_i,Gmaa_in,kind='cubic') def Gmaa_i(d): return interp_func_Gmaa(d) Amp13= open("./Two_way_coefficients/Gmaap.txt","r") Amp13s = Amp13.read() def Gmaa_p(d): return eval(Amp13s) Gmaap_in = Gmaa_p(d_i) interp_func_Gmaap = interp1d(d_i,Gmaap_in,kind='cubic') def Gmaa_pi(d): return interp_func_Gmaap(d) Aa14= open("./Two_way_coefficients/Gamm.txt","r") Aa14s = Aa14.read() def Gamm(d): return eval(Aa14s) Gamm_in = Gamm(d_i) interp_func_Gamm = interp1d(d_i,Gamm_in,kind='cubic') def Gamm_i(d): return interp_func_Gamm(d) Aap14= open("./Two_way_coefficients/Gammp.txt","r") Aap14s = Aap14.read() def Gamm_p(d): return eval(Aap14s) Gammp_in = Gamm_p(d_i) interp_func_Gammp = interp1d(d_i,Gammp_in,kind='cubic') def Gamm_pi(d): return interp_func_Gammp(d) Am14= open("./Two_way_coefficients/Gmmm.txt","r") Am14s = Am14.read() def Gmmm(d): return eval(Am14s) Gmmm_in = Gmmm(d_i) interp_func_Gmmm = interp1d(d_i,Gmmm_in,kind='cubic') def Gmmm_i(d): return interp_func_Gmmm(d) Amp14= open("./Two_way_coefficients/Gmmmp.txt","r") Amp14s = Amp14.read() def Gmmm_p(d): return eval(Amp14s) Gmmmp_in = Gmmm_p(d_i) interp_func_Gmmmp = interp1d(d_i,Gmmmp_in,kind='cubic') def Gmmm_pi(d): return interp_func_Gmmmp(d) Aa16= open("./Two_way_coefficients/Gaam.txt","r") Aa16s = Aa16.read() def Gaam(d): return eval(Aa16s) Gaam_in = Gaam(d_i) interp_func_Gaam = interp1d(d_i,Gaam_in,kind='cubic') def Gaam_i(d): return interp_func_Gaam(d) Aap16= open("./Two_way_coefficients/Gaamp.txt","r") Aap16s = Aap16.read() def Gaam_p(d): return eval(Aap16s) Gaamp_in = Gaam_p(d_i) interp_func_Gaamp = interp1d(d_i,Gaamp_in,kind='cubic') def Gaam_pi(d): return interp_func_Gaamp(d) Am16= open("./Two_way_coefficients/Gmam.txt","r") Am16s = Am16.read() def Gmam(d): return eval(Am16s) Gmam_in = Gmam(d_i) interp_func_Gmam = interp1d(d_i,Gmam_in,kind='cubic') def Gmam_i(d): return interp_func_Gmam(d) Amp16= open("./Two_way_coefficients/Gmamp.txt","r") Amp16s = Amp16.read() def Gmam_p(d): return eval(Amp16s) Gmamp_in = Gmam_p(d_i) interp_func_Gmamp = interp1d(d_i,Gmamp_in,kind='cubic') def Gmam_pi(d): return interp_func_Gmamp(d) B1= open("./Two_way_coefficients/Ja.txt","r") B1s = B1.read() def Ja(d): return eval(B1s) Ja_in = Ja(d_i) interp_func_Ja = interp1d(d_i,Ja_in,kind='cubic') def Ja_i(d): return interp_func_Ja(d) Bp1= open("./Two_way_coefficients/Jap.txt","r") Bp1s = Bp1.read() def Ja_p(d): return eval(Bp1s) Jap_in = Ja_p(d_i) interp_func_Jap = interp1d(d_i,Jap_in,kind='cubic') def Ja_pi(d): return interp_func_Jap(d) B2= open("./Two_way_coefficients/Jm.txt","r") B2s = B2.read() def Jm(d): return eval(B2s) Jm_in = Jm(d_i) interp_func_Jm = interp1d(d_i,Jm_in,kind='cubic') def Jm_i(d): return interp_func_Jm(d) Bp2= open("./Two_way_coefficients/Jmp.txt","r") Bp2s = Bp2.read() def Jm_p(d): return eval(Bp2s) Jmp_in = Jm_p(d_i) interp_func_Jmp = interp1d(d_i,Jmp_in,kind='cubic') def Jm_pi(d): return interp_func_Jmp(d) B4= open("./Two_way_coefficients/Ka.txt","r") B4s = B4.read() def Ka(d): return eval(B4s) Ka_in = Ka(d_i) interp_func_Ka = interp1d(d_i,Ka_in,kind='cubic') def Ka_i(d): return interp_func_Ka(d) Bp4= open("./Two_way_coefficients/Kap.txt","r") Bp4s = Bp4.read() def Ka_p(d): return eval(Bp4s) Kap_in = Ka_p(d_i) interp_func_Kap = interp1d(d_i,Kap_in,kind='cubic') def Ka_pi(d): return interp_func_Kap(d) B5= open("./Two_way_coefficients/Km.txt","r") B5s = B5.read() def Km(d): return eval(B5s) Km_in = Km(d_i) interp_func_Km = interp1d(d_i,Km_in,kind='cubic') def Km_i(d): return interp_func_Km(d) Bp5= open("./Two_way_coefficients/Kmp.txt","r") Bp5s = Bp5.read() def Km_p(d): return eval(Bp5s) Kmp_in = Km_p(d_i) interp_func_Kmp = interp1d(d_i,Kmp_in,kind='cubic') def Km_pi(d): return interp_func_Kmp(d) B7= open("./Two_way_coefficients/La.txt","r") B7s = B7.read() def La(d): return eval(B7s) La_in = La(d_i) interp_func_La = interp1d(d_i,La_in,kind='cubic') def La_i(d): return interp_func_La(d) Bp7= open("./Two_way_coefficients/Lap.txt","r") Bp7s = Bp7.read() def La_p(d): return eval(Bp7s) Lap_in = La_p(d_i) interp_func_Lap = interp1d(d_i,Lap_in,kind='cubic') def La_pi(d): return interp_func_Lap(d) B8= open("./Two_way_coefficients/Lm.txt","r") B8s = B8.read() def Lm(d): return eval(B8s) Lm_in = Lm(d_i) interp_func_Lm = interp1d(d_i,Lm_in,kind='cubic') def Lm_i(d): return interp_func_Lm(d) Bp8= open("./Two_way_coefficients/Lmp.txt","r") Bp8s = Bp8.read() def Lm_p(d): return eval(Bp8s) Lmp_in = Lm_p(d_i) interp_func_Lmp = interp1d(d_i,Lmp_in,kind='cubic') def Lm_pi(d): return interp_func_Lmp(d) B10= open("./Two_way_coefficients/Ma.txt","r") B10s = B10.read() def Ma(d): return eval(B10s) Ma_in = Ma(d_i) interp_func_Ma = interp1d(d_i,Ma_in,kind='cubic') def Ma_i(d): return interp_func_Ma(d) Bp10= open("./Two_way_coefficients/Map.txt","r") Bp10s = Bp10.read() def Ma_p(d): return eval(Bp10s) Map_in = Ma_p(d_i) interp_func_Map = interp1d(d_i,Map_in,kind='cubic') def Ma_pi(d): return interp_func_Map(d) B11= open("./Two_way_coefficients/Mm.txt","r") B11s = B11.read() def Mm(d): return eval(B11s) Mm_in = Mm(d_i) interp_func_Mm = interp1d(d_i,Mm_in,kind='cubic') def Mm_i(d): return interp_func_Mm(d) Bp11= open("./Two_way_coefficients/Mmp.txt","r") Bp11s = Bp11.read() def Mm_p(d): return eval(Bp11s) Mmp_in = Mm_p(d_i) interp_func_Mmp = interp1d(d_i,Mmp_in,kind='cubic') def Mm_pi(d): return interp_func_Mmp(d) # csta,cstm Da1= open("./Two_way_coefficients/Csta.txt","r") Da1s = Da1.read() def Csta(d): return eval(Da1s) Csta_s = sp.Function('Csta_s',real=True)(d) Csta_in = Csta(d_i) interp_func_Csta = interp1d(d_i,Csta_in,kind='cubic') def Csta_i(d): return interp_func_Csta(d) Dap1= open("./Two_way_coefficients/Cstap.txt","r") Dap1s = Dap1.read() def Csta_p(d): return eval(Dap1s) Cstap_s = sp.Function('Cstap_s',real=True)(d) Cstap_in = Csta_p(d_i) interp_func_Cstap = interp1d(d_i,Cstap_in,kind='cubic') def Csta_pi(d): return interp_func_Cstap(d) Ea1= open("./Two_way_coefficients/Iw_af.txt","r") Ea1s=Ea1.read() def Iw_af(d): return eval(Ea1s) Iw_af_s = sp.Function('Iw_af_s',real=True)(d) Iw_af_in = Iw_af(d_i) interp_func_Iw_af = interp1d(d_i,Iw_af_in,kind='cubic') def Iw_afi(d): return interp_func_Iw_af(d) Eap1= open("./Two_way_coefficients/Iw_apf.txt","r") Eap1s=Eap1.read() def Iw_apf(d): return eval(Eap1s) Iw_apf_in = Iw_apf(d_i) interp_func_Iw_apf = interp1d(d_i,Iw_apf_in,kind='cubic') def Iw_apfi(d): return interp_func_Iw_apf(d) Iw_apf_s = sp.Function('Iw_apf_s',real=True)(d)
In [9]:
h0 = 0.0000 ## mean curvature def K(d): return 1/d-ep**2*(D1c@(d))**2/(2*d)-ep**2*(D2c@(d))+ep*h0**6/((1-d)**9) def D1_k(d): return D1c@(K(d))
All the WRIBL terms except ∂t​
In [10]:
def Space_derivative(t,d,Qa,Qt): Qm = Qt-Qa D1_d = D1f@(d) D1_Qa = D1c@(Qa) D1_Qm = D1c@(Qm) D2_Qa = D2c@(Qa) D2_Qm = D2c@(Qm) D2_d = D1f@(D1c@(d)) return (-Ca*D1_k(d)*Iw_afi(d)+Qm+Csta_i(d)*PI_mu*Qa+\ (M1av@Ja_i(d))*Qa*(D1_d)**2+(M1av@Jm_i(d))*Qm*(D1_d)**2+\ (M1av@Ka_i(d))*D1_Qa*D1_d+(M1av@Km_i(d))*D1_Qm*D1_d+\ (M1av@La_i(d))*Qa*D2_d+(M1av@Lm_i(d))*Qm*D2_d+\ (M1av@Ma_i(d))*D2_Qa+(M1av@Mm_i(d))*D2_Qm-(PI_mu*Re_a*(M1av@Faaa_i(d)*Qa*D1_Qa+M1av@Faam_i(d)*Qa*D1_Qm)+\ M1av@Fama_i(d)*Qm*D1_Qa+M1av@Famm_i(d)*Qm*D1_Qm+\ M1av@Gaaa_i(d)*Qa**2*D1_d+M1av@Gaam_i(d)*Qa*Qm*D1_d+\ M1av@Gamm_i(d)*Qm**2*D1_d+\ Re_m*(M1av@Fmaa_i(d)*Qa*D1_Qa+M1av@Fmam_i(d)*Qa*D1_Qm+\ M1av@Fmma_i(d)*Qm*D1_Qa+M1av@Fmmm_i(d)*Qm*D1_Qm+\ M1av@Gmaa_i(d)*Qa**2*D1_d+M1av@Gmmm_i(d)*Qm**2*D1_d+\ M1av@Gmam_i(d)*Qa*Qm*D1_d))) def Space_derivative_p(t,d,Qa,Qt): Qm = Qt-Qa D1_d = D1f@(d) D1_Qa = D1c@(Qa) D1_Qm = D1c@(Qm) D2_Qa = D2c@(Qa) D2_Qm = D2c@(Qm) D2_d = D1f@(D1c@(d)) return (Ca*D1_k(d)*Iw_apfi(d)+Qm+Csta_pi(d)*PI_mu*Qa+\ M1av@Ja_pi(d)*Qa*(D1_d)**2+M1av@Jm_pi(d)*Qm*(D1_d)**2+\ M1av@Ka_pi(d)*D1_Qa*D1_d+M1av@Km_pi(d)*D1_Qm*D1_d+\ M1av@La_pi(d)*Qa*D2_d+M1av@Lm_pi(d)*Qm*D2_d+\ M1av@Ma_pi(d)*D2_Qa+M1av@Mm_pi(d)*D2_Qm-(PI_mu*Re_a*(M1av@Faaa_pi(d)*Qa*D1_Qa+M1av@Faam_pi(d)*Qa*D1_Qm+\ M1av@Fama_pi(d)*Qm*D1_Qa+M1av@Famm_pi(d)*Qm*D1_Qm+\ M1av@Gaaa_pi(d)*Qa**2*D1_d+M1av@Gamm_pi(d)*Qm**2*D1_d+\ M1av@Gaam_pi(d)*Qa*Qm*D1_d)+\ Re_m*(M1av@Fmaa_pi(d)*Qa*D1_Qa+M1av@Fmam_pi(d)*Qa*D1_Qm+\ M1av@Fmma_pi(d)*Qm*D1_Qa+M1av@Fmmm_pi(d)*Qm*D1_Qm+\ M1av@Gmaa_pi(d)*Qa**2*D1_d+M1av@Gmmm_pi(d)*Qm**2*D1_d+\ M1av@Gmam_pi(d)*Qa*Qm*D1_d)))/(Iw_apfi(d))
In [0]:
Algebric manupulations to make an explicit equation of ∂t​
In [11]:
Space_derivative_s = sp.Function('Space_derivative_s',real=True)(d) Space_derivative_p_s = sp.Function('Space_derivative_p_s',real=True)(d) dpadz_d = sp.Function('dpadz_d',real=True)(d) Full_equ= PI_mu*Re_a*(Saa_s*sp.diff(Q_a,t)+Sam_s*sp.diff(Q_m,t))+Re_m*(Sma_s*sp.diff(Q_a,t)+Smm_s*sp.diff(Q_m,t))-Space_derivative_s Full_equ_p= PI_mu*Re_a*(Saap_s*sp.diff(Q_a,t)+Samp_s*sp.diff(Q_m,t))+Re_m*(Smap_s*sp.diff(Q_a,t)+Smmp_s*sp.diff(Q_m,t))-Space_derivative_p_s Full_equ = (Full_equ.subs(Q_m,Q_t-Q_a)).doit() Full_equ_p = (Full_equ_p.subs(Q_m,Q_t-Q_a)).doit() dQadt_sol1 = sp.solve(Full_equ,sp.diff(Q_a,t)) dQadt = (sp.expand(dQadt_sol1[0])) Full_equ_p = sp.expand(Full_equ_p.subs(sp.diff(Q_a,t),dQadt)) Coeff_dqtdt=Full_equ_p.coeff(sp.diff(Q_t,t)) Res_Full_equ_p = sp.expand(sp.cancel(Full_equ_p-Coeff_dqtdt*sp.diff(Q_t,t))) coeff_s = ['Saa_s(d(z, t))','Sam_s(d(z, t))','Sma_s(d(z, t))','Smm_s(d(z, t))','Saap_s(d(z, t))','Samp_s(d(z, t))','Smap_s(d(z, t))','Smmp_s(d(z, t))','Space_derivative_s(d(z, t))','Space_derivative_p_s(d(z, t))'] coeff_f = ['Saa_i(d)','Sam_i(d)','Sma_i(d)','Smm_i(d)','Saa_pi(d)/Iw_apf(d)','Sam_pi(d)/Iw_apf(d)','Sma_pi(d)/Iw_apf(d)','Smm_pi(d)/Iw_apf(d)','Space_derivative(t,d,Qa,Qt)','Space_derivative_p(t,d,Qa,Qt)'] Coeff_dqtdt1 = str(Coeff_dqtdt) for i in range (0,len(coeff_s)): Coeff_dqtdt1 = Coeff_dqtdt1.replace(coeff_s[i],coeff_f[i]) Res_Full_equ_p1 = str(Res_Full_equ_p ) for i in range (0,len(coeff_s)): Res_Full_equ_p1 = Res_Full_equ_p1.replace(coeff_s[i],coeff_f[i])
Defining functions for the ODE solver
In [12]:
def C_dQtdt_press(d): return eval(Coeff_dqtdt1) def R_press(t,d,Qa,Qt): return eval(Res_Full_equ_p1) def odesys_1(t,d,Qa,Qt): return -trapezoid(R_press(t,d,Qa,Qt),lesz)/trapezoid(C_dQtdt_press(d),lesz) dqadt1 = str(dQadt) for i in range (0,len(coeff_s)): dqadt1 = dqadt1.replace(coeff_s[i],coeff_f[i]) dqadt2 = dqadt1.replace('Derivative(Q_t(t), t)','odesys1') def odesys_2(t,d,Qa,Qt,odesys1): return eval(dqadt2) def odesys_3(t,d,Qa): return -M1av@(1/(d))*(D1c@Qa) def odesys(t,arr): d = (arr[0:nG]) Qa = (arr[nG:2*nG]) Qt = arr[2*nG] odesys1=odesys_1(t,d,Qa,Qt) Z01 = np.concatenate([(odesys_3(t,d,Qa),odesys_2(t,d,Qa,Qt,odesys1))]) darrdt=np.append(Z01,odesys1) return darrdt
Performing simulations
In [13]:
## Initial condition d01=d0+0.001*np.cos(2*np.pi/L*lesz) Qa0 = np.zeros(nG) Qt0 = 0 Z01 = np.concatenate((d01,Qa0)) I0 = np.zeros(2*nG+1) I0[2*nG] = Qt0 I0[0:2*nG] = Z01
In [14]:
#Specifying event-handlers to stop the simulation when # the interface approaches the wall (drainout) def event1(t,arr): return (arr[0:nG]).max()-0.997 event1.terminal = True # or the interface approaches the centreline (pinchoff) def event2(t,arr): print((arr[0:nG]).min()) return (arr[0:nG]).min()-0.2 event2.terminal = True
In [15]:
# total run time T = 500000
In [22]:
## time stepping with solve_ivp sol = solve_ivp(odesys,[0.,T],I0,method='LSODA',events = [event1,event2],max_step = 1000)
Postprocessing and validation with Dietze and Ruyer-Quil (J. Fluid Mech., 762, 2015, 68)
As an example, we reproduce figure 3d of Dietze and Ruyer-Quil (2015) below. This figure plots the evolution of the minimum of the annular interface, for a case that ends in pinch-off. The results are in agreement.
In [0]:
In [17]:
# obtaining the minimum of the interface profile d(z,t) def min_data(data): min_d_data = np.zeros(data.shape[1]) for i in range(0,min_d_data.shape[0]): min_d_data[i] = data[:nG,i].min() return min_d_data def max_data(data): max_d_data = np.zeros(data.shape[1]) for i in range(0,max_d_data.shape[0]): max_d_data[i] = data[:nG,i].max() return max_d_data
In [18]:
## data extracted from Figure 9b with liquid volume 2.01*pi*R^3 of Dietze and Ruyer-Quil (2015) using OriginLab (www.originlab.com) t_min_data = np.array([0.008417508417508417, 0.08417508417508418, 0.16835016835016836, 0.2777777777777778, 0.3872053872053872, 0.5134680134680134, 0.6060606060606061, 0.6986531986531986, 0.7912457912457912, 0.9090909090909091, 1.0101010101010102, 1.1111111111111112, 1.1952861952861953, 1.2794612794612794, 1.3804713804713804, 1.5151515151515151, 1.6161616161616161, 1.7424242424242424, 1.8434343434343434, 1.9696969696969697, 2.079124579124579, 2.2306397306397305, 2.4074074074074074, 2.5673400673400675, 2.7861952861952863, 3.005050505050505, 3.274410774410774, 3.5016835016835017, 3.7794612794612794, 4.006734006734007, 4.2592592592592595, 4.486531986531986, 4.713804713804714, 4.856902356902356, 4.865319865319865, 4.873737373737374]) d_min_data = np.array([0.861904761904762, 0.861904761904762, 0.861904761904762, 0.8603809523809524, 0.8588571428571429, 0.8573333333333333, 0.8542857142857143, 0.8497142857142858, 0.8436190476190476, 0.832952380952381, 0.8192380952380952, 0.8024761904761905, 0.7826666666666666, 0.7598095238095238, 0.7293333333333333, 0.6881904761904761, 0.6577142857142857, 0.6257142857142857, 0.6013333333333333, 0.5784761904761905, 0.5586666666666666, 0.5373333333333333, 0.516, 0.5007619047619047, 0.480952380952381, 0.46571428571428575, 0.448952380952381, 0.4352380952380952, 0.4215238095238095, 0.4093333333333333, 0.39409523809523805, 0.3773333333333333, 0.3544761904761905, 0.29200000000000004, 0.24323809523809525, 0.2020952380952381])
In [23]:
plt.plot(sol.t*((R/u_char)),min_data(sol.y),'k',label = 'present work') # plt.plot(sol.t*((R/u_char)),max_data(sol.y),'k') plt.plot(t_min_data,d_min_data,'^b',label = 'Dietze and Ruyer-Quil (2015)') plt.xlabel(r'$t^{*}$') plt.ylabel(r'$d_{min}/R$') plt.legend()
Out[23]:
<matplotlib.legend.Legend at 0x7f963c15dbd0>
Note: While graphically extracting data from figure 9 of Dietze and Ruyer-Quil (2015), we have incurred errors that produce the small differences seen above.