2024年3月26日发(作者：东方莎莉)

第４５卷第１期　

西安交通大　学　学报　

Ｖｏ１．４５　Ｎｏ．１　

２０１１年１月　

ＪＯＵＲＮＡＬ　ＯＦ　ＸＩ　ＡＮ　ＪＩＡｏＴＯＮＧ　ＵＮＩＶＥＲＳＩＴＹ　

Ｊａｎ．２Ｏｌ１　

具有运动边界不可压缩黏性流动的　

ＣＢＳ有限元解法　

孙旭，张家忠　

（西安交通大学能源与动力工程学院，７１００４９，西安）　

摘要：基于特征分裂（ＣＢＳ）算法，提出了一种适用于具有运动边界的不可压缩黏性流动问题的有　

限元解法，给出了不可压缩黏性流动问题的控制方程，并推导了其无量纲形式．将传统ＣＢＳ算法与　

动网格技术相结合，提出了一种动边界流动问题的求解流程，并对所采用的基于弹簧近似的动网格　

方法进行了说明．采用该方法对静水和均匀来流中的简谐振动圆柱绕流问题进行了求解，计算所得　

圆柱受力、圆柱表面压力和涡量分布以及流场速度分布等结果均与已有数值或实验结果相吻合，初　

步验证了所提方法的可行性．　

关键词：特征分裂算法；有限元法；不可压缩流动；动边界　

中图分类号：Ｏ２４１．８２文献标志码：Ａ文章编号：０２５３—９８７Ｘ（２０１１）０１—００９９—０６　

Ａ　Ｃｈａｒａｃｔｅｒｉｓｔｉｃ　Ｂａｓｅｄ　Ｓｐｌｉｔ・－ＦＥＭ　Ｓｃｈｅｍｅ　ｆｏｒ　Ｉｎｃｏｍｐｒｅｓｓｉｂｌｅ　Ｖｉｓｃｏｕｓ　

Ｆｌｏｗ　ｗｉｔｈ　Ｍｏｖｉｎｇ　Ｂｏｕｎｄａｒｉｅｓ　

ＳＵＮ　Ｘｕ，ＺＨＡＮＧ　Ｊ　ｉａｚｈｏｎｇ　

（Ｓｃｈｏｏｌ　ｏｆ　Ｅｎｅｒｇｙ　ａｎｄ　Ｐｏｗｅｒ　Ｅｎｇｉｎｅｅｒｉｎｇ，Ｘｉ　ａｎ　Ｊｉａｏｔｏｎｇ　Ｕｎｉｖｅｒｓｉｔｙ，Ｘｉ　ａｎ　７１００４９，Ｃｈｉｎａ）　

Ａｂｓｔｒａｃｔ：Ａ　ｃｈａｒａｃｔｅｒｉｓｔｉｃ　ｂａｓｅｄ　ｓｐｌｉｔ（ＣＢＳ）一ＦＥＭ　ｓｃｈｅｍｅ　ｆｏｒ　ｔｈｅ　ｉｎｃｏｍｐｒｅｓｓｉｂｌｅ　ｖｉｓｃｏｕｓ　ｆｌｏｗ　

ｗｉｔｈ　ｍｏｖｉｎｇ　ｂｏｕｎｄａｒｉｅｓ　ｗａｓ　ｐｒｏｐｏｓｅｄ，ｉｎｃｌｕｄｉｎｇ　ｔｈｅ　ｇｏｖｅｒｎｉｎｇ　ｅｑｕａｔｉｏｎｓ　ｅｓｔａｂｌｉｓｈｅｄ　ａｎｄ　ｔｈｅ　ｃｏｒ—　

ｒｅｓｐｏｎｄｉｎｇ　ｎｏｎ－ｄｉｍｅｎｓｉｏｎａｌ　ｅｑｕａｔｉｏｎｓ　ｄｅｄｕｃｅｄ．Ｍｏｒｅｏｖｅｒ，ｃｏｍｂｉｎｅｄ　ｔｈｅ　ｃｏｎｖｅｎｔｉｏｎａｌ　ＣＢＳ　ａｌｇｏ—　

ｒｉｔｈｍ　ｗｉｔｈ　ｔｈｅ　ｍｅｓｈ　ｍｏｖｉｎｇ　ｔｅｃｈｎｉｑｕｅ，ａ　ｓｏｌｕｔｉｏｎ　ｐｒｏｃｅｄｕｒｅ　ｉｎ　Ｅｕｌｅｒｉａｎ　ｃｏｏｒｄｉｎａｔｅｓ　ｗａｓ　ｏｂｔａｉｎｅｄ　

ａｎｄ　ｔｈｅ　ｓｐｒｉｎｇ　ａｎａｌｏｇｙ－ｂａｓｅｄ　ｍｅｓｈ　ｍｏｖｉｎｇ　ｔｅｃｈｎｉｑｕｅ　ｗａｓ　ｄｉｓｃｕｓｓｅｄ．Ｉｎ　ａｄｄｉｔｉｏｎ，ｔｈｅ　ｐｒｏｐｏｓｅｄ　

ｍｅｔｈｏｄ　ｗａｓ　ａｐｐｌｉｅｄ　ｔｏ　ｔｗｏ　ｆｌｏｗ　ｐｒｏｂｌｅｍｓ，ｉ。ｅ．ｔｈｅ　ｆｌｏｗ　ｉｎｄｕｃｅｄ　ｂｙ　ｔｈｅ　ｈａｒｍｏｎｉｃ　ｏｓｃｉｌｌａｔｉｎｇ　ｃｉｒ—　

ｃｕｌａｒ　ｃｙｌｉｎｄｅｒ　ｉｎ　ｗａｔｅｒ　ａｔ　ｒｅｓｔ，ａｎｄ　ｔｈｅ　ｕｎｉｆｏｒｍ　ｆｌｏｗ　ａｒｏｕｎｄ　ａ　ｔｒａｎｓｖｅｒｓｅｌｙ　ｏｓｃｉｌｌａｔｉｎｇ　ｃｉｒｃｕｌａｒ　ｃｙｌ—　

ｉｎｄｅｒ．Ｔｈｅ　ｎｕｍｅｒｉｃａｌ　ｒｅｓｕｌｔｓ　ｓｕｃｈ　ａｓ　ｔｈｅ　ｄｒａｇ　ａｎｄ　ｌｉｆｔ　ｆｏｒｃｅｓ　ｏｆ　ｔｈｅ　ｃｙｌｉｎｄｅｒ，ｔｈｅ　ｐｒｅｓｓｕｒｅ　ａｎｄ　ｖｏｒ—　

ｔｉｃｉｔｙ　ｏｎ　ｔｈｅ　ｃｙｌｉｎｄｅｒ　ｓｕｒｆａｃｅ，ａｎｄ　ｔｈｅ　ｖｅｌｏｃｉｔｙ　ｄｉｓｔｒｉｂｕｔｉｏｎ　ｏｆ　ｔｈｅ　ｆｌｏｗ　ｆｉｅｌｄ　ａｇｒｅｅ　ｗｅｌｌ　ｗｉｔｈ　ｔｈｅ　

ｎｕｍｅｒｉｃａｌ　ａｎｄ　ｅｘｐｅｒｉｍｅｎｔａｌ　ｒｅｓｕｌｔｓ　ｉｎ　ｌｉｔｅｒａｔｕｒｅｓ．　

Ｋｅｙｗｏｒｄｓ：ｃｈａｒａｃｔｅｒｉｓｔｉｃ　ｂａｓｅｄ　ｓｐｌｉｔ；ｆｉｎｉｔｅ　ｅｌｅｍｅｎｔ　ｍｅｔｈｏｄ；ｉｎｃｏｍｐｒｅｓｓｉｂｌｅ　ｆｌｏｗ；ｍｏｖｉｎｇ　

ｂｏｕｎｄａｒｙ　

基本特征分裂（Ｃｈａｒａｃｔｅｒｉｓｔｉｃ－Ｂａｓｅｄ　Ｓｐｌｉｔ，　

ｌｏｒ—Ｇａｌｅｒｋｉｎ有限元等方法相比，其求解过程更加简　

ＣＢＳ）有限元方法［　是一种重要的针对流动问题的　

有限元解法，与Ｐｅｔｒｏｖ—Ｇａｌｅｒｋｉｎ迎风有限元、Ｔａｙ—　

单．另外，对于不可压缩黏性流动问题，ＣＢＳ算法在　

求解过程中允许速度和压力采用同阶的基函数，从　

收稿日期：２０１０—０５—１０．　作者简介：孙旭（１９８２－－），男，博士生；张家忠（联系人），男，教授，博士生导师．　基金项目：新世纪　

优秀人才支持计划资助项目（ＮＣＥＴ一０７—０６８５）．　

网络出版时间：２０１０—１０—１９　网络出版地址：ｈｔｔｐ：／／ｗｗｗ．ｃｎｋｉ．ｎｅｔ／ｋｃｍｓ／ｄｅｔａｉｌ／６１—１０６９一ｔ．２０１０１０１９．２２１１．０１９．ｈｔｍｌ　

西安交通大学学报　第４５卷　

而避开了ＢａｂｕＳｋａ—Ｂｒｅｚｚｉ（ＢＢ）条件．因此，该方法　

生扭曲甚至重叠；如果规定计算网格固结于空间坐　

