导弹与航天运载技术（中英文）

攻角$\left({}^{\circ }\right)$	来流速度$/\left({\mathrm{m}\cdot {\mathrm{s}}^{-1}}\right)$	高度 /km	密度$/\left({\mathrm{g}\cdot {\mathrm{{mm}}}^{-3}}\right)$	压强 /Pa
0	262	3	$9 \times {10}^{-7}$	$7 \times {10}^{4}$

攻角$\left({}^{\circ }\right)$	来流速度$/\left({\mathrm{m}\cdot {\mathrm{s}}^{-1}}\right)$	高度 /km	密度$/\left({\mathrm{g}\cdot {\mathrm{{mm}}}^{-3}}\right)$	压强 /Pa
0	262	3	$9 \times {10}^{-7}$	$7 \times {10}^{4}$

参数	密度 /$\left({\mathrm{{kg}}\cdot {\mathrm{{mm}}}^{-3}}\right)$	弹性模量 /GPa	泊松比	屈服强度 /MPa	延伸率
2A14	2 730	67	0.3	350	0.1

参数	密度 /$\left({\mathrm{{kg}}\cdot {\mathrm{{mm}}}^{-3}}\right)$	弹性模量 /GPa	泊松比	屈服强度 /MPa	延伸率
2A14	2 730	67	0.3	350	0.1

参数	来流速度$/\left({\mathrm{m}\cdot {\mathrm{s}}^{-1}}\right)$	高度 /mm	密度$/\left({\mathrm{g}\cdot {\mathrm{{mm}}}^{-3}}\right)$	压强 /MPa	温度 /K
值	262	$3 \times {10}^{6}$	$9 \times {10}^{-7}$	0.07	290

参数	来流速度$/\left({\mathrm{m}\cdot {\mathrm{s}}^{-1}}\right)$	高度 /mm	密度$/\left({\mathrm{g}\cdot {\mathrm{{mm}}}^{-3}}\right)$	压强 /MPa	温度 /K
值	262	$3 \times {10}^{6}$	$9 \times {10}^{-7}$	0.07	290

基于CESE/IBM算法的NACA0018翼型气动模拟及抖振分析

PDF下载

张旭童 ¹^,² , 于煜斌 ³ , 落龔寿 ³ , 郑新 ³ , 李京杰 ¹^,²

导弹与航天运载技术（中英文） | 运载器及导弹总体技术 2025,48(1): 17-28

收起

导弹与航天运载技术（中英文） | 运载器及导弹总体技术 2025, 48(1): 17-28

基于CESE/IBM算法的NACA0018翼型气动模拟及抖振分析

全屏

张旭童¹^,², 于煜斌³, 落龔寿³, 郑新³, 李京杰¹^,²

作者信息

¹大连理工大学,工业装备智能控制与优化教育部重点实验室,大连,116024

²大连理工大学,航天分离动力学与控制工程技术中心,大连,116024

³北京宇航系统工程研究所,北京,100076

张旭童(2000-),女,博士研究生,主要研究方向为分离动力学及控制。

于煜斌(1981—),男,博士,研究员,主要研究方向为飞行器总体设计、分离总体设计及分离结构设计。

落龔寿(1983—),男,博士,高级工程师,主要研究方向为分离总体设计及分离结构设计。

郑新(1984—),男,博士,高级工程师,主要研究方向为分离总体设计及分离结构设计。

李京杰(1992—),男,博士,讲师,主要研究方向为流固耦合计算和分离系统控制。

Aerodynamic Characteristics and Flutter Analysis of NACA0018 Airfoil based on CESE/IBM

Xutong ZHANG¹^,², Yubin YU³, Yanshou LUO³, Xin ZHENG³, Jingjie LI¹^,²

Affiliations

¹Key Laboratory of the Ministry of Education of Intelligent Control and Optimization of Industrial Equipment, Dalian University of Technology, Dalian ,116024

²Aerospace Separation Dynamics and Control Engineering Technology Center, Dalian University of Technology, Dalian ,116024;

³Beijing Institute of Astronautical Systems Engineering,Beijing,100076

出版时间: 2025-02-25 doi: 10.7654/j.issn.2097-1974.20250103

文章导航

摘要

收起

飞行器在飞行过程中,流体气动力和结构弹性体之间会形成耦合作用,这种相互作用可能会使弹性体发散抖动产生不同程度的破坏,造成安全隐患。针对以上问题,提出基于时-空守恒元解元和浸入式边界的流固耦合数值模拟方法,利用时-空守恒元解元方法进行流体域计算,利用浸入式贴体网格边界法进行流固耦合边界面识别。以NACA0018为例,通过仿真数据得出不同攻角下机翼外流场压力与速度云图;同时研究机翼在不同来流速度下的抖振幅度的变化及其规律。研究表明,该流固耦合求解方法,在求解具有复杂流型的高速可压缩流包括激波或爆轰与大变形问题上具有较高的精度及稳定性,可为相关研究提供借鉴。

关键词

流固耦合 / 时-空守恒元解元 / 浸入式边界 / 数值模拟 / 抖振

Abstract

收起

During the flight of aircraft, a coupling effect between the fluid aerodynamic forces and the structural elastomer will be formed, and this interaction may cause different degrees of damage to the divergence and jitter of the elastomer, resulting in safety risks. A numerical simulation method of fluid-structure coupling based on time-space conserved element solution and immersion boundary is proposed. The method of time-space conserved element solution is used to calculate the fluid domain, and the submerged body-fitted mesh boundary method is used to identify the fluid-structure coupling boundary surface. Taking NACA0018 as an example, the cloud image of wing outflow field pressure and velocity at different angles of attack is obtained through simulation data. At the same time, the buffeting amplitude of the wing under different inlet velocity and its rule are studied. The research shows that the fluid-structure coupling method has high accuracy and stability in solving high-speed compressible flows with complex flow patterns, including shock wave or detonation and large deformation problems, providing a reference research method for related research.

Key words

fluid-structure interaction / CESE / immersed boundary / numerical simulation / buffeting analysis

引用本文

张旭童, 于煜斌, 落龔寿, 郑新, 李京杰. 基于CESE/IBM算法的NACA0018翼型气动模拟及抖振分析. 导弹与航天运载技术（中英文）, 2025 , 48 (1) : 17 -28 . DOI: 10.7654/j.issn.2097-1974.20250103

