风与膜结构流固耦合作用强耦合方法的降阶模型研究
发布时间:2021年9月23日 点击数:1723
引言
风与膜结构的流固耦合作用是被研究人员广泛关注但尚未深入研究的前沿课题之一。目前研究膜结构风振中的流固耦合效应的数值方法包括弱耦合分区法、强耦合分区法和强耦合整体法[1,2]。文献[3,4]中,作者引入线弹性有限元模型推导出求解风与膜结构流固耦合作用的强耦合整体方法并在膜结构中得到了百富策略白菜网。研究表明强耦合整体法对于膜结构这种大变形结构显示出其优势:精度高、计算稳定、收敛性好。但由于强耦合整体方法是在同一时间步内求解流体和结构方程,耗费计时多,计算繁琐。降阶模型是研究非线性流固耦合问题的重要手段之一[5,6],使用降阶模型可以减少计算计时,节省计算资源,一般建立降阶模型常采用的方法是本征正交分解法(POD),它是从时域或频域计算中提取结构的动力信息[7,8,9] 。目前降阶模型在流固耦合问题中的研究还比较有限,且大都针对分区法展开[10,11]。
本文针对之前提出的强耦合整体方法,提出了基于本征正交分解(POD)的降阶模型,利用POD方法对强耦合整体方程的一系列解进行基模态分解,得到一组基向量,并百富策略白菜网Galerkin方法将各控制方程对各基向量进行投影重构,构造基于POD基的降阶模型,将所求解的流固耦合问题的自由度下降几个量级,达到减少计算时间,提高计算效率的目的。
1 强耦合整体方法简介
强耦合整体方程的推导就是对流体域、结构域以及线弹性模型方程进行时间和空间上的离散的过程。流体方程、结构方程和线弹性方程均采用Galerkin有限元法进行离散,时间上的离散采用隐式有限差分方法,具体推导过程可参见文献[3,4]。膜结构风致流固耦合问题的加权残差方程 fFSI可以写作[3],
fF=∫Ωf(ω∇⋅νf+−ζ⋅ρf(∂νf∂t+νf⋅∇νf)−∇ζ∶σf+ζ⋅ff)dΩ=0fF=∫Ωf(ω∇⋅νf+-ζ⋅ρf(∂νf∂t+νf⋅∇νf)-∇ζ∶σf+ζ⋅ff)dΩ=0 (1)
fS=∫ΩS0(∇η∶σl-η·fs)dΩ=0 (2)
fLE=∫ΩS0(∇η∶σLEllLE)dΩ=0 (3)
fFSIF=∫ΓfN∪ΓFSIζ·(σf·nf1)dΓ=0 (4)
fFSIS=∫ΓSN0∪ΓFSI0((η·ns0)∶σl)dΓ+∫ΓSN0∪ΓFSI0((η·ns0)∶σl)dΓ=0 (5)
式中, Ω是边界由ΓD,ΓN,ΓFSI构成的空间域,Ωf是流体域,ΓN和ΓFSI分别是Neumann边界和流固交界面;ω,ζ,η是检验函数。方程中的未知量包括流体速度νf,压力p、结构位移us和线弹性模型位移uES。
求解强耦合整体方程时采用Newton-Raphson方法线性化,求解方程组A(xi−1tti-1)xitti=b(xi−1t−1t-1i-1,xi−1tti-1),引入边界条件对方程进行迭代运算,得到t时刻第i次迭代的未知量xitti。
2 基于本征正交分解的降阶模型
这里将采用本征正交分解(POD)建立上述强耦合整体方程的降阶模型,其主要思路是利用POD方法对上述强耦合整体方程的一系列解进行基模态分解,得到一组基向量,并百富策略白菜网Galerkin方法将各控制方程对各基向量进行投影重构,构造基于POD基的降阶模型,将所求解的流固耦合问题的自由度下降几个量级,达到减少计算时间,提高计算效率的目的。
这里引入3个有限维度子空间Uhs⊂Us,Vhf⊂Vf,UhLE⊂ULE,其中h>0与网格尺寸相关,假定POD在这3个子空间上的基向量分别为Ψi,Ζi,Φi,i=1,…,N,且基向量数量N满足条件,∑i=1Nλi/∑i=1Mλi>a,a>0.999∑i=1Νλi/∑i=1Μλi>a,a>0.999,λi表示的是POD的第i个特征值,M是快照数[12],各子空间上的快照矩阵可以写作,
Mvf=[M1vf|⋯|MNvf],Mus=[M1us|...|MNus],MuLE=[M1uLE|...|MNuLE]Μvf=[Μvf1|⋯|ΜvfΝ],Μus=[Μus1|...|ΜusΝ],ΜuLE=[ΜuLE1|...|ΜuLEΝ] (6)
则对于t∈[0,T],未知量对应的基向量空间可以写作,
Vf=({Ψi}Nfi=1),Us=({Zi}Nsn=1),ULE=({Φi}NLEn=1)Vf=({Ψi}i=1Νf),Us=({Ζi}n=1Νs),ULE=({Φi}n=1ΝLE) (7)
在POD的基向量空间上的求解域内定义内积,
(m,n)=∫Ωm(x,y,z)n(xy,z)dΩ (8)
在上述缩减空间上采用Galerkin投影方法获得整个流固耦合问题的降阶模型,那么对于∀t∈[0,T],n=1,…,N,有
∑i=1NdαidtAt1in+∑i=1N∑j=1Nαi(t)aj(t)At2ijn+∑i=1Nαi(t)Atin+Dt1n=0∑i=1ΝdαidtA1int+∑i=1Ν∑j=1Ναi(t)aj(t)A2ijnt+∑i=1Ναi(t)Aint+D1nt=0 (9)
∑i=1Nβi(t)∑i=1Νβi(t)Btinint+Dt2n2nt=0 (10)
∑i=1Nγi(t)∑i=1Νγi(t)Ctinint=0 (11)
其中,
At1in=∫Ωfω(x,y,z,t)∇ψi⋅ψndΩ,At2ijn=∫Ωf−ζ(x,y,z,t)ρf(x,y,z,t)(∇ψi⋅ψj)⋅ψndΩ,Atin=∫Ωf∇ζ(x,y,z,t):σf(ψi)dΩ,Dt1n=∫Ωfζ(x,y,z,t)ff⋅ψn⋅ndΩ,Btin=∫Ωs0∇η(x,y,z,t):σl(Zi)dΩ,Dt2n=∫Ωs0η(x,y,z,t)fs⋅Zn⋅ndΩ,Ctin=∫Ωs0∇η(x,y,z,t):σLEl(Φi)dΩ。A1int=∫Ωfω(x,y,z,t)∇ψi⋅ψndΩ,A2ijnt=∫Ωf-ζ(x,y,z,t)ρf(x,y,z,t)(∇ψi⋅ψj)⋅ψndΩ,Aint=∫Ωf∇ζ(x,y,z,t):σf(ψi)dΩ,D1nt=∫Ωfζ(x,y,z,t)ff⋅ψn⋅ndΩ,Bint=∫Ω0s∇η(x,y,z,t):σl(Ζi)dΩ,D2nt=∫Ω0sη(x,y,z,t)fs⋅Ζn⋅ndΩ,Cint=∫Ω0s∇η(x,y,z,t):σlLE(Φi)dΩ。
vf=∑i=1Nαi(t)Vfvf=∑i=1Ναi(t)Vf,us=∑i=1Nβi(t)Usus=∑i=1Νβi(t)Us,uLE=∑i=1Nγi(t)ULEuLE=∑i=1Νγi(t)ULE (12)
当上述计算收敛时,进行下一次迭代计算,得到第(i+1)个迭代步的流固耦合系统的未知量:
ui+1s≈μfuissi, vi+1f≈μfvif , ui+1LELEi+1=μLEuiLELEi (13)
其中,
μf=[Ψ1…ΨNf],μs=[Ζ1…ΖNs],μLE=[Φ1…ΦNLE]。
图1给出了本文降阶模型的流程和主要步骤。
3 算例分析
采用降阶模型对风与某椭圆形膜结构的流固耦合问题进行研究。该椭圆形膜结构的几何形状和尺寸如图2所示,预应力为2 kN/m,薄膜厚度为1 mm,密度ρ=643.8 kg / m3,泊松比ν=0.14, 入口处风速为10.32 m/s,风向角为0°~90°。流域上壁面为滑移边界,下壁面和流固耦合边界均为无滑移边界。当相对残差满足‖b-Axn‖2≤10-4‖b‖2时停止迭代。
为了研究本文降阶模型的特性,在膜结构上取点1~4,其坐标分别为1(x=0, y=0),2(x=a/2, y=0),3(x=0, y=b/2),4(x=a/2, y=b/2)。
这里采用数值模拟的方法获得快照,基于FLUENT18.0平台,按照上述初始条件,速度场采用无滑边界条件,计算得到238组速度场,即为快照。再百富策略白菜网本征正交分解(POD),得到基向量及对应的特征值,得到的每一个特征值可以用来刻画所对应的模态所捕捉到的能量,也就是说特征值的大小可以反应该基向量在整个系统中所占的权重。表1给出了前20阶模态特征值和能量累计比例。
表1 前20阶特征值和能量比 导出到EXCEL
Table 1 First 20 eigenvalues and energy percentage
| i | 特征值 | 能量累计比例 | i | 特征值 | 能量累计比例 |
|
1 |
0.423 | 0.246 | 12 | 0.479 | 0.687 |
|
2 |
0.324 | 0.304 | 13 | 0.641 | 0.705 |
|
3 |
0.285 | 0.352 | 14 | 0.482 | 0.779 |
|
4 |
0.463 | 0.395 | 15 | 0.257 | 0.835 |
|
6 |
0.268 | 0.423 | 16 | 0.317 | 0.893 |
|
7 |
0.332 | 0.465 | 17 | 0.542 | 0.903 |
|
8 |
0.437 | 0.524 | 18 | 0.472 | 0.932 |
|
9 |
0.375 | 0.563 | 19 | 0.272 | 0.966 |
|
10 |
0.211 | 0.603 | 20 | 0.577 | 0.983 |
|
11 |
0.289 | 0.656 |
从表1数据可以看出,各场量的前20阶特征向量能够捕获系统约99%的能量,因此选择前20个特征向量作为降阶模型的POD正交基,可以替代原全阶模型。
为了说明降级模型中POD基向量对强耦合整体方程解的影响,表2给出了不同风向角下采用整体方程计算该膜结构得到的未知量残差随POD基向量数量的变化。
表2 不同风向角下整体方程未知量残差随POD基向量数量的变化
Table 2 Variation of unknown residuals with POD basis vector at different wind azimuths 下载原表
从表2可以看出,虽然风向角不同,但随着降阶模型中POD基向量N的增加,强耦合整体方程的未知量流体速度νf、压力p、结构位移us和线弹性模型位移uES的计算残差均呈减小趋势,且N越大,计算残差下降越快,说明了该降阶模型很好的还原了整体方程的未知量特性,从一个侧面证实了该降阶模型的正确性和有效性。
这里采用降阶模型计算了上述膜结构4个典型测点1~4在不同风向角(0°和90°)的位移时程响应,并与全阶模型结果进行了对比,如图3~图4所示。
由图3~图4可以看出,在不同的风向角作用下,本文的降阶模型的位移时程变化趋势和未降阶前的变化趋势一致,且结果比较一致,说明经过快照预处理的降阶模型对原模型的近似程度较高,误差的原因可能是由降阶模型中采用数值模拟生成的快照造成的。
为验证本文降阶模型计算的正确性,表3和表4分别给出了采用全阶模型和降阶模型计算,当采用不同迭代步时,点1~点4的风压系数平均值和风压系数极值对比。
表3 不同风向角下采用全阶模型和降阶模型的风压系数平均值 导出到EXCEL
Table 3 Average wind pressure coefficient using full-order model and reduced-order model under different wind direction angles
|
|
模型 |
风向角 |
||
|
0° |
45° | 90° | ||
|
点1 |
全阶模型 |
-0.43 | -0.79 | -0.57 |
|
降阶模型 |
-0.35 | -0.82 | -0.49 | |
|
点2 |
全阶模型 |
0.21 | -0.65 | -0.76 |
|
降阶模型 |
0.24 | -0.73 | -0.89 | |
|
点3 |
全阶模型 |
-0.63 | -1.24 | 0.83 |
|
降阶模型 |
-0.57 | -1.06 | 0.68 | |
|
点4 |
全阶模型 |
1.22 | -1.54 | -2.51 |
|
降阶模型 |
1.35 | -1.39 | -2.17 | |
表4 不同风向角下采用全阶模型和降阶模型的风压系数极值 导出到EXCEL
Table 4 Wind pressure coefficient extremum using full-order model and reduced-order model under different wind direction angles
|
|
模型 |
风向角 |
||
|
0° |
45° | 90° | ||
|
点1 |
全阶模型 |
-1.25 | -1.04 | -0.86 |
|
降阶模型 |
-1.07 | -1.22 | -0.95 | |
|
点2 |
全阶模型 |
0.87 | -0.91 | -1.23 |
|
降阶模型 |
0.94 | -1.04 | -1.04 | |
|
点3 |
全阶模型 |
-1.12 | -1.43 | 0.57 |
|
降阶模型 |
-0.97 | -1.26 | 0.76 | |
|
点4 |
全阶模型 |
1.85 | -2.14 | -3.05 |
|
降阶模型 |
1.67 | -1.94 | -3.57 | |
分析表3和表4可以看出,采用全阶模型和降阶模型时各测点在不同风向角下的风压系数平均值和极值分布比较一致,两种模型的平均风压系数误差平均为12.31%,极值误差平均为13.62%,说明本文提出的降阶模型是比较准确的。
表5给出了采用降阶模型和全阶模型时的计算效率对比,主要对比了耦合迭代次数和耗费计时在不同计算收敛容差下的变化。
表5 采用不同方法计算位移时程的平均耦合迭代次数和计算计时 导出到EXCEL
Table 5 Comparison of mean couplingiterations and computing time using different methods
|
计算耗时 收敛容差 |
10-2 | 10-3 | 10-4 | 10-5 | |
|
耦合迭代次数 |
全阶模型 |
23 | 35 | 48 | 63 |
|
降阶模型 |
12 | 16 | 21 | 31 | |
|
耗时(s) |
全阶模型 |
27 654 | 48 654 | 58 976 | 79 435 |
|
降阶模型 |
18 773 | 30 652 | 39 213 | 52 469 |
分析表5可以看出,达到同样的计算容差时,采用本文的降阶模型要比全阶模型的计算效率高,其中耦合迭代次数比全阶模型平均减少约50%,耗费机时比全阶模型平均减少了约36%,主要原因在于降阶模型实质上未引入新的方程或数学模型代替原强耦合整体方程,只是在求解过程中将所得解向低维子空间进行投影,所以大大降低了计算成本。
4 结论
本文针对作者之前提出的强耦合整体方程,提出了基于POD的降阶模型,将该降阶模型百富策略白菜网于膜结构与风荷载的流固耦合分析中,取得了较好效果,得到的主要结论有:
(1)随着降阶模型中POD基向量N的增加,强耦合整体方程的未知量的计算残差均呈减小趋势,且N越大,计算残差下降越快。
(2)本文的降阶模型要比全阶模型的计算效率高,其中耦合迭代次数比全阶模型平均减少约50%,耗费机时比全阶模型平均减少了约36%,大大降低了计算成本。
(3)采用POD对强耦合整体方程建立强耦合整体方法的降阶模型是可行的。