自Ｚｉｅｎｋｉｅｗｉｃｚ等人于２Ｏ世纪９Ｏ年代首次提出以　

标上，不随边界移动，这将给边界位置的捕捉带来　

来，已经广泛应用到了各种流动问题之中．　

困难．　

近年来，越来越多的学者开始采用ＣＢＳ算法来　采用动网格技术，可以令边界上的节点随边界　

求解具有运动边界的流动问题，如ＮｉｔｈｉａｒａｓｕＥ　］将　

运动，同时允许求解区域内部的节点发生独立的运　

ＣＢＳ算法与ＡＬＥ方法　相结合，对溃坝流动和水　

动，这样既可以准确描述边界的位置，也可以保证网　

面孤立波的传播问题进行了求解，Ｎｏｂａｒｉ和Ｎａｄｅ—　格的质量．本节将传统ＣＢＳ算法与动网格技术相　

ｒａｎ［４　采用Ｎｉｔｈｉａｒａｓｕ的方法求解了自由来流中振　

结合，提出了一种网格运动时方程（２）的求解方法．　

动圆柱的绕流问题，段庆林ｌ＿６ｊ将ＣＢＳ算法与ＡＬＥ　为了方便起见，本节公式在推导过程中省略了无量　

方法和无网格解法相结合，成功地求解了成型充填　

纲标示符号．　

过程中非等温非牛顿黏性流动问题．这些解法均是　２．１网格运动时的ＣＢＳ解法　

通过将ＣＢＳ与ＡＬＥ相结合实现的，求解过程比较　

根据ＣＢＳ算法求解思路＿】］，通过坐标变换可消　

复杂．另外，由于ＡＬＥ坐标下ＮＳ方程的对流项含　除控制方程（２）中动量方程的对流项，得到新坐标系　

有网格速度项，网格速度将会给计算的稳定性带来　

Ｘ　下的动量方程．就流动问题而言，随流体质点运　

新的影响．　

动的Ｌａｇｒａｎｇｅ坐标系可以作为ｘ　坐标系，其上的　

本文根据ＣＢＳ算法的特点，针对具有运动边界　

动量方程可以表示为　

的不可压缩黏性流动问题，提出了一种简单、有效的　

Ｄｕ

ｇ

——

一一　

＋　：丝＋　ｆ３）　

求解方法，并采用所提出方法对静水中和自由来流　

Ｄｔ　ａｚｌ　Ｒｅ　３ｘｉａ　１　。　

中的振动圆柱绕流问题进行了求解．　

式中第一项表示流体质点的速度变化率．　

１控制方程　

下面对式（３）中的时间导数项进行离散．图１中　

给出了时间间隔　内的流体质点和网格点的运动　

非守恒形式的黏性不可压缩ＮＳ方程为　

情况．取ｔ　时刻网格节点Ａ　处的流体质点为研究　

３ｕ；ｑ

ｌ（１）　

对象，并假设该流体质点在ｔ　时刻位于Ｂ点，Ａ　的　

空间坐标为　，Ａ和Ｂ两点间的距离为　，通过式　

（３）可得　

３ｔ　Ｕｊ　Ｏ

—

一

：ｒｊ

　十ｕ古差＋ｕ　十　Ｊ＋　ｊ　

式中：“　为流体速度；ｐ为密度；ｐ为压强；ｕ为运动　

ｆ一　Ｊ　

——　——

　！！　一　

黏性系数；　为体积力．设参考速度为Ｕ，特征长度　

Ａｔ　

为Ｌ，参考密度为　，定义　一　Ｕｉ，；一芒，　一　，　

ｃ　（一差＋去　＋　）”　＋　

一　（不可压缩流动中　一１．ｏ），　，于　＝　

臼（一ａ］？＿　ｑ＿１　３２

ｕｉ

…

ｋｆｉ）州ｉ　（４）　

而ＪｉＬ，

，

则式（１）的无量纲形式为　

式中：Ｏ≤　１，　一０时为显式格式，　＝＝＝１时为隐式　

誊　一

ａ　ｕ

　）　

格式．式（４）中Ａｔ内流体质点移动的距离　可以近　

似为　

￣＋ｕ，　一一　３

１：△　ｆ　—Ａｍ”Ｉ＋ｏ（Ａｔ　）　（５）　

ｆ

ｘｉ　Ｒｅ　　３，：．ｃｊ３

～

ｘ

　Ｉ（

ｊ１十　ｊ　

２

式中：Ｒｅ＝ＵＬ／ｕ．　

２数值求解方法　

点　

对于具有运动边界的流动问题，流场中部分边　

质点　

界的位置将随时间不断变化．采用数值方法进行求　

解时，如果在计算过程中规定动边界上的网格节点　

・　

随边界同步运动，将造成动边界附近的网格单元发　

图１　Ａｔ内流体质点和网格点的运动情况　

ｈｔｔｐ：／／ｗｗｗ．ｊｄｘｂ．ｃｎ　

第１期　孙旭，等：具有运动边界不可压缩黏性流动的ＣＢＳ有限元解法　

通过Ｔａｙｌｏｒ展开来表示　一６　处的变量，并将　

式（５）代人式（４）整理可得　

，　

Ｅｎｄｄｏ　

对于式（９）～式（１１）的有限元离散，详细过程见　

Ｊ　

“　Ｉ～ｕ７　Ｉ＝＝＝　

～

ｌ　

文献Ｅｌｉ．　

２．２　网格运动　

为了满足运动边界条件，前文给出的求解流程　

需要不断地进行网格调整．目前，针对三角形非结构　

化网格已有多种动网格生成方法，本文采用了文献　

△　（“　一　１　３２Ｕ　ｉ　ｎ＠　＋　）ｆ　＋　

／Ｎｔ２　

２　（　

一　

Ｏ２，Ｃ　ｉ＋　）　

＋（卜　）　

～　ａｚ　

（６）　

（７）

…　

３　ｆ＇＋ｅ

ａ　

．　

７　

ＵＩ－　

［７３中的边弹簧近似法．首先，将ｔ　时刻的计算网格　

等效为由相互连接的弹簧组成的系统．定义相邻节　

式中：ｕ７＋。、Ｐ　为网格节点Ａ　上的变量；ｕ７、Ｐ　和　

为Ａ处的变量．式（６）即为空间坐标　处动量方　

点ｆ和　间的弹簧平衡长度等于两节点的距离，定　

义弹簧刚度　

ｋｕ—ａｌ（（Ｘｌ—ｚＪ）　＋（　一　）　）“２　（１２）　

程的离散形式．需要说明的是，在式（６）的推导过程　

中，除了压强项之外，其他变量均采用显式格式．　

网格运动时，式（６）中的网格节点Ａ　在ｔ”时刻　

并不在Ａ处，Ａ处任意流动参数　可以由网格点Ｃ　

式中：以　、ａ　为控制参数，本文中ａ　一１．０，Ｃ／２一　

一

１．０．由于所有弹簧的实际长度均等于其平衡长　

度，所以弹簧系统处于平衡状态．在ｔ卅　时刻，运动　

边界上的网格节点发生移动，部分弹簧上的受力不　

求得，而　的Ｔａｙｌｏｒ展开式为　

　｝

　。

一　

　．

～　

ａ声”ｌ　

再为０，弹簧系统的受力平衡也将遭到破坏．此时，　

任意节点ｊ上的弹性合力为　

Ｆｓ一　ｋｕ（　～６』）　（１３）　

＋。（１　Ｉ　）（８）　

式中：　。表示在　内网格节点Ａ　移动的距离．　

对式（６）求解可以通过以下３步进行．　

（１）引人中间变量　

一　一

Ｊ＝１　

式中：Ｎ，表示与节点　相邻的节点个数；　和　，表　

示节点，和　的位移．在弹性力的作用下，网格节点　

△　（　一　１　

“　

＋　）＋　

将发生移动．根据受力平衡准则，当弹簧系统重新达　

到平衡时，任意节点上　

Ｆｓ一　ｋｕ（　，一　）一０　（１４）　

Ａｔ２

杀（Ｕｊ　＋　）”　（９）　

（２）求解ｎ＋１时刻的压强Ｐ　，即　

Ｊ：１　

蠡（　）：　蠡（“　一Ａｔ　一臼　）　

（１０）　

由式（１４）和边界上的位移便可求得每个节点上　

的位移．本文将采用预条件共轭梯度法［８　对式（１４）　

进行求解．根据每一个节点上的位移，ｔ　时刻网格　

节点所处的位置可以表示为　

一ｘ７＋　ｒ　（１５）　

（３）求解　＋１时刻的速度“　州，即　

一　一　

ｏｚ　己０　＼（差１０　／”　

＋　

Ｘｋ

（１１）　

３计算实例及分析　

为了验证算法的可行性，本文对静水中和均匀　

因此，根据上文的描述，对于具有运动边界的不　

可压缩黏性流动问题，可以按照如下流程进行求解：　

Ｄｏ　一１，总的计算时间步　

（１）根据ｔ　时刻的网格和边界运动方式，采用　

来流中的振动圆柱绕流流场进行了计算，这２个问　

题的求解过程均未考虑体积力的影响．　

３．１静水中的圆柱振动　

动网格技术得到ｔ　时刻网格和每一个节点的位　

移　２；　

图２给出了静水中简谐振动圆柱绕流的计算区　