Xutong ZHANG, Yubin YU, Yanshou LUO, Xin ZHENG, Jingjie LI. Aerodynamic Characteristics and Flutter Analysis of NACA0018 Airfoil based on CESE/IBM[J]. Missiles and Space Vehicles, 2025 , 48 (1) : 17 -28 . DOI: 10.7654/j.issn.2097-1974.20250103

正文

收起

0 引言

收起

流固耦合是指流体与结构体之间相互作用的一种普遍现象, 而流固耦合问题涉及的领域也很广泛, 如航空航天、生物工程、石油化工、船舶等。目前流固耦合仍然是工程中比较关注的问题^[1],人们对流固耦合的早期认识源于航空工程中的气动弹性问题。机翼是飞机重要结构之一, 在飞行过程中, 弹性体机翼在空气的作用下,会产生相应的弹性变形^[2],机翼结构的变形使其气动力发生变化, 而气动力又反作用于弹性体机翼, 这种相互耦合造成的作用会使机翼逐渐达到平衡状态, 否则会造成机翼不同程度的破坏, 严重时会导致机翼直接断裂,对飞机本身以及地面环境造成不可估量的损失^[3]。抖振现象是飞机结构对气流分离引起的非定常气流脉动压力的随机激振的响应, 通常是由于流场中的激波和边界层的耦合作用而产生的不稳定现象, 加载在弹性结构上的气动力呈现周期性。极限环振荡系统在非线性效应下的周期性振动与非线性因素有关, 通常表现为封闭的轨迹, 而抖振不同于极限环振荡是一种强迫振动。由于机翼具有大尺度与低刚度^[4],通过风洞试验研究抖振问题成本高、难度大,因此仿真模拟成为研究以上问题的主要方式。

1951年, Humphreys^[5]在进行风洞试验中发现抖振现象。由于抖振机理的复杂性, 各学者对其机理并没有统一的说法,其中Pearcey^[6]、Tijdeman^[7]和Lee^[8]等人对跨声速条件下的抖振特性做了大量的研究,其观点也各有不同。姚猛^[9]认为机翼抖振激励来自于引起机翼变形的分离流, 包括低速时大攻角条件下因分离流诱导的机翼抖振以及跨声速区激波-附面层相互作用诱导的机翼抖振。由于抖振风洞试验成本高^[10]以及计算流体力学的发展,对于抖振相关的数值模拟研究逐渐成为趋势。很多研究人员如张伟伟^[11]、韩冰^[12]、刘金利^[13]等对抖振特性数值模拟进行了研究, 包括选择合适湍流模型、跨声速下的抖振特性、大迎角下的结构抖振、数值模拟算法RANS、 DES、LES和DNS方法等,通过对典型翼型展开数值模拟研究, 分析了其在特定条件下抖振过程中的气动特性及流场结构。

针对以上研究, 本文提出了基于时-空守恒元解元(Space-time Conservation Element and Solution Element Method, CESE)与浸入式边界(Immersed Boundary Method, IBM)的数值模拟算法。

CESE方法是美国NASA Lewis研究中心的Chang^[14]博士提出的,在之后的发展与应用中不断被完善与改进。该方法是一种全新的守恒方程数值计算方法, 主要思想是将时间与空间统一处理, 不像大部分计算方法将两者分开处理; 在统一的时-空内, 引入守恒元(CE)与解元(SE)使得局部及整体满足守恒律；通过定义守恒元与解元便可推广至多维情况。与传统的差分格式相比, 在相同基点下可大大提高格式精度^[15]。主要应用于不定长可压缩流,例如高速复杂流动。经过十几年的发展, CESE方法也在不断地变化。Yang等^[16]提出了时空守恒元和解元方法的$\alpha -\mu$方案并证明了其收敛性,得到了一个误差限。Wang等^[17]通过改进二阶精度的二维CESE格式提出了新的$\mathrm{{CE}}$和$\mathrm{{SE}}$结构,并验证了新格式的准确性。Dwivedi等^[18]在$\alpha -\varepsilon -\alpha -\beta$格式的解中加入另一项来抑制由于不连续性而引起的摆动, 该方案能够同时捕获小扰动和尖锐间断。Shen等^[19]在定义的$\mathrm{{SE}}$上进行了改进, 通过两个不同守恒单元界面的数值通量不是直接由反时间方向的泰勒展开式导出的, 而是由一个迎风程序确定。Yu等^[20]提出了一种扩展的守恒元素/解元素方法模拟守恒定律与刚性源项。与迎风方案相比, 该方法不使用重构和黎曼求解器作为构建块, 逻辑要简单得多。在实际工程应用仿真中, 印华融等^[21]基于二维黏性CESE方法对爆轰波由小管径到大管径爆轰管的变化过程进行了数值仿真。崔树鑫等^[22]在CESE基础上结合非敏感克朗数,对激波翼型流场进行数值模拟。张永祥等^[23]在该方法的基础上采用任意四边形网格来离散区域, 进行了水位和流速的数值模拟, 结果表明, 该数学模型能够较好地模拟复杂条件下山区河流的流场特征。

浸入边界法最初是由Peskin^[24]提出的,用来模拟心脏中的血液流动问题。该方法将分布力源项加载至N-S方程中来处理复杂边界^[25],不需要生成贴体网格, 简化网格生成过程, 提高计算效率。其主要思想是由两套相互独立的网格构成, 流场采用固定直角网格, 流场中的物体采用一系列拉格朗日点来表示。仿真计算时, 流场中的物体与流场之间相互作用, 实现数据信息的传递, 对物体边界的流场信息进行处理。该法在模拟复杂物体中的巨大优势, 多用于不可压缩流领域, 近些年来在可压缩流领域的研究也逐渐广泛。Lai等^[26]提出并验证了一种形式上的二阶精度浸入边界法, 新格式具有更小的数值黏性, 模拟高雷诺数浸没边界流动的效果较好。$\mathrm{{Ye}}$等^[27]发展了一种直角坐标网格方法, 结合有限体积格式, 可模拟具有复杂浸入边界的二维非定常、黏性、不可压缩流动。 Hamid等^[28]在CESE-IBM FSI求解器中引入有限速率化学模型, 可以准确地预测氢/空气混合物在壁面前后的爆炸特性。

本研究采用基于CESE与IBM的数值模拟方法, 并与LS-DYNA结构有限元求解器耦合, 将流体压力、密度、结构形变等参数通过耦合边界进行传递与作用, 以NACA0018翼型为例进行数值模拟, 对其进行气动仿真模拟分析, 验证该算法的可行性, 并将机翼材料选用弹性材料进行数值模拟, 研究机翼在发生抖振时振动幅度与速度的变化规律。

1 数值模拟方法与验证

收起

1.1 CESE基本原理

CESE方法的主要特点是采用了时-空积分形式的控制方程, 其中时间和空间在同一框架下进行处理, 以二维可压缩流的$\mathrm{N}- \mathrm{S}$方程为例:

(1)$\frac{\partial \mathbf{U}}{\partial t}+ \frac{\partial {\mathbf{F}}_{\text{inv }}}{\partial x}+ \frac{\partial {\mathbf{G}}_{\text{inv }}}{\partial y}= \frac{\partial {\mathbf{F}}_{\text{vis }}}{\partial x}+ \frac{\partial {\mathbf{G}}_{\text{vis }}}{\partial y}$

式中$\mathbf{U}$为守恒变量向量;${\mathbf{F}}_{\text{inv }}$为$x$方向的通量向量(无黏);${\mathbf{G}}_{\mathrm{{vis}}}$为$y$方向的通量向量(黏性);${\mathbf{G}}_{\mathrm{{inv}}}$为$y$方向的通量向量(无黏);${\mathbf{F}}_{\mathrm{{vis}}}$为$x$方向的通量向量(黏性)。

(2)$\mathbf{U}= \left(\begin{matrix}\rho \\{\rho u}\\{\rho v}\\{\rho e}\end{matrix}\right)$

式中$\rho$为流体密度,$\mathrm{{kg}}/{\mathrm{m}}^{3};u, v$分别为流体流速在$x$方向和$y$方向上的分量,$\mathrm{m}/\mathrm{s}$；$e$为单位质量流体的总能量,$\mathrm{J}/\mathrm{{kg}}$。

(3)${\mathbf{F}}_{\text{inv }}= \left(\begin{matrix}{\rho u}\\\rho {u}^{2}+ p \\{\rho uv}\\\left({{\rho e}+ p}\right) u \end{matrix}\right),\;{\mathbf{G}}_{\text{inv }}= \left(\begin{matrix}{\rho v}\\{\rho uv}\\\rho {v}^{2}+ p \\\left({{\rho e}+ p}\right) v \end{matrix}\right)$

(4)${\mathbf{F}}_{\mathrm{{vis}}}= \left(\begin{matrix} 0 \\{\tau }_{xx}\\{\tau }_{xy}\\{\tau }_{xx}u +{\tau }_{xy}v -{q}_{x}\end{matrix}\right),\;{\mathbf{G}}_{\mathrm{{vis}}}= \left(\begin{matrix} 0 \\{\tau }_{xy}\\{\tau }_{yy}\\{\tau }_{xy}u +{\tau }_{yy}v -{q}_{y}\end{matrix}\right)$

式中$p$为静压,$\mathrm{{Pa}};\tau$为黏性应力张量;$q$为热通量向量; 基于牛顿黏性定律和斯托克斯假设,$\tau$可写为

(5)$\mathbf{\tau }= \mu \left({\left({\nabla \mathbf{u}}\right)+ {\left(\nabla \mathbf{u}\right)}^{\mathrm{T}}- \frac{2}{3}\left({\nabla \mathbf{u}}\right) I}\right)$

式中$\mu$为黏度,$\mathrm{{Pa}}\cdot \mathrm{s};u$为流体速度矢量;$I$为单位张量；$\mathbf{v}$为梯度算子。根据傅里叶热传导定律,热通量向量$\mathbf{q}$可写为

(6)$\mathbf{q}= -\kappa \nabla T $

式中$\kappa$为热传导率,$\mathrm{W}/\left({\mathrm{m}\cdot \mathrm{K}}\right)$;$T$为温度,$\mathrm{K}$。

当传输性质和状态方程被确定后,式(1)将变成一组封闭的方程组。为了简化式(1), 引入下式表示:

(7)$\mathbf{F}= {\mathbf{F}}_{\text{inv }}- {\mathbf{F}}_{\text{vis }},\;\mathbf{G}= {\mathbf{G}}_{\text{inv }}- {\mathbf{G}}_{\text{vis }}$

接着对时间和空间进行统一处理,定义$x$、$y$和$t$为三维空间的坐标系。式(1)的无发散形式可写为

(8)$\nabla \cdot \mathbf{h}= 0,\;\mathbf{h}= \left({\mathbf{F},\mathbf{G},\mathbf{U}}\right)$

将高斯散度定理应用于三维空间中的任意控制体积$V$中,$V$的边界面用$S\left( V\right)$表示,可将控制方程的形式写为

(9)${\oint }_{S\left( V\right)}\mathbf{h}\cdot \mathbf{n}\mathrm{d}S = 0 $

其中,$\mathbf{n}$表示边界面$S\left( V\right)$上的单位向外法向量。式(9)给出了流体流动中质量、动量和能量的时空守恒描述, 在CESE方法中作为构造数值格式的起始点。

以一维(1-D)的$\alpha -\alpha$格式为例,将CESE方法应用于这种形式的一维标量守恒型方程中可得:

(10)$\frac{\partial u}{\partial t}+ \frac{\partial f\left( u\right)}{\partial x}= 0 $

首先, CESE方法需要先将与计算相关的时空域进行离散化。离散化包括解点的划分, 其作用为选择待计算和存储的未知变量, 以及构建控制方程时-空积分形式的控制体。这些特征可以通过引入两个特殊的概念来表示: 守恒元CE和解元SE。

通过一维的物理空间进行均匀划分以及确定恒定的时间步长,可将二维的$x - t$(时间-空间)平面根据时空网格进行离散化处理, 如图1(红实线)所示。第$j$个网格节点可用空间坐标${x}_{j}$表示。每个单元中心位置可表示为${x}_{j + 1/2}= \left({{x}_{j}+ {x}_{j + 1}}\right)/2$,单元大小为${\Delta x}=$${x}_{j + 1}- {x}_{j}$。在以上规定的时间推进算法中,每步$\left\lbrack {{t}_{n - 1},{t}_{n}}\right\rbrack$由两个半步$\left\lbrack {{t}_{n - 1},{t}_{n - 1/2}}\right\rbrack$和$\left\lbrack {{t}_{n - 1/2},{t}_{n}}\right\rbrack$组成,时间步长可定义为${\Delta t}= {t}_{n}- {t}_{n - 1}$。

