膜结构是一种柔性空间结构.膜结构的动力特性研究表明, 其自振频率分布密集且相互耦合.它的轻柔特点使其在风荷载作用下的反应成为人们关注的焦点.膜结构承担面外荷载的刚度主要来源于预应力, 当结构跨度较小时, 对结构施加较大的预应力可以大大提高其刚度, 此时风荷载产生的反应不会太大, 可以采用一般的动力反应分析方法或随机振动方法进行计算;但是当结构的跨度较大时, 膜结构较轻柔, 风荷载作用下会有较大的位移、速度和加速度反应, 使得结构发生剧烈的形状改变, 而此变化又会反过来影响风场的分布, 从而引起膜结构和风场之间的耦合作用[1].

流体与膜结构的相互作用可能会导致整个系统产生不稳定的发散振动, 即气弹失稳.研究者目前仅仅提出了一些简单的气弹模型, 很难解决实际问题.此外, 可以通过风洞试验来研究膜结构的流固耦合问题, 但是成本较大.随着计算机技术的发展, 用计算流体力学和计算固体力学来解决流固耦合问题有了可能, 这就是数值风洞的概念[2].流固耦合数值模拟主要有强耦合和弱耦合2种方式, 强耦合目前仅仅有少数简单算例, 弱耦合则可以发挥已有有限元和CFD软件的优势, 仅仅需要进行耦合域层面问题的处理, 比如, 进行结构域和流体域的数据传输程序设计.另外, 土木工程结构的特征湍流问题也在研究中比较热点的大涡模拟技术已经开始百富策略白菜网于膜结构流固耦合计算中.随着计算机硬件技术的提升和湍流模型研究的深入, 直接求解基于膜结构特征湍流模型的Navier-Stokes (N-S) 方程的方法将会被愈来愈广泛百富策略白菜网.

国内外学者对膜结构的流固耦合问题做了大量研究.Uematsu等人[3]对四周封闭、底部通风的菱形双曲抛物张拉膜结构进行了大量刚性模型和弹性模型风洞试验, 比较了膜材厚度、质量、预张力不同的9组模型在不同风速、风向角、湍流强度下的位移.试验表明, 膜面平衡位置的变化引起了刚度变化, 从而导致结构振动频率随着风速提高有所偏移;各模型测点位移近似与风速平方成正比, 不受预张力影响.Minami[4]采用薄型机翼理论, 分析了平坦矩形膜结构在势流下的附加质量, 该附加质量从膜结构振动以及诱导的空气振动之间的能量关系导出.其数值计算结果表明, 附加质量和结构振动频率、振幅无关, 仅和空气密度、结构跨度有关.沈世钊、武岳[5]基于简化气弹模型方法和计算流体力学 (CFD) 数值模拟方法, 开展了相应的理论推导、风洞试验和数值模拟工作, 获得了关于风与膜结构相互作用机理的一些初步认识, 提出了一种具有操作性的简化数值分析方法.而且, 武岳基于大涡模拟编制了二维流固耦合计算程序, 在国内首次利用数值仿真技术再现风与二维柔性结构的动力耦合全过程.

尽管可以百富策略白菜网数值方法模拟风-膜结构的流固耦合过程, 但目前的CFD软件不具备成熟的流固耦合计算功能.而且, CFD软件大都基于有限体积法 (FVM) , 需要对整个流场进行单元划分, 特别在解决对计算精度要求较高的问题中, 划分出的单元量巨大.另外, 膜结构变形的几何非线性程度大, 会在计算过程中引起初始划分的单元奇异, 必须重新划分单元才可以继续计算.所以单元划分耗费了人们大量精力, 导致总体计算时间非常漫长.