域和无量纲边界条件．计算区域外边界的直径为　

１００Ｄ，其上速度为０，并设任意一点的压强为０．圆　

柱沿ｚ方向作简谐振动，其无量纲运动方程为　

３２（　）一一　ｓｉｎ（２￣ｆｔ）　（１６）　

（２）通过式（８）计算ｔ　时刻网格节点所在位置　

在ｔ　时刻的流动参数ｕ７、Ｐ”和－厂　；　

（３）通过式（９）求解　；　

（４）通过式（１Ｏ）求解　；　

（５）通过式（１１）求解ｕ７¨；　

式中：．２７＝ｘ／Ｄ表示圆柱在水平方向上的位移；　：　

ｈｔｔｐ：／／ＷＷＷ．ｊｄｘｂ．ＣＦＩ　

西安交通大学学报　第４５卷　

Ｄ为振动幅值；ｆ—Ｄｆ／　为简谐振动频率，　

为圆柱振动的最大速度，圆柱表面采用无滑移　

边界条件．定义Ｒｅｙｎｏｌｄ数和Ｋｅｕｌｅｇａｎ—Ｃａｒｐｅｎｔｅｒ　

数分别为　

Ｒｅ—

ｒ　Ｔ　

Ｕ＿堕　

ｎ　

；ＫＣ—Ｕ

ｒ　Ｔ　

—　ｍａｘ　（１７）　

ｕ　Ｕ　

式中：Ｄ为圆柱直径．对Ｒｅ＝１００、ＫＣ＝：＝５的情况进　

行了计算，通过式（１６）和式（１７）可得　＝０．７９５　８、　

厂＝０．２．　

图３给出了静水中简谐振动圆柱绕流的非结构　

计算网格，它是沿周向和径向对流场进行剖分后得　

到的．采用３种不同密度类型的网格对该圆柱绕流　

进行了计算，表ｌ给出了网格相关信息以及计算所　

得阻力系数Ｃｒ】．从表１可以看出，采用Ｍｅｓｈ一３型　

网格计算所得ＣＤ与文献［９—１０］的结果之间的偏差　

均小于１　．图４给出了采用Ｍｅｓｈ一３型网格．计算　

所得振动相位为１８０。、２１０。和３３０。时流场中部分位　

置的流速分布，图中实际流速Ｕ　、　由计算所得无量　

纲流速“、口换算得到．从图４可以看出，本文计算结　

果与文献［１Ｏ］的实验结果非常吻合．　

图４　Ｒｅ为１００、ＫＣ为５、网格数为１４０×１２０，计算、　

实验所得振动相位为１８０。、２１０。和３３０。时流场　

中部分位置的流速分布　

图２静水中振动圆柱绕流　图３静水中振动圆柱绕流　

的计算模型　的计算网络　

袁ｌ　Ｒｅ为１００、ＫＣ为５时的圆柱阻力系数ＣＤ　

图５均匀来流中振动圆柱绕流计算模型　

穿透边界条件；出口位于圆柱下游３０Ｄ处，其上压　

力为０　Ｐａ．圆柱沿垂直于来流方向作简谐振动，其　

无量纲运动方程为　

，～　～

～

～　

（￡）一　ｃｏｓ（２７【－厂ｔ）　（１８）　

式中：　一　为圆柱在垂直于来流方向上的位移；　

３．２均匀来流中的圆柱振动　

＝

Ａ／Ｄ为振幅；　＝ＤＦ／ｕ。。为无量纲振动频率，根　

对Ｒｅ＝１８５时均匀来流中振动圆柱绕流进行　

据圆柱固定时旋涡自然脱落频率　来选取．圆柱表　

了计算，图５给出了计算区域和无量纲边界条件．流　

面为无滑移边界条件．　

场进口位于圆柱上游１０Ｄ处，流动参数为自由来流　

为了获得Ｒｅ＝１８５时圆柱绕流的旋涡自然脱　

的参数；上下边界距离圆柱１０Ｄ处，边界条件为无　

落频率　，首先对圆柱固定时的流场进行计算，计　

ｈｔｔｐ：／／ＷＷＷ．ｊｄｘｂ．ｃｎ　

第１期　孙旭，等：具有运动边界不可压缩黏性流动的ＣＢＳ有限元解法　

算网格通过先将求解区域分成９部分然后再剖分得　

到，如图６所示．表２给出了不同网格下计算所得圆　

柱平均阻力系数Ｃ。、均方根升力系数ＣＬ　，以及根　

据升力系数曲线得出的旋涡脱落的斯特罗哈数Ｓ　．　

通过对比可以看出，本文计算所得升力、阻力系数比　

文献［９，ｌ１］中的结果略大，斯特罗哈数则非常吻合．　

通过Ｓ　可以求得Ｒｅ一１８５时无量纲旋涡自然脱落　

频率２ｏ一０．１９６．　

表２　Ｒｅ为１８５时均匀来流中固定圆柱绕流计算结果　

图７　圆柱运动到上方最大位置时表面压强系数分布　

ｏ／（。）　

图８　圆柱运动到上方最大位置时表面涡量的分布　

图６均匀来流中振动圆柱绕流计算网络　

采用Ｍｅｓｈ一２型网格对　一０．２、ｆ／Ｊ；一０．８，　

１．０，１．２时振动圆柱绕流流场进行了计算．图７和　

图８给出了３种振动频率下圆柱运动到上方最大位　

置时表面压强系数Ｃ　及无量纲涡量　的分布情　

况．图中０表示圆柱表面节点相对于圆柱中心的转　

角，其值在迎风面驻点处为０，并按顺时针方向增　

加．如图７所示，本文计算所得Ｃ　与文献Ｅ９］的结果　

非常吻合．图８中，在涡量绝对值较小的地方，本文　

计算结果与文献［９１ｔ￣常吻合，在涡量较大的地方二　

者略有偏差．在求解过程中，流速采用线性插值基函　

数，涡量（流速的一阶导数）在单元相邻处不连续，这　

是造成二者涡量偏差的主要原因．　

图９给出了不同振动频率下圆柱升力系数Ｃ　、　

阻力系数Ｃ。随时间的变化情况及圆柱运动到上方　

最大位置处的瞬时流场．从图９中可以看出：当圆柱　

的振动频率较低（７／－厂Ｏ一０．８，１．０）时，Ｃ。　、ＣＤ均随　

时间呈简谐形式周期变化，／的增加将改变Ｃ　、ＣＤ　

变化的幅值，但并不改变流场中流线的拓扑结构；当　

尹较高时，ｃ　、ｃ。呈现出完全不同的变化形式，流场　

中流线的拓扑结构也将发生变化．关于厂对振动圆　

柱绕流流场的影响，可参见文献［９，１１］．　

（ｃ）ｆ～／ｆ～ｏ一１．２　

图９不同振动频率下圆柱升力、阻力系数随时间的变化及　

圆柱运动到上方最大位置时的瞬时流场　

ｈｔｔｐ　ｆｆ　ＷＷＷ．　ｄｘｂ．ｃｎ　

１０４　西安交通大学学报　第４５卷　

４结论　

学与实践，２００２（２）：７－ｌ１．　

ＹＵＥ　Ｂａｏｚｅｎｇ．Ｉ　Ｉ　Ｘｉａｏｔｉａｎ．Ｓｔｕｄｙ　ｏｆ　ｔｈｅ　ＡＩ　Ｅ　ｆｉｎｉｔｅ　

本文给出了一种动网格上的ＣＢＳ计算流程，用　

于求解具有运动边界的不可压缩黏性流动问题．作　

为验证，对Ｒｅ一１００、ＫＣ＝５时静水中振动圆柱绕　

流和Ｒｅ＝１８５、

厂／／ｏ一０．８，１．０，１．２时均匀来流中　

ｅｌｅｍｅｎｔ　ｍｅｔｈｏｄ　ａｎｄ　ｉｔｓ　ａｐｐｌｉｃａｔｉｏｎｓ［Ｊ］．Ｍｅｃｈａｎｉｃｓ　ｉｎ　

Ｅｎｇｉｎｅｅｒｉｎｇ，２００２（２）：７－１１．　

［４］ＮＯＢＡＲＩ　Ｍ　Ｒ　Ｈ，ＮＡＤＥＲＡＮ　Ｈ．Ａ　ｎｕｍｅｒｉｃａｌ　ｓｔｕｄｙ　

ｏｆ　ｆｌｏｗ　ｐａｓｔ　ａ　ｃｙｌｉｎｄｅｒ　ｗｉｔｈ　ｃｒｏｓｓ　ｆｌｏｗ　ａｎｄ　ｉｎｌｉｎｅ　ｏｓｃｉｌ—　

的振动圆柱绕流问题进行了求解，所得结果与已有　

实验和数值结果非常吻合，说明了算法的可行性．与　

ｌａｔｉｏｎ［Ｊ］．Ｃｏｍｐｕｔ　Ｆｌｕｉｄｓ，２００６，３５（４）：３９３—４１５．　

［５］　ＮＯＢＡＲＩ　Ｍ　Ｒ　Ｈ，ＧＨＡＺＡＮＦＡＲＩＡＮ　Ｊ．Ａ　ｎｕｍｅｒｉｃａｌ　

