空气动力学学报

姓名
邮箱
手机号码
标题
留言内容
验证码

《空气动力学学报》2021年3期pdf合集

2021, 39(3): 1-1.

摘要(84) HTML (43) PDF(38)

摘要:

吴希明, 牟晓伟

2021, 39(3): 1-10. doi: 10.7638/kqdlxxb-2021.0042

摘要(590) HTML (243) PDF(118)

摘要:
本文针对直升机特有的技术特点，尤其是其区别于固定翼飞机的飞行原理、复杂气动和结构动力学特性、极致重量设计等，对包括气动、振动、噪声、地面共振/空中共振、抗坠毁、安全性及重量控制等方面的直升机技术难点进行剖析，对问题产生的原因、机理进行深入浅出的分析，对当前技术发展情况及解决方法进行了阐述。随后对世界直升机技术发展方向进行了分析与探讨，提出未来直升机技术将以高速、绿色、智能化、无人及安全飞行等方向为重点，并结合我国直升机技术发展现状和未来发展趋势，提出了我国直升机技术未来发展的三个重要方向—先进直升机旋翼桨叶翼型设计技术、智能旋翼技术和高速直升机技术，提出了各发展方向的关键技术的科学内涵和解决途径，为我国直升机未来发展提供思路。

高压捕获翼构型的跨流域气动特性

田鹏, 李广利, 崔凯, 李志辉, 张俊

2021, 39(3): 11-20. doi: 10.7638/kqdlxxb-2020.0144

摘要(226) HTML (58) PDF(19)

摘要:
高压捕获翼（HCW）构型是一种满足高速飞行器高容积、高升力、高升阻比的设计需求的新型气动布局。最近研究表明，HCW构型能够提高飞行器在连续流区的升力和升阻比，缓解飞行器设计中高容积率与高升阻比间的矛盾。为探究该气动布局在过渡流域（70～100 km）的气动特性，以一种楔—平板组合的高压捕获翼原理性构型作为模型，采用直接模拟Monte Carlo （DSMC）方法，详细分析了该模型在典型高超声速条件（马赫数20）下的流场结构和壁面气动力/热分布。结果表明，随着飞行高度增加，稀薄效应增强，机体压缩产生的激波厚度增加，激波边缘逐渐模糊，机体与捕获翼之间的开放通道内出现压力干扰。同时，高压捕获翼表面的摩擦系数迅速上升，气动摩擦成为制约捕获翼构型升阻比的重要因素。针对这一问题，分析了捕获翼材料表面的适应系数对飞行器的气动力/热的影响，结果表明，降低适应系数可以显著减小壁面摩擦和热流量，可通过选用适应系数较小的表面材料进一步提高该类飞行器气动性能。

黎曼边界条件在高阶精度非结构有限体积方法中的验证与应用

孔令发, 刘伟, 董义道

2021, 39(3): 21-32. doi: 10.7638/kqdlxxb-2020.0087

摘要(226) HTML (130) PDF(21)

摘要:
黎曼边界条件是一种弱施加边界条件。通过引入有限波模型，对亚声速入口、出口以及远场边界可采用精确求解黎曼问题来统一处理，有效简化了此类边界条件的施加过程，避免了基于特征关系式与黎曼不变量的推导，并已在二阶精度非结构有限体积方法中取得了较好的数值表现。为进一步验证该边界条件的实用价值，将其推广至高阶精度非结构有限体积离散。通过基于制造解方法（Method of Manufactured Solutions, MMS）的流动、亚声速无黏圆柱绕流及添加初始高斯脉冲扰动的非定常流动这三类数值算例，分别检验了黎曼边界条件在高阶精度非结构有限体积求解器中的数值表现。从计算结果来看，施加黎曼边界条件不会破坏离散格式的设计精度，同时，相比基于一维黎曼不变量的无反射边界条件，黎曼边界条件的施加过程简便，且维持了较好的出口特性，为基于非结构有限体积方法的高精度数值模拟提供了一种更加简单有效的亚声速边界处理方式。

冷却塔群塔干扰条件下脉动风压分布模式

赵林, 陈翰林, 钱永丰, 李敬生, 葛耀君

2021, 39(3): 33-43. doi: 10.7638/kqdlxxb-2020.0180

摘要(199) HTML (40) PDF(7)