本文的数值模拟基于边界单元法.边界单元法是在经典积分方程的基础上, 吸收了有限元法 (FEM) 的离散技术而发展起来的计算方法.边界元法仅仅在流场边界上离散计算区域, 大大减少了未知量, 节省了流体域求解的机时.膜结构作为一种曲面结构, 其与流体的耦合恰恰发生在流场的边界, 边界解可以直接用于膜结构变形计算, 这恰好利用了边界元的优势.从计算格式形成的过程来看, 关键问题有2个[6]:一个是问题的边界化, 即将给定区域上的定解问题化为可以只考虑边界的问题, 这一步的关键是格林公式, 也是边界元法的基石.边界化的结果使得问题降维, 如果是各维尺度相近的大型问题, 代数方程组的未知数按指数分布减少, 这无疑将大大减少计算存储量与计算时间.有限单元法则要将全部区域及边界离散, 并要求将全部节点纳入方程进行计算, 这是因为这些节点值包含在同一个封闭方程组中.边界元法的封闭方程组中只有边界节点上的未知量, 求解完方程组再根据需要有目的地求解内部节点值, 因此减少了计算的盲目性;另一个问题是边界的离散化, 离散技术与有限元法基本相同.但是边界元法仅仅对边界离散, 计算误差仅仅来自边界, 所以可以将减小误差的工作专注于边界上, 这将提高计算的准确度.而且, 区域内部点的物理量可以由解析公式计算, 这也提高了计算精度, 并保证了内部值的连续性、可微性.而现行的其他数值计算方法, 如有限差分法 (FDM) 、有限体积法 (FVM) , 都是基于区域网格的数值方法, 也都在整个域内考虑问题.

尽管边界元法有很多优势, 但目前在膜结构流固耦合问题的研究中, 边界元法的百富策略白菜网还很少Sygulski[7]利用边界元法并考虑膜结构流固耦合效应对膜结构进行稳态分析, 不过距百富策略白菜网边界元法对膜结构流固耦合全过程进行数值模拟有很大距离.国内研究者徐闻等人[8,9]将边界元法百富策略白菜网在膜结构与风耦合计算中, 假定流场为势流场, 给出了流固耦合模拟流程, 并对风-膜结构流固耦合效应进行了参数分析, 为边界元法在该领域的百富策略白菜网研究打下良好基础但在膜结构边缘、尖角处理上比较粗糙, 导致了这些区域模拟误差大;另外, 计算时间仍然很长, 没有充分发挥边界元法的优势.

针对前人研究的不足, 本文基于小扰动位势理论, 百富策略白菜网非协调边界元法 (BEM) 和有限元法 (FEM) 采取弱耦合策略对风-膜结构进行了流固耦合分析由于要计算的流场域边界非连续光滑, 有尖角和尖边存在, 本文引入了非协调边界元技术, 使得计算精度大大提高.另外, 对于系数矩阵不对称的大型边界元方程组, 本文引入了非常高效的预处理GMRES迭代算法, 使得边界元法的优势得到了充分发挥.

假定流场为无黏、不可压缩, 且为势流场.引入流体运动的连续方程[10,11]

由于流体不可压缩, 即ρ为常数, 故∂t∂ρ=0, ∇ (ρV) =ρ∇V, 代入 (1) 式得

对于势流, 可以定义速度势函数ϕ满足V=∇ϕ, 代入 (2) 式得到

此方程即是非常著名的数学物理方程——Laplace方程.

在基于小扰动理论的流场计算中, 流场边界的外法向速度都是已知的, 流场求解可以总结为下列问题:

控制方程 (3) 两端同时乘以权函数ϕ*并在流场整个区域Ω积分得

利用Green第二公式得

设ϕ*为下面方程确定的基本解

式中, δ (Q, P) 是Dirac-δ函数, Q和P分别表示场点和源点 (关于场点和源点的概念见参考文献[6]) .

将 (7) 式代入方程 (6) , 得到

已经导出的积分方程 (8) 仅仅对内部点有效.而对于边界点P, 可在P处扩充域的范围, 加一个半径为ε的部分球面 (图1) , 形成新的积分域Γ.分别在球面部分Γε和Γ-Γε, 部分进行积分, 并考虑取极限ε→0, 最后得到方程[6]

为了理解的方便, 边界积分方程 (9) 可以写为

式中, P代表源点, Q代表场点.

为了得到边界积分方程组, 首先将边界离散为NE个单元, 有

