雷诺方程推导及FDM求解

https://en.wikipedia.org/wiki/Reynolds_equation
https://www.tribonet.org/wiki/reynolds-equation/
Fluid Film Lubrication.  Andras Z. Szeri
https://baike.baidu.com/item/高阶无穷小
https://baike.baidu.com/item/雷诺方程解法
http://www.cfd2030.com/index.html

https://www.youtube.com/watch?v=QE9lDqZ9kZg

雷诺方程Reynolds equation

粘性流体的运动方程首先由纳维在1827年提出，只考虑了不可压缩流体的流动。泊松在1831年提出可压缩流体的运动方程。圣维南与斯托克斯在1845年独立提出粘性系数为一常数的形式。
雷诺方程由雷诺1886年提出，在润滑理论种描述粘性液膜（气膜）的压力分布，不要将雷诺方程与雷诺数以及雷诺平均NS方程搞混。（雷诺平均NS方程即RANS，是流体力学中一种用来描述湍流的时均NS方程，其思想是将湍流运动看作时间平均与瞬时脉动两种流动的叠加；雷诺数是1908以雷诺命名的，是液体从层流到湍流的转变的指标，是惯性力与粘性力的比值。）。
雷诺方程由连续性方程与NS方程导出。

定常：不随时间变化

NS方程的简化形式：
雷诺数趋于无穷，无粘性，简化为Euler流
雷诺数趋于0，无惯性，简化为Stokes流
仅在体积力中考虑密度的变化，简化为Boussinesq流
边界层近似用于固体边界附近和高雷诺数

高阶无穷小：变量的高阶无穷小与常量的高阶无穷小，如一个常量是$O(10^{-3})$，那么这个数比$10^{-3}$的量级还小。

journal bearing：径向滑动轴承
slope：斜率
induced pressure
no-slip boundary conditions：流体的壁面速度为0，一般大尺度的流场分析中，壁面都属于无滑移边界。
Scale analysis (or order-of-magnitude analysis)：尺度分析（量级分析）

推导

润滑膜$ϵ=(Ly/Lxz)$，量级比$10^{-3}$还小。

雷诺学习 Beauchamp Tower关于铁路轴承的报告，想到了：“the film of oil might be sufficiently thick for the unknown boundary actions to disappear, in which case the results would be deducible from the equations of hydrodynamics.”

雷诺推导方程时的假设：
满足连续性方程
满足NS方程
不可压缩（密度为常数）
粘性为常数
润滑流无涡
润滑流惯性忽略（相对于粘性力与压力，油膜的体积力太小而忽略）
$ϵR_e<1$

在径向滑动轴承中，雷诺数为1000左右，依然可以是层流。

$u,v,w$为$x,y,z$方向的速度

$\frac{\partial p}{\partial x}=\mu \frac{{\partial }^{2}u}{\partial y^2}$
$\frac{\partial p}{\partial z}=\mu \frac{{\partial }^{2}w}{\partial y^2}$
$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0$
沿着y方向的压力为常数

将前两个式子沿着$y$积两次分，$u=\frac{1}{2\mu }\frac{\partial p}{\partial x}{y}^2+Ay+B$，$w=\frac{1}{2\mu }\frac{\partial p}{\partial z}{y}^2+Cy+D$
边界条件，$y=0,w=0,u=U_1$，$y=h,w=0,u=U_2$
将边界条件带入积分后的式子，求出常数
$u=\frac{1}{2\mu }\frac{\partial p}{\partial x}(y^2-yh)+(1-\frac{y}{h})U_1+\frac{y}{h}{U}_2$
$w=\frac{1}{2\mu }\frac{\partial p}{\partial z}(y^2-yh)$

至此$u,w$中还有一个未知数，即压力梯度，即$p$未知。
确保质量守恒，将方程带入连续方程，那么未知数为$p,v$。
一个方程，两个未知数，不可解。但实际上