未知函数$u\left({x, t}\right)$由$u$在一组特定解点处的离散值表示。对于整数时间层$\left\{{{t}_{0},{t}_{1},\cdots ,{t}_{n}}\right\}$,解点则为网格节点(图1中的实心圆)。对于半整数时间$\left\{{{t}_{1/2},{t}_{3/2},\cdots ,{t}_{n + 1/2}}\right\}$,解点为单元的中心(图1中的空心圆), 因此在每个中间时间级采用的都是交错网格。在CESE算法方案中,在每个解点(j, n)处需要计算、存储的未知变量不止$u$,还有其空间导数${u}_{x}$:

(11)${u}_{j}^{n}\equiv u\left({{x}_{j},{t}_{n}}\right),{\left({u}_{x}\right)}_{j}^{n}\equiv \frac{\partial u}{\partial x}\left({{x}_{j},{t}_{n}}\right)$

图2显示了单个CESE时间步长下数据流的路径示意。如图2所示, CESE方案在空间和时间上形成了高度紧凑的样式。如果将半步视为基本迭代, 则对称模式的半宽度为${\Delta x}/2$,因为解点(j, n)处的未知数仅取决于存储在$\left({j - 1/2, n - 1/2}\right)$和$\left({j + 1/2, n - 1/2}\right)$处的数据。

在图1所示的时空网格中, 可以构造一组小单元的时-空守恒元素$\mathrm{{CE}}$。对于每个解点,都会有与之对应的一个CE。点(j, n)的$\mathrm{{CE}}$表示为${\left(\mathrm{{CE}}\right)}_{j}^{n}$,其形状为顶点是点$\left({j - 1/2, n - 1/2}\right)\text{、}\left({j + 1/2, n - 1/2}\right)\text{、}(j +$$1/2, n)$和$\left({j - 1/2, n}\right)$的矩形,由此若干$\mathrm{{CE}}$便可覆盖整个时空域, 无重叠。在两个连续的半时间步长中的$\mathrm{{CE}}$是交错分布的。在每个守恒元中,时-空积分形式由式(1)进行数值计算, 并建立未知数的离散方程。

当在$\mathrm{{CE}}$上对式(1)进行时-空积分时,一个重要的问题是如何沿$\mathrm{{CE}}$的边界计算$u$和$f$,这就需要为每个解点引入一个解元素${\left(\mathrm{{SE}}\right)}_{j}^{n}$,划分定义如图3所示,两条线段于交点处相互平分,形成$\left({j + 1/2, n}\right)$、$\left({j, n + 1/2}\right)\text{、}\left({j - 1/2, n}\right)$和$\left({j, n - 1/2}\right)$四个端点的交叉。用同样的方法, 可以定义半整数点的SE。每个${\left(\mathrm{{CE}}\right)}_{j}^{n}$元素都由3个解元素约束:${\left(\mathrm{{SE}}\right)}_{j}^{n},{\left(\mathrm{{SE}}\right)}_{j - 1/2}^{n - 1/2}$和(SE)${}_{j + 1/2}^{n - 1/2}$。SE覆盖了每一个CE的边界,无重叠。在每个$\mathrm{{SE}}$内部,若函数和定义为线性,且可以用一阶泰勒展开来近似表示$\mathrm{{SE}}$的中心。则在解元素${\left(\mathrm{{SE}}\right)}_{j}^{n}$中,$u$和$f$构造为

(12)$ u\left({x, t}\right)= {u}_{j}^{n}+ {\left({u}_{x}\right)}_{j}^{n}\left({x -{x}_{j}}\right)+ {\left({u}_{t}\right)}_{j}^{n}\left({t -{t}_{n}}\right),\left({x, t}\right)\in {\left(\mathrm{{SE}}\right)}_{j}^{n}$

(13)$ f\left({x, t}\right)= {f}_{j}^{n}+ {\left({f}_{x}\right)}_{j}^{n}\left({x -{x}_{j}}\right)+ {\left({f}_{t}\right)}_{j}^{n}\left({t -{t}_{n}}\right),\left({x, t}\right)\in {\left(\mathrm{{SE}}\right)}_{j}^{n}$

式中${u}_{t}$和${f}_{t}$分别为$u$和$f$的时间导数。

时空通量矢量定义为

(14)$\mathbf{h}= \left({f,\mu }\right)$

式中$f$和$\mu$分别为通量矢量$\mathbf{h}$在$x$方向和$t$方向上的分量。

考虑到从时间$n - 1/2$到时间$n$的半步,在解点(j, n)处两个独立的未知数${u}_{j}^{n}$和${\left({u}_{x}\right)}_{j}^{n}$需要同时计算。因此需要将积分守恒律离散化, 得到两个代数方程。以此,将守恒元${\left(\mathrm{{CE}}\right)}_{j}^{n}$分成${\left({\mathrm{{CE}}}^{- }\right)}_{j}^{n}$和${\left({\mathrm{{CE}}}^{+ }\right)}_{j}^{n}$两个子元素。如图4所示,${\left({\mathrm{{CE}}}^{- }\right)}_{j}^{n}$为矩形${ACDE},{\left({\mathrm{{CE}}}^{+ }\right)}_{j}^{n}$为矩形${CBFD}$。这些矩形的每条边都属于${\left(\mathrm{{CE}}\right)}_{j}^{n}$与相关联的3个$\mathrm{{SE}}$中的一个。

在每个子$\mathrm{{CE}}$上实现时-空积分守恒定律。控制体积$V$依次为${\left({\mathrm{{CE}}}^{- }\right)}_{j}^{n}$和${\left({\mathrm{{CE}}}^{+ }\right)}_{j}^{n}$,由此可得:

(15)${\oint }_{S\left({\mathrm{{CE}}}^{- }\right)}\mathbf{h}\cdot \mathbf{n}\mathrm{d}S = 0 $

(16)${\oint }_{S\left({\mathrm{{CE}}}^{+ }\right)}\mathbf{h}\cdot \mathbf{n}\mathrm{d}S = 0 $

式中$n$为${\left({\mathrm{{CE}}}^{- }\right)}_{j}^{n}$或${\left({\mathrm{{CE}}}^{+ }\right)}_{j}^{n}$边界上的单位向外法向量。

上式说明了每个子$\mathrm{{CE}}$的时空通量平衡。

1.2 IBM基本原理

浸入边界法网格点示意如图5所示,对于$2\mathrm{D}$物体的边界, 采用一系列拉格朗日点来表示, 其会因物体的移动、变形或流场的作用而发生改变; 而固定直角流场网格, 不会因物体的移动或变形而发生改变, 被称为欧拉网格点。

