17、使用MatLab程序求解磁流体力学方形管道流动与传热的有限差分问题

于 2025-11-24 11:24:37 发布
阅读量15 收藏
点赞数
CC 4.0 BY-SA版权
分类专栏： MATLAB实战应用精粹 文章标签： 磁流体力学 方形管道流动 传热
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.youkuaiyun.com/wine/article/details/155250563

使用MatLab程序求解磁流体力学方形管道流动与传热的有限差分问题

1. 引言

自上世纪以来，磁流体动力学（MHD）因其重要的应用而受到众多研究者的关注。例如，现代发电厂中使用的MHD蒸汽厂和MHD发电机，其基本原理是利用导电液体在垂直磁场中运动产生电能，并且MHD单元的存在提高了卡诺效率。此外，MHD泵和流量计也是重要应用，在这类泵中，电能直接转化为作用于工作流体的力。在金属铸造中，使用超导线圈进行MHD分离也是一个重要应用。

在核聚变反应堆中，锂冷却包层是一个涉及MHD的非常有用的应用。高温等离子体通过环形磁场维持在反应堆中，位于等离子体和磁绕组之间的液态锂循环回路被称为锂包层。锂具有两个功能：吸收反应释放的热能（随后用于发电）以及参与产生氚的核反应。因此，锂包层是反应堆的一个非常重要的组件。然而，包层会受到极强的磁场作用，所以需要了解适当的MHD关系来计算位于不同角度的通道或管道中液态金属的流动，并确定所需的压降、传热等。

MHD的研究可以追溯到19世纪，但直到20世纪初，该领域的广泛研究才加速发展。上世纪30年代，进行了管道和导管中MHD流动的首次理论和实验室研究。Williams发表了电解质在绝缘管中流动的实验结果，通过测量跨流的电位差来研究。Hartmann和Lazarus用汞进行了更全面的理论和实验研究，汞的电导率比电解质高100,000倍，这使得他们能够观察到比Williams实验更广泛的现象，特别是研究了磁场引起的阻力变化和湍流抑制。Hartmann获得了两平行非导电壁之间流动的精确解。

此后，许多研究者对MFM流动进行了研究。Shercliff在1956年解决了矩形管道的问题，发现对于高Hartmann数M，速度分布由均匀核心和壁附近的边界层组成，并以此近似解决了圆形管道的问题。Gold和Lykoudis在1962年获得了零壁电导率圆形管中MFM流动的解析解，Gardner和Lykoudis在1968年通过实验获得了圆形管有和无传热的结果。Al - Khawaja等人对圆形管有传热的MFM流动进行了数值研究，Al - Khawaja和Selmi用谱方法获得了均匀壁温下MFM方形管道流动的解。

同时，MFM的混合自由和强制对流管道流动也受到了很多关注。不同研究者采用了不同的方法解决相关问题，如Chang & Lundgren考虑了壁电导率的影响，Gold解析求解了零壁电导率圆形管的MHD问题，Shercliff用二阶近似求解大M时的问题，Gardner用Gold的解评估温度分布等。

2. 背景

本文考虑的问题是在均匀垂直横向磁场B0中，水平圆形管道内的强制对流问题。假设流体为均质、不可压缩、粘性、导电的流体，管道壁面具有均匀表面温度和均匀表面热通量。为了定义这个问题，做出以下假设：
- 所有流体属性是常数（流体为不可压缩），且与温度无关。
- 管道足够长，可假设流动和传热充分发展，忽略入口或出口效应，除压力和温度外，其他变量沿轴向无线性变化。
- 能量方程中的粘性和焦耳耗散贡献较小，可忽略。
- 由于施加磁场B0与主流动或二次流动相互作用产生的感应磁场，与B0相比可忽略不计，这是因为基于流动的磁雷诺数在典型应用条件下远小于1。

3. 基本守恒方程

