Kernel: Python 3 (ipykernel)
Import useful modules
In [1]:
%matplotlib inline import numpy as np from pylab import * import numpy as np import matplotlib.pyplot as plt import sympy as sym
Plot format
In [2]:
SMALL_SIZE = 10 MEDIUM_SIZE = 20 BIGGER_SIZE = 24 plt.rc('font', size=BIGGER_SIZE) # controls default text sizes plt.rc('axes', titlesize=SMALL_SIZE) # fontsize of the axes title plt.rc('axes', labelsize=BIGGER_SIZE) # fontsize of the x and y labels plt.rc('xtick', labelsize=SMALL_SIZE) # fontsize of the tick labels plt.rc('ytick', labelsize=SMALL_SIZE) # fontsize of the tick labels plt.rc('legend', fontsize=SMALL_SIZE) # legend fontsize plt.rc('axes', labelpad=30) # rotation plt.rc('text', usetex=True) mpl.rcParams['figure.dpi'] = 300
Load the reference for validation: Dabade et al., JFM, 2016
In [3]:
###Particle inertia #prolate: Figure 2.b, \xi_0 = 1.01 rod_pine = np.loadtxt('Rod_particle_ine_drift.txt') #oblate: Figure 3, \xi_0 = 1.0001 disk_pine = np.loadtxt('Disk_particle_ine_drift.txt') ###Fluid inertia #prolate: Figure 5.b, \xi_0 = 1.0001 rod_fine = np.loadtxt('Rod_fluid_ine_drift.txt') #oblate: Figure 6.b, \xi_0 = 1.0001 disk_fine = np.loadtxt('Disk_fluid_ine_drift.txt')
Define the symbols:
In [4]:
###kappa k = sym.symbols('\kappa') ###\xi_0 xi = sym.symbols('xi', positive=True) ###\overline{\xi_0} xib= sym.symbols("\overline{xi}", positive=True) ###orbit constant C C = sym.symbols('C')
Define the inverse hyperbolic cotangent function
In [5]:
###inverse hyperbolic cotangent function def invcoth(val): return 0.5 * sym.log((val+1)/(val-1))
Define the aspect ratio depending functions Fip​ ,Gip​ for particle inertia: equations 5.7−5.12
In [6]:
F1p_def = (4*xi**4-7*xi**2+3)*(xi**2*invcoth(xi)+invcoth(xi)-xi)/(160*xi*(1-2*xi**2)**2) F2p_def = xib**4*(invcoth(xi)+xi*(xi*invcoth(xi)-1))/(40*xi*(1-2*xi**2)**2) F3p_def = - xib**2*(xi**2*invcoth(xi)+invcoth(xi)-xi)/(160*xi*(1-2*xi**2)**2) F4p_def = F3p_def F5p_def = -F3p_def / 2 F6p_def = F5p_def ### G1p_def = -(3*xi**4-5*xi**2+2)*((xi**2+1)*invcoth(xi)-xi) / (40*xi*(1-2*xi**2)**2) G2p_def = (xi*(-xi**3+xi+(xi**4-1)*invcoth(xi))) / (40*(1-2*xi**2)**2) G3p_def = G2p_def / xi**2 G4p_def = - G2p_def / xi**2
Define the aspect ratio depending functions Fif​, Gif​ for fluid inertia: equations 6.1−6.8
In [7]:
F1f_def = (xi**(2)*(-648*xi**(12)+1350*xi**(10)-5571*xi**(8)+11841*xi**(6)-9269*xi**(4)+2263*xi**(2)+6) \ -27*xi**(2)*(24*xi**(8)-14*xi**(6)-19*xi**(4)+16*xi**(2)-3)*xib**(8)*invcoth(xi)**(4)+9*xi \ *(288*xi**(12)-564*xi**(10)-20*xi**(8)+799*xi**(6)-743*xi**(4)+261*xi**(2)-29)*xib**(4) \ * invcoth(xi)**(3)+xi*(2592*xi**(14)-7020*xi**(12)+13932*xi**(10)-21123*xi**(8)+14255*xi**(6) \ -577*xi**(4)-2711*xi**(2)+652)*invcoth(xi)-3*(1296*xi**(16)-4320*xi**(14)+5346*xi**(12) \ -1477*xi**(10)-4260*xi**(8)+6116*xi**(6)-3492*xi**(4)+849*xi**(2)-58)*invcoth(xi)**(2)) \ *(480*xi**(2)*(-1+2*xi**(2))**(3)*(-3*xi**(2)+3*xib**(2)*xi*invcoth(xi)+2) \ *(-3*xi**(3)+5*xi+3*xib**(4)*invcoth(xi))*((3*xi**(2)-1)*invcoth(xi)-3*xi))**(-1) F2f_def = -(xib**(2)*(-9*xi**(9)+30*xi**(7)-115*xi**(5)+90*xi**(3)-12*xi \ +9*xib**(8)*(xi**(2)+1)*xi**(2)*invcoth(xi)**(3)-3*xib**(4)*(9*xi**(6)-10*xi**(4)-17*xi**(2)+14)*xi*invcoth(xi)**(2) \ +(27*xi**(10)-87*xi**(8)+133*xi**(6)-33*xi**(4)-52*xi**(2)+12)*invcoth(xi))) \ *(40*(xi-2*xi**(3))**(2)*(-3*xi**(2)+3*xib**(2)*xi*invcoth(xi)+2) \ *(-3*xi**(3)+5*xi+3*xib**(4)*invcoth(xi)))**(-1) F3f_def = -(xi**(2)*(378*xi**(10)+801*xi**(8)-4731*xi**(6)+5551*xi**(4)-2369*xi**(2)+342) \ -27*xi**(2)*(6*xi**(6)+xi**(4)-4*xi**(2)+1)*xib**(8)*invcoth(xi)**(4)+9*xi*(12*xi**(10)+28*xi**(8)-201*xi**(6) \ +273*xi**(4)-147*xi**(2)+27)*xib**(4)*invcoth(xi)**(3)+(-972*xi**(13)-324*xi**(11)+7365*xi**(9) \ -10409*xi**(7)+5143*xi**(5)-847*xi**(3)+44*xi)*invcoth(xi)+3*(216*xi**(14)-378*xi**(12) \ +109*xi**(10)-412*xi**(8)+1204*xi**(6)-1028*xi**(4)+311*xi**(2)-22)*invcoth(xi)**(2)) \ *(480*xi**(2)*(-1+2*xi**(2))**(3)*(-3*xi**(2)+3*xib**(2)*xi*invcoth(xi)+2)*(-3*xi**(3)+5*xi \ +3*xib**(4)*invcoth(xi))*((3*xi**(2)-1)*invcoth(xi)-3*xi))**(-1) F4f_def = F3f_def F5f_def = - F3f_def / 2.0 F6f_def = - F3f_def / 2.0 ### G1f_def = (xi**(2)*(81*xi**(10)-414*xi**(8)+1074*xi**(6)-1162*xi**(4)+479*xi**(2)-54) \ +9*xi**(2)*(9*xi**(6)-7*xi**(2)+2)*xib**(8)*invcoth(xi)**(4)-3*xi*(108*xi**(10)-246*xi**(8)+69*xi**(6) \ +167*xi**(4)-129*xi**(2)+23)*xib**(4)*invcoth(xi)**(3)+(-324*xi**(13)+1566*xi**(11)-3309*xi**(9) \ +3133*xi**(7)-1023*xi**(5)-79*xi**(3)+36*xi)*invcoth(xi) \ +(486*xi**(14)-2214*xi**(12)+3819*xi**(10)-2568*xi**(8)-222*xi**(6)+1036*xi**(4)-355*xi**(2)+18) \ *invcoth(xi)**(2))*(40*(xi-2*xi**(3))**(2)*(-3*xi**(2)+3*xib**(2)*xi*invcoth(xi)+2) \ *(-3*xi**(3)+5*xi+3*xib**(4)*invcoth(xi))*((3*xi**(2)-1)*invcoth(xi)-3*xi))**(-1) G2f_def = (-xi**(2)*(27*xi**(10)-180*xi**(8)+204*xi**(6)+68*xi**(4)-133*xi**(2)+18) \ -9*xi**(4)*(3*xi**(4)+2*xi**(2)-1)*xib**(8)*invcoth(xi)**(4)+3*xi \ *(36*xi**(10)-78*xi**(8)+73*xi**(6)-69*xi**(4)+35*xi**(2)-5)*xib**(4)*invcoth(xi)**(3) \ +xi*(108*xi**(12)-630*xi**(10)+1041*xi**(8)-617*xi**(6)+115*xi**(4)-29*xi**(2)+12) \ *invcoth(xi)+(-162*xi**(14)+810*xi**(12)-1551*xi**(10)+1600*xi**(8) \ -1054*xi**(6)+448*xi**(4)-97*xi**(2)+6)*invcoth(xi)**(2)) \ *(40*(xi-2*xi**(3))**(2)*(-3*xi**(2)+3*xib**(2)*xi*invcoth(xi)+2) \ *(-3*xi**(3)+5*xi+3*xib**(4)*invcoth(xi))*((3*xi**(2)-1)*invcoth(xi)-3*xi))**(-1) G3f_def = (xi**(2)*(378*xi**(10)+801*xi**(8)-4731*xi**(6)+5551*xi**(4)-2369*xi**(2)+342) \ -27*xi**(2)*(6*xi**(6)+xi**(4)-4*xi**(2)+1)*xib**(8)*invcoth(xi)**(4) \ +9*xi*(12*xi**(10)+28*xi**(8)-201*xi**(6)+273*xi**(4)-147*xi**(2)+27)*xib**(4)*invcoth(xi)**(3) \ +(-972*xi**(13)-324*xi**(11)+7365*xi**(9)-10409*xi**(7)+5143*xi**(5)-847*xi**(3)+44*xi) \ *invcoth(xi)+3*(216*xi**(14)-378*xi**(12)+109*xi**(10)-412*xi**(8)+1204*xi**(6)-1028*xi**(4) \ +311*xi**(2)-22)*invcoth(xi)**(2))*(120*xi**(2)*(2*xi**(2)-1)**(3)*(-3*xi**(2)+3*xib**(2)*xi*invcoth(xi)+2) \ *(-3*xi**(3)+5*xi+3*xib**(4)*invcoth(xi))*((3*xi**(2)-1)*invcoth(xi)-3*xi))**(-1) G4f_def = - G3f_def
Define the particles
In [8]:
###dictionary with the suggested values of \xi_0 to match the reference figures for validation ee = {'pro_pine' : 1.01, 'obl_pine' : 1.0001, 'pro_fine' : 1.1, 'obl_fine' : 1.0001} ###definition of particle aspect ratio according to prolate/oblate particle shape rr = {'prolate' : lambda xi0 : xi0/(xi0**2-1)**0.5, 'oblate' : lambda xi0 : (xi0**2-1)**0.5/xi0}
Function that evaluates the shape dependent coefficients for both particle and fluid inertia:the prolate−oblate transformation is automatically applied.
In [9]:
def shape_coeffs(key,ee): if(key == 'prolate'):#prolates ###Particle inertia e0 = ee['pro_pine'] F1p_eff = F1p_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf() F2p_eff = F2p_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf() F3p_eff = F3p_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf() F4p_eff = F4p_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf() F5p_eff = F5p_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf() F6p_eff = F6p_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf() G1p_eff = G1p_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf() G2p_eff = G2p_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf() G3p_eff = G3p_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf() G4p_eff = G4p_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf() ###Fluid inertia e0 = ee['pro_fine'] F1f_eff = F1f_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf() F2f_eff = F2f_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf() F3f_eff = F3f_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf() F4f_eff = F4f_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf() F5f_eff = F5f_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf() F6f_eff = F6f_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf() G1f_eff = G1f_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf() G2f_eff = G2f_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf() G3f_eff = G3f_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf() G4f_eff = G4f_def.subs({xib:(xi**2-1)**0.5,xi:e0}).evalf() else:#oblates ###Particle Inertia e0 = ee['obl_pine'] F1p_tmp = ((F1p_def*xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) F1p_eff = -sym.re(F1p_tmp.subs({xi:e0}).evalf()) F2p_tmp = ((F2p_def*xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) F2p_eff = -sym.re(F2p_tmp.subs({xi:e0}).evalf()) F3p_tmp = ((F3p_def*xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) F3p_eff = -sym.re(F3p_tmp.subs({xi:e0}).evalf()) F4p_tmp = ((F4p_def*xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) F4p_eff = -sym.re(F4p_tmp.subs({xi:e0}).evalf()) F5p_tmp = ((F5p_def*xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) F5p_eff = -sym.re(F5p_tmp.subs({xi:e0}).evalf()) F6p_tmp = ((F6p_def*xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) F6p_eff = -sym.re(F6p_tmp.subs({xi:e0}).evalf()) G1p_tmp = ((G1p_def*xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) G1p_eff = -sym.re(G1p_tmp.subs({xi:e0}).evalf()) G2p_tmp = ((G2p_def*xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) G2p_eff = -sym.re(G2p_tmp.subs({xi:e0}).evalf()) G3p_tmp = ((G3p_def*xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) G3p_eff = -sym.re(G3p_tmp.subs({xi:e0}).evalf()) G4p_tmp = ((G4p_def*xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) G4p_eff = -sym.re(G4p_tmp.subs({xi:e0}).evalf()) ###Fluid inertia e0 = ee['obl_fine'] F1f_tmp = ((F1f_def * xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) F1f_eff = -sym.re(F1f_tmp.subs({xi:e0}).evalf()) F2f_tmp = ((F2f_def * xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) F2f_eff = -sym.re(F2f_tmp.subs({xi:e0}).evalf()) F3f_tmp = ((F3f_def * xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) F3f_eff = -sym.re(F3f_tmp.subs({xi:e0}).evalf()) F4f_tmp = ((F4f_def * xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) F4f_eff = -sym.re(F4f_tmp.subs({xi:e0}).evalf()) F5f_tmp = ((F5f_def * xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) F5f_eff = -sym.re(F5f_tmp.subs({xi:e0}).evalf()) F6f_tmp = ((F6f_def * xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) F6f_eff = -sym.re(F6f_tmp.subs({xi:e0}).evalf()) G1f_tmp = ((G1f_def * xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) G1f_eff = -sym.re(G1f_tmp.subs({xi:e0}).evalf()) G2f_tmp = ((G2f_def * xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) G2f_eff = -sym.re(G2f_tmp.subs({xi:e0}).evalf()) G3f_tmp = ((G3f_def * xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) G3f_eff = -sym.re(G3f_tmp.subs({xi:e0}).evalf()) G4f_tmp = ((G4f_def * xi**2).subs({xib:(xi**2-1)**0.5,xi:sym.I*(xi**2-1)**0.5}))/(xi**2) G4f_eff = -sym.re(G4f_tmp.subs({xi:e0}).evalf()) return (F1p_eff,F2p_eff,F3p_eff,F4p_eff,F5p_eff,F6p_eff,G1p_eff,G2p_eff,G3p_eff,G4p_eff, F1f_eff,F2f_eff,F3f_eff,F4f_eff,F5f_eff,F6f_eff,G1f_eff,G2f_eff,G3f_eff,G4f_eff)
Define the trigonometric integrals Ii​ and Ji​ introduced in AppendixC
In [10]:
###I1 I1=2*sym.pi ###I2 I2=2*sym.pi*(k-1)*(k+1)**(-1) ###I3 I3=2*sym.pi*(2*((C**(2)+1)*(C**(2)*k**(2)+1))**(-1/2)-1) ###I4 I4=2*sym.pi*(k-1)**(2)*(k+1)**(-2) ###I5+I6 I5I6 = -(4*sym.pi*(2*k**(2)*(3*sym.sqrt((C**(2)+1)*(C**(2)*k**(2)+1))-8*C**(2)-6) \ +4*k*sym.sqrt((C**(2)+1)*(C**(2)*k**(2)+1)) \ +sym.sqrt((C**(2)+1)*(C**(2)*k**(2)+1))+4*(4*C**(2)+1)*k**(3)*sym.sqrt((C**(2)+1)*(C**(2)*k**(2)+1)) \ +k**(4)*(sym.sqrt((C**(2)+1)*(C**(2)*k**(2)+1))-16*(C**(4)+C**(2))-2)-2)) \ *((k**(2)-1)**(2)*sym.sqrt((C**(2)+1)*(C**(2)*k**(2)+1)))**(-1) ###J1 J1=sym.pi*(k**(2)-1)*(k+1)**(-2) ###J2 J2 = -sym.pi*(-4*((C**2+1)*(C**2*k**2+1))**0.5+C**2*(k+1)**2+4)*C**(-2)*(k**2-1)**(-1) ###J3 J3=sym.pi*(k-1)**(2)*(k+1)**(-2) ###J4 J4 = -sym.pi*(k**(2)-1)*(8*C**(4)*k**(3)+C**(2)*((k+1)**(4)-8*k**(2)*sym.sqrt((C**(2)+1)*(C**(2)*k**(2)+1))) \ -4*(k**(2)+1)*(sym.sqrt((C**(2)+1)*(C**(2)*k**(2)+1))-1))*C**(-2)*(k**(2)-1)**(-3)
Evaluate the shape dependent coefficients for a prolate particle
In [11]:
F1p,F2p,F3p,F4p,F5p,F6p,G1p,G2p,G3p,G4p,F1f,F2f,F3f,F4f,F5f,F6f,G1f,G2f,G3f,G4f = shape_coeffs('prolate',ee)
Evaluate the irreversible drift of a prolate particle due to particle inertia ΔCp​, as in eq. 5.19
In [12]:
r = rr['prolate'](ee['pro_pine']) #select the correct aspect ratio to macth figure 2 PINE_pro = (C*(2*xi**2-1)/(xi*xib)*((I1*F1p+I2*F2p+I3*F3p+I4*F4p+(I5I6)*F5p) \ +(J1*G1p+J2*G2p+J3*G3p+J4*G4p))).subs({xib:(xi**2-1)**0.5,xi:ee['pro_pine'],k:r}).evalf()
Evaluate the irreversible drift of a prolate particle due to fluid inertia ΔCf​, as in eq.5.19
In [13]:
r = rr['prolate'](ee['pro_fine']) #select the correct aspect ratio to match figure 5.b FINE_pro = (C*(2*xi**2-1)/(xi*xib)*((I1*F1f+I2*F2f+I3*F3f+I4*F4f+(I5I6)*F5f) \ +(J1*G1f+J2*G2f+J3*G3f+J4*G4f))).subs({xib:(xi**2-1)**0.5,xi:ee['pro_fine'],k:r}).evalf()
Evaluate the shape dependent coefficients for an oblate particle
In [14]:
F1p,F2p,F3p,F4p,F5p,F6p,G1p,G2p,G3p,G4p,F1f,F2f,F3f,F4f,F5f,F6f,G1f,G2f,G3f,G4f = shape_coeffs('oblate',ee)
Evaluate the irreversible drift of an oblate particle due to particle inertia ΔCp​, as in eq. 5.19
In [15]:
r = rr['oblate'](ee['obl_pine']) #select the correct aspect ratio to macth figure 3 PINE_obl = (C*(2*xi**2-1)/(xi*xib)*((I1*F1p+I2*F2p+I3*F3p+I4*F4p+(I5I6)*F5p) \ +(J1*G1p+J2*G2p+J3*G3p+J4*G4p))).subs({xib:(xi**2-1)**0.5,xi:ee['obl_pine'],k:r}).evalf()
Evaluate the irreversible drift of an oblate particle due to fluid inertia ΔCf​, as in eq.5.19
In [16]:
r = rr['oblate'](ee['obl_fine']) #select the correct aspect ratio to match figure 6.b FINE_obl= (C*(2*xi**2-1)/(xi*xib)*((I1*F1f+I2*F2f+I3*F3f+I4*F4f+(I5I6)*F5f) \ +(J1*G1f+J2*G2f+J3*G3f+J4*G4f))).subs({xib:(xi**2-1)**0.5,xi:ee['obl_fine'],k:r}).evalf()
Define a custom range of orbit constants C
In [17]:
Crange = np.concatenate((np.linspace(0.00001,0.9,100),np.linspace(1,4,100))) Crange = np.concatenate((Crange,np.linspace(4,100,100)))
Calculate the drift over the C range
In [18]:
CDC_pro_pine = [] #empty list to store the result for particle inertia on prolate CDC_pro_fine = [] #empty list to store the result for fluid inertia on prolate CDC_obl_pine = [] #empty list to store the result for particle inertia on oblate CDC_obl_fine = [] #empty list to store the result for fluid inertia on oblate for j in range(len(Crange)): #loop on the C range CVAL = Crange[j] #set the orbit constant value ###prolate, particle inertia data = PINE_pro.subs({C:CVAL}) #substitute the numerical values in the symbolic equation CDC_pro_pine.append([CVAL,data]) #store the result in the list ###prolate, fluid inertia data = FINE_pro.subs({C:CVAL}) #substitute the numerical values in the symbolic equation CDC_pro_fine.append([CVAL,data]) #store the result in the list ###oblate, particle inertia data = PINE_obl.subs({C:CVAL}) #substitute the numerical values in the symbolic equation CDC_obl_pine.append([CVAL,data]) #store the result in the list ###oblate, fluid inertia data = FINE_obl.subs({C:CVAL}) #substitute the numerical values in the symbolic equation CDC_obl_fine.append([CVAL,data]) #store the result in the list #convert to numpy arrays CDC_pro_pine = np.array(CDC_pro_pine) CDC_pro_fine = np.array(CDC_pro_fine) CDC_obl_pine = np.array(CDC_obl_pine) CDC_obl_fine = np.array(CDC_obl_fine)
Validate the method: fig.2,b for the prolate particle inertia and fig. 3 for the oblate particle inertia
In [19]:
###Create the figure with 4 sub-plots fig, axs = plt.subplots(2,2,constrained_layout=True,sharex=True,figsize=(7,4)) ###Figure 2.b of Dabade et al., 2016 ftr = -1/((ee['pro_pine']-1)**(3/2)*np.log(ee['pro_pine']-1)) axs[0,0].scatter(CDC_pro_pine[:,0]/(1+CDC_pro_pine[:,0]),CDC_pro_pine[:,1]/(1+CDC_pro_pine[:,0]**2)*ftr,label=r'$\mathit{\xi_0 = 1.01}$') axs[0,0].plot(rod_pine[:,0],rod_pine[:,1],c='k',label=r'$\mathit{Dabade \ et \ al., \ 2016}$',ls='--') axs[0,0].set_ylim(-0.01,0.4) axs[0,0].set_xlim(0.0,1.0) axs[0,0].legend(fontsize=8,loc='lower center') axs[0,0].title.set_text(r'$\mathit{Prolate, \ particle \ inertia}$') axs[0,0].set_ylabel(r'$\mathit{ \frac{ {St}^{-1} \left( \Delta C_p / (C^2+1) \right) }{ \left( (\xi_0^2-1)^{3/2} \log(\xi_0-1) \right)} }$',fontsize=7,labelpad=5) ###Figure 5.b of Dabade et al., 2016 ftr = -((ee['pro_fine']-1)**(0.5)*np.log(ee['pro_fine']-1)) axs[0,1].scatter(CDC_pro_fine[:,0]/(1+CDC_pro_fine[:,0]),CDC_pro_fine[:,1]/(1+CDC_pro_fine[:,0]**2)*ftr,label=r'$\mathit{\xi_0 = 1.1}$') axs[0,1].plot(rod_fine[:,0],rod_fine[:,1],c='k',label=r'$\mathit{Dabade \ et \ al., \ 2016}$',ls='--') axs[0,1].title.set_text(r'$\mathit{Prolate, \ fluid \ inertia}$') axs[0,1].set_ylim(-0.01,0.2) axs[0,1].set_ylabel(r'$\mathit{\frac{{Re_p}^{-1} \left( \Delta C_p / (C^2+1) \right) }{ \left( (\xi_0-1)^{0.5} \log (\xi_0-1) \right)^{-1}} }$',fontsize=7,labelpad=5) axs[0,1].legend(fontsize=8,loc='upper center') ###Figure 3 of Dabade et al., 2016 ftr = ee['obl_pine']**2 axs[1,0].scatter(CDC_obl_pine[:,0]/(1+CDC_obl_pine[:,0]),CDC_obl_pine[:,1]/(1+CDC_obl_pine[:,0]**2)*ftr,label=r'$\mathit{\xi_0 = 1.0001}$') axs[1,0].plot(disk_pine[:,0],disk_pine[:,1],c='k',label=r'$\mathit{Dabade \ et \ al., \ 2016}$',ls='--') axs[1,0].title.set_text(r'$\mathit{Oblate, \ particle \ inertia}$') axs[1,0].legend(fontsize=8,loc='upper center') axs[1,0].set_ylim(-0.11,0.01) axs[1,0].set_ylabel(r'$\mathit{{St}^{-1} \left( \Delta C_p / (C^2+1) \right) \xi_0^2 }$',fontsize=7,labelpad=5) ###Figure 6.b of Dabade et al., 2016 ftr = (ee['obl_fine']-1)**0.5 axs[1,1].scatter(CDC_obl_fine[:,0]/(1+CDC_obl_fine[:,0]),CDC_obl_fine[:,1]/(1+CDC_obl_fine[:,0]**2)*ftr,label=r'$\mathit{\xi_0 = 1.0001}$') axs[1,1].plot(disk_fine[:,0],disk_fine[:,1],c='k',label=r'$\mathit{Dabade \ et \ al., \ 2016}$',ls='--') axs[1,1].title.set_text(r'$\mathit{Oblate, \ fluid \ inertia}$') axs[1,1].legend(fontsize=8,loc='upper center') axs[1,1].set_ylim(-0.16,0.02) axs[1,1].set_ylabel(r'$\mathit{{Re_p}^{-1} \left( \Delta C_f / (C^2+1) \right) (\xi_0-1)^{0.5} }$',fontsize=6.5,labelpad=5) ###Format for j in range(2): axs[1,j].set_xlabel(r'$\mathit{C/(C+1)}$',labelpad=5,fontsize=13)
Out[19]:
In [ ]: