Commit 29142e98 authored by Stefan Hiemer's avatar Stefan Hiemer
Browse files

Added truncated Power-law inlet pressure for n=1/3 and changed some plotting.

parent c9250ea3

PS_240C_exp_data.txt

0 → 100644
+20 −0
Original line number Diff line number Diff line 
U in mm/s  

F in N



0.49946	2.5771

1.04957	3.20268

1.59717	3.83596

2.12195	4.4409

2.6108	5.07469

3.11025	5.70986

3.65511	6.48184

4.22108	7.43868

4.72305	8.72071

5.2022	10.37918

6.08429	15.1482

6.87888	21.98964

7.73485	26.18548





240 °C 

Polystrol
 
 No newline at end of file

force_plots.py

+18 −8
Original line number Diff line number Diff line
@@ -13,7 +13,7 @@ R = 8.314 


def geometry_1():

    alpha= 30 # degree

    Ro= 0.002/2 # m

    Ro= 0.001 # m

    Ri= 0.0002 # m

    Lc= 0.0006 # m
@@ -24,7 +24,7 @@ def PP_material(): 
    Tm = 0# K

    rho_s = 0# kg m**(-3)

    rho_m = 0# kg m**(-3)

    h = 0# J/mol               (0 für PS)

    h = 0# J/mol

    Cs = 0# J/ mol K 

    Q = 0# J/mol

    eta0 = 0# Pas
@@ -42,8 +42,10 @@ def PS_material(): 
    eps = 2303.99 # K

    Tv = 287.03 # K

    eta0 = 0.08549 # Pas 

    gammac = 6.27 # s**(-1)

    n = 1/3

    

    return km,Tm,rho_s,rho_m,h,Cs,eps,Tv,eta0

    return km,Tm,rho_s,rho_m,h,Cs,eps,Tv,eta0,gammac,n



def PS_settings():

    Th = 513.15 # K
@@ -66,7 +68,7 @@ def vft_forces(Us, 
               debug_non=False,

               split_non=False):

    

    km,Tm,rho_s,rho_m,h,Cs,eps,Tv,eta0 = material()

    km,Tm,rho_s,rho_m,h,Cs,eps,Tv,eta0,gammac,n = material()

    Th,T0 = settings()

    alpha,Ro,Ri,Lc = geometry()
@@ -81,7 +83,7 @@ def vft_forces(Us, 
                         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,

                       Us, eta0, eps, Tv, alpha, Ro, Ri, etac, Lc,gammac,n,

                       debug=debug_non,

                       split=split_non)
@@ -122,8 +124,16 @@ if __name__=="__main__": 
    #

    fig,axs = plt.subplots(1,2)

    

    # plot experimental data

    data = np.loadtxt("PS_240C_exp_data.txt",

                      skiprows=3,

                      max_rows=13)

    

    axs[0].scatter(data[:,0],data[:,1],

                   color="r",label="exp.")

    

    # create grid of velocities

    Us = np.linspace(0,0.006,101)[1:]

    Us = np.linspace(0,0.008,101)[1:]

    

    # calculate forces

    F_iso,F_non = vft_forces(Us,
@@ -150,7 +160,7 @@ if __name__=="__main__": 
                             split_iso=True,

                             debug_non=False,

                             split_non=True)

    

    print(F_non[1])

    # plot individual force contributions

    axs[0].plot(Us*10**3,F_iso[0],color="k",

                linestyle="--",label=r"$F_{p}$")
@@ -181,7 +191,7 @@ if __name__=="__main__": 
        axs[i].legend()

        axs[i].set_xlabel(r"$U_{s} (mm/s)$")

        axs[i].set_ylabel(r"$U_{F} (N)$")

        axs[i].set_yscale("symlog")

        axs[i].set_yscale("symlog",linthresh=1e-5)

    axs[0].set_title("PS")

    axs[1].set_title("PP")

    plt.show()
 
 No newline at end of file

forces.py

+21 −7
Original line number Diff line number Diff line
@@ -150,7 +150,7 @@ def W(z,U0,alpha,A,B,delta): 
def Ft_newt_noniso(U0, eta0, delta, alpha, Ro, Ri, A, B):

    

    

    Ft = 2*np.pi*eta0*eta(delta,eta0,A,B)*(-np.cos(alpha)**(-1) *(U0*((2/3)*Ro**3 \

    Ft = 2*np.pi*eta0*eta(delta,eta0,A,B)*(np.cos(alpha)**(-1) *(U0*((2/3)*Ro**3 \

        -Ri*Ro**2 +Ri**3 /3)-W(delta,U0,alpha,A,B,delta)*\

        (Ro**2 /2 -2*Ro*Ri+Ro**2 /2)) *\

        dhdz(delta,U0,alpha,A,B,delta)/H(delta,A,B,delta) \
@@ -158,13 +158,27 @@ 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,

def Fp_newt_noniso(U0, eta0, delta, alpha, Ro, Ri, etac, Lc, A ,B, gammac, n,

                   debug, split):

    

    # inlet pressure

    # Newtonian inlet pressure 

    if n is None and gammac is None:

        inlet_p = 8*np.pi*etac*U0*Lc*(Ro ** 4)/(Ri**4)

    # Power law inlet pressure

    elif n is not None and gammac is None:

        inlet_p = ((1/n + 3)/(Ri**3) * Ro**2 * U0)**n *eta0*2*Lc/Ri 

    # truncated pwoer law inlet pressure

    elif np.isclose(n,1/3) and gammac is not None:

        prefac = np.pi * Ri**4 / (2*(1/n + 3) * eta0 * Lc) 

        a = prefac * (Ri / (eta0 * gammac * 2 * Lc))**(1/n - 1)

        b = prefac * 1/4 * (1/n - 1) * (eta0 * gammac * 2 * Lc/Ri)**4 

        Q = np.pi * Ro**2 * U0

        inlet_p = Ro**2 * np.pi * np.cbrt(Q + np.sqrt(Q**2 - 4 * a *b) / (2*a))

    else:

        raise NotImplementedError()

    

    Fp =(-np.pi/4)*eta0*H(delta,A,B,delta)**(-1) * np.cos(alpha)**(-1) \

    # pressure from stuff above

    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))
@@ -178,7 +192,7 @@ def Fp_newt_noniso(U0, eta0, delta, alpha, Ro, Ri, etac, Lc, A ,B, 
        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,

               Us, eta0, eps, Tv, alpha, Ro, Ri, etac, Lc, gammac=None, n=None,

               debug=False, split=False):

    """

    Force for the Vogel-Fulcher-Tamann fluid. eta = eta0 exp(eps/(T-Tv))
@@ -262,7 +276,7 @@ def F_newt_vft(km, Th, Tm, rho_s, rho_m, h, Cs, T0, 
    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, split)

                        gammac,n,debug, split)

    

    if debug:

        print("Ft",Ft,"\n",
Impressum   Privacy Policy   WW8 Homepage   Explore   Help