derivation_fluid_equations.sagews

%auto
typeset_mode(True)

# Defining the phase space variables
var('x, v, t')

($\displaystyle x$, $\displaystyle v$, $\displaystyle t$)

# Defining the constants
var('k, m')

($\displaystyle k$, $\displaystyle m$)

# Definitions of the functions that vary with space and time:
var('v_b, T, rho')

($\displaystyle v_{b}$, $\displaystyle T$, $\displaystyle \rho$)

# Defining collision time-scale
var('tau')

$\displaystyle \tau$

# Definition of the distribution function
distribution_function_f_0 = rho * (sqrt(m/(2*pi*k*T))) * exp(-m*(v - v_b)^2/(2*k*T))
distribution_function_f_0

$\displaystyle \sqrt{\frac{1}{2}} \rho \sqrt{\frac{m}{\pi T k}} e^{\left(-\frac{m {\left(v - v_{b}\right)}^{2}}{2 \, T k}\right)}$

# Declaring the positive quantities:
assume(m>0)
assume(k>0)
assume(T>0)
assume(rho>0)

# Redefining the variables that vary with space and time:
v_bulk = function('v_{bulk}')(x, t)
T      = function('T')(x, t)
rho    = function('rho')(x, t)

# Defining the various velocity moments:
zeroth_moment = integrate(distribution_function_f_0,      v, -infinity, infinity).subs(T=T, rho=rho, v_b=v_bulk)
first_moment  = integrate(distribution_function_f_0*v,    v, -infinity, infinity).subs(T=T, rho=rho, v_b=v_bulk)
second_moment = integrate(distribution_function_f_0*v**2, v, -infinity, infinity).subs(T=T, rho=rho, v_b=v_bulk)
third_moment  = integrate(distribution_function_f_0*v**3, v, -infinity, infinity).subs(T=T, rho=rho, v_b=v_bulk)
fourth_moment = integrate(distribution_function_f_0*v**4, v, -infinity, infinity).subs(T=T, rho=rho, v_b=v_bulk)

# Simplifying the expressions
zeroth_moment = zeroth_moment.expand().simplify()
first_moment  = first_moment.expand()
second_moment = second_moment.expand()
third_moment  = third_moment.expand()
fourth_moment = fourth_moment.expand()

# Obtaining the Euler equations
euler_continuity_eqn = zeroth_moment.diff(t) + first_moment.diff(x)
euler_momentum_eqn   = first_moment.diff(t)  + second_moment.diff(x)
euler_energy_eqn     = second_moment.diff(t) + third_moment.diff(x)

# Obtaining the NS equations

# The expression for the distribution function f_approx is f = f_0 - tau * v * df_0/dx
# Due to mass, momentum and energy conservation:
# \int     f dv = \int     f_0 dv
# \int v   f dv = \int v   f_0 dv
# \int v^2 f dv = \int v^2 f_0 dv

NS_continuity_eqn = zeroth_moment.diff(t) + first_moment.diff(x)
NS_momentum_eqn   = first_moment.diff(t) + second_moment.diff(x)
NS_energy_eqn     = second_moment.diff(t)  + (third_moment - tau * fourth_moment.diff(x)).diff(x)

# Linearizing the system:
var('rho0, delta_rho, T0, delta_T, v0, delta_v, I, k1, delta_rho_dt, delta_T_dt, delta_v_dt')

($\displaystyle \rho_{0}$, $\displaystyle \delta_{\rho}$, $\displaystyle T_{0}$, $\displaystyle \delta_{T}$, $\displaystyle v_{0}$, $\displaystyle \delta_{v}$, $\displaystyle I$, $\displaystyle k_{1}$, $\displaystyle \delta_{\rho_{\mathit{dt}}}$, $\displaystyle \delta_{T_{\mathit{dt}}}$, $\displaystyle \delta_{v_{\mathit{dt}}}$)

def linearize(term):
    return taylor(term, (delta_rho, 0), \
                        (delta_v, 0),  \
                        (delta_T, 0),   \
                        (delta_rho_dt, 0), \
                        (delta_v_dt, 0),  \
                        (delta_T_dt, 0), 1 \
                 ).simplify_full()