摘要:
冷却塔群塔为常见的群体建筑物组合形式，存在复杂的塔群干扰效应。传统的群塔干扰效应研究通常面向塔筒等效静风荷载，忽略了群塔组合对于脉动风压分布模式干扰效应的量化研究。为了研究群塔组合条件下脉动风压的干扰分布模式，选取具备多种气动干扰特点的八塔组合矩形布置方案为研究对象，对群塔进行刚性模型同步测压风洞试验，获取了不同风向角条件下塔筒表面的脉动风压分布。在“迎风来流”、“侧风来流”和“背风来流”等三种典型来流模式下，结合风向角的变化将脉动风压分布形式总结为“V”型、“M”型、半“V”型+半“M”型三类，探讨了干扰影响脉动风压分布的机理，发现可将脉动荷载的干扰效应分为“弱干扰”、“非对称强干扰”和“对称强干扰”三种类别；在量化表达脉动风压干扰分布方面提出了脉动风压分布的九项式拟合公式，并按照三类分布模式给出了脉动风压的推荐表达式。

一种基于STL文件的LBM前处理高效算法

张晓潇, 张武, 褚学森, 王良军, 杨广文

2021, 39(3): 44-50. doi: 10.7638/kqdlxxb-2020.0009

摘要(194) HTML (62) PDF(15)

摘要:
LBM是一种计算流体数值方法，计算过程中需要确定流场格点属于流体点还是固体点，通过前处理获得格点类型信息。其中一种方法是通过STL文件格式描述的物体模型判断格点在物体内外信息来确定格点类型，从而实现STL几何信息到LBM计算模型信息的重构。为了能够快速重构计算模型，本文提出了一种快速生成计算模型算法。该算法根据面三角形找到模型的边界点集，然后由边界点判定出模型内部点和外部点从而完成对模型的重建。算法从面三角形出发，大量减少了计算量，节约了时间，提高了效率。基于本文算法和直接法对不同复杂程度的圆球、NACA0012翼型、CHN-T1飞机标模开展了前处理建模比较，结果表明，随着几何复杂度增加，直接法耗时急剧增加，而本文算法一直保持低耗时（例如CHN-T1模型前处理网格生成，直接法采用120核并行运算耗时11 h，而本文算法采用单核仅耗时20 s），极大地提高了LBM针对复杂几何外形的前处理效率。使用基于本文算法前处理获得的网格开展流场计算，验证了该前处理方法的适用性。

基于改进的k-ω-γ转捩模式预测高超声速飞行器气动特性

王亮, 周玲

2021, 39(3): 51-61. doi: 10.7638/kqdlxxb-2019.0146

摘要(285) HTML (61) PDF(41)

摘要:
为验证边界层转捩对高超声速飞行器气动特性的影响，采用改进的k-ω-γ转捩模式对类X-51A高超声速飞行器进行了全机边界层转捩预测。不仅系统分析了飞行攻角和雷诺数对边界层转捩的影响规律，同时研究了边界层转捩对飞行器气动力和进气道性能的影响。发现边界层转捩对飞行器升力系数和俯仰力矩系数影响较小，对阻力系数影响较大。针对本文算例，全层流计算较转捩计算预测的阻力系数偏低可达20%～30%。此外，边界层转捩可减小高超声速飞行器前体压缩面拐角处的分离，降低喉道截面马赫数、提高增压比。以上研究结果可为高超声速飞行器的控制系统和推进系统设计提供技术参考，显示了改进的k-ω-γ转捩模式具有较大的工程应用潜力。

黏性边界层修正对前体/进气道一体化乘波布局气动性能影响分析

王旭东, 王江峰, 程克明, 陈江涛

2021, 39(3): 62-70. doi: 10.7638/kqdlxxb-2020.0021

摘要(200) HTML (153) PDF(19)