将连续方程沿着y方向积分
${\int }_0^h-\frac{\partial v}{\partial y}dy={\int }_0^h\frac{\partial u}{\partial x}dy+{\int }_0^h\frac{\partial w}{\partial z}dy$
$-v{|}_0^h={\int }_0^h\frac{\partial u}{\partial x}dy+{\int }_0^h\frac{\partial w}{\partial z}dy$
将$u,w$带入得出
$\frac{\partial }{\partial x}(\frac{h^3}{\mu }\frac{\partial p}{\partial x})+\frac{\partial }{\partial z}(\frac{h^3}{\mu }\frac{\partial p}{\partial z})=6(U_1+U_2)\frac{\partial h}{\partial x}+12\frac{\partial h}{\partial t}$

h不随时间变化，且$U_1=0$
$\frac{\partial }{\partial x}(\frac{h^3}{\mu }\frac{\partial p}{\partial x})+\frac{\partial }{\partial z}(\frac{h^3}{\mu }\frac{\partial p}{\partial z})=6U\frac{\partial h}{\partial x}$

求解

下面写的是另一种形式由youtube一个博士推导的，最终化简成的形式为稳态雷诺方程与上面的方程的最终推导结果一样，但我不清楚我之前是怎么得到这个公式的。
总之，你可以略过求解部分的公式推导，直接看代码。
通过FDM FVM FEM等方法可求解，在此运用FDM求解稳态雷诺方程。
当然刚开始写的那么复杂，最终求解的也只是一个化简后的方程。

$\frac{\partial }{\partial x}(\frac{h^3}{12\mu }\frac{\partial p}{\partial x})+\frac{\partial }{\partial z}(\frac{h^3}{12\mu }\frac{\partial p}{\partial z})=\frac{1}{2}\frac{\partial (U_2-U_1)h}{\partial x}+\frac{1}{2}\frac{\partial (W_2-W_1)h}{\partial z}+(V_1-V_2)$

等式右端为Source term
$V_1-V_2=\frac{\partial h}{\partial t}$，squeeze action
$\frac{\partial h}{\partial x},\frac{\partial h}{\partial z}$，wedge action
$\frac{\partial U}{\partial x},\frac{\partial W}{\partial z}$，stretching action

假设：
忽略惯性
忽略油膜厚度方向的压力梯度
牛顿流体（任一点上的剪应力都同剪切变形速率呈线性函数关系的流体称为牛顿流体。）
常粘性
无滑移边界
忽略坐标系的倾斜角
不可压流
沿z方向速度不变
上下面为刚性面
只有倾斜面移动

$\frac{\partial }{\partial x}(h^3\frac{\partial p}{\partial x})+\frac{\partial }{\partial z}(h^3\frac{\partial p}{\partial z})=6\mu \frac{\partial (U_2-U_1)h}{\partial x}+6\mu \frac{\partial (W_2-W_1)h}{\partial z}+12\mu \frac{\partial h}{\partial t}$

$\frac{\partial }{\partial x}(h^3\frac{\partial p}{\partial x})+\frac{\partial }{\partial z}(h^3\frac{\partial p}{\partial z})=6\mu \frac{\partial (U_2-U_1)h}{\partial x}+12\mu \frac{\partial h}{\partial t}$

$\frac{\partial }{\partial x}(h^3\frac{\partial p}{\partial x})+\frac{\partial }{\partial z}(h^3\frac{\partial p}{\partial z})=6\mu (U_2-U_1)\frac{\partial h}{\partial x}+12\mu \frac{\partial h}{\partial t}$

$\frac{\partial }{\partial x}(h^3\frac{\partial p}{\partial x})+\frac{\partial }{\partial z}(h^3\frac{\partial p}{\partial z})=6\mu U\frac{\partial h}{\partial x}+12\mu \frac{\partial h}{\partial t}$

进行无量纲化
$\bar{x}=x/X,\bar{z}=z/Z,\bar{h}=h/C$
$\bar{p}=\frac{C^{3}p}{6\eta U{X}^2},\bar{t}=\frac{tU}{C}$

$\frac{\partial }{\partial x}({\bar{h}}^3\frac{\partial \bar{p}}{\partial \bar{x}})+\frac{X^2}{Z^2}\frac{\partial }{\partial \bar{z}}({\bar{h}}^3\frac{\partial \bar{p}}{\partial \bar{z}})=\frac{C}{X}\frac{\partial \bar{h}}{\partial \bar{x}}+2\frac{\bar{h}}{\partial \bar{t}}$

如果$X=10Z$，等式左端第一项可忽略
如果$X=0.1Z$，等式左端第二项可忽略