ｉｎｖｅｓｔｉｇａｔｉｏｎ　Ｏｆ　ｆｌｕｉｄ　ｆｌｏｗ　ｏｖｅｒ　ａ　ｒｏｔａｔｉｎｇ　ｃｙｌｉｎｄｅｒ　ｗｉｔｈ　

Ｎｉｔｈｉａｒａｓｕ所提方法相比，本文计算仍在欧拉坐标　

系下完成，程序主体与传统的ＣＢＳ解法基本相同，　

求解过程更加简单．　

ｃｒｏｓｓ　ｆｌｏｗ　ｏｓｃｉｌｌａｔｉｏｎ［Ｊ］．Ｃｏｍｐｕｔ　Ｆｌｕｉｄｓ，２００９，３８　

（１０）：２０２６－２０３６．　

［６］段庆林．成型充填过程中非等温非牛顿粘性流动的　

ＡＬＥ有限元与无网格自适应耦合模拟［Ｄ］．大连：大连　

理工大学能源与动力学院，２００７．　

需要说明的是，该算法在计算时的每一步都要　

先根据边界运动求出ｔ”¨时刻的网格，因而只适用　

于动边界运动方式已知的情况．对于动边界运动未　

知的问题（如具有自由水面的流动问题），需要对已　

有的求解流程进行改进，这也是下一步需要研究的　

内容．　

参考文献：　

』ｌ　ｊ　ＺＩＥＮＫＩＥＷＩＣＺ　０　Ｃ，ＴＡＹＬＯＲ　Ｒ　Ｉ　Ｔｈｅ　ｆｉｎｉｔｅ　ｅｌｅ—　

［７］ＢＬＯＭ　Ｆ　Ｊ．Ｃｏｎｓｉｄｅｒａｔｉｏｎｓ　ｏｎ　ｔｈｅ　ｓｐｒｉｎｇ　ａｎａｌｏｇｙ［Ｊ］．　

Ｉｎｔ　Ｊ　Ｎｕｍｅｒ　Ｍｅｔｈ　Ｆｌｕｉｄｓ，２０００，３２（６）：６４７－６４８．　

Ｅ８］王勖成．有限单元法［Ｍ］．北京：清华大学出版社，　

２００３：２４５—２５０．　

［９］ＧＵＩＩ　ＭＩＮＥＡＵ　Ｅ，ＱＵＥＵＴＥＹ　Ｐ．Ａ　ｎｕｍｅｒｉｃａｌ　ｓｉｍｕ—　

ｌａｔｉｏｎ　ｏｆ　ｖｏｒｔｅｘ　ｓｈｅｄｄｉｎｇ　ｆｒｏｍ　ａｎ　ｏｓｃｉｌｌａｔｉｎｇ　ｃｉｒｃｕｌａｒ　

ｃｙｌｉｎｄｅｒ［Ｊ］．Ｊ　Ｆｌｕｉｄｓ　Ｓｔｒｕｃｔ，２００２，１６（６）：７７３—７９４．　

［ＩＯ］ＤＵＴＳＣＨ　Ｈ，ＤＵＲＳＴ　Ｆ，ＢＥＣＫＥＲ　Ｓ，ｅｔ　ａ１．Ｌｏｗ－　

Ｒｅｙｎｏｌｄｓ－ｎｕｍｂｅｒ　ｆｌｏｗ　ａｒｏｕｎｄ　ａｎ　ｏｓｃｉｌｌａｔｉｎｇ　ｃｉｒｃｕｌａｒ　

ｍｅｒｉｔ　ｍｅｔｈｏｄ：ｖｏｌｕｍｅ　３　ｆｌｕｉｄ　ｄｙｎａｍｉｃｓ（Ｍ　Ｉ．５ｔｈ　

ｅｄ．Ｏｘｆｏｒｄ，ＵＫ：Ｅｌｓｅｖｉｅｒ　Ｂｕｔｔｅｒｗｏｒｔｈ－Ｈｅｉｎｅｍａｎｎ，　

２００５．　

ｃｙｌｉｎｄｅｒ　ａｔ　ｌｏｗ　Ｋｅｕｌｅｇａｎ－Ｃａｒｐｅｎｔｅｒ　ｎｕｍｂｅｒｓ［Ｊ］．Ｊ　

Ｆｌｕｉｄ　Ｍｅｃｈ，１９９８，３６０：２４９—２７１．　

Ｅ２］　ＮＩＴＨＩＡＲＡＳＵ　Ｐ．Ａｎ　ａｒｂｉｔｒａｒｙ　Ｌａｇｒａｎｇｉａｎ　Ｅｕｌｅｒｉａｎ　

（ＡＩ　Ｅ）ｆｏｒｍｕｌａｔｉｏｎ　ｆｏｒ　ｆｒｅｅ　ＳＵｒｆａｃｅ　ｆｌｏｗｓ　ｕｓｉｎｇ　ｔｈｅ　

［１１］Ｉ．Ｕ　Ｘ　Ｙ，ＤＡＬＴＯＮ　Ｃ．Ｃａｌｃｕｌａｔｉｏｎ　ｏｆ　ｔｈｅ　ｔｉｍｉｎｇ　ｏｆ　

ｖｏｒｔｅｘ　ｆｏｒｍａｔｉｏｎ　ｆｒｏｍ　ａｎ　ｏｓｃｉｌｌａｔｉｎｇ　ｃｙｌｉｎｄｅｒ［Ｊ］．Ｊ　

Ｆ１ｕｉｄｓ　Ｓｔｒｕｃｔ，１９９６，１０（５）：５２７—５４１．　

ｃｈａｒａｃｔｅｒｉｓｔｉｃ—ｂａｓｅｄ　ｓｐｌｉｔ（（　Ｓ）ｓｃｈｅｍｅ＿Ｊ］．Ｉｎｔ　Ｊ　Ｎｕ—　

ｍｅｒ　Ｍｅｔｈ　Ｆｌｕｉｄｓ，２００５，４８（１２）：１４１５－１４２８．　

［３］岳宝增，李笑天．ＡＬＥ有限元方法研究及应用［Ｊ］．力　

（编辑苗凌）　

（上接第５７页）　

［１Ｏ］ＺＨＡＮＧ　Ｚ，ＸＩＥ　Ｙ，ＺＨＡＮＧ　Ｘ，ｅｔ　ａ１．Ａｎａｌｙｓｉｓ　ｏｆ　

ｐｉｓｔｏｎ　ｓｅｃｏｎｄａｒｙ　ｍｏｔｉｏｎ　ｃｏｎｓｉｄｅｒｉｎｇ　ｔｈｅ　ｖａｒｉａｔｉｏｎ　ｉｎ　

２Ｏ１Ｏ，４４（４）：０５５５—０５５９．　

［１２］王跃武，胡芄，陈则韶．ＨＦＣ一２２７ｅａ的ＭＢＷＲ状态方　

程［Ｊ］．西安交通大学学报，２００７，４１（１）：３７—４０．　

ＷＡＮＧ　Ｙｕｅｗｕ，ＨＵ　Ｗａｎ，ＣＨＥＮ　Ｚｅｓｈａｏ．ＭＢＷＲ　

ｔｈｅ　ｓｙｓｔｅｍ　ｉｎｅｒｔｉａ［Ｊ］．Ｐｒｏｃ　Ｉｎｓｔｎｍｅｃｈ　Ｅｎｇｒｓ：Ｐａｒｔ　Ｄ　

Ｊｏｕｒｎａｌ　ｏｆ　Ａｕｔｏｍｏｂｉｌｅ　Ｅｎｇｉｎｅｅｒｉｎｇ，２００９，２２３（４）：５４９　

５６３．　

ｅｑｕａｔｉｏｎ　ｏｆ　ｓｔａｔｅ　ｏｆ　ＨＦＣ　２２７ｅａ［Ｊ］．Ｊｏｕｒｎａｌ　ｏｆ　Ｘｉ　ａｎ　

Ｊｉａｏｔｏｎｇ　Ｕｎｉｖｅｒｓｉｔｙ，２００７，４１（１）：３７—４Ｏ．　

［１１］张执南，谢友柏．考虑连杆一曲轴变惯性影响的活塞二　

阶运动分析［Ｊ］．上海交通大学学报，２０１０，４４（４）：　

５５５—５５９．　

Ｅｌ３］盂祥慧．复杂产品时变性能设计和设计中的知识服务　

研究ＪＲ］．上海：上海交通大学机械与动力工程学　

院，２００８．　

ＺＨＡＮＧ　Ｚｈｉｎａｎ，ＸＩＥ　Ｙｏｕｂａｉ．Ｅｆｆｅｃｔｓ　ｏｆ　ｃｏｎｎｅｃｔｉｎｇ　

ｒｏｄ－ｃｒａｎｋｓｈａｆｔ　ｖａｒｉａｂｌｅ　ｉｎｅｒｔｉａ　ｏｎ　ｐｉｓｔｏｎ　ｓｅｃｏｎｄａｒｙ　ｍｏ—　

ｔｉｏｎ［Ｊ］．Ｊｏｕｒｎａｌ　ｏｆ　Ｓｈａｎｇｈａｉ　Ｊｉａｏｔｏｎｇ　Ｕｎｉｖｅｒｓｉｔｙ，　

（编辑管咏梅）　

ｈｔｔｐ：∥ｗｗｗ．ｊｄｘｂ．ｃｎ　

