当前位置：文档之家 › 方腔环流的流场计算

方腔环流的流场计算

  • 格式：pdf
  • 大小：612.02 KB
  • 文档页数：10

下载文档原格式

  / 10

合集下载

流场的计算(流体力学)

流场的计算(流体力学)

SIMPLE算法
◆SIMPLE算法的流程图 算法的流程图
SIMPLE算法
◆SIMPLE算法的讨论 算法的讨论
①在速度修正方程中，略去邻点速度修正值的影响，这一个做法并 在速度修正方程中，略去邻点速度修正值的影响， 不影响最后收敛的值， 的负担。 不影响最后收敛的值，但加重了修正压力 p 的负担。
' ' ' '
(
)
任一点上速度修正由两部分组成： 任一点上速度修正由两部分组成：一部分是与该速 度在同一方向上的相邻两节点压力修正之差， 度在同一方向上的相邻两节点压力修正之差，这是 产生速度修正的直接动力； 产生速度修正的直接动力；另一部分由相邻点速度 修正所引起， 修正所引起，这又可以视为四周压力修正位置上对 速度修正的间接或隐含影响。
SIMPLE（Semi-Implicit Method for Pressure（ Linked Equations）算法是求解压力耦合方程的半隐 ） ' ∑ anb u nb 式法。在得到速度修正方程式的过程中， 式法。在得到速度修正方程式的过程中，略去了 去掉这一项就称为“半隐” 而保留这一部分时， 项，去掉这一项就称为“半隐”，而保留这一部分时u e' ， 方程就成为一个“全隐”的代数方程。 方程就成为一个“全隐”的代数方程。
SIMPLE算法
◆压力与速度的修正
a e u e = ∑ a nb u nb + S u + p P − p E Ae
' ' ' '
(
)
a n v n = ∑ a nb v nb + S v + p P − p N An
' ' ' '

方腔环流

方腔环流

sn1 0 ，在四边

n 1 s
2 sn sn ，在两侧和底边 h2 2 sn sn h ，在顶面 h2


sn1
3.计算步骤 求解过程，可按照下列步骤进行，


1) 给定 ， 的初值。可取 i,0j 0， i,0j 0
4h 2 4h 2 1 i 1, j 2 i , j i 1, j i , j 1 2 i , j i , j 1 2 2 Re h h
引进松弛因子 方程可化为下列的迭代格式


i 1, j i , j 1 i , j 1
3 n 1 1 n n 其中， R1 max in, ~ 10 4 j i , j ， R2 max i , j i , j ，C 为给定精度，可取 C r 10
7) 4.计算结果 随着雷诺数的增加， 减少，对于一个 25 25 的网格，当 Re 10 时，可取 1 ，而当
0 ，在腔体四边
v 0 ，在 AB 和 CD 上 x
u 0 ，在 AD 和 BC 上 y
其中，Re 为雷诺数
图 1 方腔环流
2.取正方形网格（如图所示） ，网格数视情况可取 25、50、100 等，一般不建议超过 200，雷诺 数 Re 一般不超过 5000，因为超过 5000 流场已经呈现出非线性，开始出现偏差，雷诺数超过 10000 甚至根本不能算出结果，这是因为强非线性情况下描述流场行为的方程已经不能做如上简化。
图 2 方腔 1, j
2 i , j i 1, j h2


i , j 1

方腔流动

方腔流动

边界条件
流函数边界条件：根据已知条件，在四个壁上流函数均为 0，上边界平板速度 1，根据不可 滑移条件确定流函数偏导数条件。
涡量边界条件，采用（1）：
（1）可以用 Taylor 展开建立一般形式的 Thom 公式。假设某壁面切向速度 v ，沿其内法
向 n 有一节点，距离壁面距离为 h, 此点上的流函数为1 ，如果壁面上的流函数值为 0 ，那么
0
(1
0 2h2
v h)
时间导数采用向前 Euler 法，空间导数项可以用中心差分格式离散
计算步骤 1、计算 n + 1 时刻内点的涡量，需要考虑时间步长和空间步长和粘性项的关系。
2、计算 n + 1 时刻的流函数 n1 ，超松弛迭代法，松弛因子为-1.8。
3、计算 n + 1 时刻边界上的涡量
计算流体力学作业
题目
方腔流指顶部平板以恒定速度驱动规则区域内封闭的不可压流体（例如水）的流动，在方 腔流的流动中可以观察到几乎所有可能发生在不可压流体中的流动现象，如图 1 所示方腔流计 算模型图。
图 1 方腔流动示意图 流函数-涡量法以流函数和涡量为未知量，可以消去控制方程中的压力项。根据差分法编写 程序计算方腔流动。 控制方程
-0.5
0
0.5
1
1.5
3.995s 方腔速度矢量图
0
0.5
1
x0=0; x1=1; dx=(x1-x0)/(N-1); %空间步长 y0=0; y1=1; dy=(y1-y0)/(M-1); %空间步长
Thom 公式
1
0.8
0.6
0.4
0.2
0
-0.2
0
0.2
0.4

集水面积计算方法及隧道涌水量计算方法

集水面积计算方法及隧道涌水量计算方法

集水面积计算方法及隧道涌水量计算方法（1）集水面积计算方法：
多平面积计算法是一种常用的计算集水面积的方法，它具体指用地形或水文断面图上的海拔高度点绘制平面，并由断面剖面和它们所确定的轮廓线的交点构成的边线，将区域划分成多个三角形片段，采用三角函数求和法来计算这些三角形片段的面积之和，以计算该区域的面积。

（2）隧道涌水量计算方法：
1）基本量计算法：该方法是根据现测获得的隧道方面粗糙度获取其其涌水量公式中最重要的参数--水力半径，然后通过计算水力半径来确定涌水量L（m3/s），其计算公式如下：L=1.456*R*R*V，其中，R为水力半径（m），V为流速（m/s）。

2）涡折尔定律：这种方法是根据涡折尔定律来确定排水量，这个定律是指在固定的圆形管道中，流体流量Q和流速V成反比，其公式为
Q=πR2V，其中R为管径，V为流速。

3）视比算法：该算法是利用粗糙度、流速和流量的比值来计算涌水量的，通过将当前的粗糙度和流速比值与预先确定的粗糙度和流速比值进行比较，就可以从而计算出涌水量。

不可压N-S方程的高效并行直接求解

不可压N-S方程的高效并行直接求解