对于不可压缩牛顿液态金属流体和稳态条件，考虑磁场体力影响的修正Navier - Stokes方程（包括感应方程）和能量方程的向量形式分别为：
[
\begin{align}
\nabla\cdot(\rho\mathbf{V}\mathbf{V})&=-\nabla p+\mu\nabla^2\mathbf{V}+\mathbf{J}\times\mathbf{B}\
\nabla\cdot(\mu\mathbf{H})&=0\
\rho c_p\nabla\cdot(\mathbf{V}T)&=k\nabla^2T+\Phi+\sigma J^2
\end{align}
]
同时，两个向量满足无散条件：
[
\nabla\cdot\mathbf{V}=0,\quad\nabla\cdot\mathbf{H}=0
]

对于非常小的磁雷诺数RM，感应方程可由麦克斯韦方程和无散条件推导得出。能量方程右侧的最后两项分别代表粘性和焦耳耗散，与其他项相比可忽略。

经过假设充分发展流动（二维问题）且由于磁场抑制湍流，流动为层流，该流动的无量纲控制方程变为：
[
\begin{align}
-\nabla^2w^ +M\frac{\partial H^ }{\partial x^ }&=1\
-\nabla^2H^ +M\frac{\partial w^ }{\partial x^ }&=0\
\nabla^2\theta+4Nu w^ \theta&=0\
\nabla^2\theta - 4w^ &=0
\end{align}
]
其中，负无量纲压力梯度γ与w 的关系为：
[
\gamma=\frac{1}{\int_{0}^{1}\int_{0}^{1}w^ dx^ dy^ }
]
从力和能量平衡可以分别得出：$fRe = - 2\gamma$和$Nu = -\frac{1}{\theta_m}$，其中平均无量纲温度$\theta_m$为：
[
\theta_m=\frac{\int_{0}^{1}\int_{0}^{1}w^ \theta dx^ dy^ }{\int_{0}^{1}\int_{0}^{1}w^ dx^ dy^ }
]
边界条件为：$w^ = 0$（无滑移条件），$H^ = 0$（电绝缘表面），$\theta = 0$（等温表面和恒定表面热通量）。

4. 数值研究

4.1 有限差分格式

由于所考虑的偏微分方程具有稳态特性，属于椭圆型方程。为了将这些方程离散化并得到具有合理精度和稳定性的线性代数方程组，需要引入一些重要的差分格式。
- 中心差分格式 ：中心差分格式常用于椭圆型方程，如Laplace方程和Poisson方程。对于二维问题，可表示为五点公式、对角五点公式或九点公式等。然而，该格式并非适用于所有椭圆型方程，例如，使用中心差分格式通过松弛方法求解稳态Navier - Stokes方程时，如果网格雷诺数超过2，数值解可能会变得不稳定且无法收敛，这是因为在高雷诺数下，标准中心差分对对流项的处理会破坏差分方程的椭圆性，导致所得矩阵失去对角优势。
- 一阶和二阶精确迎风格式 ：一阶精确迎风格式被许多研究者用于处理高雷诺数（在本文中为Grashof数的平方根）下对流项的稳定性问题，扩散算子和源项可采用中心差分。但该格式存在两个主要缺点：一是在偏差方向引入了较大的人为扩散，导致精度显著损失；二是算法的整体精度为$O({Gr}^{1/2}h)$，在高Grashof数下，即使网格较细，误差也可能占主导地位，掩盖了物理扩散对流动的影响。Atias等人采用二阶精确迎风格式离散涡度输运方程中的对流项，该格式的整体精度为$O(Re h^2)$，基于冯·诺伊曼型稳定性分析，当网格雷诺数小于$2 + \sqrt{8}$时，二阶迎风格式的Gauss - Seidel解是稳定的。
- 三阶精确迎风格式 ：该格式由Agarwal首次引入，用于在高雷诺数下计算Navier - Stokes方程的解。他用该格式解决了多种流动问题，如二维驱动方腔流动、通道突然对称扩张流动等，并获得了与其他研究者计算结果和实验数据良好的一致性。然而，该格式也有两个主要缺点：一是边界条件的数值处理需要额外的注意，因为算法使用五点差分公式而不是标准的三点公式来计算一阶导数；二是线松弛求解需要进行五对角矩阵求逆。
- 修正三阶精确迎风格式 ：Agarwal对三阶迎风格式进行了修改，使该格式变为二阶精确，但消除了上述缺点。与二阶迎风格式相比，新算法的人为扩散较小。