浸入边界通过将固体边界上的回复力离散成体力来实现对流场的作用。对于二维不可压黏流问题, 其控制方程为

(17)$\rho \left({\frac{\partial \mathbf{u}}{\partial t}+ \mathbf{u}\cdot \nabla \mathbf{u}}\right)+ \nabla p ={\mu \Delta }\mathbf{u}+ \mathbf{f}$

(18)$\nabla \mathbf{u}= 0 $

式中$u$为欧拉点上的速度矢量;$p$为压强;$\rho$为密度;$f$为流体体力项。

1.3 基于CESE和IBM的流固耦合方法

如图6所示, CESE流固耦合采用分区策略实现流体一结构相互作用(FSI)方式^[29],流场计算采用的是CESE方法(欧拉), 而物体结构分析则使用LS-DYNA中的结构求解器FEM (拉格朗日), 两者交替计算。在物体与流场的接触界面上, 采用浸入式边界识别方法, CESE求解器从物体结构模块处得到接触界面处的运动速度、位移等数据信息, 传递给物体结构模块压力数据, 以此作为物体结构的外部压力边界条件。CESE求解器中的流体自动与其中的固体结构发生耦合作用。

流固耦合计算程序的主要步骤为:

a)生成流体和固体结构的网格, 且相互独立, 进行初始化;

b)计算从固体结构边界到单个解元$\mathrm{{SE}}$点的最短距离如图7所示。根据这个距离,解元点被分成4类: 势流体点(图7中的A点)、虚拟流体点(图7中的B点)、流体附近点(图7中的C点)和内部流体点(图7中的D点);

c)基于结构加载和从流体求解器得到的流体-结构界面边界条件, 采用LS-DYNA FEM求解器求解固体结构方程;

d)从FEM求解器获得更新的流体-结构界面位置和界面速度;

e)为邻近界面解元点更新最短距离, 获得新的内部流体点、邻近流体点和虚拟流体点;

f)采用CESE算法为内部流体点更新流体解；

g)采用直接强制IBM法或常规CESE格式为流体附近点更新流体解;

h)采用虚拟流体法来处理虚拟流体点, 这是因为某些点要用于下一时间步内部流体点求解计算；

i)将流体压力反馈给FEM求解器, 作为外部力施加在流体-结构界面, 充当边界条件;

j)如果没有达到终止时间,则回到步骤$\mathrm{b}$。

1.4 算法验证

基于CESE/IBM的流固耦合算法, 在对NACA0018机翼模型仿真时, 为验证其算法可行性, 将其与典型的湍流算法模型进行仿真对比。

建立机翼及流场区域$3\mathrm{D}$模型并创建网格,在$\mathrm{{Dy}}$- na软件中, 将机翼与流场两套独立网格进行耦合处理, 机翼部分设置为刚性材料的壳体结构; 在Fluent软件中, 将机翼壁面设置为wall边界, 创建流场网格, 选用湍流模型算法进行模拟仿真, 对比试验网格数量均为15万左右。其中, 工况条件设置如表1所示。

仿真结果如图8所示,图8a和图8c分别为CESE算法针对NACA0018翼型的压力与速度云图,图8b和图8d分别为采用$k -\varepsilon$湍流算法得到的压力与速度云图。

对比Fluent计算结果可知, 两种算法均在机翼上下端面产生了低压区, 其前缘处的流体速度迅速增大超过声速；高压区均处于机翼前缘处,当空气与机翼相接触时, 前缘处空气压力迅速增大, 流速减小, 产生激波, 由于NACA0018翼型前缘沿弦线对称, 因此机翼上下端面、前缘后缘的压力及速度分布沿弦线对称, 如图8所示, 两种算法得出的压力与流体速度的变化趋势一致, 但是在流固耦合交界面附近处, CESE相较于Fluent的梯度变化的层次更丰富细致, CESE算法在捕捉激波上以及小扰动流的计算仿真上更灵敏。 CESE算法的流速最大值为${300}\mathrm{\;m}/\mathrm{s}, k -\varepsilon$湍流模型算法的流速最大值为${312}\mathrm{\;m}/\mathrm{s}$,相对误差为${4.58}\%$；CESE算法的压强最大值为${9.145}\times {10}^{4}\mathrm{\;{Pa}}, k -\varepsilon$湍流模型算法的表压最大值为${2.356}\times {10}^{4}\mathrm{\;{Pa}}$,总压最大值为${9.356}\times {10}^{4}\mathrm{\;{Pa}}$,即总压相对误差约为${2.26}\%$,由此得证, CESE流固耦合算法对模拟NACA0018翼型仿真可行,相对误差低于10%。

2 动力学建模

收起

浸入式边界算法是由两套相互独立的网格构成的, 因此需要分别建立流场与机翼结构的独立网格来实现数据信息的计算与传递。

2.1 结构模型

选用NACA0018翼型,机翼结构如图9所示,图9a为机翼二维截面,其中弦长为${100}\mathrm{\;{mm}}$,图9b为机翼三维结构, 展弦比为5, 建立机翼模型, 内部为空心结构,将其$3\mathrm{D}$模型进行网格划分处理。

建模过程中, 为简化计算过程, 提高后续仿真计算效率, 将机翼结构设定为壳单元进行简化处理。在对其进行网格划分时, 提取机翼结构表面来代替其整体结构, 进行surface域划分。由于机翼结构厚度较薄,对SHELL单元设定厚度参数为$1\mathrm{\;{mm}}$实现机翼结构的设计。图10为机翼的网格划分图, 表面区域采用共节点网格进行划分, 共包含机翼弧形表面及机翼侧面两边界面,命名为B_1和B_2。在后续探究抖振特性时,需将B_2靠近机身一侧固定。

在研究机翼抖振特性时, 涉及到机翼结构的位移变化及变形, 壳体采用多线性本构模拟弹塑性, 铝合金材料2A14的参数如表2所示。

其中, 铝合金材料2A14多线性弹塑性本构曲线如图11所示。

2.2 气动模型