eqn_continuity = linearize(NS_continuity_eqn.subs(rho==rho0+delta_rho, rho.diff(x)==I*k1*delta_rho, rho.diff(t)==delta_rho_dt, \
                                                  T==T0+delta_T, T.diff(x)==I*k1*delta_T, T.diff(t)==delta_T_dt, \
                                                  v_bulk==v0+delta_v, v_bulk.diff(x)==I*k1*delta_v, v_bulk.diff(t)==delta_v_dt
                                                 )
                          )
eqn_momentum   = linearize(NS_momentum_eqn.subs(rho==rho0+delta_rho, rho.diff(x)==I*k1*delta_rho, rho.diff(t)==delta_rho_dt, \
                                                T==T0+delta_T, T.diff(x)==I*k1*delta_T, T.diff(t)==delta_T_dt, \
                                                v_bulk==v0+delta_v, v_bulk.diff(x)==I*k1*delta_v, v_bulk.diff(t)==delta_v_dt
                                               )
                          )
eqn_energy     = linearize(NS_energy_eqn.subs(rho==rho0+delta_rho, rho.diff(x)==I*k1*delta_rho, rho.diff(t)==delta_rho_dt, rho.diff(x, 2)==-k1**2*delta_rho,\
                                              T==T0+delta_T, T.diff(x)==I*k1*delta_T, T.diff(t)==delta_T_dt, T.diff(x, 2)==-k1**2*delta_T,\
                                              v_bulk==v0+delta_v, v_bulk.diff(x)==I*k1*delta_v, v_bulk.diff(t)==delta_v_dt, v_bulk.diff(x, 2)==-k1**2*delta_v
                                             )
                          )

eqns          = [eqn_continuity == 0, eqn_momentum == 0, eqn_energy == 0]
delta_vars    = [delta_rho, delta_v, delta_T]
delta_vars_dt = [delta_rho_dt, delta_v_dt, delta_T_dt]

solutions = solve(eqns, delta_vars_dt, solution_dict=True)

solns_delta_vars_dt = []
for dvar_dt in delta_vars_dt:
    solns_delta_vars_dt.append(solutions[0][dvar_dt])

M = jacobian(solns_delta_vars_dt, delta_vars)
M = M.apply_map(lambda x : x.simplify_full())

pretty_print("Linearized system : ")
print("\n")
pretty_print(Matrix(delta_vars_dt).transpose(), " = ", M, Matrix(delta_vars).transpose())

Linearized system :

$\displaystyle \left(\begin{array}{r} \delta_{\rho_{\mathit{dt}}} \\ \delta_{v_{\mathit{dt}}} \\ \delta_{T_{\mathit{dt}}} \end{array}\right)$ = $\displaystyle \left(\begin{array}{rrr} -I k_{1} v_{0} & -I k_{1} \rho_{0} & 0 \\ -\frac{I T_{0} k k_{1}}{m \rho_{0}} & -I k_{1} v_{0} & -\frac{I k k_{1}}{m} \\ -\frac{k_{1}^{2} m^{2} \tau v_{0}^{4} + 6 \, T_{0} k k_{1}^{2} m \tau v_{0}^{2} + 3 \, T_{0}^{2} k^{2} k_{1}^{2} \tau}{k m \rho_{0}} & -\frac{2 \, {\left(2 \, k_{1}^{2} m \tau v_{0}^{3} + 6 \, T_{0} k k_{1}^{2} \tau v_{0} + I T_{0} k k_{1}\right)}}{k} & -\frac{6 \, k_{1}^{2} m \tau v_{0}^{2} + 6 \, T_{0} k k_{1}^{2} \tau + I k_{1} m v_{0}}{m} \end{array}\right)$ $\displaystyle \left(\begin{array}{r} \delta_{\rho} \\ \delta_{v} \\ \delta_{T} \end{array}\right)$

rho0_num = 1
v0_num   = 0
T0_num   = 1
k1_num   = 2 * pi
tau_num  = 0.01

M_numerical = M.subs(rho0=rho0_num, v0=v0_num, T0=T0_num, k1=k1_num, tau = tau_num)

M_numerical = M_numerical.change_ring(CDF)
eigenvecs   = M_numerical.eigenvectors_right()