當前位置：首頁 > 编程资源 > 编程问答 >内容正文

编程问答

雷诺方程推导及FDM求解

發布時間：2023/12/20 编程问答 39 豆豆

生活随笔收集整理的這篇文章主要介紹了雷诺方程推导及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.htmlhttps://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)\epsilon = (Ly/Lxz)$ ，量級比 $10^{-3}$ 還小。

雷諾學習 Beauchamp Tower關于鐵路軸承的報告，想到了：“the ?lm of oil might be suf?ciently thick for the unknown boundary actions to disappear, in which case the results would be deducible from the equations of hydrodynamics.”

雷諾推導方程時的假設：
滿足連續性方程
滿足NS方程
不可壓縮（密度為常數）
粘性為常數
潤滑流無渦
潤滑流慣性忽略（相對于粘性力與壓力，油膜的體積力太小而忽略）
$?Re<1\epsilon R_e <1$

在徑向滑動軸承中，雷諾數為1000左右，依然可以是層流。

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

$?p?x=μ?2u?y2\frac{\partial p}{\partial x} = \mu \frac{\partial^2 u}{\partial y^2}$
$?p?z=μ?2w?y2\frac{\partial p}{\partial z} = \mu \frac{\partial^2 w}{\partial y^2}$
$?u?x+?v?y+?w?z=0\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} = 0$
沿著y方向的壓力為常數

將前兩個式子沿著 $y$ 積兩次分， $u=12μ?p?xy2+Ay+Bu=\frac{1}{2\mu}\frac{\partial p}{\partial x}y^2+Ay+B$ ， $w=12μ?p?zy2+Cy+Dw=\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=12μ?p?x(y2?yh)+(1?yh)U1+yhU2u=\frac{1}{2\mu}\frac{\partial p}{\partial x}(y^2-yh)+(1-\frac{y}{h})U_1+\frac{y}{h}U_2$
$w=12μ?p?z(y2?yh)w=\frac{1}{2\mu}\frac{\partial p}{\partial z}(y^2-yh)$

至此 $u, w$ 中還有一個未知數，即壓力梯度，即 $p$ 未知。
確保質量守恒，將方程帶入連續方程，那么未知數為 $p, v$ 。
一個方程，兩個未知數，不可解。但實際上

將連續方程沿著y方向積分
$∫0h??v?ydy=∫0h?u?xdy+∫0h?w?zdy\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∣0h=∫0h?u?xdy+∫0h?w?zdy-v|_0^h = \int_0^h \frac{\partial u}{\partial x} dy + \int_0^h \frac{\partial w}{\partial z} dy$
將 $u, w$ 帶入得出
$??x(h3μ?p?x)+??z(h3μ?p?z)=6(U1+U2)?h?x+12?h?t\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$
$??x(h3μ?p?x)+??z(h3μ?p?z)=6U?h?x\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求解穩態雷諾方程。
當然剛開始寫的那么復雜，最終求解的也只是一個化簡后的方程。

$??x(h312μ?p?x)+??z(h312μ?p?z)=12?(U2?U1)h?x+12?(W2?W1)h?z+(V1?V2)\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
$V1?V2=?h?tV_1 - V_2 = \frac{\partial h}{\partial t}$ ，squeeze action
$?h?x,?h?z\frac{\partial h}{\partial x}, \frac{\partial h}{\partial z}$ ，wedge action
$?U?x,?W?z\frac{\partial U}{\partial x}, \frac{\partial W}{\partial z}$ ，stretching action

假設：
忽略慣性
忽略油膜厚度方向的壓力梯度
牛頓流體（任一點上的剪應力都同剪切變形速率呈線性函數關系的流體稱為牛頓流體。）
常粘性
無滑移邊界
忽略坐標系的傾斜角
不可壓流
沿z方向速度不變
上下面為剛性面
只有傾斜面移動

$??x(h3?p?x)+??z(h3?p?z)=6μ?(U2?U1)h?x+6μ?(W2?W1)h?z+12μ?h?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} + 6\mu \frac{\partial (W_2 - W_1)h}{\partial z} + 12 \mu \frac{\partial h}{\partial t}$

$??x(h3?p?x)+??z(h3?p?z)=6μ?(U2?U1)h?x+12μ?h?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}$

$??x(h3?p?x)+??z(h3?p?z)=6μ(U2?U1)?h?x+12μ?h?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}$

$??x(h3?p?x)+??z(h3?p?z)=6μU?h?x+12μ?h?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}$

進行無量綱化
$xˉ=x/X,zˉ=z/Z,hˉ=h/C\bar{x}=x/X, \bar{z}=z/Z, \bar{h}=h/C$
$pˉ=C3p6ηUX2,tˉ=tUC\bar{p}=\frac{C^3p}{6\eta UX^2}, \bar{t}=\frac{tU}{C}$

$??x(hˉ3?pˉ?xˉ)+X2Z2??zˉ(hˉ3?pˉ?zˉ)=CX?hˉ?xˉ+2hˉ?tˉ\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 = 10 Z$ ，等式左端第一項可忽略
如果 $X = 0.1 Z$ ，等式左端第二項可忽略

穩態雷諾方程
$??x(hˉ3?pˉ?xˉ)+X2Z2??zˉ(hˉ3?pˉ?zˉ)=CX?hˉ?xˉ\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$ ，以壓力為例
節點值
$pˉi,j\bar{p}_{i,j}$
一階偏導
$?pˉi,j?xˉ=pˉi+1,j?pˉi?1,j2Δxˉ\frac{\partial \bar{p}_{i,j}}{\partial \bar{x}} = \frac{\bar{p}_{i+1,j} - \bar{p}_{i-1,j}}{2 \Delta \bar{x}}$
二階偏導
$?pˉi,j2?xˉ2=?pˉ?xˉi+0.5,j??pˉ?xˉi?0.5,j\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 # xbarh = 2/3*(2 - xbar) # h_{i,j} hbar開始為4/3 末端為2/3 hbar沿xbar線性變化 xbar總長為1; 由徑向間隔可以推出 h開始為0.04mm 末端為0.02mmhm = 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**3cubhm = hm**3cubhp = hp**3# 常數2const2 = cubhp + cubhm + 2*const1*cubhA = const1*cubh/const2B = AC = cubhp/const2D = cubhm/const2E = -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_sumif 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()

之前的粘性為常數，考慮粘性隨壓力變化， $η=η0eap\eta = \eta_0 e^{ap}$ ，再求解一次。

總結

以上是生活随笔為你收集整理的雷诺方程推导及FDM求解的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。

上一篇： Fantasia （Tarjan+树
下一篇： TEXT宏