稳态雷诺方程
$\frac{\partial }{\partial x}({\bar{h}}^3\frac{\partial \bar{p}}{\partial \bar{x}})+\frac{X^2}{Z^2}\frac{\partial }{\partial \bar{z}}({\bar{h}}^3\frac{\partial \bar{p}}{\partial \bar{z}})=\frac{C}{X}\frac{\partial \bar{h}}{\partial \bar{x}}$

变量有$p,h$，以压力为例
节点值
${\bar{p}}_{i,j}$
一阶偏导
$\frac{\partial {\bar{p}}_{i,j}}{\partial \bar{x}}=\frac{{\bar{p}}_{i+1,j}-{\bar{p}}_{i-1,j}}{2\Delta \bar{x}}$
二阶偏导
$\frac{\partial {\bar{p}}_{i,j}^2}{\partial {\bar{x}}^2}={\frac{\partial \bar{p}}{\partial \bar{x}}}_{i+0.5,j}-{\frac{\partial \bar{p}}{\partial \bar{x}}}_{i-0.5,j}$
二阶偏导取的是一阶偏导计算完成后中点上的值，这是为了保证精度。

将高度，压力梯度，高度梯度都用差分表示，带入方程中，可得到差分方程，编程求解。

# -*- coding: utf-8 -*-
"""
@time: 2021-08-29 下午 05:19
@author: leslie lee
"""
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm


# ----------------------------------------
# ****************FDM*********************
# ----------------------------------------
# 单位 mm
# x方向划分网格数
N = 50
# x方向网格间隔
delxbar = 1/N
# z方向划分网格数
M = 50
# z方向网格间隔
delzbar = 1/M
# 径向间隔
c = 0.03
# x方向的尺寸
X1 = 100
# z方向的尺寸
Z1 = 25
# 常数1
const1 = X1**2/Z1**2
# 迭代次数
ITER = 1000
# 预定义矩阵
p = np.zeros([N,M])
# 对p求和
current_sum = np.sum(p)

# 迭代
for K in range(ITER):
    # 遍历节点x方向的索引
    for I in range(1,N-1):
        xbar = delxbar*I # xbar
        h = 2/3*(2 - xbar) # h_{i,j} hbar开始为4/3 末端为2/3  hbar沿xbar线性变化 xbar总长为1; 由径向间隔可以推出 h开始为0.04mm 末端为0.02mm
        hm = 2/3*(2 - xbar - 0.5*delxbar) # h_{i+0.5,j}
        hp = 2/3*(2 - xbar + 0.5*delxbar) # h_{i+0.5,j}
        hm1 = 2/3*(2 - xbar - delxbar) # h_{i-1,j}
        hp1 = 2/3*(2 - xbar + delxbar) # h_{i+1,j}
        cubh = h**3
        cubhm = hm**3
        cubhp = hp**3
        # 常数2
        const2 = cubhp + cubhm + 2*const1*cubh
        A = const1*cubh/const2
        B = A
        C = cubhp/const2
        D = cubhm/const2
        E = -0.5*delxbar*(hp1-hm1)*c/(const2*X1)
        # 遍历节点z方向的索引
        for J in range(1,M-1):
            p[I,J] = A*p[I,J+1]+B*p[I,J-1]+C*p[I+1,J]+D*p[I-1,J]-E
    # 如果满足条件 则停止迭代
    before_sum = np.abs(current_sum)
    current_sum = np.abs(np.sum(p))
    percentage = (current_sum - before_sum)/current_sum
    if percentage < 0.0001:
        break

# ----------------------------------------
# ****************绘图*********************
# ----------------------------------------
"""
https://www.cnblogs.com/wuwen19940508/p/8638266.html
这种写法已经过时
fig = plt.figure()
ax = fig.gca(projection='3d')
"""
ax = plt.axes(projection='3d')
# Plot the surface.
x = np.arange(0,1,delxbar)
z = np.arange(0,1,delzbar)
x, z = np.meshgrid(x,z)
surf = ax.plot_surface(x, z, p, cmap=plt.get_cmap('rainbow'))
plt.show()

之前的粘性为常数，考虑粘性随压力变化，$\eta ={\eta }_0{e}^{ap}$，再求解一次。