Fortran语言涡量流函数法中心差分格式的二维方腔顶盖驱动计算.docx

涡量流函数法中心差分格式的二维方腔顶盖驱动计算本题是关于粘性流体方腔顶盖驱动的问题。采用涡量流函数法，在均分网格下，用中心差分格式进行计算，结果与文献中所采用的其他方法和格式进行比较，认为中心差分格式符合计算的精度，但是其明显缺点是计算过程不稳定。1、参数的无量纲化令其中，由涡量的定义式：将速度的导数项无量纲化，并将上述定义式代入，即令则涡量的无量纲形式可以写为：再由流函数的定义式：同理，令可以得到流函数的无量纲表达式为：流函数边界条件涡量的边界条件（Thom公式）+2、方程的无量纲处理过程2.1、非守恒型方程的处理（1）将以上假设的各式代入到非守恒型方程组：得到以下无量纲形式的方程：①涡量控制方程②流函数控制方程（2）采用有限差分法离散非守恒型涡量的无量纲控制方程：整理得到：简化为：写成一般形式是: 其中；；；2.2、守恒型方程的处理（1）将假设的各个无量纲量代入守恒型涡量方程得到其无量纲形式为：（2）采用有限容积法离散守恒型涡量的无量纲控制方程：对于扩散项，有对于对流项，有对于界面上的取值，采用如下形式：；；；；其中，界面上涡量的插值采用中心差分格式，令；；；所以，关于界面上速度的插值，可以用上一层次的流函数来表示，形式如下：由以上各式，整理得到守恒型涡量控制方程无量纲形式的离散方程为：=综上，可整理得到（同位的均分网格）：写成一般形式为：其中，2.3、对于流函数方程，离散可得到：整理并写成一般形式是：其中，；；；；并且，采用有限容积法离散的结果与采用这里的方法得到的结果是相同的。2.4、将压力的泊松方程进行无量纲化并离散令代入压力泊松方程，得到其无量纲形式为：边界上的压力满足N-S方程，同样将其无量纲化，并与上式组成压力的封闭方程组，X方向上： Y方向上：将以上三式所组成的方程组离散，得到：内点压力计算式为：左边界压力为：右边界压力为：上边界压力为：下边界压力为：由此可得到压力的封闭方程组，通过求得的速度场来计算压力场。3、程序编写思路：4、计算结果4.1、两种离散方式的计算时间的比较：（EPS=1E-6, i 代表松弛因子）RE有限容积法与有限差分法计算时间的比较/s网格数100(i=0.3)1000(i=0.1)有限容积法FVM有限差分法FDM有限容积法FVM有限差分法FDM20×2037.0287.481973.202132.5640×40509.17517.584415.804713.2080×806889.286906.0914359.8015536.744.2、雷诺数为100时各物理量场的分布图形（1）有限差分法的情况下：A 流函数等值线（图1a/b/c）N=20×20N=40×40N=80×80B 涡量等值线（图2a/b/c）N=20×20N=40×40N=80×80C 压力等值线（图3a/b/c）N=20×20N=40×40 N=80×80D 中心线的速度分布（图4a/b）（2）有限容积法的情况下：A 流函数等值线（图5a/b/c）N=20×20N=40×40 N=80×80B 涡量等值线（图6a/b/c）N=20×20N=40×40 N=80×80C 压力等值线（图7a/b/c）N=20×20 N=40×40 N=80×80D 中心线的速度分布（图8a/b）4.3、雷诺数为1000时各物理量场的分布图形（1）有限差分法的情况下：A 流函数等值线（图9a/b/c）N=20×20 N=40×40 N=80×80B 涡量等值线（图10a/b/c）N=20×20 N=40×40 N=80×80C 压力等值线（图11a/b/c）N=20×20 N=40×40 N=80×80D 中心线的速度分布（图12a/b/c）（2）有限容积法的情况下：A 流函数等值线（图13a/b/c）N=20×20 N=40×40 N=80×80B 涡量等值线（图14a/b/c）N=20×20 N=40×40 N=80×80C 压力等值线（图15a/b/c）N=20×20 N=40×40 N=80×80D 中心线的速度分布（图16a/b/c）4.4、文献[1]上的结果（1）雷诺数100时候的流函数和涡量图形（图17a/b）（2）雷诺数1000时的流函数和涡量图形（图1