In [1]:
# A simple solver of the Couette flow of pure atomic argon for slip conditions (with no temperature jump)
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)
omint_p = omegaint.omega(T + d_step, 2, 2, d)
omint_m = omegaint.omega(T - d_step, 2, 2, d)
this_eta = eta(T, omint)
this_eta_deriv = (eta(T + d_step, omint_p) - eta(T - d_step, omint_m)) * T_wall / (2. * d_step)
this_lamb = lamb(T, omint)
this_lamb_deriv = (lamb(T + d_step, omint_p) - lamb(T - d_step, omint_m)) * 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)
return [dx1, dx2, rhs1, -(1. / this_lamb) * rhs2]
In [8]:
def give_solution_dimensionless(dx1, dx2, lower_wall_vel):
"""
Returns a solution of the Couette ODE system in dimensionless variables,
given the derivatives of the dimensionless velocity and temperature
"""
return odeint(sysvec, [lower_wall_vel, T_lower_wall / T_wall, dx1, dx2], y_span)
In [9]:
def dimensionless_obj(param_vec, *args):
"""
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
"""
upper_wal_vel, lower_wall_vel = args
res = give_solution_dimensionless(param_vec[0], param_vec[1], lower_wall_vel)
return [res[-1, 0] - upper_wal_vel, res[-1, 1] - 1.]
In [10]:
def full_solve_dimensionless(upper_wal_vel=1., lower_wall_vel=0.):
"""
return a solution for given values of the upper and lower wall velocities
"""
start_cond_T = (T_wall - T_lower_wall) / T_wall
dvx_0, dT_0 = fsolve(dimensionless_obj, [1., start_cond_T], (upper_wal_vel, lower_wall_vel, ))
return give_solution_dimensionless(dvx_0, dT_0, lower_wall_vel)
In [11]:
def simple_slip(loc_T, loc_eta, loc_dv):
"""
Calculate the value of the slip velocity given local temperature (loc_T), local value of viscosity (loc_eta),
and local value of velocity gradient (loc_dv)
"""
loc_v = (loc_T * 2 * constants.k / Ar.mass) ** 0.5
# print(loc_v * sigma_p * (loc_eta * loc_dv / (n * constants.k * loc_T)),
# sigma_t * (loc_eta * loc_dT / (loc_T * rho)))
return sigma_p * loc_v * (loc_eta * loc_dv / (n * constants.k * loc_T))
# return sigma_p * loc_v * (loc_eta * loc_dv / (n * constants.k * loc_T)) + sigma_t * (loc_eta * loc_dT / (loc_T * rho))
In [12]:
def full_solve_dimensionless_with_slip():
counter = 0
this_upper_v = 1.
this_lower_v = 0.
err = 1.
while err > tol and counter < N_steps_max:
counter += 1
tmp_sol = full_solve_dimensionless(this_upper_v, this_lower_v)
T_l, T_u = tmp_sol[0, 1] * T_wall, tmp_sol[-1, 1] * T_wall
om_l, om_u = omegaint.omega(T_l, 2, 2, d, dimensional=True,
model=om_mod), omegaint.omega(T_u, 2, 2, d, dimensional=True, model=om_mod)
tmp_eta_l, tmp_eta_u = eta(T_l, om_l), eta(T_u, om_u)
this_lower_v = simple_slip(T_l, tmp_eta_l, tmp_sol[0, 2] * v_wall) / v_wall
this_upper_v = 1. - simple_slip(T_u, tmp_eta_u, tmp_sol[-1, 2] * v_wall) / v_wall
err = abs((this_lower_v - tmp_sol[0, 0]) / this_lower_v) + abs((this_upper_v - tmp_sol[-1, 0]) / this_upper_v)
print('Step:', counter, 'error:', err, this_lower_v, this_upper_v)
return tmp_sol
In [13]:
def convert_to_dimensional(sol_vx, sol_T):
return [sol_vx * v_wall, sol_T * T_wall]
In [65]:
v_wall = 300. # the velocity of the moving wall
T_wall = 273. # the temperature of the moving wall
T_lower_wall = 273. # the temperature of the stationary wall
Y = 1. # the distance between the walls
om_mod = 'VSS' # the elastic collision integral model
# sigma_p = 1.18 # experimental value
sigma_p = 0.8862 # basic fully diffuse scattering
# sigma_t = 0.75 # based on Ref.14, Felix Sharipov (Data on Velocity Slip)
tol = 0.00005
N_steps_max = 50
n = 1.4 * (10. ** 20)
rho = n * Ar.mass
In [66]:
points = 100.
y_span = np.linspace(0, 1., points)
In [67]:
sol_dimless = full_solve_dimensionless()
s2 = convert_to_dimensional(sol_dimless[:,0], sol_dimless[:,1])
In [68]:
sol_dimless.shape
Out[68]:
In [69]:
sol_slip_dimless = full_solve_dimensionless_with_slip()
s3 = convert_to_dimensional(sol_slip_dimless[:,0], sol_slip_dimless[:,1])
In [70]:
plt.rcParams['figure.figsize'] = (12.0, 12.0)
plt.plot(y_span, (sol_dimless[:,0] - y_span) * v_wall, 'r')
plt.plot(y_span, (sol_slip_dimless[:, 0] - y_span) * v_wall, 'b')
plt.xlabel(r'$y / Y$', fontsize=38)
plt.ylabel(r'$v - \frac{y}{Y} v_{wall}$', fontsize=38)
plt.legend(['No slip', 'Slip'], loc='upper right', fontsize=20)
plt.tick_params(labelsize=20)
# plt.plot(y_span_dimensional, sol_dimensional[:, 0], 'k*')
In [71]:
plt.rcParams['figure.figsize'] = (12.0, 12.0)
plt.plot(y_span, sol_dimless[:,1], 'r')
plt.plot(y_span, sol_slip_dimless[:,1], 'b')
plt.xlabel(r'$y / Y$', fontsize=38)
plt.ylabel(r'$T / T_{wall,2}$', fontsize=38)
plt.legend(['No slip', 'Slip'], loc='upper right', fontsize=20)
plt.tick_params(labelsize=20)
In [72]:
eta_vals = np.zeros(sol_dimless[:,1].shape[0])
eta_vals_slip = 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)
T = sol_slip_dimless[j, 1] * T_wall
omint = omegaint.omega(T, 2, 2, d)
eta_vals_slip[j] = eta(T, omint)
In [73]:
plt.rcParams['figure.figsize'] = (12.0, 12.0)
plt.plot(y_span, eta_vals * (10. ** 6), 'r')
plt.plot(y_span, eta_vals_slip * (10. ** 6), 'b')
# 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.legend(['No slip', 'Slip'], loc='upper right', fontsize=20)
plt.tick_params(labelsize=20)
In [ ]: