【CFD理论】扩散项-01

咔咔酷酷冲冲

已于 2022-08-14 18:44:20 修改

阅读量865

点赞数 2

分类专栏： of初阶文章标签：算法机器学习人工智能

于 2022-07-31 23:38:26 首次发布

本文链接：https://blog.csdn.net/zq93538196/article/details/126065123

版权

of初阶专栏收录该内容

16 篇文章 7 订阅

订阅专栏

【CFD理论】扩散项-01

分别描述扩散项和对流项

$\nabla \cdot \boldsymbol U=0 \\ \frac{\partial \boldsymbol U}{\partial t}+\nabla \cdot (\boldsymbol U \boldsymbol U)=-\frac{\nabla p}{\rho}+\nabla \cdot (\nu \nabla \boldsymbol U)+S$
二维扩散方程
$\begin{aligned} & \nabla \cdot (D\nabla\phi)+Q=0\\ & D(\nabla ^2 \phi)\boldsymbol S=D( \Delta \phi)+Q=0\\ & D(\frac{\partial ^2 \phi}{\partial x^2}+\frac{\partial ^2 \phi}{\partial y^2})+Q=0\\ & \lmoustache_V\nabla \cdot(D\nabla \phi)dV+ \lmoustache_V Q dV=0\\ & \sum_f\lmoustache_f(D\nabla \phi)\cdot d\boldsymbol S+Q_pV_p=0\\ & \sum_f(D\nabla \phi)\cdot \boldsymbol S+Q_pV_p=0 \end{aligned}$
$\begin{aligned} & (\nabla \phi)_f:\\ & 1.\frac{\phi_N-\phi_P}{\Delta PN}\\ & 2.^a(\nabla \phi)_P+^b(\nabla\phi)_N\\ & a=(1-\gamma),b=\gamma\\ \ & \gamma=\frac{pf}{pn}=\frac{V_P}{V_P+V_N} \end{aligned}$

在这里插入图片描述
推导：

下游减上游，w p e从上游到下游
结果为neighbor-owner
$\begin{aligned} & e \\ & D(\nabla\phi)_e\cdot \boldsymbol S_e\\ & =D(\frac{\partial \phi}{\partial x}\boldsymbol i+\frac{\partial \phi}{\partial y}\boldsymbol j)\cdot(\Delta y\boldsymbol i)\\ & =\Delta y D(\frac{\partial \phi}{\partial x})_p\\ & =\Delta y D \frac{\phi_E-\phi_P}{\delta x_{PE}}\\ &=a_E(\phi_E-\phi_P)\\ & w\\ & D(\nabla\phi)_w\cdot \boldsymbol S_w\\ & =D(\frac{\partial \phi}{\partial x}\boldsymbol i+\frac{\partial \phi}{\partial y}\boldsymbol j)\cdot(-\Delta y\boldsymbol i)\\ & =-\Delta y D(\frac{\partial \phi}{\partial x})_p\\ & =-\Delta y D \frac{\phi_P-\phi_W}{\delta x_{PW}}\\ &=a_W(\phi_W-\phi_P)\\ & n \\ & D(\nabla\phi)_n\cdot \boldsymbol S_n\\ & =D(\frac{\partial \phi}{\partial x}\boldsymbol i+\frac{\partial \phi}{\partial y}\boldsymbol j)\cdot(\Delta y\boldsymbol j)\\ & =\Delta x D(\frac{\partial \phi}{\partial y})_p\\ & =\Delta x D \frac{\phi_N-\phi_P}{\delta y_{PN}}\\ &=a_N(\phi_N-\phi_P)\\ & s\\ & D(\nabla\phi)_s\cdot \boldsymbol S_s\\ & =D(\frac{\partial \phi}{\partial x}\boldsymbol i+\frac{\partial \phi}{\partial y}\boldsymbol j)\cdot(-\Delta y\boldsymbol j)\\ & =-\Delta x D(\frac{\partial \phi}{\partial y})_p\\ & =-\Delta x D \frac{\phi_P-\phi_S}{\delta y_{PS}}\\ &=a_S(\phi_P-\phi_S)\\ \\ & \sum_f(D\nabla \phi)\cdot \boldsymbol S+Q_pV_p=0\\ & Q_pV_P=Q_u+Q_s\phi_p\\ & a_P=-(a_S+a_N+a_W+a_E)+Q_s\\ & a_P\phi_P=a_E\phi_E+a_W\phi_W+a_S\phi_S+a_N\phi_N+Q_u\\ \end{aligned}$
对角占优矩阵， $\left|a_{ii}\right|\ge \sum_{j\neq i} a_{ij}, for \ all \ i$
注意边界条件（内外网格）

基本边界条件

Dirichlet

$\begin{aligned} & b\\ & D \cdot (\nabla \phi)_b \cdot \boldsymbol S_b\\ & = D\frac{\phi_b-\phi_p}{\delta x_{pb}}\cdot \Delta y\\ & =a_b\phi_b-a_p\phi_p\\ & a_p\phi_p+a_N\phi_N+a_W\phi_W+a_S\phi_S+a_b\phi_b=0\\ & a_p\gt \left| a_N\right|+ \left| a_W\right|+ \left| a_S\right|\\ & if \ D \cdot (\nabla \phi)_b \cdot \boldsymbol {\hat{i}}=q\\ & D \cdot (\nabla \phi)_b \cdot \boldsymbol S_b=q\Delta y\\ & a_p\phi_p+a_N\phi_N+a_W\phi_W+a_S\phi_S+q\Delta y=0\\ \end{aligned}$
$zeroGradient\Rightarrow \phi_p=\phi_b$

Neumann

Robin

symmetry

interface diffusivity

$D_e=(1-\gamma_e)D_p+\gamma_eD_E\\ \gamma_e=\frac{d_{pe}}{d_{pe}+d_{pE}}$

example

二维热传导控制方程
在这里插入图片描述
$\begin{aligned} &a_1T_1&+a_2T_2 & &+a_4T_4&&&&&=b_1\\ &a_1T_1&+a_2T_2 &+a_3T_3 &&+a_5T_5&&&&=b_2\\ &&&+a_3T_3 &+a_4T_4&&+a_6T_6&&&=b_3\\ &a_1T_1&&&+a_4T_4&+a_5T_5&&+a_7T_7&&=b_4\\ &&a_2T_2&&+a_4T_4&+a_5T_5&a_6T_6&&+a_8T_8&=b_5\\ &&&a_3T3&&+a_5T_5&&&&+a_9T_9=b_6\\ &&&&a_4T_4&&&+a_7T_7&+a_8T8&=b_7\\ &&&&&a_5T_5&&+a_7T_7&+a_8T8&+a_9T_9=b_8\\ &&&&&&a_6T_6&&+a_8T_8&+a_9T_9=b_9\\ \end{aligned}$
node坐标
$\begin{aligned} & x_{21} = x_{20}=x_{19}=0 & y_{10} = y_{11}=y_{12}=0\\ & x_{10}=x_{1}=x_7=x_{17}=0.05 & y_{21}=y_1=y_2=y_3=y_{13}=0.05\\ & x_{11}=x_2=x_5=x_8=x_{17}=0.2 & y_{20}=y_4=y_5=y_6=y_{14}=0.15\\ & x_{12}=x_3=x_6=x_{16}=0.45 &y_{19}=y_7=y_8=y_9=y_{15}=0.25\\ & x_{13}=x_{14}=x_{15}=0.6 & y_{18}=y_{17}=y_{16}=0.3 \end{aligned}$
node之间的距离
$\begin{aligned} & \delta x_{21-1}=\delta x_{20-4}=\delta x_{19-7}=0.05 &\delta y_{10-1}=\delta y_{11-2}=\delta y_{12-3}=0.05\\ & \delta x_{1-2}=\delta x_{4-5}=\delta x_{7-8}=0.15 &\delta y_{1-4}=\delta y_{2-5}=\delta y_{3-6}=0.1\\ & \delta x_{2-3}=\delta x_{5-6}=\delta x_{8-9}=0.25 &\delta y_{4-7}=\delta y_{5-8}=\delta y_{6-9}=0.1\\ & \delta x_{3-13}=\delta x_{6-14}=\delta x_{9-15}=0.15 &\delta y_{7-18}=\delta y_{7-18}=\delta y_{9-16}=0.05\\ \end{aligned}$
element的大小
$\begin{aligned} & \Delta x_1=\Delta x_4=\Delta x_7=0.1 &\Delta y_1=\Delta y_4=\Delta y_7=0.1\\ & \Delta x_2=\Delta x_5=\Delta x_8=0.2 &\Delta y_2=\Delta y_5=\Delta y_8=0.1\\ & \Delta x_3=\Delta x_6=\Delta x_9=0.3 &\Delta y_3=\Delta y_6=\Delta y_9=0.1\\ \end{aligned}$
cell体积
$\begin{aligned} & V_1=V_4=V_7=0.01\\ & V_2=V_5=V_8=0.02\\ & V_3=V_6=V_9=0.03\\ \end{aligned}$
插值因子
$\begin{aligned} & (ge)_1=(ge)_4=(ge)_7=\frac{V_1}{V_1+V_2}=0.333\\ & (ge)_2=(ge)_5=(ge)_8=\frac{V_2}{V_2+V_3}=0.4\\ & (gn)_1=(gn)_2=(gn)_3=\frac{V_1}{V_1+V_4}=0.5\\ & (gn)_4=(gn)_5=(gn)_6=\frac{V_4}{V_{4}+V_{7}}=0.5\\ \end{aligned}$
热传导系数值
$\begin{aligned} &k_1=k_4=k_7=k_{10}=k_{21}=k_{19}=k_{18}=10^{-3}\\ &k_2=k_3=k_5=k_{6}=k_{8}=k_{9}=k_{11}=k_{12}=k_{13}=k_{14}=k_{15}=k_{16}=k_{17}=10^{2}\\ &face\\ & k_{1-2}=\frac{k_1k_2}{[1-(ge)_1]k_1+(ge)_1k_2}\approx 3\times 10^{-3}\\ & k_{1-2}=k_{4-5}=k_{7-8}=3\times 10^{-3}\\ & k_{1-4}=k_{4-7}=10^{-3}\\ & k_{2-5}=k_{3-6}=k_{6-9}=10^{2}\\ \end{aligned}$
element 1
$\begin{aligned} &e:\frac{\Delta y_1}{\delta x_{1-2}}=\frac{0.1}{0.15}=0.667 ,&\ k_e=k_{1-2}=3\times 10^{-3},& a_2=-\frac{\Delta y_1}{\delta x_{1-2}}k_{1-2}=-0.002\\ &w:\frac{\Delta y_1}{\delta x_{21-1}}=\frac{0.1}{0.05}=2 &,\ k_w=k_{21}=10^{-3},& a_w=-\frac{\Delta y_1}{\delta x_{21-1}}k_{21-1}=-0.002\\ &n:\frac{\Delta x_1}{\delta y_{1-4}}=\frac{0.1}{0.1}=1&,\ k_n=k_{1-4}=10^{-3},& a_4=-\frac{\Delta x_1}{\delta y_{1-4}}k_{1-4}=-0.001 \\ & a_1=-(a_2+a_4+a_w)=0.005\\ &a_1T_1+a_2T_2 +a_4T_4=b_1\\ &\\ &0.005T_1-0.002T_2 -0.001T_4=b1=0+a_wT_w=0.002\times 320=0.64\\ \end{aligned}$
element 2
$\begin{aligned} &e:\frac{\Delta y_2}{\delta x_{2-3}}=\frac{0.1}{0.25}=0.4 ,&\ k_e=k_{2-3}=10^{2},&\ a_3=-\frac{\Delta y_2}{\delta x_{2-3}}k_{2-3}=-40\\ &w:\frac{\Delta y_2}{\delta x_{1-2}}=\frac{0.1}{0.15}=0.667, &\ k_w=k_{1-2}=3\times10^{-3},& \ a_1=-\frac{\Delta y_2}{\delta x_{1-2}}k_{1-2}=-0.002\\ &n:\frac{\Delta x_2}{\delta y_{3-6}}=\frac{0.2}{0.1}=2, &\ k_n=k_{2-5}=10^{2},& \ a_5=-\frac{\Delta x_2}{\delta y_{2-5}}k_{2-5}=-200 \\ & a_2=-(a_1+a_3+a_5)=240.002\\ &a_1T_1+a_2T_2+a_3T_3 +a_5T_5=b_1\\ &\\ &-0.002T_1+240.002T_2-40T_3 -200T_5=b2=q_b \boldsymbol j\cdot \boldsymbol S_1=100\times 0.2=20\\ \end{aligned}$
element 3
$\begin{aligned} &w:\frac{\Delta y_3}{\delta x_{2-3}}=\frac{0.1}{0.25}=0.4, &\ k_w=k_{2-3}=10^{2},& \ a_2=-\frac{\Delta y_3}{\delta x_{2-3}}k_{2-3}=-40\\ &n:\frac{\Delta x_3}{\delta y_{3-6}}=\frac{0.3}{0.1}=3, &\ k_n=k_{3-6}=10^{2},& \ a_5=-\frac{\Delta x_3}{\delta y_{3-6}}k_{3-6}=-300 \\ & a_3=-(a_2+a_5)=340\\ &a_2T_2+a_3T_3 +a_5T_5=b_3\\ &\\ &-40T_2+340T_3 -300T_5=b3=q_b \boldsymbol j\cdot \boldsymbol S_1=100\times 0.3=30\\ \end{aligned}$
element 4
$\begin{aligned} &e:\frac{\Delta y_4}{\delta x_{4-5}}=\frac{0.1}{0.15}=0.667, &\ k_e=k_{4-5}=3\times10^{-3},& \ a_e=-\frac{\Delta y_4}{\delta x_{4-5}}k_{4-5}=-0.002\\ &w:\frac{\Delta y_4}{\delta x_{20-4}}=\frac{0.1}{0.05}=2, &\ k_w=k_{20-4}=10^{-3},& \ a_5=-\frac{\Delta y_4}{\delta x_{20-4}}k_{20-4}=-0.002\\ &n:\frac{\Delta x_4}{\delta y_{4-7}}=\frac{0.1}{0.1}=1, &\ k_n=k_{4-7}=10^{-3},& \ a_7=-\frac{\Delta x_4}{\delta y_{4-7}}k_{4-7}=-0.001 \\ &s:\frac{\Delta x_4}{\delta y_{1-4}}=\frac{0.1}{0.1}=1, &\ k_s=k_{1-4}=10^{-3},& \ a_1=-\frac{\Delta x_4}{\delta y_{1-4}}k_{1-4}=-0.001 \\ & a_4=-(a_1+a_5+a_7+a_e)=0.006\\ &a_1T_1+a_4T_4+a_5T_5+a_7T_7 =b_4\\ &\\ &-0.001T_1+0.006T_4-0.002T_5-0.001T_7 =b_4=a_wT_w=0.002\times 320=0.64\\ \end{aligned}$
element 5
$\begin{aligned} &-200T_2-0.002T_4+440.002T_5-40T_6-200T_8 =0\\ \end{aligned}$
element 6
$\begin{aligned} &-300T_3-40T_5+640T_6-300T_9 =0\\ \end{aligned}$
element 7
$\begin{aligned} &-0.001T_4+0.006998T_7-0.002T_8 =1.2394\\ \end{aligned}$
element 8
$\begin{aligned} &-200T_5-0.002T_7+243.9624T_8-40T_9 =1188.12\\ \end{aligned}$
element 9
$\begin{aligned} &-300T_6-40T_8+345.9406T_9 =1782.18\\ \end{aligned}$
用Matlab使用Gauss-shield方法求解
高斯-赛德尔（Gauss-Seidel）解线性方程组的Matlab实现](https://blog.csdn.net/weixin_45102840/article/details/105779775)

clc;
clear;

% //  输入量：
% A: 线性方程组的系数矩阵（n*n，非奇异）
% b: 方程组右边的常数项列向量
% n: 方程组维数
% x0: 初始值
% tol: 精度上限值
% N:  最大迭代次数

A=zeros(9,9);
A(1,:)=[0.005 -0.002 0  -0.001 0 0 0 0 0];
A(2,:)=[-0.002 240.002 -40 0 -200 0 0 0 0];
A(3,:)=[0  -40  340 0  0 -300 0 0 0];
A(4,:)=[-0.001 0 0 0.006 -0.002 0 -0.001 0 0];
A(5,:)=[0 -200 0 -0.002 440.002 -40 0 -200 0];
A(6,:)=[0 0 -300 0 -40 640 0 0 -300];
A(7,:)=[0 0 0 -0.001 0 0 0.006998 -0.002 0];
A(8,:)=[0 0 0 0 -200 0  -0.002 243.9624 -40];
A(9,:)=[0 0 0 0 0 -300 0 -40  345.9406];

b=[0.64;20;30;0.64;0;0;1.2394;1188.12;1782.18];

n=9;

x0=[300;300;300;300;300;300;300;300;300];

tol=1e-9;

N=10000;

x=Gauss_Seidel_fun(A,b,n,x0,tol,N);


function x=Gauss_Seidel_fun(A,b,n,x0,tol,N)
x=zeros(n,1);  % 给x赋值
k=1;
while k<N
    for i=1:n
        if i==1
            x(1)=(b(1)-A(1,2:n)*x0(2:n))/A(1,1);
        elseif i==n
            x(n)=(b(n)-A(n,1:n-1)*x(1:n-1))/A(n,n);
        else
            x(i)=(b(i)-A(i,1:i-1)*x(1:i-1)-A(i,i+1:n)*x0(i+1:n))/A(i,i);       
        end  
    end
    if norm(x-x0)<tol
        break;
    end
    
    x0=x;
    k=k+1;
    
    disp(['when k=',num2str(k)])
    disp('x=');
    disp(x);                       %输出中间结果
end

if k==N
    disp('迭代次数已到达上限!');
end
disp(['迭代次数 k=',num2str(k)])

end

在这里插入图片描述

沿着下边界的T值
$\begin{aligned} &-k_b\frac{\partial T}{\partial y}=q_b\Rightarrow-k_b\frac{T_C-T_b}{y_C-y_b}=q_b\Rightarrow T_b=T_C+\frac{q_b}{k_b}(y_C-y_b)\\ & T_{10}=T_1=312.490\\ & T_{11}=T_2+\frac{q_b}{k_{11}}(y_2-y_{11})=305.254+\frac{100}{100}(0.05-0)=305.304\\ & T_{12}=T_3+\frac{q_b}{k_{12}}(y_3-y_{12})=305.254+\frac{100}{100}(0.05-0)=305.304 \end{aligned}$
沿着右边边界的T值
$\begin{aligned} & -k_b\frac{\partial T}{\partial y}=0 \Rightarrow T_b=T_C\\ & T_{13}=T_3=305.254\\ & T_{14}=T_6=305.154\\ & T_{15}=T_9=305.054\\ \end{aligned}$
沿着顶部边界的T值
$\begin{aligned} & T_{18}=\frac{h_bT_\infty+(k_{18}/\delta y_{7-18})T_7}{h_b+(k_{18}/\delta y_{7-18})}=\frac{20\times 300+(10^{-3}/0.05)308.867}{20+(10^{-3}/0.05)}=300.009\\ & T_{17}=305.0040\\ & T_{16}=305.0040\\ \end{aligned}$
左边界全部热传递
$\begin{aligned} Q_{left}&=q_{21}\Delta y_{21}+q_{20}\Delta y_{20}+q_{19}\Delta y_{19}\\ &=k_{21}\frac{T_1-T_{21}}{\delta x_{21-1}}\Delta y_{21}+k_{20}\frac{T_4-T_{20}}{\delta x_{20-4}}\Delta y_{20}+k_{19}\frac{T_7-T_{19}}{\delta x_{7-19}}\Delta y_{19}\\ &=\frac{10^{-3}\times0.1}{0.05}[312.49+311.944+308.867-3\times 320]\\ &=-0.053398W \end{aligned}$
顶部边界全部热传递
$\begin{aligned} Q_{top}&=q_{18}\Delta y_{18}+q_{17}\Delta y_{17}+q_{16}\Delta y_{16}\\ &=h_\infty [\Delta x_{18}(T_{18}-T_\infty)+\Delta x_{17}(T_{17}-T_{\infty})+\Delta x_{16}(T_{16}-T_{\infty})]\\ &=20[0.1(300.009-300)+0.2(305.004-300)+0.3(305.004-300)]\\ &=50.058W \end{aligned}$
底部边界全部热传递
$\begin{aligned} Q_{bottom}&=-q_{b}\Delta y\\ &=-100\times 0.5\\ &=-50W \end{aligned}$
能量守恒
$Q_{top}+Q_{left}+Q_{bottim}=-0.053398+50.058-50=0.004602\approx 0$