设单元节点数为M.对第e个单元, 以ϕα, qα表示此单元第α个节点的速度势、法向速度, 有

式中, Nα (α=1, ⋅⋅⋅, M) 为插值函数.

将 (12) 式代入边界积分方程 (11) 得到

(13) 式即源点为边界点P时的积分方程的离散形式假设有n个边界点, 势函数在计算边界域的n个边界点取值分别为ϕ1, ϕ2, ⋅⋅⋅, ϕn, 势函数的法向导数在n个边界点取值分别为q1, q 2, ⋅⋅⋅, qn, 可以得到n个形如 (13) 式的方程.令

代入 (14) 式整理后, 可以得到n阶线性方程组

式中, c=diag{c 1, c 2, ⋅⋅⋅, c n}, ϕ= (ϕ1, ϕ2, ⋅⋅⋅, ϕn) T.

利用法向速度边界条件, 可以通过方程 (16) 求解出所有边界点的速度势.

对于三维问题, 当场点Q和源点P位于同一个单元时, 这时积分方程式 (13) 中核函数ϕ*和q*分别含有因子r-1和r-2, 相应的积分分别为弱奇异和强奇异积分[12].关于奇异积分的数值求解, 较多学者做了相关研究[13~18].

对于弱奇异积分, 可以采取子分单元法 (subdivision technique) [19~22].以八节点平面单元为例, 当源点P处于单元的角点, 将原单元分割为2个子单元 (参见图2) .

如图3, 直接将图2中分割出来的子单元A映像为局部坐标空间 (ξ, η) 里面的八节点单元, 而此变换的雅可比行列式正比于r, 恰好消除了核函数ϕ*中因子r-1带来的积分奇异性, 使之变为正常积分.假设原单元被分割为Ns个子单元, 奇异积分将有如下的形式

式中, Js (ξ, η) 是从全局坐标空间映像到局部坐标空间 (ξ, η) 的雅可比变换行列式.

边界积分方程 (13) 的强奇异积分项 (在Cauchy主值意义下) 及自由项c构成了边界积分方程组 (16) 中矩阵H的对角项.尽管利用某些坐标变换[23]可以对强奇异积分项进行精确计算, 但过程较为复杂, 故本文利用H矩阵各行元素之和等于0的特性来计算对角项.

当场点靠近源点时, 在场点所在单元内进行的高斯积分误差较大[24].高斯积分是在局部坐标系 (ξ, η) 中进行的, 表达式为

式中, (ξi, ηj) 是高斯积分点局部坐标, wiξwjη是相应于积分点 (ξi, ηj) 的权重, 1m和2m分别对应坐标ξ和η方向的高斯积分阶次, 1E和2E分别对应坐标ξ和η方向的积分误差.

文献[19]给出的高斯积分阶次表达式为

式中, ek=Ek/I (k=1, 2) 为高斯积分相对误差, p是被积函数的奇异阶次, kL (k=1, 2) 为单元长度, r是源点到被积分单元的最小距离, p′的定义为

由 (19) 式可以得到

由 (20) 式可知, 高斯积分误差与Lk/r成正指数阶相关.当场点靠近源点时, 积分误差较大.为了减少数值积分误差, 可以将原单元进一步划分为一系列更小单元以减小误差.为了百富策略白菜网以上误差标准进行积分计算, 需要确定几何参数r和Lk.

由微分几何知识知, 微元向量表达式为

单元长度为

式中, Nα (ξ, η) 是形函数, M是单元节点数目, xjα是单元中的第α个节点的全局坐标.

源点到积分单元最小距离r不可以直接确定, 本文百富策略白菜网Newton-Raphson迭代[25]进行求解.此方法的实质是通过迭代使得计算点不断逼近源点, 从而找到被积单元边界上距离源点最近场点的局部坐标.

本文取允许误差限ei=0.00000001, 最大允许高斯积分阶次mmax=10, 利用 (19) 式得到积分阶次mi如果mi≤mmax直接在单元上进行高斯积分 ( (18) 式) ;如果im>mmax, 令im=mmax, 计算子单元长度L is (i=1, 2) ( (20) 式) 后进一步划分原单元[26,27].