4.2 求解线性代数方程组的迭代方法

求解线性代数方程组的方法分为直接法和迭代法。直接法通过复杂算法在有限且可预测的操作次数内给出精确解（如果不存在舍入误差），而迭代法通过重复应用简单算法，只有在迭代过程收敛时，才能在有限但通常不可预测的操作次数内接近精确解。因此，本文采用迭代法，这类方法有时也称为松弛法，可分为点（或显式）迭代法和块（或隐式）迭代法。
- 点 - Gauss - Seidel迭代 ：该方法是显式的，其应用于一般代数方程组的步骤如下：
1. 对所有未知量进行初始猜测。
2. 求解每个方程中系数绝对值最大的未知量，初始使用猜测值，后续使用最近计算的值。
3. 重复上述步骤，直到未知量的变化变得很小。
- Gauss - Seidel过程收敛的充分条件 ：点 - Gauss - Seidel迭代法简单，但只有在系数矩阵具有对角优势的特定条件下才收敛。幸运的是，许多稳态守恒方程的差分格式提供了这种对角优势。该方法应用于代数方程组收敛的充分条件是系统不可约（不能通过求解少于m个方程来确定部分未知量），且差分方程得到的系数矩阵具有对角优势。
- 逐次超松弛（SOR） ：逐次超松弛是一种可用于加速任何迭代过程的技术，本文主要将其作为Gauss - Seidel方法的改进。当某点的收敛呈现振荡模式并趋于超过最终解时，逐次欠松弛（SUR）似乎更为合适。超松弛通常适用于具有Dirichlet边界条件的Laplace方程的数值解，而欠松弛在椭圆型非线性问题中有时是必要的。一般来说，没有具体公式确定最佳松弛因子，有时可通过数值实验确定。
- 线迭代松弛法 ：线迭代松弛算法有时也称为块迭代法，由于其具有隐式性质，也被称为隐式迭代法。该方法可以与几乎任何迭代算法结合使用，通常在Gauss - Seidel方法的框架下结合SOR或SUR使用，并且有多种应用SOR或SUR的方式。

5. 求解过程

本文针对方形管道的两种传热极限情况（恒定温度和恒定热通量边界条件），使用MatLab软件和Gauss - Seidel迭代法进行数值研究。将具有均匀温度条件（能量方程为式(7)）和均匀热通量条件（能量方程为式(8)）的修正无量纲Navier - Stokes方程转换为有限差分方程（使用中心差分格式）：
[
\begin{align}
w_{i,j}^ &=\frac{1}{4}(w_{i + 1,j}^ +w_{i - 1,j}^ +w_{i,j + 1}^ +w_{i,j - 1}^ )-\frac{\Delta x M}{8}(H_{i + 1,j}^ -H_{i - 1,j}^ )-\frac{(\Delta x)^2}{4}\
H_{i,j}^ &=\frac{1}{4}(H_{i + 1,j}^ +H_{i - 1,j}^ +H_{i,j + 1}^ +H_{i,j - 1}^ )-\frac{\Delta x M}{8}(w_{i + 1,j}^ -w_{i - 1,j}^ )\
\theta_{i,j}&=\frac{1}{4}(\theta_{i + 1,j}+\theta_{i - 1,j}+\theta_{i,j + 1}+\theta_{i,j - 1})-(\Delta x)^2Nu w_{i,j}\
\theta_{i,j}&=\frac{1}{4}(\theta_{i + 1,j}+\theta_{i - 1,j}+\theta_{i,j + 1}+\theta_{i,j - 1})-(\Delta x)^2w_{i,j}
\end{align}
]
同时，无量纲压力梯度γ和平均无量纲温度$\theta_m$的定义分别为：
[
\gamma=\frac{1}{\sum_{i = 0}^{I}\sum_{j = 0}^{J}(\Delta x^ )^2w_{i,j}^ },\quad\theta_m=\frac{\sum_{i = 0}^{I}\sum_{j = 0}^{J}(\Delta x^ )^2\theta_{i,j}w_{i,j}}{\sum_{i = 0}^{I}\sum_{j = 0}^{J}(\Delta x^ )^2w_{i,j}}
]
边界条件为$H^ = 0$（电绝缘壁），$w^ = 0$（无滑移条件），$\theta = 0$（无量纲温度定义）。

