求解Euler方程的区域分解方法与并行算法.doc

[文章编号] (20000420360207 求解Euler 方程的区域分解方法与并行算法 吕晓斌1,　兰黔章2,　朱自强3 (11中国科学院力学研究所,北京　100080;　21上海交通大学动力系,上海　200030; 31北京航空航天大学流体所,北京　100083 [摘　要]　将复杂形状区域划分成多块子区域,研究发展了一种多块区域之间迎风守恒型的内边界耦合方法,实现相邻子区域解的光滑过渡,使多区耦合得到总体流场的数值解。对二维翼型跨音速流动和圆弧形隆起物超音速流动等进行了分区数值计算,并将计算结果与单区计算结果和实验结果作了比较。并行分区计算引入“先进先出”的同步控制等待机制,实现了高效率并行计算,还分析了影响并行效率的主要因素。 [关键词]　计算流体力学;区域分解方法;并行处理;Euler 方程[中图分类号]　O35;V211　　　　[文献标识码]　A [收稿日期]1999204220 [基金项目]国家自然科学基金资助项目 [作者简介]吕晓斌(1965~,男,山东,博士后,从事计算空气动力学并行计算方面的研究.0　引言复杂外形无粘性可压缩流场的数值计算,主要涉及到偏微分方程组的求解。该类问题的特征波速是有限的,如果把复杂几何区域分解成若干块子区域,子区域之间通过一定的耦合条件进行信息交换,以实现相邻子区域解的光滑过渡,复杂区域总体流场可实现并行求解。因而分区方法是一个有效地并行处理复杂区域流场的数值方法基础[1]。 本文结合矢通量分裂数值计算方法和信息传播的特征理论,采用完全对接网格,研究发展了一种迎风型内边界通量守恒的耦合条件,应用于二维跨音速翼型绕流、圆弧形隆起物超音速流动激波反射现象的分区并行计算。 1　控制方程和数值方法 在直角坐标系(x ,y 中,二维非定常Euler 方程组的微分形式为5U t +5F (U x +5G (U y =0(1式中,U 是守恒矢变量,F (U 、G (U 是矢通量 U =ρρu ρv E ,　F (U =ρu ρu 2+p ρuv (E +p u ,　G (U =ρv ρuv ρv 2+p (E +p v (2第17卷第4期 2000年7月计　　算　　物　　理　CHIN ESE J OU RNAL OF COMPU TA TIONAL PH YSICS 　Vol.17,No.4J ul.,2000 　 理想气体状态方程为p =(γ-1[E -12ρ(u 2+v 2](3 其中,ρ为密度,u 、v 为速度分量,p 为压力,γ为比热比,空气γ=114,E 为单位体积的总能。 用有限体积方法在每一个网格单元控制体上进行积分,计算变量定义在网格中心点上,在边界上计算通量,得到一组空间离散方程V i ,j d U i ,j =-d t F i +1/2,j S i +1/2,j - F i -1/2,j S i -1/2,j + F i ,j +1/2S i ,j +1/2- F i ,j -1/2S i ,j -1/2(4 式中 F =n x F (U +n y G (U (5 其中,S 是控制体V i ,j 的边界面,n 为边界面的外法向单位矢量。如果网格是足够光滑的,边界物理量由二阶插值公式计算,则方程式(4是二阶的。 本文按Van Leer 矢通量分裂方法[3]计算控制体边界面对流矢通量。对二维流动情况,根据边界面的法向流动速度进行矢通量分裂[4]。公式(4中 F i +1/2,j , F i ,j +1/2等不能直接应用Van Leer 通量分裂公式得到 F ±i +1/2,j , F ±i ,j +1/2等,必须先用一个变换T 得到局部坐标系中法向通量,然后按一维分裂方法对法向通量作分裂,分裂后再用一个反变换T -1,将通量变换回来才能得到 F ±i +1/2,j 等。 亚音速流动(|M n |1时,边界面法向矢通量分裂公式为 H n =±ρc 4 (1±M n 2 h ±1c[(γ-1M n ±2]/γ h ±1cM t h ±1(E +p /ρ (6式中h ±1表示第一项,即h ±1=±(ρc/4(1±M n 2,M n 为边界面的法向马赫数,M t 表示与界 面法向垂直的马赫数分量。超音速流动时,|M n |Ε1,此时法向特征速度(q n +c ,q n ,q n -c 全为正或全为负。由信息传播理论,边界面中点物理量全部由上游一侧插值计算,固壁物面边界上满足壁面无穿透条件,即物面法向速度等于零。远场根据Riemann 不变量分析,采用无反射边界条件。时间积分方法采用分区隐式因式分解(AF 方法[5]。2　内边界耦合条件 图1　(a 重叠网格　　　　(b 对接网格Fig.1　(a overlapping grid 　　 (b patched grid