如图4, 如果要在某个单元上进行积分, 依照上文高斯积分误差标准进行分析后, 该单元需要被进一步分割为Ns个子单元.在具体的积分过程中, 需要进行2次坐标系之间的雅可比变换.第1个变换是将坐标空间 (x, y, z) 中的该单元整体映像在局部坐标 (ξ, η) 中, J (ξ, η) 是相应的雅可比变换行列式;第2个变换是将坐标空间 (ξ, η) 中的子单元 (例如图4中的子单元S) 映像于局部坐标 (ξ′, η′) 中, Js (ξ′, η′) 是相应的雅可比变换行列式.

(23) 式中, (ξi′, η′j) 是高斯积分点坐标, wiξ′wjη′是相应于积分点 (ξi′, η′j) 的权重, m1和m2分别对应坐标ξ′和η′方向的高斯积分阶次.

这2个变换过程在编程实现时, 只需要编写一段变换函数代码即可以.原因是2种雅可比变换过程是相似的.

广义极小残余法 (generalized minimum residual) , 简称为GMRES, 是一种求解大型非对称线性方程组的Krylov子空间投影法, 1986年由Saad和Schultz提出[28].该方法源自基于Galerkin原理的Arnoldi算法.

对线性方程组

上式中, A∈Rn×n为一非奇异矩阵, b∈Rn为一指定向量.设x0∈Rn为任一向量, 令x=x0+z, 0r=bAx0, 则求解 (24) 式等价于求解

设所要求解的方程组为 (25) 式, 取

设zm∈Km, 利用Galerkin原理, 需让残余向量rm=r0-Azm与Lm中的所有向量正交, 从而求得 (25) 式的近似解zm.而对于这组特定的Lm和Km取法, 求解zm相当于在Krylov子空间Km中极小化残余向量的2-范数[29,30], 故此方法被命名为广义极小残余法.本文采用下面的循环型GMRES (m) 算法[30].

(a) 选择x0∈Rn, 计算r0=b-Ax0, V1=r0/||r0||, 记β=||r0||;

(b) 选择适当大小的m, 开始Arnoldi过程:

2) 对k=1, 2, ⋅⋅⋅, m, 计算

Hm被称为上海森伯 (Heissenberg) 矩阵, em={0, 0, ⋅⋅⋅, T0, 1}∈.mR

(e) 如果残余向量范数|rm|2=||r0-Azm||2<ε, 近似解为x*=xm, 停止计算;

(f) x0=xm, r0=xm, v1=rm/||rm||, 转到 (b) .

对线性方程组 (24) , 如果其系数矩阵条件数较大则方程组是病态的, 为了发挥GMRES算法高速收敛的特性, 需对方程组 (24) 进行预处理, 以改善其系数矩阵的条件数.本文采用左预处理, 设预处理矩阵为ML (本文采用 (29) 式) , 在方程组 (24) 的左右两端分别乘以预处理矩阵的逆矩阵ML-1, 得到

用GMRES法求解方程组 (28) 即可.

本文的预处理采用简单易操作的对角线法, 不需要对系数矩阵进行繁琐的数学处理.设n阶系数矩阵A= (aij) , 左预处理矩阵为n阶阵ML= (mij) , 有

实际问题中流场域边界非光滑, 有尖角和尖边存在.而尖角尖边处节点的法向速度q难以确定, 因为此处节点同时处于不同的面, 各个面q值未必一致即法向速度在此处突变.在边界元求解时, 要求合理给出这些点的q值, 否则对整个方程组求解结果的精度有较大影响.处理这个问题的手段主要有增加单元法、混合单元法和非协调单元法.本文采用非协调单元法, 它将所有节点均取在光滑边界处, 从而避免了法向速度的不连续现象.和普通单元节点取在单元端点不同[31], 非协调单元的节点取在单元的内部 (如图5和6) .