对于低到中等Hartmann数（M = 0到100），使用均匀的101×101网格；对于高Hartmann数（M = 200），使用均匀的201×201网格。通过根均方残差R（在程序中定义）来测试解的收敛性，当R小于一个很小的数（如$10^{-7}$）时，认为解是可接受的。

从图中可以看出，归一化轴向速度残差的收敛性随着Hartmann数的增加而提高。例如，在M = 200时，残差在500次迭代后达到$10^{-3}$，而无磁场时，残差在5800次迭代后达到$10^{-3}$。归一化磁场残差也有类似的行为。然而，无量纲温度的情况有所不同，对于均匀温度情况，残差收敛较慢，特别是在M = 200时，需要超过29500次迭代才能使残差达到$10^{-3}$；对于均匀热通量情况，在M = 200时，需要33800次迭代使残差收敛到$10^{-3}$，而在M = 0时，需要超过10600次迭代。

以下是整个求解过程的mermaid流程图： 

graph TD;
    A[开始] --> B[设置参数和初始条件];
    B --> C[选择网格大小（根据M）];
    C --> D[初始化变量（ws, Hs, t）];
    D --> E[设置边界条件];
    E --> F[迭代求解w*和H*];
    F --> G{是否收敛};
    G -- 否 --> F;
    G -- 是 --> H[计算γ, w, H, f];
    H --> I{边界条件类型};
    I -- 均匀热通量 --> J[迭代求解θ];
    I -- 均匀温度 --> K[假设Nu，迭代求解θ];
    J --> L{θ是否收敛};
    K --> L;
    L -- 否 --> J;
    L -- 是 --> M[计算thetam和Nu];
    M --> N{是否满足误差要求};
    N -- 否 --> K;
    N -- 是 --> O[输出结果];
    O --> P[结束];

下面是不同有限差分格式的特点总结表格：
| 差分格式 | 优点 | 缺点 |
| ---- | ---- | ---- |
| 中心差分格式 | 常用于椭圆型方程 | 高雷诺数下可能不稳定，破坏椭圆性 |
| 一阶精确迎风格式 | 处理高雷诺数对流项稳定性好 | 引入人为扩散，精度损失大 |
| 二阶精确迎风格式 | 精度为$O(Re h^2)$，稳定性较好 | - |
| 三阶精确迎风格式 | 适用于高雷诺数问题，结果一致性好 | 边界处理复杂，需五对角矩阵求逆 |
| 修正三阶精确迎风格式 | 二阶精确，人为扩散小 | - |

6. 结果分析

6.1 温度分布结果

对于MFM方形管道流动，在均匀温度和均匀热通量边界条件下得到了一些显著的传热结果。轴向速度由于磁场的存在而变得更加平坦，摩擦因数随磁场增强而增大。在中平面（沿或垂直于磁场方向），负无量纲温度分布总是随着Hartmann数M的增加而减小。这是因为磁场开启后，速度分布变得更加平坦，使得温度分布更加均匀。并且，对于任何Hartmann数，沿磁场方向和垂直于磁场方向的温度分布几乎相同。

