In [1]:
from kineticlib import particles, omegaint, loaddata
from scipy import constants
from scipy.integrate import odeint
from scipy.optimize import fsolve
import numpy as np
import matplotlib.pyplot as plt
In [2]:
%matplotlib inline
In [3]:
Ar = particles.Atom('Ar')
d = loaddata.load_elastic_parameters(Ar, Ar)
In [4]:
d_step = 2.
In [5]:
def eta(T, om):
return (5. * constants.k * T) / (8. * om)
In [6]:
def lamb(T, om):
return (75. * constants.k * constants.k * T) / (32. * Ar.mass * om)
In [7]:
def sysvec(y, t):
"""
Right-hand side of the Couette ODE system in dimensionless variables
"""
x1, x2, dx1, dx2 = y
T = x2 * T_wall
omint = omegaint.omega(T, 2, 2, d)
# print(T, omint)
omdiv = (omegaint.omega(T + d_step, 2, 2, d, dimensional=True)
- omegaint.omega(T - d_step, 2, 2, d, dimensional=True)) * T_wall / (2. * d_step)
this_eta = eta(T, omint)
this_eta_deriv = (eta(T + d_step, omint) - eta(T - d_step, omint)) * T_wall / (2. * d_step)
this_lamb = lamb(T, omint)
this_lamb_deriv = (lamb(T + d_step, omint) - lamb(T - d_step, omint)) * T_wall / (2. * d_step)
rhs1 = - dx1 * dx2 * this_eta_deriv / this_eta
rhs2 = (this_eta * ((v_wall * dx1) ** 2) / T_wall + (dx2 ** 2) * this_lamb_deriv)
# print(Y * dx1, Y * dx2, rhs1, rhs2)
# return [Y * dx1, Y * dx2, Y * rhs1, -(Y / this_lamb) * rhs2]
return [dx1, dx2, rhs1, -(1. / this_lamb) * rhs2]
In [8]:
def sysvec_dimensional(y, t):
"""
Right-hand side of the Couette ODE system in dimensional variables
"""
vx, T, dvx, dT = y # velocity, temperature, derivative of velocity, derivative of temperature
omint = omegaint.omega(T, 2, 2, d)
omdiv = (omegaint.omega(T + d_step, 2, 2, d, dimensional=True)
- omegaint.omega(T - d_step, 2, 2, d, dimensional=True)) / (2. * d_step)
this_eta = eta(T, omint)
this_eta_deriv = (eta(T + d_step, omint) - eta(T - d_step, omint)) / (2. * d_step)
this_lamb = lamb(T, omint)
this_lamb_deriv = (lamb(T + d_step, omint) - lamb(T - d_step, omint)) / (2. * d_step)
rhs1 = - dT * dvx * this_eta_deriv / this_eta
rhs2 = (this_eta / this_lamb) * (dvx ** 2) + (this_lamb_deriv / this_lamb) * (dT ** 2)
return [dvx, dT, rhs1, -rhs2]
In [9]:
def give_solution_dimensionless(dx1, dx2):
"""
Returns a solution of the Couette ODE system in dimensionless variables,
given the derivatives of the dimensionless velocity and temperature
"""
return odeint(sysvec, [0., T_lower_wall / T_wall, dx1, dx2], y_span)
In [10]:
def give_solution_dimensional(dvx, dT):
"""
Returns a solution of the Couette ODE system, given the derivatives of the velocity and temperature
"""
return odeint(sysvec_dimensional, [0., T_lower_wall, dvx, dT], y_span_dimensional)
In [11]:
def dimensionless_obj(param_vec):
"""
The objective function for which we need to find the roots to solve the Couette ODE system, dimensionless version
param_vec is a list containing an approximation of the velocity and temperature derivatives
"""
res = give_solution_dimensionless(param_vec[0], param_vec[1])
return [res[-1, 0] - 1., res[-1, 1] - 1.]
In [12]:
def dimensional_obj(param_vec):
"""
The objective function for which we need to find the roots to solve the Couette ODE system, dimensional version
param_vec is a list containing an approximation of the velocity and temperature derivatives
"""
res = give_solution_dimensional(param_vec[0], param_vec[1])
return [res[-1, 0] - v_wall, res[-1, 1] - T_wall]
In [13]:
def full_solve_dimensionless():
omint = omegaint.omega(T_lower_wall, 2, 2, d)
eta_lw = eta(T_lower_wall, omint)
lamb_lw = lamb(T_lower_wall, omint)
# start_cond_T = -eta_lw * (v_wall ** 2) / (2. * lamb_lw)
start_cond_T = (T_wall - T_lower_wall) / T_wall
dvx_0, dT_0 = fsolve(dimensionless_obj, [1., start_cond_T])
return give_solution_dimensionless(dvx_0, dT_0)
In [14]:
def full_solve_dimensional():
dvx_0, dT_0 = fsolve(dimensional_obj, [v_wall / Y, (T_wall - T_lower_wall) / Y])
return give_solution_dimensional(dvx_0, dT_0)
In [15]:
def convert_to_dimensional(sol_vx, sol_T):
return [sol_vx * v_wall, sol_T * T_wall]
In [16]:
v_wall = 700. # the velocity of the moving wall
T_wall = 1000. # the temperature of the moving wall
T_lower_wall = 1000. # the temperature of the stationary wall
Y = 1. # the distance between the walls
In [17]:
points = 100.
points_per_meter = 100
y_span = np.linspace(0, 1., points)
# y_span_dimensional = np.linspace(0., Y, points_per_meter * Y)
In [18]:
sol_dimless = full_solve_dimensionless()
s2 = convert_to_dimensional(sol_dimless[:,0], sol_dimless[:,1])
In [19]:
# sol_dimensional = full_solve_dimensional()
In [20]:
s2[0].shape
Out[20]:
In [24]:
plt.rcParams['figure.figsize'] = (15.0, 15.0)
plt.plot(y_span, (sol_dimless[:,0] - y_span) * v_wall, 'r')
# plt.plot(y_span, y_span, 'b--')
plt.xlabel(r'$y / Y$', fontsize=38)
plt.ylabel(r'$v - \frac{y}{Y} v_{wall}$', fontsize=38)
plt.tick_params(labelsize=20)
# plt.plot(y_span_dimensional, sol_dimensional[:, 0], 'k*')
In [26]:
plt.rcParams['figure.figsize'] = (15.0, 15.0)
plt.plot(y_span, sol_dimless[:,1], 'r')
# plt.plot(y_span, y_span, 'b--')
plt.xlabel(r'$y / Y$', fontsize=38)
plt.ylabel(r'$T / T_{wall,2}$', fontsize=38)
plt.tick_params(labelsize=20)
In [27]:
eta_vals = np.zeros(sol_dimless[:,1].shape[0])
for j, i in enumerate(sol_dimless[:,1]):
T = i * T_wall
omint = omegaint.omega(T, 2, 2, d)
eta_vals[j] = eta(T, omint)
In [29]:
plt.rcParams['figure.figsize'] = (15.0, 15.0)
plt.plot(y_span, eta_vals * (10. ** 6), 'r')
# plt.plot(y_span, y_span, 'b--')
plt.xlabel(r'$y / Y$', fontsize=38)
plt.ylabel(r'$\eta$, $\mu$Pa$\cdot$s', fontsize=38)
plt.tick_params(labelsize=20)
In [ ]: