Changed force plots and fixed some inconsistencies between VFT and Eyring Fluid.

from forces import F_newt_isoth,F_newt_vft,F_newt_eyr

if __name__=="__main__":

def geometry_1():

    alpha= 30 # degree

    Ro= 0.002/2 # m

    Ri= 0.0002 # m

    Lc= 0.0006 # m 

    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 PS_material():

    km = 0.26 # m s /K

    Tm = 378.15 # K

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

    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 

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

def PS_settings():

    Th = 513.15 # K

    T0 = 293.15 # K

    return Th,T0

def vft_forces(Us,

               material,

               geometry,

               settings,

               debug_iso=False,

               debug_non=False):

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

    Th,T0 = settings()

    alpha,Ro,Ri,Lc = geometry()

    # calculate viscosity in capillary

    etac = eta0 * np.exp(eps/(Th-Tv))

    # 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)), 

                         alpha=alpha, Ro=Ro, Ri=Ri, etac=etac, Lc=Lc,

                         debug=debug_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)

    return F_iso,F_non

def eyr_forces(Us,

               material,

               geometry,

               settings,

               debug_iso=False,

               debug_non=False):

    # extract material data, geometry data and process settings 

    km,Th,Tm,rho_s,rho_m,h,Cs,Q,eta0 = material()

    Th,T0 = settings()

    alpha,Ro,Ri,Lc = geometry()

    # calculate viscosity in capillary

    etac = eta0 * np.exp(Q/(R*Th))

    # 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)), 

                         alpha=alpha, Ro=Ro, Ri=Ri, etac=etac, Lc=Lc,

                         debug=debug_iso)

    # calculate nonisothermal

    F_non = F_newt_eyr(km=km, Th=Th, Tm=Tm, rho_s=rho_s, rho_m=rho_m, h=h, 

                       Cs=Cs, T0=T0, Us=Us, eta0=eta0, Q=Q, R=R, alpha=alpha, 

                       Ro=Ro, Ri=Ri, etac=etac, Lc=Lc, 

                       debug=debug_non)

    return F_iso,F_non

if __name__=="__main__":

    km = 0.163 # m s /K

    Th = 473 # K

    Tm = 424 # K

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

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

    h = 93600 # J/mol

    Cs = 1800 # J/ mol K 

    T0 = 293 # K

    alpha=56.31 # degree

    Ro= 0.00197/2 # m

    Ri= 0.0002 # m

    Lc= 0.001 # m 

    Q = 80.574 # J/mol

    R = 8.314 # J /(Kmol)

    eps = Q/R # K

    Tv = 0 # K

    eta0 = 3.546 * 10**3 # Pas 

    etac = eta0 * np.exp(Q/(R*Th))

    #

    fig,ax = plt.subplots(1,1)

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

    # create grid of velocities

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

    # calculate isothermal

    F = 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)), 

                     alpha=alpha, Ro=Ro, Ri=Ri, etac=etac, Lc=Lc,

                     debug=False)

    # calculate forces

    F_iso,F_non = vft_forces(Us,

                             PS_material,

                             geometry_1,

                             PS_settings)

    # plot isoth. plot

    ax.plot(Us*10**3,F,label="isoth.")

    axs[0].plot(Us*10**3,F_iso,label="isoth.")

    # plot nonisoth. plot

    axs[0].plot(Us*10**3,F_iso,label="nonisoth.")

    # arrhenius

    F = F_newt_eyr(km=km, Th=Th, Tm=Tm, rho_s=rho_s, rho_m=rho_m, h=h, Cs=Cs, 

                   T0=T0, Us=Us, eta0=eta0, Q=Q, R=R, alpha=alpha, Ro=Ro, Ri=Ri, 

                   etac=etac, Lc=Lc, debug=False)

    # plot arrh.

    ax.plot(Us*10**3,F,label="arrh.")

    # vogel fulcher tamann 

    F = F_newt_vft(km, Th, Tm, rho_s, rho_m, h, Cs, T0,

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

                   debug=False)

    # plot vft.

    ax.plot(Us*10**3,F,label="vft.")

    ax.legend()

    ax.set_xlabel(r"$U_{s} (m/s)$")

    ax.set_ylabel(r"$U_{F} (N)$")

    ax.set_yscale("log")

    # calculate forces

    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.")

    # plot nonisoth. plot

    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")

    plt.show()

 No newline at end of file