粘弹性圆柱绕流场的稀疏化动态模态分解方法.pdf

2024年3月第41卷第2 期枣庄学院学报JOURNALOFZAOZHUANGUNIVERSITYMar.2024Vol.41 NO.2粘弹性圆柱绕流场的稀疏化动态模态分解方法李璇，苏进（西安工程大学理学院，陕西西安7 10 0 48）【摘要】为研究粘弹性圆柱绕流涡团时空演化特征及状态，给出动态模态分解和稀疏促进优化相结合的流动关键模态分析方法（sparsity promoting dynamic modedecomposition，SP-D M D）。粘弹性圆柱绕流的数值模拟结果表明：利用SP-DMD方法所得到的模态不仅可以将原始流场的整体流动形态重现，而且能够较好地刻画出较大尺度的局部流场结构。与传统动态模态分解方法相比，SP-DMD能够识别和提取出具有稳定状态的粘弹性流场结构信息。该方法可为研究弹性诱导的复杂粘弹性流体流场特性提供良好参考。【关键词】粘弹性圆柱绕流；动态模态分解；稀疏化；特征提取【中图分类号0 351.2文献标识码A文章编号10 0 4-7 0 7 7（2 0 2 4）0 2-0 0 33-100引言粘弹性流体是一种具有粘性和弹性的特殊非时变性非牛顿流体，在生物工程等多种领域有着广泛应用。针对粘弹性圆柱绕流问题，不仅要考虑圆柱壁面（边界）的流动特点，还需考虑粘弹性圆柱绕流所具有的减阻特性。因而，分析弹性诱导的涡团结构特征对粘弹性流体的工程应用具有重要价值1-2。近年来，学者们对流体流动结构进行了较多研究：Thomases等3针对四辊轧机几何中的粘弹性流体，利用本征正交分解（proper orthogonal decomposition，PO D）提取了弹性模态，并根据其能量贡献分析了流体动力学特征；Henshaw等4采用POD量化流体动力学系统的模型结构；Shakeri 等5利用POD 提取了聚合物的弹性湍流结构并发现反向旋转涡对流动的总动能有显著贡献。最近在流体动力学领域内，动态模态分解6】（dynamic mode decomposition，D M D）广泛用于分析非定常流场7 的流动特征和构建低阶的流场动力学模型8-9。与POD 形成的完全基于空间相关性和能量含量的模态层次结构相比，DMD方法分解后的每个模态由空间的相关结构组成，在时间上具有相同的线性行为10。DMD可以看作是将奇异值分解（singular value decomposition，SVD）在空间降维方面的优点与快速傅里叶变换（fastfouriertransform，FFT）在时间频率识别方面的优点相结合。然而，在高速绕流条件下，粘弹性圆柱绕流场流体粘性减小阻力较低，随着速度的增加减阻效果逐渐减弱，这导致粘弹性圆柱绕流流场结构中有很多复杂的模态特征。DMD虽然在牛顿流体流场数据降维方面起到了一定的作用11-12，但是对于粘弹性圆柱绕流数据的减阻现象往往不能起到很好的降维效果。针对粘弹性圆柱绕流问题，为实现粘弹流高效的结构分析和动态预测，引人一种基于稀疏促进优化的DMD方法（SP-DMD），与传统的DMD相比，该方法可以剔除粘弹性流体中的非关键模态。在粘弹性流场演化中，SP-DMD使粘弹性流场中关键模态的数量和降维近似残差之间处于平衡状态，收稿日期】2 0 2 3-11-2 0基金项目】陕西省自然科学基础研究计划项目（2 0 2 3-JC-YB-063）。【作者简介李璇（19 9 9 一），女，云南昆明人，西安工程大学理学院数学专业2 0 2 1级研究生，主要从事数据驱动动力学的方法研究。苏进（198 3一），男，陕西榆林人，西安工程大学理学院副教授，理学博士，主要从事复杂流体的可计算建模与计算方面的研究。33枣庄学院学报从而更利于分析和发现那些能够对粘弹性圆柱绕流场发展起关键性作用的因素。此外，当模态数量增加时，SP-DMD重构的流场可以得到更多的局部流动细节，并且与原始流场进行对比，局部较小尺度的涡团形态以及空间分布的特点都能够很好地重现出来。1 DMD 算法从数学角度看，DMD是对Koopman算子模态的近似，从谐波的角度反映流场的特性。DMD方法的运用需要从数值模拟或物理实验中收集一系列快照x1,x2：，m+1，并形成一个数据矩阵：其中矩阵的每一列表示每一时刻的空间状态。此外，假定数据在时间上是等距采样的，时间步长为At。通常情况下，每一个快照x；=(it）是一个有n分量的复向量，即x；C”。从快照序列中构造两个快照矩阵：其中，X=C,X=C。并假设快照是由离散时间线性时变系统生成的，该系统表示如下：线性算子A记录着快照序列的动态特性。在粘弹性流体问题中，矩阵A中包含大量数据，且一般情况下为复数。每个快照x；中分量的数量通常远大于快照的数量，即n。根据两个连续时间步长之间的线性关系xk+1=Ax，k=1,2.,m，将数据矩阵X和X通过线性算子A连接起来，则矩阵X有如下表示：X=x2x.xm+=A x A x.因此，求出线性算子A是DMD的主要目标。由于快照x；的维数n较大，所以很难直接对线性算子A进行分析。而DMD算法绕过了A的直接特征分解，用基于POD模态矩阵U的投影矩阵A来表示A的降维近似矩阵。然后通过最小二乘法来确定矩阵A的形式。具体步骤如下。（1）对矩阵XCnxm进行奇异值分解（SVD）：X=UZV*=U Ure其中,Ue R,Ze R*,V R,Ue R*,e R,V R*,rm表示剩余的m-r个奇异值，*表示复共轭转置，并且U和满足假设数据总维数是m，上述SVD通过选择适当的截断值r来减少数据矩阵的维数，被消除的允余项就用rem（r e d u n d a n c y t e r m）表示，其维数是m。（2）在X的POD模态U的基础上，线性算子A可近似为A=CxrA=UAU*,A=UAU=UXV-l（3）矩阵A可以用快照X和X来确定，利用最小二乘法，即通过A=UAU*和X=UZV*来最小化X和AX差值之间的Frobenius范数：其中，矩阵Q的Frobenius范数由下式来决定：342024年第2 期(xX2.Xm+1)。X=xiX2.Xm门，X=X12:,xm+1。X+1=Ax,k=1,2,.,m。AxmJ=AX。U*U=I,V*V=Imin IIX-UAZV*II,Q,=trace(QQ)=trace(QQ*)。李璇，苏进粘弹性圆柱绕流场的稀疏化动态模态分解方法（4）矩阵AC在X的POD模态所跨域的子空间上确定了快照间映射AC的最优低维表示。因此，这个r维子空间的动力学系统为：且POD模态矩阵U可将映射到高维空间C”上，XUk。如果A有线性无关的特征向量（www），对应的特征值为（入，入,，则可以将其转化为对角坐标形式：A=wW2小k+I=Ak,入入2W(1)Z122*Z入其中,令W=wW2入2w,1,D二Z2*Z且12 2.：z,是水A*的特征向量，对应的特征值为1，.：动力学方程（1）所生成的每一组数据为：其中，;=z;表示初始条件为的第i个模态的贡献。（5）使用DMD模态的线性组合来近似实验或数值快照，即；=U w;，且每个;是对应DMD模态的振幅。同样，在矩阵形式中有：入21入2x2.m$12上式表明动态模态演化特性由范德蒙矩阵来控制，且范德蒙矩阵由矩阵A的个复特征值决定，其特征值包含相关时间频率和生长/衰减速率的信息。入2,D.=diag(),Vand令$=9 P2.,1,=则有（6）未知振幅向量D。=2入，将它们适当缩放以满足如下双正交条件：1,i=jlo,i牛j小=WD,Zi=2w入z4xxU=Z9/.he(1,2,m),1X dD.Vando,T 可通过下式来求解：11入，入入入”入入211入，入”35(2)枣庄学院学报使用X的截断奇异值分解=UV*和矩阵的定义=UW，将问题（2）转化为以下形式：(3)J）可以等价表示为：J)=p-q-q+s,式中：p=（*）（a n d V），q=d i a g(Va n d Vz*），s=t r a c e(Z*Z)，是两个矩阵中元素逐元相乘的乘积运算。（7）通过最小化关于的二次函数（4）可以得到该优化问题（3）的DMD振幅的最优矢量为：D=plq=(YY)(Vand Vna)-diag(VandVY)。2稀疏化动态模态分解（SP-DMD）方法DMD算法虽然能够求出振幅最优向量和提取出关键模态，但由于在维数较高的粘弹流中还存在部分非关键模态未剔除，因此引人稀疏化动态模态分解（SP-DMD）方法来剔除这些非关键模态。2.1SP-DMD核心思想SP-DMD需要在所有DMD模态中提取出若干关键的具有非零幅值的模态，然后通过算法来调整其幅值，具体方法如下：具有非零幅值的DMD模态可以通过求解如下凸优化问题而得到：式（5）中，矩阵E反映了振幅向量的稀疏结构信息。E的列均是单位向量，其非零元素对应于具有零幅值的DMD模态。获得具有非零幅值的DMD模态后，通过改变其幅值实现降维近似。将优化问题（3）进行修改，通过给式（3）中的目标函数J（）增加一个正则化项card(）来解决稀疏性的问题，这就惩罚了未知振幅向量D。中非零元素的数量：在修改后的优化问题（6）中，是正则化参数，它反映了“稀疏化”振幅向量D。的权重性。越大，说明越重视向量D。中非零元素的数量。因此，鼓励用（6）这个更稀疏的解决方案。一般来说，找到问题（6）的解决方案相当于组合搜索，对于任何问题，组合搜索会变得很棘手。为了绕过这个问题，引人问题（6）的更宽泛的表述，通过用向量D。的L1范数来代替基函数：SP-DMD的核心思想就是找到式（7）的解，这是一个凸优化问题，解决该问题使用了交替方向乘子法（alternatingdirectionmethod of multipliers，A D M M），在实验或数值快照降维近似质量和DMD模态数量之间达到预期的平衡后，从而确定最终优化后的稀疏化振幅向量D。2.2交替方向乘子法交替方向乘子法（ADMM）是一种用于求解优化问题的计算框架，适用于求解分布式凸优化问题。ADMM算法将大的全局问题分解为多个较小且较容易求解的局部子问题，并通过协调子问题的解而得到大的全局问题的解。ADMM通常用于求解两个优化变量且只含等式约束条件的优化问题，其一般表示形式为：其中，和z是优化变量，）和g（）都是凸函数。362024年第2 期min Il X-D.Vana Il,min.J()=/V*-WDaVmall,min J(),s.t.ET=O。min J()+card(z)。min J()+ll ll 1min f(x)+g(z),s.t.Ax+Bz=Co(4)(5)(6)(7)李璇，苏进粘弹性圆柱绕流场的稀疏化动态模态分解方法引人增广拉格朗日函数：L,(x,z,)=x)+g(z)+y(Ax+Bz-c)+(p/2)I Ax+Bz-c ll2,其中：入为拉格朗日乘子;p为惩罚系数且p0。由于ADMM算法是基于增广拉格朗日函数最小化的迭代算法，迭代方式如下所示：从初始点开始进行迭代，直到满足如下优化条件和停止准则：.k+1=argminL,(x,z,入),+=argmin,(x*,z,*),Z*+=*+p(Ax*I+Bz*+-c)。Z其中，epim，e d a 表示停止准则的阈值。3数值算例3.1粘弹性圆柱绕流模型参数设置为了模拟粘弹性圆柱绕流在高Re数下的湍流减阻现象，设置流场中的平均流速为U，圆柱直径的特征长度尺寸为d，几何特征停留时间为d/U。流体弹性特征的参数为Wi=入U/d，表征流体黏性特征的参数为 Re=pUd/m。选取如图1所示的几何示意图进行二维粘弹性圆柱绕流流减阻的模拟。8dualoL=30dL-10dLi-10dL2-20d图1粘弹性圆柱绕流的几何模型采用双分布LBM-IBM耦合算法模拟粘弹圆柱绕流，设置Re=1000，粘度比=n/m=0.01。数值模拟参数如表1所示。表1数值模拟参数1Ln10d220d30d10d宏观方程的处理采用IBM耦合D2Q9格子模型，D2Q9格子模型中碰撞迁移过程的边界处理为非平衡外推格式，区域上下边界采用无滑移条件，即duu=0,u=00。dydy出口处按充分发展条件处理，出口处的水平速度的法向导数为0，其中n是单位法向量，即37dx枣庄学院学报这里入口处的初始水平速度u设置成抛物状，初始垂直速度=0，人口的初始速度其中，Umax是通道水平中线处的最大速度且U=2Umax/3对于本构方程采用Oldroyd-B方程的格子Boltzmann模型，应力张量的分布函数上下壁面的边界条件应用非平衡态外推方法，计算区域的上下边界的应力分量为2024年第2 期0。dnu(0,y)=3Umax(10-y)/50,Tx=2Wi(1-)(duQuTxy(1dyTyy=0。人口和出口边界的应力分量采用充分发展边界条件，即aTxaTY00dxx3.2粘弹性圆柱绕流数值结果不同的Re（Re y n o l d s）数对应的Wi（W e i s s e n b e r g）不同，且Wi是Re的函数。根据3.1节所述，模拟出Re=1000时的粘弹性圆柱绕流流场的数据，并且对流体弹性特征数Wi=1的数据进行研究。每一时刻的流场数据构成一个30 16 0 1的网格。在研究粘弹性圆柱绕流问题的过程中，首先要将二维数据转化成一维，由此可知每一时刻所对应的空间维数高达18 0 90 1维。由于各种运动形态和不同结构的尾流流态会表现出不同的能量和频率，仅通过对瞬态流场的观察，很难对粘弹性圆柱绕流流场的特性进行深入理解。因此研究过程中，选择了连续10 0 0 个时刻的二维粘弹性流场数据进行DMD算法处理，可以得到30 2 个模态。得到的30 2 个DMD模态的特征值在复数坐标系内的分布情况如图2 所示，这些特征值进行对数化后的分布情况如图3所示。1.00.80.60.40.2OTxY0-0.2-0.4-0.6-0.8-1图2 DMD未对数化的特征值分布由图2 可知，大部分的DMD模态捕捉到的流场结构强度近似处于稳定的状态。此外，还可以发现有个别特征值落在单位圆内部，并且距离单位圆较远，此现象表示这部分特征值所对应的流场结构强度的衰减率比较大，流场结构的状态不稳定，将在随后的圆柱绕流流场演化过程中逐渐消38-10-12-140.80.60.4 0.2特征值虚部10.20.40.60.81.0-50-4030图3DMD对数化的特征值分布-20.10特征值虚部。1020304050李璇，苏进粘弹性圆柱绕流场的稀疏化动态模态分解方法散。由图3可知，每个特征值的实部（即衰减率）都基本接近于0，这更直观地表明了大部分DMD模态所对应的流场结构强度是近似稳定的。采用SP-DMD和标准DMD算法得到的粘弹性流场非对数化特征值的对比如图4（a）所示。由于该图不能清晰地看出特征值的分布情况，将其局部放大，如图4（b）所示。1.00.80.60.40.20-0.2-0.4-0.6-0.8-1在采用SP-DMD时，分别保留了30 个模态和10 0 个模态进行对比。SP-DMD关注的是DMD模态在整个粘弹性圆柱绕流流场的演化过程中对流场的贡献，通过SP-DMD能够更加方便地识别并提取在整个流场的演化过程中起到比较重要影响作用的模态信息。由图4可知，利用标准DMD算法所得到的具有较大衰减率的模态在 SP-DMD中都被有效剔除了，从而保留了有稳定状态的DMD模态。幅值与频率之间的关系可以定义为DMD频谱，在研究粘弹性流场问题的过程中，DMD频谱能够将流场结构所包含的能量情况反映出来，因此又将其称之为能谱。采用两种方法所获得的模态幅值与特征值虚部（频率）的关系如图5所示，对其进行局部放大如图6 所示。利用两种方法所得到的模态幅值与特值实部（衰减率）的关系如图7 所示，其局部放大如图8 所示。5000标准DMD算法?100个模态稀疏化DMD算法30个模态稀疏化DMD算法40001.011.000.99标准DMD算法100个模态稀统化DMD算法30个模态稀疏化DMD算法-0.8-0.6-0.4特征值虚部（a）特征值的全局分布图4DMD和SP-DMD未对数化的特征值分布0.980.970.960.95-0.20.20.40.6标准DMD算法100个模态稀统化DMD林法30个模态稀疏化DMD算法0.8-0.2515001000-0.2（b）特征值的局部分布-0.15-0.1特征值虚部-0.050.050.10.1530005002000标准DMD算法100个模态稀疏化DMD算法30个模态稀疏化DMD算法1000-1000-1.0（a）模态幅值实部与特征值虚部的关系图5采用DMD与SP-DMD得到的模态幅值与特征值虚部（频率）的关系-500-1000-1500-0.50特征值虚部0.51.0-1（b）模态幅值虚部与特征值虚部的关系-0.50特征值虚部0.539枣庄学院学报2001501002024年第2 期200标准DMD算法100个模态稀疏化DMD算法15030个模态稀蔬化DMD算法1005050-50-100-150-200-0.8-0.6-0.4-0.200.20.4 0.60.8特征值虚部（a）模态幅值实部与特征值虚部的关系图6 采用DMD与SP-DMD得到的模态幅值与特征值虚部（频率）的关系的局部放大5000标准DMD算法100个模态稀疏化DMD算法30个模态稀疏蔬化DMD算法400030002000-50-100标准DMD算法-150100个模态稀疏化DMD算法30个模态稀疏化DMD算法-200-0.8-0.6-0.4-0.2特征值虚部（b）模态幅值虚部与特征值虚部的关系1500标准DMD算法100个模态稀疏化DMD算法30个模态稀疏化DMD算法1000500000.20.40.60.81000-5000-1000-1000-1（a）模态幅值实部与特征值实部的关系图7 采用DMD与SP-DMD得到的模态幅值与特征值实部（衰减率）的关系10080604020-1500-0.50特征值实部0.5-1（b）模态幅值虚部与特征值实部的关系100标准DMD算法80100个模态稀疏化DMD算法30个模态稀疏化DMD算法：604020D-0.5特征值实部0.51-20-40-60-80-1000.7（a）模态幅值实部与特征值实部的关系图8 采用DMD与SP-DMD得到的模态幅值与特征值实部（衰减率）的关系的局部放大40-20-40-60标准DMD算法100个模态稀疏化DMD算法30个模态稀疏化DMD算法0.750.8特征值实部-80-1000.850.90.9510.7（b）模态幅值虚部与特征值实部的关系0.750.8特征值实部0.850.90.951李璇，苏进粘弹性圆柱绕流场的稀疏化动态模态分解方法由图5、6、7、8 可知，部分具有较大衰减率的非稳定的DMD模态所对应的流场结构承载着较多的能量，其对流场初期演化的影响较强，但在后期的演化过程中，影响力明显下降。SP-DMD将能够最大程度反映流场信息的稳定模态保留下来，并对保留的模态幅值进行调整，从而达到对原始粘弹性流场降维的效果。利用SP-DMD所获得的10 0 个模态与30 个模态对流场进行重构后的粘弹性圆柱绕流流场结构的空间分布如图9（）和图9（d）所示。为了验证SP-DMD有着较好的优化降维效果，采用标准DMD所获得的30 2 个模态对流场进行重构后的粘弹性圆柱绕流流场结构的空间分布如图9（b）所示。为了进行比较，原始的瞬态流场如图9（a）所示。(a)原始流场（c）SP-D M D 重构流场（10 0 个模态）（b）D M D 重构流场（30 2 个模态）由图9可知，用30 个模态对圆柱绕流流场进行重构就可以将原始流场的整体流动形态重现，并且能够较好地刻画出局部较大尺度的流场结构。当模态的数量增加到10 0 时，可以看到重构流场中显示出更多的局部细节。然而，标准DMD算法由30 2 个模态才能较好地刻画出原始流场的整体形态。由此可见，SP-DMD能够对将流场的非关键信息进行剔除，从而保留并提取具有稳定状态的关键流场结构。4结论针对粘弹性圆柱绕流涡团结构的研究，引入了DMD与稀疏促进相结合的特征分析方法。该方法在提取粘弹性流场关键模态时，不只考虑幅值，而且通过计算模态对流场贡献率来剔除非关键模态，从而留下产生影响较大的关键模态，使重构后的粘弹性流场仍然可以保留原始粘弹流流场中的重要动态信息。因此，利用SP-DMD可以更好地解析粘弹性圆柱绕流流场的特性，为粘弹性圆柱绕流流场的特性分析研究提供了新的视角。在未来的发展中，能否将包含控制输人的模型预测控制与SP-DMD相结合来解决实际问题，还有待进一步测试和改进。参考文献1 J ALVES M A,OLIVEIRA P J,PINHO F T.Numerical methods for viscoelastic fluid flows J.Annual Review of FluidMechanics,2021,53（1):50 9 -541.2 WU Y,LI P,TAO R,et al.Improvement of mode selection criterion of dynamic mode decomposition in a hydrofoil cavi-tation multiphase flow case J.Ocean Engineering,2022,265:112579.41（d）SP-D M D 重构流场（30 个模态）图9原始流场与DMD和SP-DND重构流场的对比枣庄学院学报3 THOMASES B,GUTIERREZ C P.Proper orthogonal decomposition(POD)of the flow dynamics for a viscoelastic fluidin a four-roll mill geometry at the Stokes limit J.Journal of Non-Newtonian Fluid Mechanics,2019,264:48-61.4 HENSHAW R J,MARTIN O G,GUASTO J S.Dynamic mode structure of active turbulence J.Physical Review Flu-ids,2023,8（2):0 2 310 1.5 J SHAKERI P,JUNG M,SEEMANN R.Characterizing purely elastic turbulent flow of a semi-dilute entangled polymer so-lution in a serpentine channel J.Physics of Fluids,2022,34（7):0 7 3112.6 KOU J Q,ZHANG W W.Dynamic mode decomposition and its applications in fluid dynamics JJ.Acta Aerodynami-ca Sinica,2018,36（2):16 3-17 9.7 LI B.Dynamic mode decomposition based algorithms as feature detection techniques:application to large flow databasesD.Madrid:Institute of Advanced Technology in Aerospace Engineering,Polytechnic University of Madrid,2021.8 SCHMID PJ.Dynamic mode decomposition and its variants J.Annual Review of Fluid Mechanics,2022,54(1):225-254.9J LEE K,TRASK N,STINIS P.Structure-preserving sparse identification of nonlinear dynamics for data-driven modeling CJ/Mathematical and Scientific Machine Learning.New York:Procedings of Machine Learning Research,2022(190):65-80.10 BADDOO PJ,HERRMANN B,MCKEON B J,et al.Kernel learning for robust dynamic mode decomposition:linearand nonlinear disambiguation optimization JJ.Proceedings of the Royal Society A,2022,478(2260):20210830.1IJ LI CY,CHEN Z,TSE TK T,et al.A parametric and feasibility study for data sampling of the dynamic mode decom-position:Range,resolution,and universal convergence states JJ.Nonlinear Dynamics,2022,107(4):3683-3707.12 秦文瑾，汪涛，齐观超，等。稀疏化动态模态分解算法在发动机缸内流场研究中的应用J内燃机学报，2020,38(5):426-432.【责任编辑：刘广德Sparsity Promoting Dynamic Mode Decomposition for Viscoelastic Flow Around Circular CylinderLI Xuan,SU Jin(School of Science,Xi an Polytechnic University,Xi an 710048,China)Abstract:In order to study the temporal-spatial evolution characteristics and states of vortexes around viscoelastic cylinders,afluid critical mode analysis method combining dynamic mode decomposition and sparse promotion optimization(SP-DMD)isproposed.The numerical simulation results of the viscoelastic flow around circular cylinder show that the modes obtained by theSP-DMD method can not only reproduce the overall flow pattern of the original flow field,but also describe the local flow fieldstructure of a larger scale.Compared with the traditional dynamic mode decomposition,SP-DMD can identify and extract theviscoelastic flow field structure information with stable state.This method can provide a good reference for studying the flow fieldcharacteristics of complex viscoelastic fluid induced by elasticity.Key words:viscoelastic flow around cylinder;dynamic mode decomposition;sparsity promoting;feature extraction2024年第2 期42