非协调单元插值函数与普通单元不同.对于图5中节点1和2所在单元, 节点1和2的局部坐标为ξ1=-1/2, ξ2=1, 对应的插值函数分别为N1 (ξ) , N2 (ξ) , 单元中任一点的法向速度为

利用插值函数性质

得到

而二维单元可用类似方法求得单元的插值函数.

是否采用非协调单元要由边界点处边界的光滑性来判断, 这一过程在计算程序中自动实现.对于三维模型, 可以在边界元程序中判断相邻边界面单元的法向向量的夹角.如果夹角接近0°, 那么可以近似认为此处表面光滑;否则, 应采取非协调单元 (图6) .而对于二维模型, 则通过判断相邻线单元的切线的夹角来判断光滑性, 在尖角处采取非协调单元 (图5) .增加非协调单元后, 边界点数量增加, 要利用新的节点及相应坐标进行边界元计算.

在耦合过程中, 有一个关键问题就是由三维空间膜面上的节点速度势得到其节点力.假定在膜面上势函数是连续函数, 且膜面光滑, 所以势函数可微假定边界元单元的节点数为M.坐标系定义如图7, 其中 (x, y, z) 为全局笛卡儿坐标系; (ξ, η) 为单元面内的局部坐标系, 坐标原点位于单元中心O (ξ=η=0) , Q (ξ, η) 为单元内任一点.另外, 向量rξ方向与过Q点的曲线L1 (ηL1=η) 的切向一致, 向量rη方向与过Q点的曲线L2 (ξL2=ξ) 的切向一致.还需要定义另外一个笛卡尔坐标系 (x′, y′, z′) , 其坐标原点位于单元内点Q (ξ, η) , x′轴方向与向量rξ方向一致, z′轴方向与向量rξ×rη方向一致.

局部坐标系 (ξ, η) 中

(x′, y′) 平面坐标系中, 膜面单元内任一点 (ξ, η) 的速度为

式中, J为坐标系 (x′, y′) 映像到局部坐标系 (ξ, η) 的雅可比变换矩阵.关于J的得出, 利用

可以得到

另外, 设矩阵α是坐标系 (x′, y′, z′) 和坐标系 (x, y, z) 间的坐标变换矩阵, 且ei (i=1, 2, 3) 和fj (j=1, 2, 3) 分别是坐标系 (x′, y′, z′) 和 (x, y, z) 的基向量, 那么

式中, f1=i, f2=j, f3=k.

将 (37) 式代入 (36) 式得

要计算出αij, 才可以确定雅可比变换矩阵J.由微分几何知识[32]可知 (如图7)

因为r=xi+yj+zk, 所以

式中

记n*=rξ×rη, n=n*/|n*|, (x′, y′, z′) 坐标系的基向量为

将Q点的局部坐标值 (ξ, η) 代入 (39) ~ (42) 式, 计算出e1, e2, e3, 再利用 (37) 式可得坐标变换矩阵α.将α代入 (34) ~ (36) 式, 可得单元内任一点Q处的流体速度

流体速度求出后, 考虑势流的Bernoulli方程

式中, V是流速, g代表重力加速度, z为流体质点处相对高度, p为压强, ρ为流体密度.

利用方程 (43) 可以求出膜面处流场压强, 进而得出膜面单元节点受到的风压力.

本文采取的是弱耦合策略.耦合过程中主要用到有限元软件ANSYS10.0和边界元计算程序以及两者之间的数据传递程序.本文设计的弱耦合流程如下.

(a) 对流场整体模型的边界面进行划分;

(b) 运行边界元程序, 得到节点速度势值.提取膜面处的节点速度势, 运行程序NodeForce.exe得到节点力 (程序编制基于 (33) ~ (43) 式) ;

(c) 建立膜面的有限元计算模型, 将步骤 (b) 得到节点力转化到有限元模型中相应节点上 (程序为BEM_NF_FEM.exe) ;

(d) 用ANSYS计算得到膜面新的节点坐标;

(e) 利用 (d) 得到节点坐标, 更新相应边界元模型中节点的坐标 (程序为ExchangeCoor.exe) ;