摘要:
针对基于流线追踪生成的乘波进气道在内流部分存在较大黏性误差的问题，以构建具有消波隔离段的内外流一体化乘波气动布局为目的，发展了一种多级波系锥导或吻切锥乘波体的边界层黏性修正方法。在基于特征线法的多级前体与截短Busemann进气道一体化乘波布局设计方法基础上，耦合了高超声速可压缩轴对称流动的冯•卡门动量方程积分方法，通过修正边界层位移厚度来提高黏性情况下与设计预期的吻合度。对修正前后的一体化乘波构型进行了设计状态下的流场对比验证以及非设计状态下给定马赫数和迎角范围下的气动性能影响分析。计算结果显示，设计状态下黏性修正构型的前体两级激波和入口反射激波位置与设计预期基本一致，消除了隔离段内多次反射激波，进气道流量比增加4.51%，出口总压损失减少8.23%；非设计状态下修正构型同样具备更强的来流捕获能力和更小的出口冲压损失。所发展的方法能够提升内外流一体化乘波布局设计方法的精度，适用于吸气式高超声速飞行器前体/进气道的精确设计。

考虑滑流影响的分布式螺旋桨布局优化设计

杨伟, 范召林, 吴文华, 余雷

2021, 39(3): 71-79. doi: 10.7638/kqdlxxb-2020.0189

摘要(120) HTML (46) PDF(26)

摘要:
相对于传统螺旋桨飞机，前缘分布式螺旋桨飞机的滑流对整个机翼周围的流场特性影响更大。在分布式螺旋桨飞机动力布局设计过程中，应当考虑螺旋桨布局对滑流和机翼绕流的影响，使得机翼在滑流作用下表现出更好的气动特性。首先，基于雷诺平均Navier-Stokes方程，用等效盘代替真实螺旋桨，并采用重叠网格方法进行网格布局，优化过程中快速自动生成流场计算网格，获得机翼在滑流影响下的气动数据。然后，基于量子粒子群寻优算法建立分布式螺旋桨布局优化设计系统，使用该系统对采用NACA4412翼型的平直机翼和5个前缘分布式螺旋桨组成的构型进行了分布式螺旋桨布局优化设计。结果表明：等效盘模型能够对分布式螺旋桨滑流进行高效的数值模拟；优化后的分布式螺旋桨动力布局能使构型在起飞状态下升力系数提高5.6%，同时阻力系数减小13.9%。因此，建立的考虑滑流影响下机翼气动特性的分布式螺旋桨动力布局优化设计方法是可行且有工程应用价值的。

舰船纵摇突变对舰面流场的影响

李通, 王逸斌, 赵宁

2021, 39(3): 80-89. doi: 10.7638/kqdlxxb-2019.0136

摘要(135) HTML (39) PDF(9)

摘要:
舰船在海上航行时会受到海浪的影响而产生横摇、纵摇等不同方向上的运动，同时海浪的随机性也会造成舰船摇摆参数发生变化，扰乱舰面流场，增加舰载直升机在甲板上方作业的危险性。为了保障舰载直升机在舰船上的起降和其他甲板作业的安全性，现基于两栖攻击舰的简化模型，通过计算流体力学方法，在纵摇状态下对舰面流场进行数值模拟，将不规则的纵摇运动简化为具有周期性的简谐运动，采用动网格方法来实现舰船的摇摆。首先定性分析了甲板上方涡结构的发展以及舰船纵摇周期和振幅的突变对流场结构带来的影响，其次通过垂向速度的变化对动态舰面流场做了进一步的定量分析。结果表明：当舰船突然开始纵摇或纵摇周期减半后，上层建筑和舰尾附近的位置要滞后于舰首一段时间才能建立起周期性流动，同时纵摇的加速导致了流场的恶化。当纵摇振幅减半后，相当于舰船的“静息期”，流场变化平缓，垂向速度变化范围明显减小，有利于舰载直升机的降落。对于不同的纵摇状态，甲板上方各观测点垂向速度的波动都随高度的降低而增强，特别是舰首靠近甲板区域的波动最强。舰载机飞行员在进行甲板作业时应针对纵摇突变带来的影响，及时采取相应措施以保障作业安全。

S1223翼型低雷诺数下分离泡及气动特性

王皓田, 朱杨柱, 车学科, 李修乾

2021, 39(3): 90-98. doi: 10.7638/kqdlxxb-2020.0075

摘要(427) HTML (170) PDF(63)