不可压N-S方程的高效并行直接求解包芸;叶丰;张义招【摘要】对不可压N-S方程的数值计算,当计算规模增大时,不论是采用湍流模型计算还是直接数值模拟(Direct Numerical Simulation,DNS),大规模的并行计算都难以实现.该问题的关键是求解全场联立的压力泊松方程的并行计算技术.利用并行近似解求解方案,创建高效大规模并行计算的不可压N-S方程的直接求解方法.三维窄方腔热对流的DNS计算结果表明,该直接求解并行计算方法具有很好的并行效率,并且计算的三维湍流热对流的特性是合理的.【期刊名称】《计算机辅助工程》【年(卷),期】2016(025)003【总页数】5页(P19-23)【关键词】泊松方程;并行计算;不可压流动;湍流热对流;直接数值模拟【作者】包芸;叶丰;张义招【作者单位】中山大学力学系,广州 510275;中山大学力学系,广州 510275;中山大学力学系,广州 510275【正文语种】中文【中图分类】O357.5在航空航天等高科技工程的推动下，计算流体力学在可压缩流动的数值模拟计算技术领域进步非凡.不可压流动的数值模拟技术也在不断进步.超级计算机硬件技术的快速发展为计算流体力学数值模拟的进一步发展提供技术支持，高效并行计算技术的发展为进一步扩大不可压N-S方程的数值计算规模提供新的平台，并使计算结果数据能更好地反映流体的流动特性.热对流问题广泛存在于自然界和工业界中，研究其对全球气候变化、海洋环流、反应堆设计、工业冷却设计等问题的影响具有重要意义.[1-2]在Boussinesq假设下，湍流热对流的描述方程为不可压N-S方程联立温度对流扩散方程，因此热对流问题也是典型的不可压流动问题.高瑞利数Ra的湍流Rayleigh-Bénard(RB)热对流的直接数值模拟(Direct Numerical Simulation，DNS)面临重大挑战.[3]随着Ra的提高，热对流进入湍流状态，DNS模拟的规模不断增大导致计算所需要的成本巨大，数值计算难以实现.目前，湍流热对流的DNS只达到Ra=1012水平.[4]大规模可并行的高Ra湍流热对流DNS及其海量数据结果分析已成为热对流研究工作者们特别关注的问题.在不可压流动N-S方程的数值模拟计算中，不论采用何种计算模型或是DNS，其压力泊松方程的求解都是大规模并行计算的难题：直接求解方法不易于并行.迭代求解压力泊松方程通常采用局部算法从而较容易实现并行，但迭代计算过程占用大量的计算时间，所以很难实现高效率的大规模计算.这使得不可压流动的大规模数值模拟难以在有效时间内满足需求，因此妨碍大规模不可压流动N-S方程的湍流模型计算或DNS的进一步发展.一种新的泊松方程块三对角近似求解方案[5-6]可解决在超级计算机上的高效并行直接求解问题.这使得建立不可压流动N-S方程模拟的大规模高效并行计算方法成为可能，并可在超级计算机上实现三维高Ra湍流热对流特性研究的DNS.无量纲化后的三维不可压N-S方程为式中：V为速度矢量；p为压力；Re为雷诺数.计算边界条件为无滑移边界.投影法是数值求解不可压N-S方程组的有效方法之一.[7]实际上，不论采用何种湍流模型或DNS，以及采用何种思路求解不可压N-S方程的速度压力法，一般的动量方程都可以采用显式格式连续方程推导出压力泊松方程并进行迭代求解.大规模并行计算的主要困难为压力泊松方程必须全场联立求解.本文主要针对矩形计算域的特定情况，结合有效的泊松方程高效并行的直接求解方案，创建不可压流动DNS的可高效并行求解计算方法.2.1 网格布置和离散格式计算区域取矩形，见图1.网格数为nx×ny×nz.在设计大规模并行计算时，消息传递接口(Message Passing Interface，MPI)的计算区域沿xOy平面对z方向分割，即图中x方向较粗的线.在该面上并行计算需要数据通信；区域内部用OpenMP并行，无须数据通信.由于直接求解压力泊松方程要用到二维快速傅里叶离散余弦变换[8]，因此在x和y方向必须采用等距网格且最好网格数是2k,在z方向上可根据计算的需要采用非等距网格.本文采用不可压流动计算时常用的交错网格，时间方向采用一阶精度离散，空间采用二阶精度离散格式.2.2 不可压N-S方程的并行求解方案在不可压N-S方程的数值求解过程中，采用投影法将计算分步骤进行.动量方程中的速度计算采用显式格式，在求解中很容易实现并行.由连续方程推导出的压力泊松方程在求解时需要全场联立，是求解过程中计算工作量最大的部分.同时，联立求解也给并行造成困难，是大规模高效并行计算的难点.高效合理并可大规模并行计算的压力泊松方程求解方案，是解决不可压N-S方程大规模并行计算的关键技术.建立三维泊松方程的直接求解方法.计算方法只用于矩形计算区域，x方向采用等距网格.具体做法为在xOy平面上使用二维快速傅里叶变换将空间3个方向上都要求联立求解的压力泊松方程解耦，使泊松方程变换为只在z方向上的三对角方程.将三对角方程变换为块三角方程，设计高效且可并行计算求解方法，求解后再使用对应的反变换得到原来压力泊松方程的全场解.压力泊松方程为x和y方向使用等距网格，z方向使用变距网格.二阶精度中心差分的压力泊松方程的离散形式为使用二维离散余弦傅里叶变换将全场联立的泊松方程在x和y方向上解耦.余弦傅里叶变换能使压力泊松方程自动满足压力梯度为0的边界条件.变换通过FFTW软件包[8]实现.二维离散余弦傅里叶变换公式为将式(4)代入压力泊松方程，并令展开式两边对应系数相等，可以得到一组只沿z方向联立的三对角方程，使压力泊松方程在x和y方向上解耦，求解过程变简单. 在以往的二维热对流DNS中，利用追赶法求解三对角方程的泊松方程直接解法[9]，与采用迭代求解方法在计算机上单线程计算相比有效得多，但三对角方程的追赶法很难进行大规模并行计算.数学和计算机的研究者们尝试建立块三对角方程的大规模高效并行求解方案[5-6]，将变换得来的三对角方程写成块三对角的形式为式中为块三对角矩阵,和为M×N阶矩阵；和为M×N维列向量，T.为已知方程的右端项，通过式(4)求在计算得到块三对角方程的解后，通过对应的二维离散余弦傅里叶的反变换公式利用以上高效并行三维压力泊松方程直接求解方法，联合其他方便并行的动量方程等计算，创建三维不可压N-S方程高效并行直接求解计算方法，使得在一些特定情况下，大规模高效并行的不可压流动N-S方程湍流模型计算或者DNS成为可能.3.1 三维湍流热对流方程RB热对流是研究流体对流传热的典型物理模型.在封闭的盒子内，下底板加热而上底板冷却后形成对流传热的研究系统.在Boussinesq假设下，无量纲化后的三维热对流的描述方程为式中：Ra为瑞利数；Pr为普朗特数；θ为相对温度，下底板为加热，θ=0.5，上底板为冷却，θ=-0.5.通过热对流方程组可以看出，整个计算过程实际上就是数值求解不可压N-S方程组联立温度的对流扩散方程.本文对三维湍流热对流进行DNS.3.2 并行计算效率采用本文建立的三维不可压流动的直接求解并行计算方法，在超级计算机“天河二号”上进行加速比测试.每个计算节点包含24个计算物理核心.测试算例的计算网格数和物理计算核心数见表1.不同计算核数时直接求解并行计算方法的加速比见图2，其中加速比以计算节点24核心数为基数.计算中z方向上网格数1 536是24核心数的64倍.可知当使用32节点即32×24=768核心数时，具有约81.7%的计算效率，当使用64节点即64×24=1 536核心数时，具有约67%的计算效率.加速比随计算机核数的进一步增加仍有较好的可增长性.由于采用快速傅里叶变换解耦压力泊松方程的需要, 并行区域只在z方向上划分，由此本文创建的直接求解方法在大规模并行计算时z方向网格数与计算核心数之间有一定关联.总体上讲，本文的不可压流动直接求解方法在大规模并行计算上已得到较好的计算效率，可以进行较大规模的三维湍流热对流的DNS.新的直接求解计算方法与本课题组原有的迭代计算方法相比，节省1/2以上的总计算时间，计算技术的进步为系列三维热对流的DNS及物理特性的研究提供有力工具.3.3 三维湍流热对流计算结果选取窄方腔，宽高比Γ=1/4，Pr=4.3，对不同Ra的三维热对流进行计算.低Ra计算采用512×64×768网格，高Ra计算采用1 024×128×1 200网格.首先计算流动的初场，根据传热特性计算时平均场统计的需要，当热对流中出现较稳定的大尺度环流流动后，继续迭代计算至少100万时间步.根据不同的Ra及流动特性，计算时间的迭代至少达到300万时间步.RB热对流主要探讨由浮力产生的对流运动对流体传热特性的影响.Ra=1010的平均场温度等值面分布见图3，下底板加热为高温，上底板冷却为低温.温度作为热对流中的重要物理量，其平均场的分布反映传热过程中对流流动带动温度运动的情况.图3中方腔的中间部分没有温度的分布，表明中间没有流动也没有传热作用.图3显示由大尺度环流带动的温度分布形态，RB热对流中对流传热的热量主要沿着侧壁和棱边向上传输.RB热对流研究的核心问题之一是对流传热效率，表征传热的物理参数是努塞尔数Nu，表示流体对流传热与热传导传热的比值.Nu与Ra之间存在一定的标度关系[3]，因此需要一系列Ra数的不同数值模拟结果进行研究.一组三维方腔热对流的Nu/Ra0.3随Ra的变化见图4.图4中同时也给出计算得到的二维热对流Nu数的变化结果[10]，以及KACZOROWSKI等[11]的三维计算结果对比.由此可见，本文的三维湍流热对流的DNS模拟结果是合理的.图4中可以看到，二维和三维热对流的计算Nu随Ra的变化都存在很好的标度率.计算结果与理论预测和大量试验结果得到的结论一致.[3]Γ=1/4三维方腔的Nu总体偏大(向上平移).在试验和计算中也发现同样的Nu平移现象，并且不同的Γ变化会引起不同的Nu向上平移量.[12]在传热Nu对Ra的标度率中，三维和二维流动计算结果产生差异的物理原因，还有待更深入的研究.在三维不可压流动特性的研究过程中，尤其是到湍流阶段，超大规模的数值模拟计算十分必要.依靠超级计算机技术，探索高效并行的计算方法和计算技术，进行大规模的高自由度三维不可压N-S方程的数值计算，对更深入研究不可压流动的物理和流动特性具有很重要的意义.大规模并行数值模拟计算三维不可压N-S方程，最难的计算方法和计算技术问题是压力泊松方程的高效并行求解.利用块三对角方程的大规模高效并行近似解求解方案，创建大规模高效并行计算三维不可压N-S 方程的直接求解方法.DNS是研究高Ra湍流热对流的重要手段之一.热对流问题是典型的不可压流动问题.利用本文创建的高效并行不可压流动的直接求解方法，对高Ra三维湍流热对流进行DNS.通过并行效率的验证计算,证明本文建立的直接求解方法具有较好的并行效率.一系列不同Ra的三维窄方腔热对流计算得到的传热Nu具有合理的标度率. 本文创建的大规模高效并行计算的直接求解方法，为高Ra的湍流热对流大规模高效并行计算和数值模拟研究提供有价值的计算参考.ZHANG W, ZHANG H. Scalable parallel algorithm of block tridiagonal systems for initial boundary value problem[J]. Journal of Shanghai University(Natural Science), 2007, 13(5): 497-503.XU W, BAO Y. An efficient solution for 2D Rayleigh-Bénard convection using FFT[J]. Chinese Journal of Theoretical and Applied Mechanics, 2013, 45(4): 1-6. DOI: 10.6052/0459-1879-12-334.。