机翼的外流场区域设置为长度为${3500}\mathrm{\;{mm}}$、宽度为${800}\mathrm{\;{mm}}$、高度为${2700}\mathrm{\;{mm}}$的长方体区域。由于浸入式边界方法两套网格互相独立, 流场区域内部不受机翼结构影响, 但建模需与机翼结构统一坐标系。根据流体区域的结构与作用特征, 其边界分为速度入口、速度出口及流场壁面三个边界面,图12为流体域的网格划分及边界面设置, 采用六面体进行体网格划分, 将流体域内的机翼位置及机翼四周网格进行细化, 随着距机翼位置变远, 流场网格尺寸增大。流体域网格内部细化网格与粗略网格共节点处理。

3 仿真模拟与结果分析

收起

本节基于CESE与浸入式边界方法对搭建的机翼结构模型及流体域模型进行仿真计算。对NACA0018翼型在不同攻角下的气动特性进行分析, 总结在不同攻角下流体压力及速度的变化规律, 通过改变来流速度, 研究机翼在不同工况下的位移变化并进行分析, 总结其抖振幅度及速率变化规律。

3.1 气动仿真

为提高计算效率, 简化计算模型, 研究机翼不同攻角下气动特性时, 将机翼壳单元设置为刚性材料进行模拟计算。本节研究了NACA0018翼型分别在攻角为$-{9}^{\circ }\text{、}- {4}^{\circ }\text{、}{0}^{\circ }\text{、}{5}^{\circ }\text{、}{8}^{\circ }\text{、}{12}^{\circ }$时的外流场的压力与速度分布, 工况条件参数设置如表3所示, 选取高度$3\mathrm{\;{km}}$、当地声速${328}\mathrm{\;m}/\mathrm{s}$、行速度${262}\mathrm{\;m}/\mathrm{s}$为初始条件, 在关键字*CESE_CONTROL_SOLVER中可设置流固耦合算法, 选用浸入式边界算法。

仿真模拟结果如图13所示, 图中是攻角分别为$-{9}^{\circ }\text{、}- {4}^{\circ }\text{、}{0}^{\circ }\text{、}{5}^{\circ }\text{、}{8}^{\circ }\text{、}{12}^{\circ }$时的压力与速度云图。

在负攻角飞行时,由图13a、图13c可观察到, 在机翼的上翼面靠近前缘的位置处大气中有一个压力的高峰值区, 该区域的压力高峰值是由气流在接触到机翼瞬间的高速流动所造成的；随着飞行方向与弦线的夹角变大,高压区域逐渐接近机翼前缘,机翼下表面低压区域逐渐向后缘移动, 气流最大压强值也会随之变大；如图13b、图13d所示,对应机翼上下翼面气流速度变化较大, 下翼面流速明显高于上翼面, 压强大对应流速小。

当机翼处于正攻角飞行时,图13e、图13g、图13i、图13k为正攻角时外流场压力云图,对比负攻角飞行压力云图,上翼面由高压区域转变为低压区域, 下翼面对应转变为高压区域, 较负攻角飞行时升力变大; 随攻角变大, 下翼面的高压区域逐渐沿下翼面向后缘移动,同时气流压强最大值也由攻角为${5}^{\circ }$时的${0.0806}\mathrm{{MPa}}$增大到${12}^{\circ }$时的${0.1024}\mathrm{{MPa}}$；但是攻角为${0}^{\circ }$时气流压强最大值为${0.9146}\mathrm{{MPa}}$,这是因为前缘的曲率半径最大, 当气流接触到机翼前缘时便产生高压区域；由速度云图13h、图13j、图131可观察到流速最大值也在随着攻角的增大而增大,攻角为${12}^{\circ }$时,最大流速达到${325.4}\mathrm{\;m}/\mathrm{s}$,符合机翼流场随攻角变化的变化规律。

3.2 抖振数值模拟

本节研究NACA0018翼型在不同来流速度下的抖振特性, 主要分析其幅值及速率变化过程。机翼选用2.1节所述铝合金材料2A14,具体参数如表1所示, 同时,在机翼上翼面处设置6个参考点(A1~A6), 如图14所示,并将机翼边界$\mathrm{B}\_ 2$靠近$\mathrm{A}1$端面固定, 采用表2工况条件所示的压强、密度、温度、高度参数数值, 进行模拟仿真。其中, 自变量为来流速度, 分别设置为${65}\mathrm{\;m}/\mathrm{s}$、${131}\mathrm{\;m}/\mathrm{s}$、${164}\mathrm{\;m}/\mathrm{s}$、${262}\mathrm{\;m}/\mathrm{s}$、${328}\mathrm{\;m}/\mathrm{s}$、${394}\mathrm{\;m}/\mathrm{s}$。

流体质点的运行速度$V$定义为

(19)$ V ={Ma}\cdot c $

式中${Ma}$为马赫数;$c$为当地声速。$3\mathrm{\;{km}}$高空处,当地声速为${328}\mathrm{\;m}/\mathrm{s}$。由此设置的来流速度分别对应马赫数0.2、0.4、0.5、0.8、1、1.2。

通过仿真得到的${262}\mathrm{\;m}/\mathrm{s}$速度下最外侧观测点位移变化最大时刻的幅值如图15所示,图16为${262}\mathrm{\;m}/\mathrm{s}$速度下外侧位移最大时刻的应力变化。其中应力最大区域为与固定面相近的上翼面区域, 最大值为${27.72}\mathrm{\;N}$。由铝合金材料$2\mathrm{\;A}{14}$在${Ma}= {0.8}$下的飞行状态, 可得出该材料的最大应力情况, 符合变化规律。该数据可用于后期对机翼进行结构优化设计, 通过对比分析出最佳结构,同时结合不同速度下机翼的抖动频率及幅值变化,也可达到有效的抑制与控制设计, 体现数据实用性。

不同来流速度下,观测点沿$Y$轴方向的幅值变化如图17所示。由图17可以观察到, 随着与固定面的位移增大,最大幅值也在变大,整体最外侧观测点A6位移变化最大, 由图15可见, 机翼整体幅值变化由固定面向最外侧自由面近似线性增大。由图17a至图${17}\mathrm{\;d}$可知,观测点的大幅抖振周期约为${16}\mathrm{\;{ms}}$,大幅振动后随之会有约为$1 \sim 2\mathrm{{ms}}$的细微抖动,振幅整体呈现周期性变化。当速度较小时,每一周期内的振幅趋势大致相同, 当速度由亚声速增大至跨声速时, 抖振周期增至${19}\mathrm{{ms}}$,$Y$轴抖振垂向位移的最大值也呈增大趋势。第1周期的幅值略高, 并随时间增大, 幅值缓慢降低, 直到各周期幅值呈现稳定相同的变化趋势。并且当跨声速飞行时, 第1周期的幅值与稳定周期的幅值差距也变大, 小幅抖动也较明显, 这是由于跨声速下飞行情况复杂造成的。