摘要:
采用表面测压技术，测量了低雷诺数下(Re = 6.0×10⁴、1.0×10⁵、2.0×10⁵) S1223翼型的表面压力分布，通过时均化处理及瞬态处理方法，分别获得了翼型稳态和瞬态压力系数、升力系数，分析了流场结构随雷诺数及攻角的变化规律，研究了雷诺数及攻角对翼型升力的影响机理。结果表明，从时均升力系数随攻角的变化规律来看，S1223翼型在低雷诺数下存在“静态滞回”效应。攻角由负逐渐增大至0°时，下翼面由完全分离转变为出现层流分离泡，随后分离泡逐渐减小直至消失，导致升力系数斜率呈现随攻角逐渐增大的非线性现象。当攻角超过临界攻角后，不同雷诺数下翼型流场结构随攻角的变化规律存在本质不同，Re = 6.0×10⁴和1.0×10⁵时，翼型周围流场迅速发生大范围流动分离，升力系数迅速减小；而Re = 2.0×10⁵时，上翼面周期性生成短泡，引发低频振荡现象，升力系数呈现准周期性变化，α = 16°时上翼面时均流场呈现40%弦长的长泡结构。

近壁湍流中微小非球形颗粒取向行为研究综述

崔智文, 赵立豪

2021, 39(3): 99-108. doi: 10.7638/kqdlxxb-2021.0045

摘要(142) HTML (60) PDF(21)

摘要:
本文回顾了近十年关于微小非球形颗粒在壁湍流中取向行为的数值研究进展，重点介绍了数值方法和物理机理方面的发现。由于壁面的存在，在近壁面区域存在较强的平均剪切以及复杂的湍流结构，使得壁面的湍流具有较强的各向异性特点。前期研究发现杆状颗粒在壁面附近会倾向性地朝着流向方向，而碟状颗粒会倾向于朝着壁面法向，且这种倾向性取向行为会随着颗粒形状偏离球形的程度的增加而增强。同时，研究发现颗粒的倾向性取向行为与流体的拉格朗日的拉伸与压缩方向具有较强相关性。但是，颗粒的取向与流体的拉格朗日的拉伸与压缩方向在壁面附近的相关性要弱于其在远离壁面的槽道中部或均匀各向同性湍流的相关性，且存在较明显的形状敏感性。前期工作以杆状颗粒为例分析了其相对于拉格朗日坐标系(由左柯西格林张量三个主轴方向确定)的概率密度分布的特点，并进一步阐释颗粒行为与拉格朗日的拉伸与压缩方向的差异。同时，基于壁面湍流的特点，建立了以平均剪切与速度梯度脉动之比的颗粒转动周期预测的二维简化模型，并进一步揭示了平均剪切与速度梯度脉动在颗粒取向行为趋近于拉格朗日拉伸与压缩方向过程中的重要作用。

微纳颗粒两相流系统中湍流脉动效应研究综述

杨海林, 林建忠

2021, 39(3): 109-120. doi: 10.7638/kqdlxxb-2021.0030

摘要(265) HTML (80) PDF(16)

摘要:
给出了微纳颗粒两相湍流中颗粒一般动力学方程，分析了该方程包含的因湍流脉动导致的未封闭源项；回顾了湍流脉动对颗粒成核、生长、凝并以及破碎作用的研究，包括从源项核函数的构造发展到对方程进行雷诺平均或亚格子滤波；叙述了用直接数值模拟（DNS）方法、蒙特卡洛法、大涡模拟方法（LES）以及雷诺平均方法（RANS）对方程进行求解的研究成果；介绍了对颗粒一般动力学方程进行封闭求解的研究途径，包括引入相关函数建立未封闭源项与颗粒性质、湍流场特性的关系以及直接将概率密度函数（PDF）代入未封闭源项，再通过统计学方法获得源项的平均效应；提出了将来值得进一步探讨的重要问题，包括湍流脉动对颗粒凝并的影响、实验仪器和实验技术以及对颗粒一般动力学方程中未知项的模化开展有针对性的实验研究。

单气泡池沸腾传热中的重力效应数值模拟

赵建福, 张良, 杜王芳, 李会雄

2021, 39(3): 121-129. doi: 10.7638/kqdlxxb-2021.0039

摘要(198) HTML (158) PDF(28)