2024年3月26日发(作者：东方莎莉)

第４５卷第１期　

西安交通大　学　学报　

Ｖｏ１．４５　Ｎｏ．１　

２０１１年１月　

ＪＯＵＲＮＡＬ　ＯＦ　ＸＩ　ＡＮ　ＪＩＡｏＴＯＮＧ　ＵＮＩＶＥＲＳＩＴＹ　

Ｊａｎ．２Ｏｌ１　

具有运动边界不可压缩黏性流动的　

ＣＢＳ有限元解法　

孙旭，张家忠　

（西安交通大学能源与动力工程学院，７１００４９，西安）　

摘要：基于特征分裂（ＣＢＳ）算法，提出了一种适用于具有运动边界的不可压缩黏性流动问题的有　

限元解法，给出了不可压缩黏性流动问题的控制方程，并推导了其无量纲形式．将传统ＣＢＳ算法与　

动网格技术相结合，提出了一种动边界流动问题的求解流程，并对所采用的基于弹簧近似的动网格　

方法进行了说明．采用该方法对静水和均匀来流中的简谐振动圆柱绕流问题进行了求解，计算所得　

圆柱受力、圆柱表面压力和涡量分布以及流场速度分布等结果均与已有数值或实验结果相吻合，初　

步验证了所提方法的可行性．　

关键词：特征分裂算法；有限元法；不可压缩流动；动边界　

中图分类号：Ｏ２４１．８２文献标志码：Ａ文章编号：０２５３—９８７Ｘ（２０１１）０１—００９９—０６　

Ａ　Ｃｈａｒａｃｔｅｒｉｓｔｉｃ　Ｂａｓｅｄ　Ｓｐｌｉｔ・－ＦＥＭ　Ｓｃｈｅｍｅ　ｆｏｒ　Ｉｎｃｏｍｐｒｅｓｓｉｂｌｅ　Ｖｉｓｃｏｕｓ　

Ｆｌｏｗ　ｗｉｔｈ　Ｍｏｖｉｎｇ　Ｂｏｕｎｄａｒｉｅｓ　

ＳＵＮ　Ｘｕ，ＺＨＡＮＧ　Ｊ　ｉａｚｈｏｎｇ　

（Ｓｃｈｏｏｌ　ｏｆ　Ｅｎｅｒｇｙ　ａｎｄ　Ｐｏｗｅｒ　Ｅｎｇｉｎｅｅｒｉｎｇ，Ｘｉ　ａｎ　Ｊｉａｏｔｏｎｇ　Ｕｎｉｖｅｒｓｉｔｙ，Ｘｉ　ａｎ　７１００４９，Ｃｈｉｎａ）　

Ａｂｓｔｒａｃｔ：Ａ　ｃｈａｒａｃｔｅｒｉｓｔｉｃ　ｂａｓｅｄ　ｓｐｌｉｔ（ＣＢＳ）一ＦＥＭ　ｓｃｈｅｍｅ　ｆｏｒ　ｔｈｅ　ｉｎｃｏｍｐｒｅｓｓｉｂｌｅ　ｖｉｓｃｏｕｓ　ｆｌｏｗ　

ｗｉｔｈ　ｍｏｖｉｎｇ　ｂｏｕｎｄａｒｉｅｓ　ｗａｓ　ｐｒｏｐｏｓｅｄ，ｉｎｃｌｕｄｉｎｇ　ｔｈｅ　ｇｏｖｅｒｎｉｎｇ　ｅｑｕａｔｉｏｎｓ　ｅｓｔａｂｌｉｓｈｅｄ　ａｎｄ　ｔｈｅ　ｃｏｒ—　

ｒｅｓｐｏｎｄｉｎｇ　ｎｏｎ－ｄｉｍｅｎｓｉｏｎａｌ　ｅｑｕａｔｉｏｎｓ　ｄｅｄｕｃｅｄ．Ｍｏｒｅｏｖｅｒ，ｃｏｍｂｉｎｅｄ　ｔｈｅ　ｃｏｎｖｅｎｔｉｏｎａｌ　ＣＢＳ　ａｌｇｏ—　

ｒｉｔｈｍ　ｗｉｔｈ　ｔｈｅ　ｍｅｓｈ　ｍｏｖｉｎｇ　ｔｅｃｈｎｉｑｕｅ，ａ　ｓｏｌｕｔｉｏｎ　ｐｒｏｃｅｄｕｒｅ　ｉｎ　Ｅｕｌｅｒｉａｎ　ｃｏｏｒｄｉｎａｔｅｓ　ｗａｓ　ｏｂｔａｉｎｅｄ　

ａｎｄ　ｔｈｅ　ｓｐｒｉｎｇ　ａｎａｌｏｇｙ－ｂａｓｅｄ　ｍｅｓｈ　ｍｏｖｉｎｇ　ｔｅｃｈｎｉｑｕｅ　ｗａｓ　ｄｉｓｃｕｓｓｅｄ．Ｉｎ　ａｄｄｉｔｉｏｎ，ｔｈｅ　ｐｒｏｐｏｓｅｄ　

ｍｅｔｈｏｄ　ｗａｓ　ａｐｐｌｉｅｄ　ｔｏ　ｔｗｏ　ｆｌｏｗ　ｐｒｏｂｌｅｍｓ，ｉ。ｅ．ｔｈｅ　ｆｌｏｗ　ｉｎｄｕｃｅｄ　ｂｙ　ｔｈｅ　ｈａｒｍｏｎｉｃ　ｏｓｃｉｌｌａｔｉｎｇ　ｃｉｒ—　

ｃｕｌａｒ　ｃｙｌｉｎｄｅｒ　ｉｎ　ｗａｔｅｒ　ａｔ　ｒｅｓｔ，ａｎｄ　ｔｈｅ　ｕｎｉｆｏｒｍ　ｆｌｏｗ　ａｒｏｕｎｄ　ａ　ｔｒａｎｓｖｅｒｓｅｌｙ　ｏｓｃｉｌｌａｔｉｎｇ　ｃｉｒｃｕｌａｒ　ｃｙｌ—　

ｉｎｄｅｒ．Ｔｈｅ　ｎｕｍｅｒｉｃａｌ　ｒｅｓｕｌｔｓ　ｓｕｃｈ　ａｓ　ｔｈｅ　ｄｒａｇ　ａｎｄ　ｌｉｆｔ　ｆｏｒｃｅｓ　ｏｆ　ｔｈｅ　ｃｙｌｉｎｄｅｒ，ｔｈｅ　ｐｒｅｓｓｕｒｅ　ａｎｄ　ｖｏｒ—　

ｔｉｃｉｔｙ　ｏｎ　ｔｈｅ　ｃｙｌｉｎｄｅｒ　ｓｕｒｆａｃｅ，ａｎｄ　ｔｈｅ　ｖｅｌｏｃｉｔｙ　ｄｉｓｔｒｉｂｕｔｉｏｎ　ｏｆ　ｔｈｅ　ｆｌｏｗ　ｆｉｅｌｄ　ａｇｒｅｅ　ｗｅｌｌ　ｗｉｔｈ　ｔｈｅ　

ｎｕｍｅｒｉｃａｌ　ａｎｄ　ｅｘｐｅｒｉｍｅｎｔａｌ　ｒｅｓｕｌｔｓ　ｉｎ　ｌｉｔｅｒａｔｕｒｅｓ．　

Ｋｅｙｗｏｒｄｓ：ｃｈａｒａｃｔｅｒｉｓｔｉｃ　ｂａｓｅｄ　ｓｐｌｉｔ；ｆｉｎｉｔｅ　ｅｌｅｍｅｎｔ　ｍｅｔｈｏｄ；ｉｎｃｏｍｐｒｅｓｓｉｂｌｅ　ｆｌｏｗ；ｍｏｖｉｎｇ　

ｂｏｕｎｄａｒｙ　

基本特征分裂（Ｃｈａｒａｃｔｅｒｉｓｔｉｃ－Ｂａｓｅｄ　Ｓｐｌｉｔ，　

ｌｏｒ—Ｇａｌｅｒｋｉｎ有限元等方法相比，其求解过程更加简　

ＣＢＳ）有限元方法［　是一种重要的针对流动问题的　

有限元解法，与Ｐｅｔｒｏｖ—Ｇａｌｅｒｋｉｎ迎风有限元、Ｔａｙ—　

单．另外，对于不可压缩黏性流动问题，ＣＢＳ算法在　

求解过程中允许速度和压力采用同阶的基函数，从　

收稿日期：２０１０—０５—１０．　作者简介：孙旭（１９８２－－），男，博士生；张家忠（联系人），男，教授，博士生导师．　基金项目：新世纪　

优秀人才支持计划资助项目（ＮＣＥＴ一０７—０６８５）．　

网络出版时间：２０１０—１０—１９　网络出版地址：ｈｔｔｐ：／／ｗｗｗ．ｃｎｋｉ．ｎｅｔ／ｋｃｍｓ／ｄｅｔａｉｌ／６１—１０６９一ｔ．２０１０１０１９．２２１１．０１９．ｈｔｍｌ　

西安交通大学学报　第４５卷　

而避开了ＢａｂｕＳｋａ—Ｂｒｅｚｚｉ（ＢＢ）条件．因此，该方法　

生扭曲甚至重叠；如果规定计算网格固结于空间坐　