不同边界条件和Hartmann数下的温度分布情况如下表所示：
| 边界条件 | Hartmann数M | 负无量纲温度分布变化 | 收敛所需迭代次数（残差到$10^{-3}$） |
| ---- | ---- | ---- | ---- |
| 均匀温度 | M = 0 | 较高 | 约7700 |
| 均匀温度 | M = 200 | 较低 | 超过29500 |
| 均匀热通量 | M = 0 | 较高 | 超过10600 |
| 均匀热通量 | M = 200 | 较低 | 33800 |

从这些结果还可以看出，均匀热通量边界条件下，管道内温度的均匀性更好，这也解释了为什么在任何Hartmann数下，该条件下的Nusselt数更高。

6.2 Nusselt数结果

正如预期的那样，对于恒定壁温和恒定热通量边界条件，Nusselt数Nu都随着Hartmann数M的增加而增加。从传统流动（M = 0）开始，本文得到的Nu值与解析方法得到的值吻合良好。例如，对于恒定热通量和恒定温度边界条件，本文得到的Nu值分别为3.606和2.977，而解析值分别为3.61和2.98。

不同流动几何形状下Nusselt数与Hartmann数的关系如下：
- 对于任何Hartmann数M，圆形管恒定热通量条件下的Nusselt数最高。
- 方形管均匀温度条件下的Nusselt数最低。

以下是Nusselt数随Hartmann数变化的mermaid流程图，展示了不同流动几何形状的对比： 

graph LR;
    A[Hartmann数M增加] --> B[圆形管恒定热通量Nu增加];
    A --> C[方形管恒定热通量Nu增加];
    A --> D[方形管均匀温度Nu增加];
    B --> E[最高Nu值];
    D --> F[最低Nu值];

7. 结论与展望

7.1 结论总结

本文研究了在方形管道中施加均匀横向磁场时，导电流体的流动和传热问题，考虑了两种传热极限情况（均匀热通量和均匀温度边界条件）。由于磁场会抑制流动中的湍流，因此层流假设在MFM流动中大多是有效的。

在之前的研究中，该问题的流体力学部分已通过谱方法进行了广泛研究，且仅展示了均匀温度边界条件下的传热结果。而本文使用迭代Gauss - Seidel方法和MatLab软件进行数值计算，对于恒定温度条件下的结果与之前的研究吻合良好。

7.2 未来研究方向

未来可以从以下两个方面扩展这项工作：
- 考虑纵横比 ：将纵横比（即一般矩形横截面）作为无量纲独立参数加入问题中，与Hartmann数M一起作为独立变量，使问题更具一般性。
- 引入自然对流 ：考虑自然对流，研究横向磁场中的混合强制和自由对流流动问题。这将使问题变得高度非线性，需要采用准确且稳定的算法来解决，并且控制方程中会引入独立的Grashof数，增加问题的复杂性。

8. 代码实现

8.1 均匀表面热通量代码 

%Solution for MHD flow inside square duct for const. heat flux B.C.'s
% a<x<b , c<y<d
% M = Hartmann number
% n = number of subintervals for x
% m = number of subintervals for y
% h = delta x*
% k = delta y*
% Note: In this program delta x = delta y
a=0; b=1; c=0; d=1; num_iter=20000; M=100;
n=100; m=100; h=(b-a)/n; k=(d-c)/m;
c1=1/4;    c2=(h*M)/8;  c3=(h^2)/4;
% ws = negative normalized axial velocity (w*)
% Hs = normalized induced axial magnetic field (H*)
ws=zeros(n+1, m+1);
Hs=zeros(n+1, m+1);
%B.C.'s at the four corners for w*(no slip conditions) & H* (electrically insulated surface)
    ws(1,1)=0;
    ws(n+1,1)=0;
    ws(1,m+1)=0;
    ws(n+1,m+1)=0;
    Hs(1,1)=0;
    Hs(n+1,1)=0;
    Hs(1,m+1)=0;
    Hs(n+1,m+1)=0;