环流定理,涡度方程和散度方程

环流定理,涡度方程和散度方程

Ca C Ce C Ca Ce 绝对环流=相对环流+牵连环流：
故相对环流定理形如：
dCa dC dCa dCe ——（*） ，其中， dt 刚已讨论，那么 dt dt dt
Ce
？
○L A dr A d ，有： 由曲线-曲面积分转换（Stokes ）定理：
N区上升，L区下沉，近地面北风，高空南风。实际上引入地转效应后， 不应是单圈环流，而是三圈环流。这就是Hadley 等环流。 当然也可用其解释一些局地风，如海陆风，山谷风等。
RT p0
总之：斜压作用是大气运动中的一个重要因子。
6
§6.2 相对环流定理
已知，绝对速度为相对速度与牵连速度之矢量和：V V r a 两端对环线L积分： ○ LVa dr ○ LV dr ○ L ( r ) dr ，可见：
算子只对Ω运算，故 可互换，且省写下标
( r ) 2 ，代之入牵连环流的表达式（6.12），有：
Ce 2 d 2 d 2 ——（6.14）


~ 在赤道面上的投影，即其法线方向与 一致。 其中，
8
（6.14）代回到（*），有相对环流定理（Bjerknes环流定理）：
由于大中尺度运动是准水平的，故水平运动引起的垂直涡度较重要，
●
故有时又称
●
v u 为涡度 ， x y
v u ) 2 sin f x y
Ωsinφ Ω j Ω
φ
k
而绝对涡度～
a
(
——（6.27）
φ Ωcosφ
பைடு நூலகம்11
§6.4 绝对涡度矢量方程，Taylor-Proudman定理

方腔环流的流场计算

方腔环流的流场计算

中山大学工学院、理论与应用力学刘广编制
实验编号及题目：实验六
方腔环流
《计算流体力学试验》课程实验报告纸
院系：工学院 姓名：刘广 学号：11309018 日期：2014 年 05 月 23 号
方腔环流示意图
我们取网格做如下说明，取正方形网格，网格数量由用户输入，但是建议不超过 200，因 为超过 200 计算量会急剧增大，雷诺数也是由用户输入，建议不超过 5000，因为超过 5000 流 场已经呈现出非线性，开始出现偏差，雷诺数超过 10000 甚至根本不能算出结果，这是因为强 非线性情况下描述流场行为的方程已经不能做如上简化， 同样的， 在计算过程中我们采用超松 弛迭代法，下面对离散形式做一下说明： 我们用二阶精度的差商代替微商做以下说明，得
姓名：刘广
学号：11309018
日期：2014 年 05 月 23 号
程序执行时，我输入的是10E-3作为计算精度，对于程序的说明程序注释已经说的很清楚， 这里不再赘述。最终会在目录下生成OUTPUT.txt 。 四：实验结果 最终会在目录下生成OUTPUT.txt 。然后我们可以用Matlab软件对计算出来的数据进一步整 理，画成云图和画出流线，我们做一下说明。

Re n n [( in, j 1 in, j 1 )( in in1, j in1, j )( in 1, j i 1, j ) ( , j 1 i , j 1 )]} 4 (1 ) in ,j
边界值的点需要特殊处理，然后由内点差分格式顺序扫描，一次次迭代直到满足用户输入的 精度为止。边界条件如下：
( i 1, j 2 i , j i 1, j ) / h 2 ( i , j 1 2 i , j i , j 1 ) / h 2 i , j ( i , j 1 i , j 1 )(

二维方腔环流计算

二维方腔环流计算

���������������,������+���+11
=
���������������������������,���������
−
������������∆������ 4ℎ2
(���������������,���������+1
−
���������������,���������−1)
(���������������+��� 1,������
cos(���⃗���,
���̂���)
4.用 fortran 编程（附件 1）进行计算，求出不同雷诺数（400，,1000,3200 等）下流场、
1
涡量场和压力场的数值解（调整ρ控制收敛速度）；再用 matlab 程序进行绘图（附件 3），进行分析和比较。 （三）使用 ADI 方法和亚松弛迭代法对非定常问题进行求解
+
���������������−���+11,������
+
���������������,������+���−11
+
������������ 4
[(���������������,���������+1
−
���������������,������+���−11)(���������������+��� 1,������
二维方腔环流计算
一、实验目的
二维方形腔体，腔体内部充满了流体。当顶板沿着水平方向被匀速（u=1）拉动时，腔体内 部的流体被带动而作环状运动。这种环流导致了腔体底边的二次涡。 方程与边界条件
用流函数涡量法求解方腔环流， 和 满足下列无量纲定常方程与边界条件

第六章 流场的计算

第六章 流场的计算

第六章 流场的计算6．1 制订一个特殊程序的必要性6．1—1 主要的困难在第五章中，我们推导了在已知流场的情况下求解对Φ的通用微分方程的公式．但是除了某些特殊的情况之外，流场是不可能已知的；相反，我们必须由相应的控制方程计算局部的速度场与密度场。

x x V B u u tu ++=+∂∂)grad div()div()(μρρu （2.11） xp u u t u ∂∂-=+∂∂)grad div()div()(μρρu （2.11）速度分量受动量方程控制，而动量方程则是对Φ的通用微分方程的特殊情况(取Φ＝u ，Г＝μ，等等)．于是我们便想得出结论：我们已经研究出了求解动量方程、因而得到速度场的方法．事实上，问题并不如想象的那么简单困难在哪儿呢?如果说动量方程的非线性是一个困难的话，我们只要提醒一下我们自己，在处理热传导的问题时，我们是如何用迭代法来处理非线性问题的就够了．特别要指出的是，作为动量方程因变量u 的一个函数的对流系数ρ u 与作为温度T 的一个函数的导热系数，这两者之间是没有什么不同的．从一个估计的速度场开始，我们可以迭代求解动量方程，从而得到速度分量的收敛解．计算速度场的真正困难在于未知压力场．压力梯度构成动量方程中源项的一部分．然而还没有一个可以用来求得压力场的明显方程．对一个给定的压力场，求解动量方程确实并没有什么特别因难．可是确定压力场的途径则似乎是相当迷惘的．()0=∂∂+∂∂j ju x t ρρ （5.1） 0)()()()(=∂∂+∂∂+∂∂+∂∂zw y v x u u t ρρρρρ （6.20） 当正确的压力场代入动量方程时，所得到的速度场满足连续性方程． 直接求解由动量方程与连续性方程所推得的整个离散化方程组． 因为即便是单个因变量的离散化方程组，我们也都乐意采用迭代法来求解，因而直接求解速度分量以及压力的整套方程似乎已经超出了我们所考虑的问题范围了。

6．1—2 以涡量为基础的方法由确定压力所带来的困难已经导致人们提出若干从控制方程中消去压力的方法．于是在二维的问题中，通过交叉微分从两个动量方程中消去压力，就可以导出一个涡量输运方程（这一推导的有关要点将在习题6、1中概述）．把涡量的定义与二维稳态情况下的流函数的定义相结合就是众所周知的“流函数／涡量法”，除了其它许多人之外要特别指出的是迈克斯(1963)、弗罗姆(Fromm)和哈洛(Harlow)(1963)、皮尔逊(Pearson)(1965)、巴拉卡特和克拉克(1966)以及朗切尔和沃尔夫茨坦(1969)等人都描述了这一方法，随后，戈斯曼、潘、朗切尔、斯波尔丁和沃尔夫茨坦(1969)等人所写的书把这一方法讲得更为清楚，并使该方法变得更易于接受了．流因数／涡量法具有某些吸引人的优点．其中，压力消失了，代替处理连续性方程和两个动量方程的问题，我们只需要解两个方程以求得流函数和涡量。

方腔顶盖驱动流动

方腔顶盖驱动流动

一、问题描述方腔顶盖驱动流动如图1所示的一个简化两维方腔（高，宽都等于L)，内部充满水分。

上表面为移动墙，非维化速度为ｕ／ｕ0 =1。

其他三面为固定墙。

试求方腔内水分流动状态。

u=1, v=0u=0,v=0u=0, v=0图1常微分方程理论只能求解极少一类常微分方程；实际中给定的问题不一定是解析表达式，而是函数表，无法用解析解法.二、离散格式数值解法:求解所有的常微分方程 计算解函数 y(x) 在一系列节点a = x 0< x 1<…<x n = b 处的近似值节点间距为步长，通常采用等距节点，即取 hi = h (常数)。

步进式：根据已知的或已求出的节点上的函数值计算当前节点上的函数值，一步一步向前推进。

因此只需建立由已知的或已求出的节点上的函数值求当前节点函数值的递推公式即可。

),...,1()(n i x y y i i =≈欧拉方法几何意义在假设 y n = y (x n )，即第 n 步计算是精确的前提下，考虑公式或方法本身带来的误差: R n = y (x n +1) - y n +1 , 称为局1(,) 0,1,...n n n n y y h f x y n +=+=部截断误差.显式欧拉公式一阶向前差商近似一阶导数223111232()[()()()()][ (,)] ()()h n n n n n n n n n h n R y x y y x hy x y x O h y hf x y y x O h +++'''=-=+++-+''=+推导如下：隐式欧拉公式x n +1点向后差商近似导数推导如下：1()()()n n n y x y x y x h+-'≈111()()() ()()(,)n n n n nn n n n n y x y x hy x y x y y x y y h f x y +++'≈+↑≈≈=+11()()()n n n y x y x y x h++-'≈11()()()()n n n n ny x y x hy x y x y ++'≈+↑≈几何意义设已知曲线上一点 P n (x n , y n ),过该点作弦线，斜率为(x n +1 , y n +1 ) 点的方向场f (x ,y )方向,若步长h 充分小，可用弦线和垂线x =x n +1的交点近似曲线与垂线的交点。

方腔环流实验报告

方腔环流实验报告

实验名称：方腔环流实验实验日期：2023年3月15日实验地点：流体力学实验室实验人员：张三、李四、王五一、实验目的1. 理解和掌握方腔环流的基本原理。

2. 观察和测量方腔环流的速度分布和压力分布。

3. 分析方腔环流在不同条件下的流动特性。

二、实验原理方腔环流是一种典型的二维流动问题，其基本原理是利用流体力学中的Navier-Stokes方程来描述。

在方腔环流中，流体在腔体内受到重力、压力梯度和边界条件的影响，形成复杂的流动结构。

本实验通过改变方腔的几何形状和边界条件，研究方腔环流的流动特性。

三、实验仪器与材料1. 方腔环流实验装置2. 高精度流量计3. 压力传感器4. 数据采集系统5. 计算机及分析软件四、实验步骤1. 安装实验装置，确保其稳定性。

2. 连接数据采集系统，对流量计和压力传感器进行校准。

3. 设置方腔的几何形状和边界条件，例如方腔的尺寸、入口和出口的流量等。

4. 启动实验装置，开始方腔环流的流动过程。

5. 在流动过程中，实时采集流量计和压力传感器的数据。

6. 实验结束后，关闭实验装置，整理实验数据。

五、实验结果与分析1. 速度分布通过对实验数据的分析，可以得到方腔环流的速度分布图。

从图中可以看出，在方腔中心区域，速度较大，而在靠近边界区域，速度较小。

这是因为中心区域的流体受到的压力梯度较大，导致流体流速加快。

2. 压力分布通过对实验数据的分析，可以得到方腔环流的压力分布图。

从图中可以看出，在方腔中心区域，压力较小，而在靠近边界区域，压力较大。

这是因为中心区域的流体流速较快，导致压力降低。

3. 不同条件下的流动特性（1）改变方腔尺寸当方腔尺寸增大时，中心区域的速度和压力都增大，边界区域的速度和压力都减小。

这是因为方腔尺寸增大，流体在腔体内的流动空间增大，导致流速和压力分布发生变化。

（2）改变入口和出口的流量当入口和出口的流量增大时，中心区域的速度和压力都增大，边界区域的速度和压力都减小。

中山大学-计算流体力学--方腔环流实验报告--陈建华

中山大学-计算流体力学--方腔环流实验报告--陈建华

图 23：松弛因子 1.9
图 24：松弛因子-0.1，出现了无穷大。 Re=3000 时 5000 步
图 25： 图 26：
图 27：5000 步
图 28： 50000 步
图 29： 图 30：
图 31：
四、结果分析
1. 在定常与非定常问题中，松弛因子都在（0，2）之间。当我把松弛因子取 2.1 带入非定常问题之后，Fortran 程序输出结果非常非常慢（就和我用手打字差不 多的速度），这样的的速度一直输出，可能要等到期末结束才能验证其正确性。 而松弛因子是负数 2. 在定常问题里面，网格划分的密集程度对流线略有影响（见图 9 和图 11 对比）， 但是大致轮廓是有的。 3. 定常问题在雷诺数 2000 以前，计算都没算出来二次涡，但是其他同学却能算 出来，也没发现程序哪里错了。2000 以后，有二次涡，5000 雷诺数的时候，和 实际情况是差不多的。但是雷诺数 10000 以后的二次涡却占了整个输出图像的半 壁江山。可能是雷诺数太大了，非线性效应增强，导致算得的偏差比较大？ 4. 非定常问题也和定常问题算出的结果相似，而且对于步数，在其他条件相同
图 2：网格划分示意图 用二阶精度的差商代替上式中的微商，得
引进松弛因子ρ方程可化为下列的迭代格式
这里和书上 P238 的公式略有不同，这里用已知的下一层的迭代数据代替部分上 一层的数据，但思路是一样的。 为了得到涡量的边界条件，令方程在腔体四边也成立，利用 Taylor 公式，有如 下边界条件：
的情况下，一般 5000 步以后就基本上不会有大的结果改变了。
5. 可能是程序在哪里出了点问题，二次涡要么不出现，要么就非常非常大（占
了 1/2）。雷诺数比较低的时候（图 31，Re=500）还是比较符合的。

河南省考研水利工程高频公式总结

河南省考研水利工程高频公式总结

河南省考研水利工程高频公式总结随着社会经济的不断发展，水利工程在城市建设、农田灌溉、水资源管理等方面发挥着重要的作用。

对于要参加河南省考研水利工程专业的考生来说，熟悉高频公式是备考的重要一环。

本文将针对河南省考研水利工程的高频公式进行总结，帮助考生更好地备考。

一、水力学部分公式1. 流速公式在水力学中，流速是指液体通过管道、河道或其他流体工程设施时的速度。

流速公式包括两种，即平均流速和最大流速。

平均流速公式：V = Q / A最大流速公式：Vmax = K × sqrt (R × S)其中，Q表示流量，A表示过流面积，K表示系数，R表示液体的密度，S表示水流的横断面积。

2. 断面流量公式断面流量公式用于计算流体在特定横断面上通过的水量。

断面流量公式：Q = V × A其中，Q表示流量，V表示流速，A表示过流面积。

3. 曼宁公式曼宁公式是描述水流在河道中的流速与河道的横截面积、湿周和水力半径之间关系的经验公式。

曼宁公式：V = (1 / n) × R^(2/3) × S^(1/2)其中，V表示流速，n表示曼宁摩擦系数，R表示水流的水力半径，S表示水流坡度。

4. 水力半径公式水力半径是指泄水管道或河道中液体流动过程中液体与流动空间的相互作用状态的平均直径。

水力半径公式：R = A / P其中，R表示水力半径，A表示过流面积，P表示湿周。

5. 孔流公式孔流公式用于计算孔或缝隙中水通过的流量。

孔流公式：Q = C × A × sqrt (2gh)其中，Q表示流量，C表示流量系数，A表示缝隙面积，g表示重力加速度，h表示水头。

二、水环境工程部分公式1. 简化沉降深度公式简化沉降深度公式用于计算土壤的沉降深度。

简化沉降深度公式：Δh = (C × H) / (1 + e_0)其中，Δh表示沉降深度，C表示可压缩性指数，H表示填筑高度，e_0表示初始土壤压缩率。

(完整版)不同水深流速分布及推力计算

(完整版)不同水深流速分布及推力计算

一、流速分布及计算自然界中的水流大部分是湍流。

湍流是一种高度复杂的非线性流体运动，在空间中不规则、时间上无秩序，具有在运动过程中液体质点不断混掺的运动特性。

实际中流速计算一般根据实测数据进行推导，具有代表性的是“六点测流法”，2014年之后，声学多普勒流速剖面仪开始被采用，随后有部分学者提出了相应的“多点法测速计算”。

水流由于受到层间切应力的作用，其流速沿水深而变化，河底流速小，水面流速大，河底流速受河床的粘滞作用，基本为零。

理论上水流流速由下往上可分成直线层、过渡层、对数区和外层区，其相应的计算公式如下：(一) 直线层水流为层流（层流是流体的一种流动状态，它作层状的流动。

流体在管内低速流动时呈现为层流，其质点沿着与管轴平行的方向作平滑直线运动。

流体的流速在管中心处最大，其近壁处最小。

管内流体的平均流速与最大流速之比等于0.5。

），只受粘滞切应力，此时流速可按下式计算：μy=√ghJJ：水力坡度；0≤y<0.5%。

水力坡度，又称比降，是指河流水面单位距离的落差，常用百分比、千分比、万分比表示。

(二) 过渡层水流由层流向紊流过度，既受粘滞切应力，又受紊动切应力。

计算方法：近似按照直线层或者对数层公式计算。

（三）对数区水流为紊流，主要受紊动切应力影响，流速分布呈对数曲线规律,一般计算公式如下：uμy=A?lgy+B其中A和B是系数，与床面粗糙情况有关，通过实际资料确定，y为计算点至河床的距离。

爱因斯坦提出的具体计算公式如下:μμy =5.75lg(30.2yk sx)其中k s为床面粗糙高度，可取床沙代表粒径；x为反映对流速分布实际影响的系数，与k sδ值有关；δ:为近壁层流层的厚度。

直线层、过度层、对数区合称为内层区，区内流速分布主要受床面的影响。

（四）外层区水流为紊流，其流速分布除受床面的影响外，还要受到上游来流条件和上部边界条件的影响，因而其分布规律偏离对数曲线而有一流速增值，计算公式的一般计算形式为：μμ?=A?lgy+B+πk?ω(yh)式中，π为尾迹强度系数；k为卡门常数；ω为函数符号；π和k通过实测资料确定。

矩形断面有压管道汇流口流场数值计算

矩形断面有压管道汇流口流场数值计算

矩形断面有压管道汇流口流场数值计算随着科技的发展，水力学的研究也在快速发展，汇流口的数值计算是水力学中的重要内容。

在现实的工程中，管道汇流是一种极其常见的现象，对管道汇流口的流场进行有效的数值计算是很有必要的。

本文采用有限体积法（FVM）基于流体动力学方程，结合矩形断面有压管道汇流口的流场，进行数值计算，为管道汇流口流场的实际应用提供决策参数和研究内容。

管道汇流口是指两条或多条河流、池塘或渠道在某一点发生汇合的水体现象，它是水体流动过程中的重要组成部分，其发生的位置和形式都会发生变化。

当管道汇流口的流场作用力施加在水体上时，会产生一系列的运动现象，影响着水体的流动状况，包含了水流的流向、速度分布及水体的混合状况。

因此，研究和分析管道汇流口的流场对于解决水体运动问题具有重要的意义。

由于水体运动受到水流和地质赋存条件的限制，一般地认为，管道汇流口主要分为覆盖汇流口和管道汇流口两种形式。

在覆盖汇流口中，汇流口的管道主要以抛物面形式出现，即管道的下缘在汇流口的下游，管道的上缘在汇流口的上游，形成一种水体覆盖汇流口的状态；而管道汇流口中，汇流口的管道主要以水平形式排列，即管道断面水平，形成两两相连的水流流入管道汇流口，形成一种排水层状态。

本文主要研究的是矩形断面有压管道汇流口的流场，该流场在实际工程中也被广泛应用。

矩形断面有压管道汇流口的流场数值计算，可以从水力学方程组出发，分析条件；或者采用有限体积法（FVM），针对矩形断面有压管道汇流口的流场，进行连续体（水）的求解。

FVM以水体作为连续体，根据水力学方程，把它分割成一系列的有限小的控制体，以求解流体在每个控制体上的运动状态。

在管道汇流口的流场中，特别是矩形断面有压管道汇流口中，如果将控制体拆分成等尺寸的矩形，就可以计算汇流口的流量、流速及其流场的分布情况。

结合实际管道汇流口流场的问题，针对矩形断面有压管道汇流口，搭建基于有限体积法（FVM）的数值计算系统，模拟管道汇流口的流场，以便更好地研究和分析管道汇流口的流场状态。

计算次级环流范文

计算次级环流范文

计算次级环流范文次级环流是指地球大气环流系统中的局地风系，常常受到地形条件和海陆分布的影响，形成气候和天气的变化。

以下是关于次级环流的详细计算。

次级环流的形成和演变通常与地形因素密切相关。

由于地球的自转和地球表面的不均匀性，形成了高气压和低气压系统，从而产生了次级环流。

地区与地区之间的温度差异、气压差异、风速和风向都是导致次级环流的因素。

次级环流可分为大尺度的环流和小尺度的环流。

大尺度环流是指以大范围和连续性为特点的环流形式，如副高和副低等系统。

小尺度环流是指尺度相对较小、范围较短的环流系统，如局地的山谷风、河谷风等。

次级环流的计算通常涉及到气象学中的物理参数，如温度、气压、湿度、风速等。

计算次级环流的方法主要有以下几种：1.静力平衡方程法：根据地球表面的气压分布和地球上空的气压分布，可以利用静力平衡方程计算次级环流的运动。

这种方法常常适用于大尺度环流系统。

2.地形解析法：通过分析地形的变化和对流体力学方程的应用，可以计算次级环流在地形上的影响。

这种方法常常适用于小尺度环流系统。

3.数值模拟法：通过建立气象模型，利用计算机进行数值模拟，可以模拟次级环流的形成和演变过程。

这种方法常常适用于大尺度和小尺度环流系统。

4.统计分析法：通过对大量实测资料的分析和统计，可以计算次级环流的特征和变化规律。

这种方法常常适用于长时间序列的研究。

以上方法各有优缺点，适用于不同的研究对象和研究目的。

在实际计算中，常常需要综合运用多种方法，并结合实测资料进行验证和调整。

在计算次级环流时，需要考虑的因素较多，包括地形条件、海陆分布、大气压力和温度分布、湿度和降水等。

这些因素相互作用，共同决定了次级环流的形成和演变过程。

最后，需要指出的是，次级环流是一个复杂而多变的系统，受到多种因素的影响，包括地球自转、地形、海洋热力等。

在计算和研究次级环流时，需要综合考虑这些因素，并进行合理的模型建立和数据分析。

第七章1-有关流场的几个基本概念

第七章1-有关流场的几个基本概念

三. 流管和流量
流管 流线
• 在流场中，取一条
不与流线重合的封闭 曲线L，在同一时刻 过 L上每一点作流 线，由这些流线围成 的 管 状 曲 面 称 为 流管。
L
•根据流管的定义易知，在
对应瞬时，流体不可能通 过流管表面流出或流入。
• 与流线一样，
流管是瞬时概念。
••过水断面 过水断面
与流动方向正交的流管的横断面
是 否 接 近 均 匀 流？
是
渐 变 流
流线虽不平行，但夹角较小； 流线虽有弯曲，但曲率较小。
否
急 变 流
流线间夹角较大； 流线弯曲的曲率较大。
•• 渐变流和急变流是工程意义上对流动是否符合均匀流条件的 渐变流和急变流是工程意义上对流动是否符合均匀流条件的
划分，两者之间没有明显的、确定的界限，需要根据实际情况 划分，两者之间没有明显的、确定的界限，需要根据实际情况 来判定 来判定
• 二维流动 • 二维流动
流场与某一空间 坐标变量无关，且沿该坐标 方向无速度分量的流动。
大展弦比机翼绕流
直角系中的平面流 直角系中的平面流 动： 动： ⎧ux = ux ( x , y , t ) uz = 0 ⎪ ∂ ⎨uy = uy ( x , y , t ) =0 ⎪u = 0 ∂z为相互平行的直线，而非均匀流的 均匀流的流线必为相互平行的直线，而非均匀流的 流线要么是曲线，要么是不相平行的直线。
• 例如，以下的流动
是均匀流：
y a o z ux x
u x = u x ( y) , u y = uz = 0
• 应注意将均匀流与完全不随空间位置而变的等速直线流动
急变流示意图
五. 流动按空间维数的分类
一维流动 二维流动 三维流动

水力学第十章流场理论基础赵

水力学第十章流场理论基础赵

角变形和旋转运动
dθ 1 ∂u x ∂u z ⎫ θy = = ( + )⎪ dt 2 ∂z ∂x ⎪ ⎪ 1 ∂u z ∂u y θx = ( + ) ⎬ ∂z 2 ∂y ⎪ ⎪ 1 ∂u y ∂u x θz = ( + ) ⎪ ∂y 2 ∂x ⎭
有涡流与无涡流
当流体微团存在绕自身轴旋转的运动，称为 有涡流；当流体微团不存在绕自身轴旋转的 运动，称为无涡流。 v 对于无涡流 ω = 0 ，即 ω x = ω y = ω z = 0
∂ρ dx ∂u x dx (ρ − )(u x − )dtdydz ∂x 2 ∂x 2 Эuz d z
Эz 2
ux Эux d x Эx 2 G
C dz E
O
B uz Эuz d z Эz 2 F dx
o x
A
dy
∂ ( ρu x ) − dxdydzdt ∂x
yy uy Эuy d Эy Nhomakorabea2∂ρ dx ∂u x dx (ρ + )(u x + )dtdydz ∂x 2 ∂x 2
液体微团运动的基本形式
设想在运动的液体中取一微分平行六面体 ABCDEFJH。由于此微团上各点的速度不 同，当经历微分时段dt后运动至新的位置 时，该微团的形状和大小都将发生变化，变 ′ 为 A′B′C′D′E′F ′J ′H，若对运动进行分解，以上 这一变化过程，可以看成是在dt微分时段经 历以下四种基本运动的结果，即平移、线变 形、角变形和旋转。
z
∂u z uz + dz ∂z ∂u x ux + dz ∂z ∂u z ∂u z dz + dx uz + ∂z ∂x
∂u x ∂u x ux + dx + dz ∂x ∂z

河道二维非恒定流场的一种计算方法

河道二维非恒定流场的一种计算方法

河道二维非恒定流场的一种计算方法
张莉
【期刊名称】《广西轻工业》
【年(卷),期】2007(023)002
【摘要】本文提出了模拟天然河道二维非恒定流场的一种通用的数学模型,该方法具有收敛快、稳定性好、运算量小、内存少的特点.用该方法对长江北支天生港到白茆河口段的流场进行了模拟,将其计算水位值与实测水位值进行比较,发现结果吻合良好.
【总页数】3页(P65-66,78)
【作者】张莉
【作者单位】河海大学理学院,江苏,南京,210098
【正文语种】中文
【中图分类】TV133.2
【相关文献】
1.河道二维非恒定流的数值模拟
2.二维非恒定流模型在大辽河河道现状行洪能力分析中的应用
3.非棱柱形河道恒定非均匀流水面线的一种数值计算方法
4.河道二维非恒定流场计算方法研究
5.基于二维非恒定流分析法在河道采砂中的应用
因版权原因，仅展示原文概要，查看原文内容请购买。

  1. 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

2 2 x 2 x 2 1 2 y x x y Re
（1）
（2）
0 ， 在腔体四边
v 0 ，在 AD 和 BC 上 x
u 0 ，在 AB 和 DC 上 y
其中，Re 为雷诺数。
Page 2 of 10
姓名：刘广
学号：11309018
日期：2014 年 05 月 23 号
Page 6 of 10
中山大学工学院、理论与应用力学刘广编制
实验编号及题目：实验六
方腔环流
《计算流体力学试验》课程实验报告纸
院系：工学院
E2=ABS(TEMPM1(I,1)-M1(I,1)) E1=ABS(TEMPM(I,1)-M(I,1)) IF(MAX(E1,E2)>E0)THEN E0=MAX(E1,E2) END IF END DO DO I=2,X M1(I,X+1)=-2*(TEMPM(I,X)-TEMPM(I,X+1)+H)/(H*H) E2=ABS(TEMPM1(I,X+1)-M1(I,X+1)) E1=ABS(TEMPM(I,X+1)-M(I,X+1)) IF(MAX(E1,E2)>E0)THEN E0=MAX(E1,E2) END IF END DO !层之间交换 TEMPM(:,:)=M(:,:) TEMPM1(:,:)=M1(:,:) END DO DO I=1,X+1,1 DO J=1,X+1,1 WRITE(8,*)I,',',J,',',M(I,J) END DO END DO DO I=1,X+1,1 DO J=1,X+1,1 WRITE(10,*)I,',',J,',',M1(I,J) END DO END DO END SUBROUTINE
i , j 1 )
Page 3 of 10
中山大学工学院、理论与应用力学刘广编制
实验编号及题目：实验六
方腔环流
《计算流体力学试验》课程实验报告纸
Байду номын сангаас院系：工学院
1 in ,j
姓名：刘广
学号：11309018
日期：2014 年 05 月 23 号

4
n n n { in 1, j i 1, j i , j 1 i , j 1
姓名：刘广
学号：11309018
日期：2014 年 05 月 23 号
程序执行时，我输入的是10E-3作为计算精度，对于程序的说明程序注释已经说的很清楚， 这里不再赘述。最终会在目录下生成OUTPUT.txt 。 四：实验结果 最终会在目录下生成OUTPUT.txt 。然后我们可以用Matlab软件对计算出来的数据进一步整 理，画成云图和画出流线，我们做一下说明。
中山大学工学院计算流体力学实验报告
实验名称：方腔环流的流场计算 姓 名： 刘广 刘广 11309018 詹杰民
参与组员： 学 号：
任课教师：
学科专业：工学院 理论与应用力学
中山大学 2014 年 05 月 23 日
中山大学工学院、理论与应用力学刘广编制
实验编号及题目：实验六
方腔环流
《计算流体力学试验》课程实验报告纸
中山大学工学院、理论与应用力学刘广编制
实验编号及题目：实验六
方腔环流
《计算流体力学试验》课程实验报告纸
院系：工学院 姓名：刘广 学号：11309018 日期：2014 年 05 月 23 号
方腔环流示意图
我们取网格做如下说明，取正方形网格，网格数量由用户输入，但是建议不超过 200，因 为超过 200 计算量会急剧增大，雷诺数也是由用户输入，建议不超过 5000，因为超过 5000 流 场已经呈现出非线性，开始出现偏差，雷诺数超过 10000 甚至根本不能算出结果，这是因为强 非线性情况下描述流场行为的方程已经不能做如上简化， 同样的， 在计算过程中我们采用超松 弛迭代法，下面对离散形式做一下说明： 我们用二阶精度的差商代替微商做以下说明，得
( i 1, j 2 i , j i 1, j ) / h 2 ( i , j 1 2 i , j i , j 1 ) / h 2 i , j ( i , j 1 i , j 1 )(
2 i 1, j
i 1, j )
4h 4h 2 1 i 1, j 2 i , j i 1, j i , j 1 2 i , j i , j 1 ( ) Re h2 h2
sn 1 0, 在四边
2( sn * sn ) , 在两侧和底边 h2 2( sn * sn h) n 1 s , 在顶边 h2
n 1 s
三：实验步骤: 学会对MS-FORTRAN的基本操作。 用Fortran编写程序计算方腔环流流场计算源代码，我们可 以写如下Fortran90程序： 主程序如下：
! Author : GUANG_LIU * owenyaa@ * ! Created Time : X+14-04-22 19:50 ! Last Revised : GUANG_LIU ,X+14-04-22 ! Remark : 以下实验结果中均为网格19×19，步长0.001，松弛因子0.5，迭代精度0.000001。 !Re=100 PROGRAM MAIN IMPLICIT NONE INTEGER::I WRITE(*,*)'请选择计算方法,流函数涡量法输入1,速度压力法输入2' READ(*,*)I IF (I==1)THEN CALL StreamFunctionVorticity() ELSE IF (I==2)THEN CALL PressureVelocity() Re=1000
院系：工学院 实验人姓名：刘广 日期：2014 年 05 月 23 号 专业：理论与应用力学 参加人姓名：刘广 温度：18° C 学号：11309018 年级：2011 级
实验六：方腔环流 一、 实验目的 1、在正方形水腔当中布满液体，水腔三方为墙壁，上方为盖板，盖板以速度 1 向右匀速 运动，当水腔内部液体达到平衡时，计算水箱内的流场。 2、计算全区域的流函数分布，并绘制流场图像。 二：实验仪器与设备: ① 装有 Fortran 的计算机 三、实验原理 1.方程与边界条件 用流函数涡量法， 和 满足下列无量纲定常方程与边界条件 1台

Re n n [( in, j 1 in, j 1 )( in in1, j in1, j )( in 1, j i 1, j ) ( , j 1 i , j 1 )]} 4 (1 ) in ,j
边界值的点需要特殊处理，然后由内点差分格式顺序扫描，一次次迭代直到满足用户输入的 精度为止。边界条件如下：
Page 4 of 10
中山大学工学院、理论与应用力学刘广编制
实验编号及题目：实验六
方腔环流
《计算流体力学试验》课程实验报告纸
院系：工学院
END IF WRITE(*,*)'请在目录下查看结果' PAUSE STOP END
姓名：刘广
学号：11309018
日期：2014 年 05 月 23 号
以下为流函数涡量法的子程序，
松弛因子 代入， i , j
n 1 2 n n n n ( in in, j , j h i 1, j i 1, j i , j 1 i , j 1 ) / 4 (1 )

( i 1, j i 1, j )(
i , j 1
Page 5 of 10
中山大学工学院、理论与应用力学刘广编制
实验编号及题目：实验六
方腔环流
《计算流体力学试验》课程实验报告纸
院系：工学院
E0=0 !循环中进行判断，但是这种方法不好，可以采用-1，0，1标记法，循环区间为大区间，计算之前判断，-1 为外点，CYCLE掉，1为内点，用迭代公式，0为边界点，值不变 DO I=2,X,1 DO J=2,X,1 M(I,J)=(1-W)*TEMPM(I,J)+(W/4)*(M(I-1,J)+TEMPM(I+1,J)+M(I,J1)+TEMPM(I,J+1)+H*H*TEMPM1(I,J)) M1(I,J)=(W/4)*(TEMPM1(I+1,J)+M1(I-1,J)+TEMPM1(I,J+1)+M1(I,J-1)+(RE/4)*((TEMPM(I,J+1)M(I,J-1))*(TEMPM1(I+1,J)-M1(I-1,J))-(TEMPM(I+1,J)-M(I-1,J))*(TEMPM1(I,J+1)-M1(I,J-1))))+(1W)*TEMPM1(I,J) E1=ABS(TEMPM(I,J)-M(I,J)) E2=ABS(TEMPM1(I,J)-M1(I,J)) IF(MAX(E1,E2)>E0)THEN E0=MAX(E1,E2) END IF END DO END DO DO J=1,X+1!左边 M1(1,J)=-2*(TEMPM(2,J)-TEMPM(1,J))/(H*H) E2=ABS(TEMPM1(1,J)-M1(1,J)) E1=ABS(TEMPM(1,J)-M(1,J)) IF(MAX(E1,E2)>E0)THEN E0=MAX(E1,E2) END IF END DO DO J=1,X+1!右边 M1(X+1,J)=-2*(TEMPM(X,J)-TEMPM(X+1,J))/(H*H) E2=ABS(TEMPM1(X+1,J)-M1(X+1,J)) E1=ABS(TEMPM(X+1,J)-M(X+1,J)) IF(MAX(E1,E2)>E0)THEN E0=MAX(E1,E2) END IF END DO DO I=2,X!底边 M1(I,1)=-2*(TEMPM(I,2)-TEMPM(I,1))/(H*H)