随着来流速度的增大,可以观察到$\mathrm{A}6$点(即$Y$方向抖振位移最大点)的最大位移也逐渐增大,图18为随速度变化的最大位移变化。可见, 低声速到亚声速飞行阶段, 最大幅值的变化增加缓慢, 亚声速到超声速飞行时, 位移变化随速度增大呈直线上升趋势。

以来流速度${262}\mathrm{\;m}/\mathrm{s}$为例,其$X$轴与$Z$轴方向的抖振幅值变化如图19所示。$X$轴与$Z$轴方向上的位移变化较$Y$轴较小,且整体呈现周期性变化。

图20为机翼在$X$、 $Y$、 $Z$方向上的抖振速度。 $\mathrm{A}6$单元在各方向的抖振速度均高于其他单元, 且随着距固定面的位移减小而减小, 固定面处的A1单元抖振位移与速度变化几乎为零。由图20b可见, $Y$方向上的速度变化波形近似呈周期性正弦形式, 但速度变化曲线的抖动较为明显,速度大小变化也均高于$X$、 $Z$方向。由于机翼抖振通常由流场自身的不稳定性引起, 当在特定条件下运动时, 流场中的涡旋、激波、分离等现象会导致机翼产生抖动,机翼的抖振变化主要集中于$Y$方向上的变化,受到气动的干扰较大,故曲线也不如$X$、 $Z$方向曲线平滑。

4 结论

收起

本文以NACA0018机翼翼型为例, 形成了完整的时-空守恒元解元CESE及浸入式边界IBM的流固耦合算法模型, 通过气动仿真与抖振的数值模拟, 得到结论如下:

a)与Fluent仿真结果对比,流速与压力云图符合规律,流速最大值的相对误差为${4.58}\%$,总压最大值的相对误差大约为$7\%$,由此得证, CESE流固耦合算法在模拟NACA0018翼型仿真上可行, 相对误差低于10%。

b)以攻角为自变量,研究其在来流速度为${262}\mathrm{\;m}/\mathrm{s}$时的流速与压力云图,规律变化结果符合相关文献的结果, 负攻角飞行时的升力低于正攻角飞行的升力,并且随着攻角变大,速度最大值可达${325.4}\mathrm{\;m}/\mathrm{s}$。由于浸入式边界方法对计算机及网格的要求较高, 故在刚性边界处的细节并不是很明显。

c)以来流速度为自变量,研究${0}^{\circ }$攻角下的抖振特性,得出机翼的抖振幅值在沿$X\text{、}Y\text{、}Z$方向呈现周期性变化, 且随着来流速度变化, 低声速飞行至亚声速阶段, 最大幅值的变化增加缓慢, 但从亚声速至超声速阶段,抖振位移变化呈直线上升趋势。

基于CESE与IBM的流固耦合求解方法中, CESE方法能够提供较快的计算速度和较高的精度, 尤其是在处理具有复杂边界条件的问题时; IBM方法在处理流固耦合问题时, 虽然在处理复杂几何形状方面具有优势, 但可能需要更多的计算资源来确保精度。因此, 在求解具有复杂流型的高速可压缩流包括激波或爆轰与大变形问题上该方法具有较高精度及稳定性, 为相关研究提供了可借鉴的研究方法。

基金

收起

辽宁省重点研发计划(2023JH2)
辽宁省重点研发计划(101800025)

参考文献

收起

文献

收起

参考文献引证文献

排序方式：

[1]

邵学博. 基于ANSYS的机翼的流固耦合分析[D]. 哈尔滨: 哈尔滨工程大学, 2012.

SHAO

Xuebo

. The fluid-solid interaction analysis of an airfoil based on ANSYS[D]. Harbin: Harbin Engineering University, 2012.

[2]

岳承宇. 基于板模型的柔性机翼气动弹性分析与控制[D]. 南京: 南京航空航天大学, 2014.

YUE

Chengyu