%B.C.'s at the four sides for w*(no slip conditions) & H* (electrically insulated surface)
    for i=2:n
        ws(i,1)=0;
        ws(i,m+1)=0;
        Hs(i,1)=0;
        Hs(i,m+1)=0;
    end
    for j=2:m
        ws(1,j)=0;
        ws(n+1,j)=0;
        Hs(1,j)=0;
        Hs(n+1,j)=0;         
    end    
    for it=1:num_iter
    wsave=w;
    Hsave=H;  
    Rws=0;
    RHs=0;
   % Solution for w* & H*    
    for i=2:n
        for j=2:m
       ws(i,j)= c1*(ws(i,j+1)+ws(i,j-1)+ws(i+1,j)+ws(i-1,j))-c2*(Hs(i+1,j)-Hs(i-1,j))-c3;
       Hs(i,j)= c1*(Hs(i,j+1)+Hs(i,j-1)+Hs(i+1,j)+Hs(i-1,j))-c2*(ws(i+1,j)-ws(i-1,j));      
            % Rws = Root mean square residuals for w*
            % RHs = Root mean square residuals for H*             
            Rws=Rws+sqrt((w(i,j)-wsave(i,j))^2);
            RHs=RHs+sqrt((H(i,j)-Hsave(i,j))^2);
        end
    end
   Rwss(it,1)=Rws;
   RHss(it,1)=RHs;
    if (RHs<1e-8 & Rws<1e-8)
        break
    end
end
% gamma = Non-dimensional pressure gradient
% w = Dimensionless axial velocity
% H = Dimensionless induced axial magnetic field
% f = friction factor
gamma=1/sum(sum(ws*h^2));
w=ws*gamma;
H=Hs*gamma;
f=-2*gamma;
% t = Dimensionless temperature (theta)
t=zeros(n+1,m+1);
%B.C.'s at the four corners for theta (uniform surface heat flux)
    t(1,1)=0;
    t(1,m+1)=0;
    t(n+1,1)=0;
    t(n+1,m+1)=0;
%B.C.'s at the four sides for theta (uniform surface heat flux)
    for i=2:n
         t(i,m+1)=0;
         t(i,1)=0;
     end
    for j=2:m
         t(1,j)=0;
         t(n+1,j)=0;
     end
for itt=1:num_iter
     tsave=t;
     Rt=0;
     % Solution for theta
     for i=2:n
         for j=2:m
             t(i,j)=1/4*(t(i-1,j)+t(i+1,j)+t(i,j-1)+t(i,j+1))-h^2*w(i,j);
             % Rt = Root mean square residuals for theta              
             Rt=Rt+sqrt((t(i,j)-tsave(i,j))^2);
         end
     end
     Rtt(itt,1)=Rt;
         if Rt<1e-8
        break
    end
 end      