自Ｚｉｅｎｋｉｅｗｉｃｚ等人于２Ｏ世纪９Ｏ年代首次提出以　

标上，不随边界移动，这将给边界位置的捕捉带来　

来，已经广泛应用到了各种流动问题之中．　

困难．　

近年来，越来越多的学者开始采用ＣＢＳ算法来　采用动网格技术，可以令边界上的节点随边界　

求解具有运动边界的流动问题，如ＮｉｔｈｉａｒａｓｕＥ　］将　

运动，同时允许求解区域内部的节点发生独立的运　

ＣＢＳ算法与ＡＬＥ方法　相结合，对溃坝流动和水　

动，这样既可以准确描述边界的位置，也可以保证网　

面孤立波的传播问题进行了求解，Ｎｏｂａｒｉ和Ｎａｄｅ—　格的质量．本节将传统ＣＢＳ算法与动网格技术相　

ｒａｎ［４　采用Ｎｉｔｈｉａｒａｓｕ的方法求解了自由来流中振　

结合，提出了一种网格运动时方程（２）的求解方法．　

动圆柱的绕流问题，段庆林ｌ＿６ｊ将ＣＢＳ算法与ＡＬＥ　为了方便起见，本节公式在推导过程中省略了无量　

方法和无网格解法相结合，成功地求解了成型充填　

纲标示符号．　

过程中非等温非牛顿黏性流动问题．这些解法均是　２．１网格运动时的ＣＢＳ解法　

通过将ＣＢＳ与ＡＬＥ相结合实现的，求解过程比较　

根据ＣＢＳ算法求解思路＿】］，通过坐标变换可消　

复杂．另外，由于ＡＬＥ坐标下ＮＳ方程的对流项含　除控制方程（２）中动量方程的对流项，得到新坐标系　

有网格速度项，网格速度将会给计算的稳定性带来　

Ｘ　下的动量方程．就流动问题而言，随流体质点运　

新的影响．　

动的Ｌａｇｒａｎｇｅ坐标系可以作为ｘ　坐标系，其上的　

本文根据ＣＢＳ算法的特点，针对具有运动边界　

动量方程可以表示为　

的不可压缩黏性流动问题，提出了一种简单、有效的　

Ｄｕ

ｇ

——

一一　

＋　：丝＋　ｆ３）　

求解方法，并采用所提出方法对静水中和自由来流　

Ｄｔ　ａｚｌ　Ｒｅ　３ｘｉａ　１　。　

中的振动圆柱绕流问题进行了求解．　

式中第一项表示流体质点的速度变化率．　

１控制方程　

下面对式（３）中的时间导数项进行离散．图１中　

给出了时间间隔　内的流体质点和网格点的运动　

非守恒形式的黏性不可压缩ＮＳ方程为　

情况．取ｔ　时刻网格节点Ａ　处的流体质点为研究　

３ｕ；ｑ

ｌ（１）　

对象，并假设该流体质点在ｔ　时刻位于Ｂ点，Ａ　的　

空间坐标为　，Ａ和Ｂ两点间的距离为　，通过式　

（３）可得　

３ｔ　Ｕｊ　Ｏ

—

一

：ｒｊ

　十ｕ古差＋ｕ　十　Ｊ＋　ｊ　

式中：“　为流体速度；ｐ为密度；ｐ为压强；ｕ为运动　

ｆ一　Ｊ　

——　——

　！！　一　

黏性系数；　为体积力．设参考速度为Ｕ，特征长度　

Ａｔ　

为Ｌ，参考密度为　，定义　一　Ｕｉ，；一芒，　一　，　

ｃ　（一差＋去　＋　）”　＋　

一　（不可压缩流动中　一１．ｏ），　，于　＝　

臼（一ａ］？＿　ｑ＿１　３２

ｕｉ

…

ｋｆｉ）州ｉ　（４）　

而ＪｉＬ，

，

则式（１）的无量纲形式为　

式中：Ｏ≤　１，　一０时为显式格式，　＝＝＝１时为隐式　

誊　一

ａ　ｕ

　）　

格式．式（４）中Ａｔ内流体质点移动的距离　可以近　

似为　

￣＋ｕ，　一一　３

１：△　ｆ　—Ａｍ”Ｉ＋ｏ（Ａｔ　）　（５）　

ｆ

ｘｉ　Ｒｅ　　３，：．ｃｊ３

～

ｘ

　Ｉ（

ｊ１十　ｊ　

２

式中：Ｒｅ＝ＵＬ／ｕ．　

２数值求解方法　

点　

对于具有运动边界的流动问题，流场中部分边　

质点　

界的位置将随时间不断变化．采用数值方法进行求　

解时，如果在计算过程中规定动边界上的网格节点　

・　

随边界同步运动，将造成动边界附近的网格单元发　

图１　Ａｔ内流体质点和网格点的运动情况　

ｈｔｔｐ：／／ｗｗｗ．ｊｄｘｂ．ｃｎ　

第１期　孙旭，等：具有运动边界不可压缩黏性流动的ＣＢＳ有限元解法　

通过Ｔａｙｌｏｒ展开来表示　一６　处的变量，并将　

式（５）代人式（４）整理可得　

，　

Ｅｎｄｄｏ　

对于式（９）～式（１１）的有限元离散，详细过程见　

Ｊ　

“　Ｉ～ｕ７　Ｉ＝＝＝　

～

ｌ　

文献Ｅｌｉ．　

２．２　网格运动　

为了满足运动边界条件，前文给出的求解流程　

需要不断地进行网格调整．目前，针对三角形非结构　

化网格已有多种动网格生成方法，本文采用了文献　

△　（“　一　１　３２Ｕ　ｉ　ｎ＠　＋　）ｆ　＋　

／Ｎｔ２　

２　（　

一　

Ｏ２，Ｃ　ｉ＋　）　

＋（卜　）　

～　ａｚ　

（６）　

（７）

…　

３　ｆ＇＋ｅ

ａ　

．　

７　

ＵＩ－　

［７３中的边弹簧近似法．首先，将ｔ　时刻的计算网格　

等效为由相互连接的弹簧组成的系统．定义相邻节　

式中：ｕ７＋。、Ｐ　为网格节点Ａ　上的变量；ｕ７、Ｐ　和　

为Ａ处的变量．式（６）即为空间坐标　处动量方　

点ｆ和　间的弹簧平衡长度等于两节点的距离，定　

义弹簧刚度　

ｋｕ—ａｌ（（Ｘｌ—ｚＪ）　＋（　一　）　）“２　（１２）　

程的离散形式．需要说明的是，在式（６）的推导过程　

中，除了压强项之外，其他变量均采用显式格式．　

网格运动时，式（６）中的网格节点Ａ　在ｔ”时刻　

并不在Ａ处，Ａ处任意流动参数　可以由网格点Ｃ　

式中：以　、ａ　为控制参数，本文中ａ　一１．０，Ｃ／２一　

一

１．０．由于所有弹簧的实际长度均等于其平衡长　

度，所以弹簧系统处于平衡状态．在ｔ卅　时刻，运动　

边界上的网格节点发生移动，部分弹簧上的受力不　

求得，而　的Ｔａｙｌｏｒ展开式为　

　｝

　。

一　

　．

～　

ａ声”ｌ　

再为０，弹簧系统的受力平衡也将遭到破坏．此时，　

任意节点ｊ上的弹性合力为　

Ｆｓ一　ｋｕ（　～６』）　（１３）　

＋。（１　Ｉ　）（８）　

式中：　。表示在　内网格节点Ａ　移动的距离．　

对式（６）求解可以通过以下３步进行．　

（１）引人中间变量　

一　一

Ｊ＝１　

式中：Ｎ，表示与节点　相邻的节点个数；　和　，表　

示节点，和　的位移．在弹性力的作用下，网格节点　

△　（　一　１　

“　

＋　）＋　

将发生移动．根据受力平衡准则，当弹簧系统重新达　

到平衡时，任意节点上　

Ｆｓ一　ｋｕ（　，一　）一０　（１４）　

Ａｔ２

杀（Ｕｊ　＋　）”　（９）　

（２）求解ｎ＋１时刻的压强Ｐ　，即　

Ｊ：１　

蠡（　）：　蠡（“　一Ａｔ　一臼　）　

（１０）　

由式（１４）和边界上的位移便可求得每个节点上　

的位移．本文将采用预条件共轭梯度法［８　对式（１４）　

进行求解．根据每一个节点上的位移，ｔ　时刻网格　

节点所处的位置可以表示为　

一ｘ７＋　ｒ　（１５）　

（３）求解　＋１时刻的速度“　州，即　

一　一　

ｏｚ　己０　＼（差１０　／”　

＋　

Ｘｋ

（１１）　

３计算实例及分析　

为了验证算法的可行性，本文对静水中和均匀　

因此，根据上文的描述，对于具有运动边界的不　

可压缩黏性流动问题，可以按照如下流程进行求解：　

Ｄｏ　一１，总的计算时间步　

（１）根据ｔ　时刻的网格和边界运动方式，采用　

来流中的振动圆柱绕流流场进行了计算，这２个问　

题的求解过程均未考虑体积力的影响．　

３．１静水中的圆柱振动　

动网格技术得到ｔ　时刻网格和每一个节点的位　

移　２；　

图２给出了静水中简谐振动圆柱绕流的计算区　