摘要:
本文数值模拟了不同重力条件下饱和状态的FC-72在厚度5 mm的SiO₂固壁上的单气泡池沸腾传热现象及相应的气泡动力学和传热性能。固壁底面给定均匀过热度10 K，其瞬态热响应被考虑在内。通过多个气泡周期的计算，得到了准稳态的沸腾过程。结果表明，在小热流密度下，气泡脱落直径D_b反比于重力的1/2次方，且与Fritz模型一致；气泡脱落频率f则正比于重力。传热特性存在2个明显不同的区域，即重力相关区和重力无关区，其分界位置约为0.03g₀：高于该临界重力时，热流密度与重力相关，恒定过热度下热流密度与重力呈确定的指数函数关系；小于该临界重力时，热流密度与重力无关。临界重力值对应于基于Laplace长度的无量纲加热器长度约2～3之间，与Raj-Kim-McQuillen重力标度模型建议的2.1相近。重力相关区标度指数随着过热度增大而单调增大，但明显大于Raj-Kim-McQuillen重力标度模型预测结果。研究结果还表明固壁瞬态热响应对重力无关区传热性能的影响比重力相关区更为明显。

水平T型管两相流流致振动实验研究

宿恒, 瞿叶高, 周劭翀, 周莹, 彭志科

2021, 39(3): 130-137. doi: 10.7638/kqdlxxb-2021.0033

摘要(159) HTML (80) PDF(20)

摘要:
以水平T型管为研究对象，开展了两相流输运管道流致振动实验，研究了T型管在气液两相流作用下的振动响应，重点考察了不同气液表观流速下段塞流对流体脉动压力和管道振动响应的影响规律。研究结果表明：T型管相分离现象主要与流体所受惯性力有关，径向及轴向分支管路相分配随气相表观流速增加趋向均匀；段塞流引起的脉动压力频段范围集中在0～15 Hz以内，压力幅值随气相表观流速的增加而增大；段塞流引起的结构振动以垂向振动为主，且振动幅值随气相表观流速增大而增大；T型管的振动与管内气塞和液塞的交替运动有关，两相流作为宽频激励主要激起管道的运行模态并以一阶振动模态为主；主管路两相流平均动量高于分支管路，主管路振动响应幅值最高，峰值频率最小；轴向分支管路由于相分离现象导致液相占比增大，流速减小，易出现因液柱堵塞引起的流体冲击现象，轴向分支管路振动响应幅值高于径向分支管路。

多相流系统的离散玻尔兹曼研究进展

许爱国, 陈杰, 宋家辉, 陈大伟, 陈志华

2021, 39(3): 138-169. doi: 10.7638/kqdlxxb-2021.0021

摘要(690) HTML (213) PDF(176)

摘要:
针对多相复杂流体系统模拟研究，简要介绍从格子气模型到离散玻尔兹曼方法的发展历程。从统计物理学基本原理出发，通过粗粒化建模思路，给出玻尔兹曼方程；分析Chapman-Enskog多尺度展开方法所蕴含的测量逐步细化的物理图像，给出离散玻尔兹曼建模的基本原则和主要步骤。简要介绍离散玻尔兹曼在相分离、燃烧、流体不稳定性等系统中的应用。对于多相复杂流体系统的动理学建模，技术关键是分子间作用力和化学反应贡献的引入。不同颜色的示踪粒子的引入，使得在单流体理论框架下即可实现混合过程中物质粒子来源的确定；示踪粒子在其速度相空间的分布所形成的结构蕴含丰富的流场信息，为复杂流场研究张开一个全新的视角。在多介质情形，离散玻尔兹曼建模与动理学宏观建模的对应关系是一对多。随着系统非平衡程度加深，相对于动理学宏观建模与模拟思路，离散玻尔兹曼建模与模拟的复杂度上升速度较慢。作为系统行为粗粒化描述的一种物理模型构建方法，离散玻尔兹曼根据研究需求，选取一个视角，研究系统的一组动理学性质，因而要求描述这组性质的动理学矩在模型简化过程中保值；是动理学直接建模方法的一种，为连续介质建模失效或物理功能不足、而分子动力学方法因适用尺度受限而无能为力的介尺度情形提供了一条方便、有效的研究途径。

风吹雪多相流运动及其在寒区雪水文中的应用

李广, 于鸿翔, 张洁, 黄宁

2021, 39(3): 170-181. doi: 10.7638/kqdlxxb-2021.0036

摘要(387) HTML (185) PDF(30)