(f) 重复步骤 (b) ~ (e) , 直到前后2次迭代的结果差值满足设定要求为止.

电脑硬件条件为:AMD Athon (tm) 64 X2;DualCore Processor Tk-55:1.79 GHz, 896 MB内存.流体无黏、不可压缩且无旋.流场长200 m, 宽200 m, 高100 m, 刚性半球直径40 m, 球心位于流场下底面中心 (图9) .风从左向右吹过, 模型入口面流速V∞恒为1 m/s.计算模型初始划分992个平面四边形单元、994个节点 (如图10) , 其中半球面部分321个单元.

运行边界元程序, 程序自动判断模型尖角尖边位置, 并自动设置非协调单元并增加新节点, 最后的节点数为1138, 图11为增加节点后的单元划分情况, 模型尖角尖边已经被非连续化.计算时间为9 min.计算结果如图12和13, 由计算结果发现:最大速度势为100.91 m2 s-1, 最小速度势为-100.909 m2 s-1, 分别位于出流面和入流面;速度势差值为201.819 m2 s-1, 比无扰动情况下的理论速度势差值200 m2 s-1稍大 (因为模型对流场产生了微小扰动) .

为验证边界元程序计算的准确性, 本文将数值解和解析解进行对比.样点分布见图14.

而理想势流绕圆球流动时势函数解析表达式为

式中, V∞为来流速度, r0为圆球半径, r, θ分别为流场中某点的极径和极角 (如图15) .

对比结果见表1.表1列出了未使用和使用非协调元2种情况下的计算结果.结果表明, 使用非协调元后, 边界元程序得到的数值解与解析解间的绝对误差要小于未使用非协调元的情况, 而且计算误差非常接近于0, 可见边界元计算程序是正确可靠的.

流体为无黏、不可压缩无旋流.流场长160 m, 宽160 m, 高80 m (图16) ;膜结构模型四周封闭, 顶面为菱形马鞍膜面, 模型下底面中心坐标 (10, 10, 0) , 膜结构模型见图17, 图17中HP表示高点, LP表示低点.风从左向右吹过, 模型入口面流速V∞恒为20 m/s同时膜结构的2个低点连线与风向一致.膜结构跨度L=28 m, 低点高度h=3.667 m, 矢跨比f/L=1/12.膜面初始张力T=2.5 kN/m, 薄膜厚度t=1.0 mm, 单位面积质量m=1.25 kg/m2, 张拉刚度Et=2550 N/cm, Gt=800 N/cm, 泊松比v=0.3.

计算模型共划分4888个平面四边形单元, 4890个节点 (图18) , 其中膜面部分900个单元 (图19) .采用非协调元后, 节点数目增加, 最后的计算节点数为5362 (图20) .终止条件是前后2次迭代得到各节点变形差值绝对值的均值小于0.02 m且方差小于0.003 m2;迭代8次达到收敛.计算机配置如上例, 每次迭代中的流体计算时间仅85 min, 而本文作者原来百富策略白菜网高斯-若当消去法求解则需要10 h左右.由于计算所用电脑配置较差, 如果利用高性能电脑计算, 速度还将大大提高.可见, GMRES算法在求解大型线性方程组方面有非常大的优势, 在工程计算中有很高的百富策略白菜网价值.

在流体力学中, 一般将压强用无量纲的数压强系数Cp来表示, 表达式为

式中, p, p∞, p0分别是某点压强、来流压强 (静压) 、驻点压强 (总压) .

本算例中, 静压取1个标准大气压, 即∞p=101325 Pa, V∞=20 m/s, 图21是利用本文方法得到的膜面风压系数分布.取膜面最后的稳定形态建立流场模型, 在Fluent软件设定流体为层流, 计算得到风压系数分布如图22.