域和无量纲边界条件．计算区域外边界的直径为　

１００Ｄ，其上速度为０，并设任意一点的压强为０．圆　

柱沿ｚ方向作简谐振动，其无量纲运动方程为　

３２（　）一一　ｓｉｎ（２￣ｆｔ）　（１６）　

（２）通过式（８）计算ｔ　时刻网格节点所在位置　

在ｔ　时刻的流动参数ｕ７、Ｐ”和－厂　；　

（３）通过式（９）求解　；　

（４）通过式（１Ｏ）求解　；　

（５）通过式（１１）求解ｕ７¨；　

式中：．２７＝ｘ／Ｄ表示圆柱在水平方向上的位移；　：　

ｈｔｔｐ：／／ＷＷＷ．ｊｄｘｂ．ＣＦＩ　

西安交通大学学报　第４５卷　

Ｄ为振动幅值；ｆ—Ｄｆ／　为简谐振动频率，　

为圆柱振动的最大速度，圆柱表面采用无滑移　

边界条件．定义Ｒｅｙｎｏｌｄ数和Ｋｅｕｌｅｇａｎ—Ｃａｒｐｅｎｔｅｒ　

数分别为　

Ｒｅ—

ｒ　Ｔ　

Ｕ＿堕　

ｎ　

；ＫＣ—Ｕ

ｒ　Ｔ　

—　ｍａｘ　（１７）　

ｕ　Ｕ　

式中：Ｄ为圆柱直径．对Ｒｅ＝１００、ＫＣ＝：＝５的情况进　

行了计算，通过式（１６）和式（１７）可得　＝０．７９５　８、　

厂＝０．２．　

图３给出了静水中简谐振动圆柱绕流的非结构　

计算网格，它是沿周向和径向对流场进行剖分后得　

到的．采用３种不同密度类型的网格对该圆柱绕流　

进行了计算，表ｌ给出了网格相关信息以及计算所　

得阻力系数Ｃｒ】．从表１可以看出，采用Ｍｅｓｈ一３型　

网格计算所得ＣＤ与文献［９—１０］的结果之间的偏差　

均小于１　．图４给出了采用Ｍｅｓｈ一３型网格．计算　

所得振动相位为１８０。、２１０。和３３０。时流场中部分位　

置的流速分布，图中实际流速Ｕ　、　由计算所得无量　

纲流速“、口换算得到．从图４可以看出，本文计算结　

果与文献［１Ｏ］的实验结果非常吻合．　

图４　Ｒｅ为１００、ＫＣ为５、网格数为１４０×１２０，计算、　

实验所得振动相位为１８０。、２１０。和３３０。时流场　

中部分位置的流速分布　

图２静水中振动圆柱绕流　图３静水中振动圆柱绕流　

的计算模型　的计算网络　

袁ｌ　Ｒｅ为１００、ＫＣ为５时的圆柱阻力系数ＣＤ　

图５均匀来流中振动圆柱绕流计算模型　

穿透边界条件；出口位于圆柱下游３０Ｄ处，其上压　

力为０　Ｐａ．圆柱沿垂直于来流方向作简谐振动，其　

无量纲运动方程为　

，～　～

～

～　

（￡）一　ｃｏｓ（２７【－厂ｔ）　（１８）　

式中：　一　为圆柱在垂直于来流方向上的位移；　

３．２均匀来流中的圆柱振动　

＝

Ａ／Ｄ为振幅；　＝ＤＦ／ｕ。。为无量纲振动频率，根　

对Ｒｅ＝１８５时均匀来流中振动圆柱绕流进行　

据圆柱固定时旋涡自然脱落频率　来选取．圆柱表　

了计算，图５给出了计算区域和无量纲边界条件．流　

面为无滑移边界条件．　

场进口位于圆柱上游１０Ｄ处，流动参数为自由来流　

为了获得Ｒｅ＝１８５时圆柱绕流的旋涡自然脱　

的参数；上下边界距离圆柱１０Ｄ处，边界条件为无　

落频率　，首先对圆柱固定时的流场进行计算，计　

ｈｔｔｐ：／／ＷＷＷ．ｊｄｘｂ．ｃｎ　

第１期　孙旭，等：具有运动边界不可压缩黏性流动的ＣＢＳ有限元解法　

算网格通过先将求解区域分成９部分然后再剖分得　

到，如图６所示．表２给出了不同网格下计算所得圆　

柱平均阻力系数Ｃ。、均方根升力系数ＣＬ　，以及根　

据升力系数曲线得出的旋涡脱落的斯特罗哈数Ｓ　．　

通过对比可以看出，本文计算所得升力、阻力系数比　

文献［９，ｌ１］中的结果略大，斯特罗哈数则非常吻合．　

通过Ｓ　可以求得Ｒｅ一１８５时无量纲旋涡自然脱落　

频率２ｏ一０．１９６．　

表２　Ｒｅ为１８５时均匀来流中固定圆柱绕流计算结果　

图７　圆柱运动到上方最大位置时表面压强系数分布　

ｏ／（。）　

图８　圆柱运动到上方最大位置时表面涡量的分布　

图６均匀来流中振动圆柱绕流计算网络　

采用Ｍｅｓｈ一２型网格对　一０．２、ｆ／Ｊ；一０．８，　

１．０，１．２时振动圆柱绕流流场进行了计算．图７和　

图８给出了３种振动频率下圆柱运动到上方最大位　

置时表面压强系数Ｃ　及无量纲涡量　的分布情　

况．图中０表示圆柱表面节点相对于圆柱中心的转　

角，其值在迎风面驻点处为０，并按顺时针方向增　

加．如图７所示，本文计算所得Ｃ　与文献Ｅ９］的结果　

非常吻合．图８中，在涡量绝对值较小的地方，本文　

计算结果与文献［９１ｔ￣常吻合，在涡量较大的地方二　

者略有偏差．在求解过程中，流速采用线性插值基函　

数，涡量（流速的一阶导数）在单元相邻处不连续，这　

是造成二者涡量偏差的主要原因．　

图９给出了不同振动频率下圆柱升力系数Ｃ　、　

阻力系数Ｃ。随时间的变化情况及圆柱运动到上方　

最大位置处的瞬时流场．从图９中可以看出：当圆柱　

的振动频率较低（７／－厂Ｏ一０．８，１．０）时，Ｃ。　、ＣＤ均随　

时间呈简谐形式周期变化，／的增加将改变Ｃ　、ＣＤ　

变化的幅值，但并不改变流场中流线的拓扑结构；当　

尹较高时，ｃ　、ｃ。呈现出完全不同的变化形式，流场　

中流线的拓扑结构也将发生变化．关于厂对振动圆　

柱绕流流场的影响，可参见文献［９，１１］．　

（ｃ）ｆ～／ｆ～ｏ一１．２　

图９不同振动频率下圆柱升力、阻力系数随时间的变化及　

圆柱运动到上方最大位置时的瞬时流场　

ｈｔｔｐ　ｆｆ　ＷＷＷ．　ｄｘｂ．ｃｎ　

１０４　西安交通大学学报　第４５卷　

４结论　

学与实践，２００２（２）：７－ｌ１．　

ＹＵＥ　Ｂａｏｚｅｎｇ．Ｉ　Ｉ　Ｘｉａｏｔｉａｎ．Ｓｔｕｄｙ　ｏｆ　ｔｈｅ　ＡＩ　Ｅ　ｆｉｎｉｔｅ　

本文给出了一种动网格上的ＣＢＳ计算流程，用　

于求解具有运动边界的不可压缩黏性流动问题．作　

为验证，对Ｒｅ一１００、ＫＣ＝５时静水中振动圆柱绕　

流和Ｒｅ＝１８５、

厂／／ｏ一０．８，１．０，１．２时均匀来流中　

ｅｌｅｍｅｎｔ　ｍｅｔｈｏｄ　ａｎｄ　ｉｔｓ　ａｐｐｌｉｃａｔｉｏｎｓ［Ｊ］．Ｍｅｃｈａｎｉｃｓ　ｉｎ　

Ｅｎｇｉｎｅｅｒｉｎｇ，２００２（２）：７－１１．　

［４］ＮＯＢＡＲＩ　Ｍ　Ｒ　Ｈ，ＮＡＤＥＲＡＮ　Ｈ．Ａ　ｎｕｍｅｒｉｃａｌ　ｓｔｕｄｙ　

ｏｆ　ｆｌｏｗ　ｐａｓｔ　ａ　ｃｙｌｉｎｄｅｒ　ｗｉｔｈ　ｃｒｏｓｓ　ｆｌｏｗ　ａｎｄ　ｉｎｌｉｎｅ　ｏｓｃｉｌ—　

的振动圆柱绕流问题进行了求解，所得结果与已有　

实验和数值结果非常吻合，说明了算法的可行性．与　

ｌａｔｉｏｎ［Ｊ］．Ｃｏｍｐｕｔ　Ｆｌｕｉｄｓ，２００６，３５（４）：３９３—４１５．　

［５］　ＮＯＢＡＲＩ　Ｍ　Ｒ　Ｈ，ＧＨＡＺＡＮＦＡＲＩＡＮ　Ｊ．Ａ　ｎｕｍｅｒｉｃａｌ　