. Aeroelastic analysis and control of flexible wings based on plate model[D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2014.

[3]

郭镤泽. 亚声速下机翼振动的流固耦合模拟[D]. 杭州: 浙江大学, 2013.

GUO

Puze

. Simulation of fluid-structure interaction about airfoil vibration at subsonic speed[D]. Hangzhou: Zhejiang University, 2013.

[4]

冒森, 张斌, 肖良华, 等.太阳能无人机机翼颤振动力学建模与分析[J]. 计算机测量与控制, 2023,31(9): 183-189.

MAO

Sen

, ZHANG

Bin

, XIAO

Lianghua

, et al.Dynamic modeling and analysis of solar UAV wing flutter[J]. Computer Measurement & Control, 2023,31(9): 183-189.

[5]

Pressure pulsations on rigid airfoils at transonic speeds[R]. NACA RM, 1951.

[6]

A method for the prediction of the onset of buffeting and other separation effects from wind tunnel tests on rigid models[R]. Reston: AGARD, 1958.

[7]

TIJDEMAN

. Investigations of the transonic flow around oscillating airfoils[D]. Delft: Technische Hoge School Delft, 1977.

[8]

LEE

BH K

, MURTY

, JIANG

. Role of Kutta waves on oscillatory shock motion on an airfoil[J]. AIAA Journal, 1994,32(4): 789-796.

[9]

姚猛. 机翼跨声速抖振特性数值模拟研究[D]. 南京: 南京航空航天大学, 2017.

YAO

Meng

. Numerical study of airfoil transonic buffet characteristics[D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2017.

[10]

董璐. 飞机抖振特性N-S方程计算[D]. 南京: 南京航空航天大学, 2007.

DONG

. Numerical studies of aircraft buffet with navier-stokes equations[D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2007.

[11]

张伟伟, 高传强, 叶正寅. 机翼跨声速抖振研究进展[J]. 航空学报, 2015,36(4): 1056-1075.

ZHANG

Weiwei

, GAO

Chuanqiang

, YE

Zhengyin

. Research advances of wing/airfoil transonic buffet[J]. Acta Aeronautica Sinica, 2015,36(4): 1056-1075.

[12]

韩冰. 三角翼大迎角动态气动特性及结构抖振数值模拟[D]. 西安: 西北工业大学, 2014.

HAN

Bing

. Numerical investigation of unsteady aerody-namics of delta wing and structural buffeting[D]. Xi'an: Northwestern Polytechnical University, 2014.

[13]

刘金利. 翼型跨声速激波振荡机理与抖振载荷减缓研究[D]. 西安: 西北工业大学, 2016.

LIU

Jinli

. Study on mechanism of transonic shock oscillation and alleviation of buffet load on the air-foil[D]. Xi'an: Northwestern Polytechnical University, 2016.

[14]

CHANG

S C

. The method of space-time conservation element and solution element: a new approach for solving the Navier-Stokes and Euler equations[J]. Journal of Computational Physics, 1995,119(2): 295-324.

[15]

刘凯欣, 王景焘, 王刚, 等.时-空守恒元解元(CE/SE)方法综述[J]. 力学进展, 2011,41(4): 447-461.

LIU

Kaixin

, WANG

Jingtao

, WANG

Gang

, et al.A review on CE/SE method[J]. Advances in Mechanics, 2011,41(4): 447-461.

[16]

YANG

, YU

, ZHAO

. Convergence and error bound analysis for the space-time CESE method[J]. Numerical Methods for Partial Differential Equations: An International Journal, 2001,17(1): 64-78.

[17]

WANG

Gang

, ZHANG

Deliang

, LIU

Kaixin

. An improved CE/SE scheme and its application to detonation propagation[J]. Chinese Physics Letters, 2007,24(12): 3563-3566.

[18]

DWIVEDI

, SUPRADEEPAN

, PARTHASARADHI

R G

. A review of space-time conservation element solution element (CESE) schemes[J]. Journal of Physics: Conference Series, 20191276:012007.

[19]

SHEN

, WEN

C Y

. A characteristic space-time conservation element and solution element method for conservation laws II. multidimensional extension[J]. Journal of Computational Physics, 2016305:775-792.

[20]

S T

, CHANG

S C

, YU

S T

, et al.Treatments of stiff source terms in conservation laws by the method of space-time conservation element/solution element[C]// 35th Aerospace Sciences Meeting and Exhibit. Reno: 1997.

[21]

印华融, 翁春生. 爆轰波在变截面突扩管道中的传播特性研究[J]. 航空兵器, 2016(6): 66-72.

YIN

Huarong

, WENG

Chunsheng

. Characteristics of detonation wave in variable propagation cross-section sudden expansion pipe[J]. Aviation Weapons, 2016(6): 66-72.

[22]

崔树鑫, 高琳, 高歌. 基于CE/SE法的翼型流场激波捕捉计算[J]. 北京航空航天大学学报, 2012,38(12): 1591-1595.

CUI

Shuxin

, GAO

Lin

, GAO

. Numerical simulation based on CE/SE method for airfoil flow with shock wave[J]. Journal of Beijing University of Aeronautics and Astronautics, 2012,38(12): 1591-1595.

[23]

张永祥, 陈景秋, 文岑, 等.时空守恒元和解元法在山区河流模拟中的应用[J]. 重庆建筑大学学报, 2007(6): 49-52.

ZHANG

Yongxiang

, CHEN

Jingqiu

, WEN

Cen

, et al.Numerical simulation of 2-D flow in the mountain rivers by using space-time conservation element and solution element scheme[J]. Journal of Chongqing Jianzhu University, 2007(6): 49-52.

[24]

PESKIN

C S

. Flow patterns around heart valves: a numerical method[J]. Journal of Computational Physics, 1972,10(2): 252-271.

[25]

王力, 田方宝. 浸入边界法及其在可压缩流动中的应用和进展[J]. 中国科学: 物理学力学天文学, 2018,48(9): 196-204.

WANG

, TIAN

Fangbao

. Recent progress of immersed boundary method and its applications in compressible fluid flow[J]. SCIENTIA SINICA Physica, Mechanica & Astronomica, 2018,48(9): 196-204.

[26]

LAI

, PESKIN

C S

. An immersed boundary method with formal second-order accuracy and reduced numerical viscosity[J]. Journal of computational Physics, 2000,160(2): 705-719.

[27]

, MITTAL

, UDAYKUMAR

H S

, et al.An accurate Cartesian grid method for viscous incompressible flows with complex immersed boundaries[J]. Journal of computational physics, 1999,156(2): 209-240.

[28]

HAMID

, MIRZABABAIE

T M

. Tracking the explosion characteristics of the hydrogen-air mixture near a concrete barrier wall using CESE IBM FSI solver in LS-DYNA incorporating the reduced chemical kinetic model[J]. International Journal of Impact Engineering, 2023172:104401.

[29]

ZHANG

, COOK

JR G

, IM

. Overview of the CESE compressible fluid and FSI solvers[C]// 16th International LS-DYNA Users Conference. Online: 2020.

2025年第48卷第1期

PDF下载

398

176

引用本文

BibTeX

文章信息

doi: 10.7654/j.issn.2097-1974.20250103

接收时间：2024-03-14
首发时间：2025-07-21
出版时间：2025-02-25

补充材料

相关文章

文章信息

作者

出版历史

收稿日期：2024-03-14
修回日期：2024-12-08

基金

辽宁省重点研发计划(2023JH2)

辽宁省重点研发计划(101800025)

作者信息

¹大连理工大学,工业装备智能控制与优化教育部重点实验室,大连,116024

²大连理工大学,航天分离动力学与控制工程技术中心,大连,116024

³北京宇航系统工程研究所,北京,100076

参考文献

分享链接

https://castjournals.cast.org.cn/joweb/ddyht/CN/10.7654/j.issn.2097-1974.20250103

分享至

全文二维码

扫描看全文

引用本文

BibTeX

本文的引用情况

2种不同金属材料的力学参数

科 Family	属数 Number of genus	种数 Number of species	占总种数比例 Percentage of total species (%)	属 Genus	种数 Number of species	占总种数比例 Percentage of total species (%)
鹅膏菌科Amanitaceae	2	11	5.26	鹅膏菌属 Amanita	10	4.78
小菇科 Mycenaceae	2	12	5.74	丝盖伞属 Inocybe	5	2.39
多孔菌科 Polyporaceae	8	14	6.70	蜡蘑属 Laccaria	5	2.39
红菇科 Russulaceae	3	23	11.00	小皮伞属 Marasmius	6	2.87
				小菇属 Mycena	11	5.26
				光柄菇属 Pluteus	5	2.39
				红菇属 Russula	17	8.13
				栓菌属 Trametes	5	2.39

关闭全屏

BibTeX
EndNote
RefWorks
TxT