摘要:
积雪作为冰冻圈的重要组成部分和地球表面最为活跃的自然要素之一，其分布及时空演化对全球气候变化具有重要的指示和反馈作用。同时，积雪也是中纬度山区主要的淡水资源，其融雪水是北半球地表径流的主要补给之一，对春季融雪径流也具有重要的影响。风吹雪是积雪在风力的作用下沿流向输运的一种常见自然现象，是地表积雪搬运的主要途径之一，是山区积雪再分布的关键驱动力。风吹雪涉及大气湍流边界层与雪颗粒相互作用，具有多相流动、多场耦合和多尺度作用的特点，是一个十分复杂的动力学过程，包含诸多科学热点和难点问题。本文从积雪质能平衡方程出发，介绍了与风吹雪相关的各个物理过程，包括降雪的优先沉降、风吹雪及其升华等对积雪质能平衡过程的影响以及其在寒区水文学中的应用等方面的研究现状及进展，并指出目前我国在雪水文学研究中存在的问题和未来的发展趋势，供读者参考。

惯性颗粒和壁湍流的相互作用：湍流调制及颗粒分布

王贵全

2021, 39(3): 182-191. doi: 10.7638/kqdlxxb-2021.0018

摘要(172) HTML (93) PDF(35)

摘要:
直接数值模拟壁面湍流与惯性点粒子方法被用来研究粒子对湍流的调制机理以及粒子在湍流场中的分布。由于在壁面湍流中，内层（约为y⁺ < 100，y⁺为黏性尺度下的壁面垂直方向坐标）与外层（约为y⁺ > 100）中存在着两种特征结构，即内层中的大尺度结构（large-scale structures, LSMs）与外层中的超大尺度结构（very-large-scale structures, VLSMs），这两种特征结构控制着壁面湍流的动力学过程。因此研究方向主要是探究粒子与LSMs以及VLSMs之间的相互作用：即通过研究粒子惯性效应对LSMs主导的湍流内层自维持过程来解释不同惯性的粒子导致湍流增强或减弱的原因；通过研究粒子惯性效应对VLSMs主导的湍流外层的调制，提出粒子的惯性效应导致VLSMs增强的直接与间接机理。同时，通过人为控制LSMs和VLSMs与粒子分别耦合来研究这两种特征结构的存在导致的粒子聚集与群聚行为。最后，推导了惯性粒子沉积速度的控制方程，用于分析控制粒子沉积速度的关键因素。

驻点流中的颗粒动力学

李青, 陈坚强, 毕林, 袁先旭

2021, 39(3): 192-200. doi: 10.7638/kqdlxxb-2021.0015

摘要(132) HTML (76) PDF(10)

摘要:
驻点流中的颗粒动力学可见于化工管道输运、航空发动机叶片等，相应基础研究的开展是为了建立可预测的两相流模型。针对化工管道中的驻点流颗粒问题，采用颗粒解析的直接数值模拟方法，对颗粒动力学与流体动力学的耦合机理进行了研究。研究发现，单个有限尺寸中性颗粒在驻点流对称轴上的动力学非常反常：远离壁面的时候，无论颗粒的有限尺寸或惯性大小，颗粒动力学呈现示踪粒子行为；靠近壁面，示踪粒子行为不再成立。开发了可以自适应颗粒惯性并且可以反应真实物理材料特性的弹性碰撞模型。研究发现，颗粒相对携带流体的惯性是驻点流颗粒动力学的关键。加速/减速流场导致的无黏环境压力是产生反常的流体力驱动的颗粒非接触反弹现象的根本原因。

高马赫数低雷诺数条件下圆球绕流曳力系数

李瑞元, 陈飞国, 葛蔚, 张永民

2021, 39(3): 201-208. doi: 10.7638/kqdlxxb-2021.0051

摘要(232) HTML (92) PDF(33)

摘要:
圆球绕流曳力系数通常表示为雷诺数的函数。而随着马赫数增大，曳力系数与马赫数的相关性逐渐增强。本文应用硬球-拟颗粒模型方法模拟研究1≤Re≤20和0.1≤Ma≤3条件下的圆球绕流过程。分析了模拟区域设置和拟颗粒-圆球作用的热适应系数对曳力系数结果的影响，得到合适的区域设置条件和热适应系数，获得高马赫数低雷诺数条件下圆球绕流曳力系数，并依此初步建立1≤Re≤20和0.1≤Ma≤3范围内马赫数对曳力系数的修正模型。

留言板

2021年第39卷第3期

留言板

2021年 第39卷 第3期

2021年第39卷第3期