ｉｎｖｅｓｔｉｇａｔｉｏｎ　Ｏｆ　ｆｌｕｉｄ　ｆｌｏｗ　ｏｖｅｒ　ａ　ｒｏｔａｔｉｎｇ　ｃｙｌｉｎｄｅｒ　ｗｉｔｈ　

Ｎｉｔｈｉａｒａｓｕ所提方法相比，本文计算仍在欧拉坐标　

系下完成，程序主体与传统的ＣＢＳ解法基本相同，　

求解过程更加简单．　

ｃｒｏｓｓ　ｆｌｏｗ　ｏｓｃｉｌｌａｔｉｏｎ［Ｊ］．Ｃｏｍｐｕｔ　Ｆｌｕｉｄｓ，２００９，３８　

（１０）：２０２６－２０３６．　

［６］段庆林．成型充填过程中非等温非牛顿粘性流动的　

ＡＬＥ有限元与无网格自适应耦合模拟［Ｄ］．大连：大连　

理工大学能源与动力学院，２００７．　

需要说明的是，该算法在计算时的每一步都要　

先根据边界运动求出ｔ”¨时刻的网格，因而只适用　

于动边界运动方式已知的情况．对于动边界运动未　

知的问题（如具有自由水面的流动问题），需要对已　

有的求解流程进行改进，这也是下一步需要研究的　

内容．　

参考文献：　

』ｌ　ｊ　ＺＩＥＮＫＩＥＷＩＣＺ　０　Ｃ，ＴＡＹＬＯＲ　Ｒ　Ｉ　Ｔｈｅ　ｆｉｎｉｔｅ　ｅｌｅ—　

［７］ＢＬＯＭ　Ｆ　Ｊ．Ｃｏｎｓｉｄｅｒａｔｉｏｎｓ　ｏｎ　ｔｈｅ　ｓｐｒｉｎｇ　ａｎａｌｏｇｙ［Ｊ］．　

Ｉｎｔ　Ｊ　Ｎｕｍｅｒ　Ｍｅｔｈ　Ｆｌｕｉｄｓ，２０００，３２（６）：６４７－６４８．　

Ｅ８］王勖成．有限单元法［Ｍ］．北京：清华大学出版社，　

２００３：２４５—２５０．　

［９］ＧＵＩＩ　ＭＩＮＥＡＵ　Ｅ，ＱＵＥＵＴＥＹ　Ｐ．Ａ　ｎｕｍｅｒｉｃａｌ　ｓｉｍｕ—　

ｌａｔｉｏｎ　ｏｆ　ｖｏｒｔｅｘ　ｓｈｅｄｄｉｎｇ　ｆｒｏｍ　ａｎ　ｏｓｃｉｌｌａｔｉｎｇ　ｃｉｒｃｕｌａｒ　

ｃｙｌｉｎｄｅｒ［Ｊ］．Ｊ　Ｆｌｕｉｄｓ　Ｓｔｒｕｃｔ，２００２，１６（６）：７７３—７９４．　

［ＩＯ］ＤＵＴＳＣＨ　Ｈ，ＤＵＲＳＴ　Ｆ，ＢＥＣＫＥＲ　Ｓ，ｅｔ　ａ１．Ｌｏｗ－　

Ｒｅｙｎｏｌｄｓ－ｎｕｍｂｅｒ　ｆｌｏｗ　ａｒｏｕｎｄ　ａｎ　ｏｓｃｉｌｌａｔｉｎｇ　ｃｉｒｃｕｌａｒ　

ｍｅｒｉｔ　ｍｅｔｈｏｄ：ｖｏｌｕｍｅ　３　ｆｌｕｉｄ　ｄｙｎａｍｉｃｓ（Ｍ　Ｉ．５ｔｈ　

ｅｄ．Ｏｘｆｏｒｄ，ＵＫ：Ｅｌｓｅｖｉｅｒ　Ｂｕｔｔｅｒｗｏｒｔｈ－Ｈｅｉｎｅｍａｎｎ，　

２００５．　

ｃｙｌｉｎｄｅｒ　ａｔ　ｌｏｗ　Ｋｅｕｌｅｇａｎ－Ｃａｒｐｅｎｔｅｒ　ｎｕｍｂｅｒｓ［Ｊ］．Ｊ　

Ｆｌｕｉｄ　Ｍｅｃｈ，１９９８，３６０：２４９—２７１．　

Ｅ２］　ＮＩＴＨＩＡＲＡＳＵ　Ｐ．Ａｎ　ａｒｂｉｔｒａｒｙ　Ｌａｇｒａｎｇｉａｎ　Ｅｕｌｅｒｉａｎ　

（ＡＩ　Ｅ）ｆｏｒｍｕｌａｔｉｏｎ　ｆｏｒ　ｆｒｅｅ　ＳＵｒｆａｃｅ　ｆｌｏｗｓ　ｕｓｉｎｇ　ｔｈｅ　

［１１］Ｉ．Ｕ　Ｘ　Ｙ，ＤＡＬＴＯＮ　Ｃ．Ｃａｌｃｕｌａｔｉｏｎ　ｏｆ　ｔｈｅ　ｔｉｍｉｎｇ　ｏｆ　

ｖｏｒｔｅｘ　ｆｏｒｍａｔｉｏｎ　ｆｒｏｍ　ａｎ　ｏｓｃｉｌｌａｔｉｎｇ　ｃｙｌｉｎｄｅｒ［Ｊ］．Ｊ　

Ｆ１ｕｉｄｓ　Ｓｔｒｕｃｔ，１９９６，１０（５）：５２７—５４１．　

ｃｈａｒａｃｔｅｒｉｓｔｉｃ—ｂａｓｅｄ　ｓｐｌｉｔ（（　Ｓ）ｓｃｈｅｍｅ＿Ｊ］．Ｉｎｔ　Ｊ　Ｎｕ—　

ｍｅｒ　Ｍｅｔｈ　Ｆｌｕｉｄｓ，２００５，４８（１２）：１４１５－１４２８．　

［３］岳宝增，李笑天．ＡＬＥ有限元方法研究及应用［Ｊ］．力　

（编辑苗凌）　

（上接第５７页）　

［１Ｏ］ＺＨＡＮＧ　Ｚ，ＸＩＥ　Ｙ，ＺＨＡＮＧ　Ｘ，ｅｔ　ａ１．Ａｎａｌｙｓｉｓ　ｏｆ　

ｐｉｓｔｏｎ　ｓｅｃｏｎｄａｒｙ　ｍｏｔｉｏｎ　ｃｏｎｓｉｄｅｒｉｎｇ　ｔｈｅ　ｖａｒｉａｔｉｏｎ　ｉｎ　

２Ｏ１Ｏ，４４（４）：０５５５—０５５９．　

［１２］王跃武，胡芄，陈则韶．ＨＦＣ一２２７ｅａ的ＭＢＷＲ状态方　

程［Ｊ］．西安交通大学学报，２００７，４１（１）：３７—４０．　

ＷＡＮＧ　Ｙｕｅｗｕ，ＨＵ　Ｗａｎ，ＣＨＥＮ　Ｚｅｓｈａｏ．ＭＢＷＲ　

ｔｈｅ　ｓｙｓｔｅｍ　ｉｎｅｒｔｉａ［Ｊ］．Ｐｒｏｃ　Ｉｎｓｔｎｍｅｃｈ　Ｅｎｇｒｓ：Ｐａｒｔ　Ｄ　

Ｊｏｕｒｎａｌ　ｏｆ　Ａｕｔｏｍｏｂｉｌｅ　Ｅｎｇｉｎｅｅｒｉｎｇ，２００９，２２３（４）：５４９　

５６３．　

ｅｑｕａｔｉｏｎ　ｏｆ　ｓｔａｔｅ　ｏｆ　ＨＦＣ　２２７ｅａ［Ｊ］．Ｊｏｕｒｎａｌ　ｏｆ　Ｘｉ　ａｎ　

Ｊｉａｏｔｏｎｇ　Ｕｎｉｖｅｒｓｉｔｙ，２００７，４１（１）：３７—４Ｏ．　

［１１］张执南，谢友柏．考虑连杆一曲轴变惯性影响的活塞二　

阶运动分析［Ｊ］．上海交通大学学报，２０１０，４４（４）：　

５５５—５５９．　

Ｅｌ３］盂祥慧．复杂产品时变性能设计和设计中的知识服务　

研究ＪＲ］．上海：上海交通大学机械与动力工程学　

院，２００８．　

ＺＨＡＮＧ　Ｚｈｉｎａｎ，ＸＩＥ　Ｙｏｕｂａｉ．Ｅｆｆｅｃｔｓ　ｏｆ　ｃｏｎｎｅｃｔｉｎｇ　

ｒｏｄ－ｃｒａｎｋｓｈａｆｔ　ｖａｒｉａｂｌｅ　ｉｎｅｒｔｉａ　ｏｎ　ｐｉｓｔｏｎ　ｓｅｃｏｎｄａｒｙ　ｍｏ—　

ｔｉｏｎ［Ｊ］．Ｊｏｕｒｎａｌ　ｏｆ　Ｓｈａｎｇｈａｉ　Ｊｉａｏｔｏｎｇ　Ｕｎｉｖｅｒｓｉｔｙ，　

（编辑管咏梅）　

ｈｔｔｐ：∥ｗｗｗ．ｊｄｘｂ．ｃｎ