Initial commit

#!/usr/bin/env python3

# -*- coding: utf-8 -*-

"""

@author: Stefan Hiemer

"""

import numpy as np

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

    prefac = 2 * np.pi * (eta0 * 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))

    term2 = - np.sin(alpha) * (Ro**2 - Ri**2) / 2

    return prefac*(term1 + term2)

def Fp_newt_isoth(U0, eta0, delta, alpha, Ro, Ri, etac, Lc):

    prefac = 3 * np.pi * (U0 * eta0) / (delta**2*np.cos(alpha)**3)

    term1 = ((2*np.log(Ro/Ri) - 3/2) * Ro**4 + 2*Ri**2 *Ro**2 + Ri**4 /2)/delta

    term1 += np.sin(alpha) * ((2*np.log(Ro/Ri) - 7/3)*Ro**3 + 2*Ri*Ro**2 +\

             Ri**2*Ro - 2/3 *Ri**3)

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

    return prefac * term1 + term2

def F_newt_isoth(km, Th, Tm, rho_s, rho_m, h, Cs, T0,

                 Us, eta0, delta, alpha, Ro, Ri, etac, Lc):

    """

    Parameters

    ----------

    km : float

        melt heat conductivity.

    Th : float

        heater temperature in Kelvin.

    Tm : float

        melting temperature in K.

    rho_s : float

        density solid

    rho_m : float

        density melt.

    h : float

        heat of fusion.

    Cs : float

        heat capacity solid.

    T0 : float

        room temperature.

    Us : float

        feeder velocity.

    eta0 : float

        Newtonian viscosity.

    delta : float

        film thickness.

    alpha : float

        angle.

    Ro : float

        outer radius.

    Ri : float

        inner radius.

    etac : float

        viscosity in capillary.

    Lc : float

        capillary length.

    Returns

    -------

    F: float

        feeder force.

    """

    U0= rho_s/rho_m * Us

    delta = get_delta(km, Th, Tm, rho_s, Us, h, Cs, T0)

    Ft = Ft_newt_isoth(U0, eta0, delta, alpha, Ro, Ri)

    Fp = Fp_newt_isoth(U0, eta0, delta, alpha, Ro, Ri, etac, Lc)

    return np.cos(alpha) * Fp - np.sin(alpha) * Ft

if __name__=="__main__":

    pass

#!/usr/bin/env python3

# -*- coding: utf-8 -*-

"""

@author: Stefan Hiemer

"""

import numpy as np 

from scipy.special import expi 

def int_eta_dx(z,A,B):

    return A*(expi(-A/(B+z))-expi(-A/(B))) + (z+B)*np.exp(-A/(B+z)) \

           - B*np.exp(-A/B)

def int_etax_dx(z,A,B):

    term1 = A*(A+2*B)*(expi(-A/B)-expi(-A/(B+z)))

    term2 = B*(A+B)*np.exp(-A/B) - (B+z)*(A+B-z)*np.exp(-A/(B+z))

    return 0.5*(term1+term2)

def intint_eta_dxdx(z,A,B):

    term1 = A*(A + 2*(B+z))*expi(-A/(B+z))

    term2 = (B+z)*(A+B+z)*np.exp(-A/(B+z)) 

    c1 = A*expi(-A/B) + B*np.exp(-A/B) 

    c2 = A*(A+ 2*B)*expi(-A/B)+B*(A+B)*np.exp(-A/B)

    return 0.5 * (term1 + term2 - 2*c1*z - c2) 

def intint_etax_dxdx(z,A,B):

    term1 = A*(A**2 + 6*A*B + 6*B**2 + 3*(A+2*B)*z)*expi(-A/(B+z))

    term2 = (B+z)*(A**2 + 5*A*B + 2*B**2 + (B + 2*A)*z - z**2)*np.exp(-A/(B+z))

    c1 =  A*(A+2*B)*expi(-A/B) + B*(A+B)*np.exp(-A/B)

    c2 = A*(A**2 + 6*A*B + 6*B**2)*expi(-A/(B))+\

         B*(A**2 + 5*A*B + 2*B**2)*np.exp(-A/(B))

    return (c2-term1 + term2)/6 + c1*z 

def eta(z, A,B):

    return np.exp(-A/(B+z))

def etax(z, A,B):

    return z*np.exp(-A/(B+z))

def test_integrals():

    from scipy.integrate import quad

    A = 2

    B = 5

    z = 10 

    if not np.isclose(int_eta_dx(z,A,B),quad(eta,a=0,b=z,args=(A,B))[0]):

        raise ValueError("Viscosity integral wrong.")

    if not np.isclose(int_etax_dx(z,A,B),quad(etax,a=0,b=z,args=(A,B))[0]):

        raise ValueError("Locally weighted viscosity integral wrong.")

    if not np.isclose(intint_eta_dxdx(z,A,B),quad(int_eta_dx,0,z, args=(A,B))[0]):

        raise ValueError("Viscosity double integral wrong.")

    if not np.isclose(intint_etax_dxdx(z,A,B),quad(int_etax_dx,a=0,b=z,args=(A,B))[0]):

        raise ValueError("Locally weighted viscosity double integral wrong.")

    return 

if __name__=="__main__":

    test_integrals()