From c9250ea36ea8f23133e40d30fe6c9d8e55f8cb8b Mon Sep 17 00:00:00 2001 From: shiemer Date: Thu, 9 Nov 2023 15:52:33 +0100 Subject: [PATCH] Added splitting the forces into their contribution. There seems to be a sign mistake in F_tau and an overflow for Fp of the nonisothermal model. --- force_plots.py | 106 ++++++++++++++++++++++++++++++++----------------- forces.py | 88 ++++++++++++++++++++++++++++------------ 2 files changed, 132 insertions(+), 62 deletions(-) diff --git a/force_plots.py b/force_plots.py index afbe4fc..2320917 100644 --- a/force_plots.py +++ b/force_plots.py @@ -9,6 +9,8 @@ import matplotlib.pyplot as plt from forces import F_newt_isoth,F_newt_vft,F_newt_eyr +R = 8.314 + def geometry_1(): alpha= 30 # degree Ro= 0.002/2 # m @@ -17,19 +19,17 @@ def geometry_1(): return alpha,Ro,Ri,Lc -def PP_values(): - - km = # m s /K - Th = # K - Tm = # K - rho_s = # kg m**(-3) - rho_m = # kg m**(-3) - h = # J/mol (0 für PS) - Cs = # J/ mol K - Q = # J/mol - eta0 = # Pas - - return km,Th,Tm,rho_s,rho_m,h,Cs,Q,eta0 +def PP_material(): + km = 0# m s /K + Tm = 0# K + rho_s = 0# kg m**(-3) + rho_m = 0# kg m**(-3) + h = 0# J/mol (0 für PS) + Cs = 0# J/ mol K + Q = 0# J/mol + eta0 = 0# Pas + + return km,Tm,rho_s,rho_m,h,Cs,Q,eta0 def PS_material(): @@ -39,11 +39,11 @@ def PS_material(): rho_m = 917 # kg m**(-3) h = 0 # J/mol Cs = 1630 # J/ mol K - eps = # K - Tv = 308.1 # K - eta0 = 2.27 * 10**3 # Pas + eps = 2303.99 # K + Tv = 287.03 # K + eta0 = 0.08549 # Pas - return km,Th,Tm,rho_s,rho_m,h,Cs,T0,eps,Tv,eta0 + return km,Tm,rho_s,rho_m,h,Cs,eps,Tv,eta0 def PS_settings(): Th = 513.15 # K @@ -52,7 +52,7 @@ def PS_settings(): return Th,T0 def PP_settings(): - Th = # K + Th = 0 # K T0 = 293.15 # K return Th,T0 @@ -62,7 +62,9 @@ def vft_forces(Us, geometry, settings, debug_iso=False, - debug_non=False): + split_iso=False, + debug_non=False, + split_non=False): km,Tm,rho_s,rho_m,h,Cs,eps,Tv,eta0 = material() Th,T0 = settings() @@ -74,12 +76,14 @@ def vft_forces(Us, # calculate isothermal F_iso = F_newt_isoth(km=km, Th=Th, Tm=Tm, rho_s=rho_s, rho_m=rho_m, h=h, Cs=Cs, T0=T0, Us=Us, - eta0=eta0*np.exp(eps/((Th+Tm)/2 -Tv)), + eta_eff=eta0*np.exp(eps/((Th+Tm)/2 -Tv)), alpha=alpha, Ro=Ro, Ri=Ri, etac=etac, Lc=Lc, - debug=debug_iso) + debug=debug_iso, + split=split_iso) F_non = F_newt_vft(km, Th, Tm, rho_s, rho_m, h, Cs, T0, Us, eta0, eps, Tv, alpha, Ro, Ri, etac, Lc, - debug=debug_non) + debug=debug_non, + split=split_non) return F_iso,F_non @@ -91,7 +95,7 @@ def eyr_forces(Us, debug_non=False): # extract material data, geometry data and process settings - km,Th,Tm,rho_s,rho_m,h,Cs,Q,eta0 = material() + km,Tm,rho_s,rho_m,h,Cs,Q,eta0 = material() Th,T0 = settings() alpha,Ro,Ri,Lc = geometry() @@ -101,7 +105,7 @@ def eyr_forces(Us, # calculate isothermal F_iso = F_newt_isoth(km=km, Th=Th, Tm=Tm, rho_s=rho_s, rho_m=rho_m, h=h, Cs=Cs, T0=T0, Us=Us, - eta0=eta0*np.exp(Q/(R*(Th+Tm)/2)), + eta_eff=eta0*np.exp(Q/(R*(Th+Tm)/2)), alpha=alpha, Ro=Ro, Ri=Ri, etac=etac, Lc=Lc, debug=debug_iso) @@ -115,9 +119,6 @@ def eyr_forces(Us, if __name__=="__main__": - - etac = eta0 * np.exp(Q/(R*Th)) - # fig,axs = plt.subplots(1,2) @@ -128,28 +129,59 @@ if __name__=="__main__": F_iso,F_non = vft_forces(Us, PS_material, geometry_1, - PS_settings) + PS_settings, + debug_iso=False, + split_iso=False, + debug_non=False, + split_non=False) # plot isoth. plot - axs[0].plot(Us*10**3,F_iso,label="isoth.") + axs[0].plot(Us*10**3,F_iso,color="k",label="isoth.") + # plot nonisoth. plot - axs[0].plot(Us*10**3,F_iso,label="nonisoth.") + axs[0].plot(Us*10**3,F_non,color="b",label="nonisoth.") + + # calculate force contributions + F_iso,F_non = vft_forces(Us, + PS_material, + geometry_1, + PS_settings, + debug_iso=False, + split_iso=True, + debug_non=False, + split_non=True) + + # plot individual force contributions + axs[0].plot(Us*10**3,F_iso[0],color="k", + linestyle="--",label=r"$F_{p}$") + axs[0].plot(Us*10**3,F_iso[1],color="k", + linestyle="-.",label=r"$F_{pi}$") + axs[0].plot(Us*10**3,F_iso[2],color="k", + linestyle=":",label=r"$F_{\tau}$") + axs[0].plot(Us*10**3,F_non[0],color="b", + linestyle="--",label=r"$F_{p}$") + axs[0].plot(Us*10**3,F_non[1],color="b", + linestyle="-.",label=r"$F_{pi}$") + axs[0].plot(Us*10**3,F_non[2],color="b", + linestyle=":",label=r"$F_{\tau}$") # calculate forces - F_iso,F_non = eyr_forces(Us, - PP_material, - geometry_1, - PP_settings) + #F_iso,F_non = eyr_forces(Us, + # PP_material, + # geometry_1, + # PP_settings) # plot isoth. plot - axs[1].plot(Us*10**3,F_iso,label="isoth.") + #axs[1].plot(Us*10**3,F_iso,label="isoth.") # plot nonisoth. plot - axs[1].plot(Us*10**3,F_iso,label="nonisoth.") + #axs[1].plot(Us*10**3,F_iso,label="nonisoth.") for i in range(2): axs[i].legend() axs[i].set_xlabel(r"$U_{s} (mm/s)$") axs[i].set_ylabel(r"$U_{F} (N)$") - axs[i].set_yscale("log") + axs[i].set_yscale("symlog") + axs[0].set_title("PS") + axs[1].set_title("PP") plt.show() \ No newline at end of file diff --git a/forces.py b/forces.py index 6f96d9a..68bcd1a 100644 --- a/forces.py +++ b/forces.py @@ -11,10 +11,10 @@ from integrals import eta_1,int_eta_dz,int_etaz_dz,intint_eta_dzdz,intint_etaz_d def get_delta(km, Th, Tm, rho_s, Us, h, Cs, T0): return (km * (Th - Tm)) / (rho_s * Us * (h + Cs * (Tm - T0))) -def Ft_newt_isoth(U0, eta0, delta, alpha, Ro, Ri): +def Ft_newt_isoth(U0, eta_eff, delta, alpha, Ro, Ri): - prefac = 2 * np.pi * (eta0 * U0) / delta + prefac = 2 * np.pi * (eta_eff * U0) / delta term1 = 3 * ((1/np.cos(alpha)) * (delta**(-1)) * \ ((2/3) * Ro**3 - Ri * Ro**2 + (1/3) * Ri**3)\ + np.tan(alpha) * (Ro**2 /2 - 2*Ro*Ri + Ro**2 /2)) @@ -23,24 +23,27 @@ def Ft_newt_isoth(U0, eta0, delta, alpha, Ro, Ri): return prefac*(term1 + term2) -def Fp_newt_isoth(U0, eta0, delta, alpha, Ro, Ri, etac, Lc, - debug): +def Fp_newt_isoth(U0, eta_eff, delta, alpha, Ro, Ri, etac, Lc, + debug, split): # inlet_p = 8*np.pi*etac*U0*Lc*(Ro / Ri)**4 - Fp=3*np.pi*(U0*eta0)/(delta**2 *np.cos(alpha)**3) * \ + Fp=3*np.pi*(U0*eta_eff)/(delta**2 *np.cos(alpha)**3) * \ (((2*np.log(Ro/Ri) - 3/2) * Ro**4 + 2*Ri**2 *Ro**2 + Ri**4 /2)/delta \ + np.sin(alpha) * ((2*np.log(Ro/Ri) - 7/3)*Ro**3 + 2*Ri*Ro**2 +\ - Ri**2*Ro - 2/3 *Ri**3)) + inlet_p + Ri**2*Ro - 2/3 *Ri**3)) if debug: print("Pi",inlet_p) - - return Fp + if not split: + return Fp + inlet_p + else: + return Fp, inlet_p def F_newt_isoth(km, Th, Tm, rho_s, rho_m, h, Cs, T0, - Us, eta0, alpha, Ro, Ri, etac, Lc, - debug=False): + Us, eta_eff, alpha, Ro, Ri, etac, Lc, + debug=False, + split=False): """ @@ -64,7 +67,7 @@ def F_newt_isoth(km, Th, Tm, rho_s, rho_m, h, Cs, T0, room temperature. Us : float feeder velocity. - eta0 : float + eta_eff : float Newtonian viscosity. alpha : float angle. @@ -76,11 +79,24 @@ def F_newt_isoth(km, Th, Tm, rho_s, rho_m, h, Cs, T0, viscosity in capillary. Lc : float capillary length. + debug : boolean + if True, print out several intermediate results for debugging. + split : boolean + if True, return different force contributions. Returns ------- - F: float - feeder force. + F: np.array float + feeder force. Only provided if split not True + Fp: np.array float + pressure part of the force not stemming from the end capillary. + Only provided if split is True + Fpi: np.array float + pressure part of the force stemming from the end capillary. + Only provided if split is True + Ftau: np.array float + deviatoric stress part of the force. + Only provided if split is True """ @@ -99,13 +115,16 @@ def F_newt_isoth(km, Th, Tm, rho_s, rho_m, h, Cs, T0, "U0",U0, "delta",delta) - Ft = Ft_newt_isoth(U0, eta0, delta, alpha, Ro, Ri) - Fp = Fp_newt_isoth(U0, eta0, delta, alpha, Ro, Ri, etac, Lc, debug) + Ft = Ft_newt_isoth(U0, eta_eff, delta, alpha, Ro, Ri) + Fp = Fp_newt_isoth(U0, eta_eff, delta, alpha, Ro, Ri, etac, Lc, + debug, split) if debug: print("Ft",Ft) print("Fp",Fp) - - return np.cos(alpha) * Fp - np.sin(alpha) * Ft + if not split: + return np.cos(alpha) * Fp - np.sin(alpha) * Ft + else: + return np.cos(alpha) * Fp[0], np.cos(alpha) * Fp[1], - np.sin(alpha) * Ft def eta(z,eta0,A,B): return eta0 * np.exp(A/(B+z)) @@ -140,7 +159,7 @@ def Ft_newt_noniso(U0, eta0, delta, alpha, Ro, Ri, A, B): return Ft def Fp_newt_noniso(U0, eta0, delta, alpha, Ro, Ri, etac, Lc, A ,B, - debug): + debug, split): # inlet pressure inlet_p = 8*np.pi*etac*U0*Lc*(Ro ** 4)/(Ri**4) @@ -148,16 +167,19 @@ def Fp_newt_noniso(U0, eta0, delta, alpha, Ro, Ri, etac, Lc, A ,B, Fp =(-np.pi/4)*eta0*H(delta,A,B,delta)**(-1) * np.cos(alpha)**(-1) \ *(U0*(Ro**4 * (2*np.log(Ro/Ri) - 3/2) + 2*Ri**2 *Ro**2 - Ri**4 /2) \ - 2*W(delta,U0,alpha,A,B,delta)*((2*np.log(Ro/Ri) - 7/3)*Ro**3 + \ - 2*Ri*Ro**2 + Ri**2 *Ro - 2/3 *Ri**3)) + inlet_p + 2*Ri*Ro**2 + Ri**2 *Ro - 2/3 *Ri**3)) if debug: print("Pi",inlet_p) - return Fp + if not split: + return Fp + inlet_p + else: + return Fp, inlet_p def F_newt_vft(km, Th, Tm, rho_s, rho_m, h, Cs, T0, Us, eta0, eps, Tv, alpha, Ro, Ri, etac, Lc, - debug=False): + debug=False, split=False): """ Force for the Vogel-Fulcher-Tamann fluid. eta = eta0 exp(eps/(T-Tv)) @@ -197,12 +219,25 @@ def F_newt_vft(km, Th, Tm, rho_s, rho_m, h, Cs, T0, viscosity in capillary. Lc : float capillary length. + debug : boolean + if True, print out several intermediate results for debugging. + split : boolean + if True, return different force contributions. Returns ------- F: float feeder force. - + Fp: np.array float + pressure part of the force not stemming from the end capillary. + Only provided if split is True + Fpi: np.array float + pressure part of the force stemming from the end capillary. + Only provided if split is True + Ftau: np.array float + deviatoric stress part of the force. + Only provided if split is True + """ if debug: @@ -226,13 +261,16 @@ def F_newt_vft(km, Th, Tm, rho_s, rho_m, h, Cs, T0, "B",B) Ft = Ft_newt_noniso(U0, eta0, delta, alpha, Ro, Ri, A, B) - Fp = Fp_newt_noniso(U0, eta0, delta, alpha, Ro, Ri, etac, Lc, A, B, debug) + Fp = Fp_newt_noniso(U0, eta0, delta, alpha, Ro, Ri, etac, Lc, A, B, + debug, split) if debug: print("Ft",Ft,"\n", "Fp",Fp) - - return np.cos(alpha) * Fp - np.sin(alpha) * Ft + if not split: + return np.cos(alpha) * Fp - np.sin(alpha) * Ft + else: + return np.cos(alpha) * Fp[0], np.cos(alpha) * Fp[1], - np.sin(alpha) * Ft def F_newt_eyr(km, Th, Tm, rho_s, rho_m, h, Cs, T0, Us, eta0, Q, R, alpha, Ro, Ri, etac, Lc, -- GitLab