% thetam = Mean dimensionless temperature
% nu = Nusselt number (Nu)
thetam=(sum(sum(t.*w*h^2)))/(sum(sum(w*h^2)));      
nu=-1/thetam;    
% Display solution with x from left to right  
ws=[ws'];  
w=[w'];
Hs=[Hs'];  
H=[H'];
t=[t'];
x= a:h:b; y=c:k:d;

8.2 均匀表面温度代码 

%Solution for MHD flow inside square duct for const. temperature B.C.'s
% a<x<b , c<y<d
% M = Hartmann number
% n = number of subintervals for x
% m = number of subintervals for y
% h = delta x*
% k = delta y*
% Note: In this program delta x = delta y
a=0; b=1; c=0; d=1; num_iter=20000; M=100;
n=100; m=100; h=(b-a)/n; k=(d-c)/m;
c1=1/4;    c2=(h*M)/8;  c3=(h^2)/4;
% ws = negative normalized axial velocity (w*)
% Hs = normalized induced axial magnetic field (H*)
ws=zeros(n+1, m+1);
Hs=zeros(n+1, m+1);
%B.C.'s at the four corners for w*(no slip conditions) & H* (electrically insulated surface)
    ws(1,1)=0;
    ws(n+1,1)=0;
    ws(1,m+1)=0;
    ws(n+1,m+1)=0;
    Hs(1,1)=0;
    Hs(n+1,1)=0;
    Hs(1,m+1)=0;
    Hs(n+1,m+1)=0;
    %B.C.'s at the four sides for w*(no slip conditions) & H* (electrically insulated surface)
    for i=2:n
        ws(i,1)=0;
        ws(i,m+1)=0;
        Hs(i,1)=0;
        Hs(i,m+1)=0;
    end
    for j=2:m
        ws(1,j)=0;
        ws(n+1,j)=0;         
        Hs(1,j)=0;
        Hs(n+1,j)=0;
    end
    for it=1:num_iter
    wsave=w;
    Hsave=H;  
    Rws=0;
    RHs=0;
   % Solution for w* & H*
    for i=2:n
        for j=2:m
       ws(i,j)= c1*(ws(i,j+1)+ws(i,j-1)+ws(i+1,j)+ws(i-1,j))-c2*(Hs(i+1,j)-Hs(i-1,j))-c3;
       Hs(i,j)= c1*(Hs(i,j+1)+Hs(i,j-1)+Hs(i+1,j)+Hs(i-1,j))-c2*(ws(i+1,j)-ws(i-1,j));      
            % Rws = Root mean square residuals for w*
            % RHs = Root mean square residuals for H*             
            Rws=Rws+sqrt((w(i,j)-wsave(i,j))^2);
            RHs=RHs+sqrt((H(i,j)-Hsave(i,j))^2);
        end
    end
   Rwss(it,1)=Rws;
   RHss(it,1)=RHs;
    if (RHs<1e-8 & Rws<1e-8)
        break
    end
end
% gamma = Non-dimensional pressure gradient
% w = Dimensionless axial velocity
% H = Dimensionless induced axial magnetic field
% f = friction factor
gamma=1/sum(sum(ws*h^2));
w=ws*gamma;
H=Hs*gamma;
f=-2*gamma;
% t = Dimensionless temperature (theta)
t=zeros(n+1,m+1);
%B.C.'s at the four corners for theta (uniform surface temperature)
    t(1,1)=0;
    t(1,m+1)=0;
    t(n+1,1)=0;
    t(n+1,m+1)=0;
%B.C.'s at the four sides for theta (uniform surface temperature)
    for i=2:n
         t(i,m+1)=0;
         t(i,1)=0;
     end
    for j=2:m
         t(1,j)=0;
         t(n+1,j)=0;
     end
for itt=1:num_iter
     tsave=t;
     Rt=0;
     % Solution for theta
     for i=2:n
         for j=2:m
             t(i,j)=(t(i-1,j)+t(i+1,j)+t(i,j-1)+t(i,j+1))/(4-4*h^2*nu*w(i,j));
             % Rt = Root mean square residuals for theta
             Rt=Rt+sqrt((t(i,j)-tsave(i,j))^2);
         end
     end
     % thetam = Mean dimensionless temperature
     % nu = Nusselt number (Nu)
     thetam=(sum(sum(t.*w*h^2)))/(sum(sum(w*h^2)));      
     nu=-1/thetam;         
     Rtt(itt,1)=Rt;
         if Rt<1e-8
        break
    end
 end      
% Display solution with x from left to right  
ws=[ws'];  
w=[w'];
Hs=[Hs'];  
H=[H'];
t=[t'];
x= a:h:b; y=c:k:d;

8.3 能量方程推导说明

简化能量方程的详细推导在相关文献中有给出。其关键在于忽略能量方程中代表粘性和焦耳耗散的最后两项，并应用能量平衡：流体从表面获得的热量等于流体通过对流带走的能量。在应用上述假设、省略轴向传导项并定义适当的无量纲变量后，即可得到简化的无量纲能量方程。

通过以上的研究和分析，我们对MFM方形管道流动的传热问题有了更深入的理解，并且为未来的研究提供了有价值的参考方向。

评论
成就一亿技术人!
拼手气红包6.0元
还能输入1000个字符  | 博主筛选后可见
 
红包 添加红包
表情包 插入表情
表情包 代码片
 条评论被折叠 查看
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值