观察图21和22发现, 边界元程序计算结果 (图21) 所示的风压呈现对称分布, 符合位势流的特点.与Fluent软件计算结果 (图22) 比较发现:膜面大部分与Fluent模拟结果很相近;仅在结构的边缘部分有差异, 这是因为此边界元程序是基于位势流的, 没有考虑旋涡的影响.而与国内哈工大的刚性模型结果 (图23) 对比发现:在马鞍面中间凹处, 风洞试验得到的风压系数绝对值偏小0.1.主要是因为本文边界元计算基于势流, 而为了方便数值方法间的对比 (尽量与边界元的势流假定接近) , Fluent的计算采取了层流, 没有引入湍流模型计算;另外风洞试验是刚性模型试验, 没有反应出流固耦合效应 (哈工大的风洞试验研究表明气弹模型试验结果与刚性模型试验相比, 平均风吸力约增大10%左右, 亦即刚性模型的风压系数绝对值偏小) , 再加上试验设备、测量误差、试验条件的影响导致了数值解和试验解的偏差.实际上, 结构前缘会卷起旋涡, 后缘和尾流区也常常有脱落的旋涡, 所以本文成果对于外形不规则、易激起大面积旋涡区域的结构应慎用.对于此类结构, 下一步的研究工作是拟通过涡方法模拟或者修正旋涡卷起区、脱落区及尾流区, 进而得到可以用于指导工程的风压值.比如, 在航空工程中常用的离散涡方法, 就是在势流场的局部区域嵌入有限数目的离散涡, 用来代表局部有旋区域的连续分布的涡量, 通过计算离散旋涡的成长、卷起及脱落的复杂演变过程来实现对整个流场的数值模拟.

边界元法仅仅在流场边界上离散计算区域, 大大减少了未知量, 节省了流体求解时间.而且, 人们可以更加专注于区域边界来解决数值计算误差问题, 使边界元计算精确不断提高.而膜结构作为一种曲面结构, 其与流体的耦合恰恰发生在流场的边界, 边界元解可以直接用于膜结构变形计算, 这恰好利用了边界元的优势.

边界元计算模型通常是非光滑、有尖角尖边存在的, 导致了边界处法向物理量非连续, 进而导致尖角尖边处节点的边界条件难以确定, 如果不对尖角尖边进行特殊处理将使边界元计算精度大大下降, 特别是在上述区域, 故本文基于小扰动位势理论, 将非协调元引入三维曲面模型以解决尖角和尖边问题.通过算例7.1 (绕刚性半球三维势流流场计算) 可以发现, 引入非协调元后的边界元计算精度比未引入非协调元时大大提高, 非常接近解析解, 说明非协调元可以很好解决边界元计算中的尖角尖边问题, 也验证了本文边界元程序的可靠性.

为了发挥边界元法的优势, 本文引入了预条件GMRES法求解边界元方程组.GMRES法是求解大型系数矩阵非对称线性方程组的有效方法, 特别是经过预条件处理后, 计算达到收敛所需迭代步数少、精度高.通过算例7.2 (菱形马鞍面膜结构流固耦合计算) 可以发现, 百富策略白菜网预条件GMRES法后, 每次流体计算时间仅为85 min, 而百富策略白菜网高斯-若当消去法则需要花费10 h左右, 可见预条件GMRES法使得边界元法的优势得到了进一步的发挥.

本文流场计算基于小扰动位势理论, 流体为理想无旋流, 与实际有差异.这种差异主要体现在旋涡脱落区域和尾流区域, 实际风场中这些区域为湍流流动, 导致了风压分布的差异主要集中在结构边缘, 所以对于形状不规则的结构应慎用位势理论.尽管有差异, 却可以通过引入空气动力学领域成熟的方法来模拟或修正旋涡卷起区、脱落区及尾流区来得到可用于指导工程设计的风压值, 这将是本文下一步的工作.20世纪, 基于位势流的机翼理论的正确性使人们肯定了位势理论指导工程设计的重大意义;而当前, 潜艇设计以及许多低马赫数空气动力问题中仍在百富策略白菜网位势理论.目前, 在人们对膜结构特征湍流研究尚不深入、甚至流体力学研究者还没有对湍流机理有准确认识的情况下, 采取适当空气动力学方法 (比如涡方法) 对位势解进行修正应该是解决各类外形结构流固耦合问题的有效途径.