EI、Scopus 收录
中文核心期刊

沟槽方向对湍流边界层流动结构影响的实验研究

崔光耀, 潘翀, 高琪, 李鹿辉, 王晋军

崔光耀, 潘翀, 高琪, 李鹿辉, 王晋军. 沟槽方向对湍流边界层流动结构影响的实验研究[J]. 力学学报, 2017, 49(6): 1201-1212. DOI: 10.6052/0459-1879-17-252
引用本文: 崔光耀, 潘翀, 高琪, 李鹿辉, 王晋军. 沟槽方向对湍流边界层流动结构影响的实验研究[J]. 力学学报, 2017, 49(6): 1201-1212. DOI: 10.6052/0459-1879-17-252
Cui Guangyao, Pan Chong, Gao Qi, Akira Rinoshika, Wang Jinjun. FLOW STRUCTURE IN THE TURBULENT BOUNDARY LAYER OVER DIRECTIONAL RIBLETS SURFACES[J]. Chinese Journal of Theoretical and Applied Mechanics, 2017, 49(6): 1201-1212. DOI: 10.6052/0459-1879-17-252
Citation: Cui Guangyao, Pan Chong, Gao Qi, Akira Rinoshika, Wang Jinjun. FLOW STRUCTURE IN THE TURBULENT BOUNDARY LAYER OVER DIRECTIONAL RIBLETS SURFACES[J]. Chinese Journal of Theoretical and Applied Mechanics, 2017, 49(6): 1201-1212. DOI: 10.6052/0459-1879-17-252
崔光耀, 潘翀, 高琪, 李鹿辉, 王晋军. 沟槽方向对湍流边界层流动结构影响的实验研究[J]. 力学学报, 2017, 49(6): 1201-1212. CSTR: 32045.14.0459-1879-17-252
引用本文: 崔光耀, 潘翀, 高琪, 李鹿辉, 王晋军. 沟槽方向对湍流边界层流动结构影响的实验研究[J]. 力学学报, 2017, 49(6): 1201-1212. CSTR: 32045.14.0459-1879-17-252
Cui Guangyao, Pan Chong, Gao Qi, Akira Rinoshika, Wang Jinjun. FLOW STRUCTURE IN THE TURBULENT BOUNDARY LAYER OVER DIRECTIONAL RIBLETS SURFACES[J]. Chinese Journal of Theoretical and Applied Mechanics, 2017, 49(6): 1201-1212. CSTR: 32045.14.0459-1879-17-252
Citation: Cui Guangyao, Pan Chong, Gao Qi, Akira Rinoshika, Wang Jinjun. FLOW STRUCTURE IN THE TURBULENT BOUNDARY LAYER OVER DIRECTIONAL RIBLETS SURFACES[J]. Chinese Journal of Theoretical and Applied Mechanics, 2017, 49(6): 1201-1212. CSTR: 32045.14.0459-1879-17-252

沟槽方向对湍流边界层流动结构影响的实验研究

基金项目: 

国家自然科学基金资助项目 11672020

国家自然科学基金资助项目 11721202

详细信息
    通讯作者:

    王晋军, 教授, 主要研究方向:湍流拟序结构、流动控制、飞行器空气动力学等.E-mail:jjwang@buaa.edu.cn

  • 中图分类号: O357.5

FLOW STRUCTURE IN THE TURBULENT BOUNDARY LAYER OVER DIRECTIONAL RIBLETS SURFACES

  • 摘要: 本文采用时间解析的二维粒子图像测速技术,对零压力梯度光滑以及汇聚和发散沟槽表面平板湍流边界层统计特性和流动结构进行了研究.结果表明在垂直于汇聚和发散沟槽表面的对称平面内,相对于光滑壁面,发散沟槽壁面使当地边界层厚度、壁面摩擦阻力、湍流脉动、雷诺应力等明显减小;而汇聚沟槽壁面对湍流边界层特性和流动结构的影响正好相反,汇聚沟槽使壁面流体有远离壁面向上运动的趋势,因而导致边界层厚度增加了约43%;同时,在汇聚沟槽表面情况下流向大尺度相干结构更容易形成,这对减阻是不利的.此外,顺向涡数量在湍流边界层的对数区均存在一个极大值,发散沟槽表面所对应的极大值位置更靠近沟槽壁面,而在汇聚沟槽表面则有远离壁面的趋势,由顺向涡诱导产生的较强的喷射和扫掠运动会在湍流边界层中产生较强的剪切作用,顺向涡数量的减少是发散沟槽壁面当地摩擦阻力降低的主要原因.
    Abstract: The statistical characteristics and flow structures for turbulent flow over smooth, convergent and divergent riblets flat plate with zero pressure gradient are investigated with two-dimensional time-resolved particle image velocimetry (TRPIV). It is shown in the wall-normal planes of the convergent and divergent riblets flat plate that, compared to the smooth flat plate, the local boundary layer thickness, wall friction velocity, turbulent fluctuation and Reynolds stress are evidently reduced over the divergent surface. Furthermore, the effect of convergent riblets flat plate on turbulent boundary layer flow is different from the divergent one, which causes the near wall fluid move away and results in an increment of about 43% for turbulent boundary layer thickness. Meanwhile, the large scale coherent structures are more likely to be formed for flow over convergent riblets surface, this is not benefit for drag reduction. Besides, the population of the prograde vortices reaches a maximum value in the log region of turbulent boundary layer, and which appears much closer to the divergent riblets surface than the convergent one. The ejection and sweep induced by the prograde vortices make a great contribution to the mean shear in turbulent boundary layer, and it is the decrease of the prograde vortices which results in the reduction of the wall friction for the diverqent riblets surface.
  • 早在20世纪70-80年代,人们就开始了对沟槽表面湍流边界层的研究,重点关注的是沟槽表面的减阻特性及其机理[1-4]. Bechert等[5]通过改变沟槽的形状并调整沟槽尺寸,以期获得最佳的减阻效果.实验研究及直接数值模拟均表明,在沟槽无量纲间距$s^+$为16$\sim$20, 高度$h/s$为0.5$\sim$0.8时沟槽表面能够达到最佳的减阻效果,其最大减阻率可达10% [4-5]. Bacher和Smith[6]认为沟槽减阻是由其诱导产生的二次涡与流向涡相互作用,使低速流体保持在沟槽内,减小了高低速流体间动量的交换,从而使阻力减小. Bechert等[5]则指出, 沟槽的存在会抑制展向动量交换的效率,进而使阻力减小.常跃峰和姜楠[7]通过比较光滑和沟槽表面边界层内展向涡量与相干结构猝发事件各脉动速度的分布特性,进一步指出沟槽表面可以抑制相干结构的猝发过程,缩短相干结构喷射和扫掠的时间尺度,从而达到减阻效果.胡海豹等[8]在对沟槽边界层的研究中得到了类似的结果,指出沟槽可以通过干扰湍流的猝发实现阻力的减小. Klumpp等[9]指出沟槽同样会对转捩的过程产生影响,这种影响既可以是二维线性的过程,又可能是三维的过程. Wang等[10]指出沟槽表面不仅使层流边界层区域增大, 而且使层流转捩为湍流的雷诺数约增大4倍, 从而使平板边界层的阻力降低.不同沟槽或表面形式对流动存在较大影响,与沟槽表面类似,Zhang等[11]指出超疏水表面可使靠近壁面的涡结构向远离壁面的方向移动,涡结构对壁面阻力的贡献减小,从而达到减阻效果.

    大多数关于沟槽表面的研究主要针对顺流向布置的沟槽,Koeltzsch等[12]在圆管流动中首次研究了汇聚和发散型沟槽的特性,指出汇聚和发散沟槽表面边界层内平均速度和速度脉动均会出现显著的改变.相对于光滑表面,沟槽表面边界层平均速度型的对数区会呈现明显的平移[13-14],其中发散表面边界层对数区相对于光滑表面向上平移,这与顺流向沟槽表面变化趋势一致,而汇聚表面对数区则呈现显著的下移,对数区的上下平移与摩擦速度$u_\tau$的改变密切相关. Nugroho等[13]的热线测速结果表明,汇聚和发散沟槽表面在其各自边界层内诱导产生的沿法向的二次流动引起了平均速度和湍流度的显著变化[15],并且这种沿法向的运动与湍流边界层中存在的大尺度相干结构产生相互作用,导致流动结构发生明显的改变,沟槽汇聚和发散的倾斜角度和沟槽的无量纲尺寸是引起这种变化的关键因素.众所周知,热线测速技术只能给出测量点处的速度信息,流场中的涡结构及大尺度相干结构并不能直接予以展示.

    因此,本实验拟采用时间解析的二维粒子图像测速技术对光滑和沟槽表面湍流边界层进行测量,对汇聚和发散沟槽引起的流动结构的变化进行深入探究.期望发现一些新的流动现象,揭示沟槽排列引起流场结构改变的机理,为流动控制和湍流减阻提出新的方法.本文第一节为实验仪器设备、实验模型及测速系统的介绍;第二节重点对测量所得速度场进行分析,包括时均速度和湍流特性、展向顺向涡分布规律、相干结构流向尺度、顺向涡对剪切应力的贡献等;第三节为论文的结论.

    本实验在北京航空航天大学的低速循环水槽中进行,实验段几何尺寸为3 000 mm $\times$ 600 mm $\times$ 700 mm (长$\times$宽$\times$高).水槽实验段的侧壁和底壁由光滑的有机玻璃构成,均可看作是水力光滑的.本实验自由来流流速为195 mm/s,来流湍流度小于1%.

    本文采用的实验模型有3种,分别为光滑表面、汇聚沟槽表面和发散沟槽表面平板.实验模型水平放置于水槽实验段的底壁上,几何尺寸均为2 400 mm $\times$ 600 mm $\times$ 20 mm (长$\times$宽$\times$厚).所有模型均覆盖实验段全部展向范围,以消除模型边缘对测量平面的影响.为简洁起见,图 1(a)仅展示了汇聚沟槽表面的示意图.汇聚沟槽平板的沟槽与来流的倾斜角$\alpha=15^\circ$,与之对应的发散沟槽表面为$\alpha=-15^\circ$.图 1(b)表明沟槽截面形状为锯齿形,其中槽间距$s=1.932$ mm,沟槽的高度$h=1.5$ mm,顶角为60$^\circ$.沟槽表面由三轴数控机床加工所得,槽底部0.2 mm的平台为60$^\circ$铣刀的退刀槽,模型表面足够光滑,表面粗糙度的影响可以忽略.本文中$x/y/z$方向分别代表流向、法向和展向,相应的各个方向的速度则由$u, v, w$表示,如图 1(a)所示.

    图  1  实验模型
    Figure  1.  Experimental model

    为加速转捩,在实验段入口处放置3 mm的绊线,视野范围中心位于绊线下游1 850 mm处,确保三种模型的测量平面位于相同的流向位置.汇聚和发散沟槽模型的测量平面位于模型汇聚和发散交界处的正上方,也即实验模型和水槽实验段展向的中心对称面.测量位置($x=1 850$ mm)处水槽侧壁边界层厚度小于30 mm,较水槽的半宽度(300 mm)小一个量级,故可忽略水槽侧壁边界层对测量平面流动的影响.

    本实验采用二维粒子图像测速系统对不同表面形状平板模型进行测量,示踪粒子为直径5$\sim$20 ${\rm{\mu }}$m的空心玻璃微珠,其密度为1.05 g/cm$^{3}$.实验采用的激光片光厚度约为1 mm,实验模型表面进行喷漆处理以消除壁面反光的影响.激光与高速CMOS相机通过同步器连接以保证采样的同步性.测量的视野范围为98 mm $\times$ 98 mm,CMOS相机的空间分辨率为2 048 $\times$ 2 048像素.实验中相机采样频率为500 Hz,相邻两幅粒子图像的主流区粒子位移超过8个像素.由于大部分粒子直径大于2个像素,因此峰值锁定现象可以忽略[16].速度场计算采用经典的MILK算法[17],查询窗口采用32 $\times$ 32像素,重叠率为75%.由于CMOS相机内存的限制,每个采样周期可得到5 456幅粒子图像.本次实验每个工况重复采样5次,全部采样时长约为1 min,为相关计算提供足够的样本,以保证计算收敛.

    对于光滑平板模型,摩擦速度可以通过Clauser方法拟合得到[18].然而对于沟槽表面模型,由于零点位置的不确定导致经典的Clauser方法失效,故如何得到沟槽表面的真实零点位置及摩擦速度是首先需要解决的问题.在沟槽表面湍流边界层速度型拟合中常采用修正的Clauser方法,该方法指出沟槽表面边界层对数区满足下式

    $$ u^+=\dfrac 1 \kappa {\rm ln} \hat y^+ +B+\Delta u^+ $$ (1)

    式中$\hat y$为真实的法向高度($\hat y=y+y_v$), 即测量点到沟槽顶部的距离($y$)与由沟槽引起的零点偏移($y_v$)之和,$\Delta u^+$为沟槽表面边界层对数区速度型与光滑表面速度型的偏移量, $\kappa$为卡门常数,$B$为对数速度型体截距. Nugroho等[13]利用该方法对汇聚和发散表面边界层摩擦速度和零点位置进行拟合,将上式(1)对$y$进行微分可得

    $$ \dfrac{\Delta u}{\Delta y}=\dfrac{u_\tau}{\kappa} \dfrac 1{y+y_v} $$ (2)

    式中$u_\tau$为沟槽表面摩擦速度.采用迭代算法计算得到摩擦速度和沟槽表面的零点位置,对于沟槽表面模型,$u_\tau$和$y_v$均未知,在迭代过程中将修正的Clauser方法与Townsend外区相似性假说相结合以最终确定零点位置及摩擦速度.

    Choi[14]在对沟槽表面的研究中指出,无论何种粗糙形式,在边界层的重叠区域内($0.002 < yu_\tau / \delta^* U_\infty < 0.15$, $\delta^*$为边界层位移厚度, $U_\infty$为自由来流速度),当无量纲的法向高度用$yu_\tau/\delta^*U_\infty$表示时,粗糙表面湍流边界层无量纲速度亏损函数$(U_\infty-u)/u_\tau$与光滑表面边界层相吻合,因此可以通过这一特性确定沟槽表面的摩擦速度和零点位置.本文中,首先将光滑表面内速度亏损函数在重叠区域内随无量纲高度的变化进行四阶多项式拟合,然后通过调整摩擦速度和零点位置初值对沟槽表面的速度亏损函数进行最小二乘拟合迭代,不断修正迭代初值最终可得到沟槽表面摩擦速度$u_\tau$和零点位置$y_v$.经比较发现由Nugroho等[13]、Choi[14]两种方法得到的结果吻合得很好(图 2(b)).

    图  2  三种模型平均速度型及速度亏损曲线对比
    Figure  2.  Comparison of mean velocity and velocity defect profiles over three surfaces

    表 1给出了光滑表面和汇聚、发散表面边界层的各个参数,$U_\infty$为自由来流流速,$\alpha$为沟槽倾斜角度,$u_\tau$为壁面摩擦速度,$\delta$为边界层厚度,$s^+$和$h^+$分别为无量纲沟槽间距和高度,$Re_\tau$和$Re_\theta$分别为摩擦雷诺数和动量雷诺数,$H$为形状因子.可见汇聚沟槽表面边界层厚度较光滑表面增大43%,而发散沟槽表面的边界层厚度则略小于光滑表面.另外,壁面摩擦速度$u_\tau$在3种表面内存在明显不同,汇聚沟槽表面$u_\tau$较之于光滑表面偏大,而发散表面则正好相反.

    表  1  3种平板模型湍流边界层参数列表
    Table  1.  Parameter lists in turbulent boundary layers over three surfaces
    下载: 导出CSV 
    | 显示表格

    图 2(a)给出了光滑表面、汇聚和发散沟槽表面边界层的平均速度分布,可见相较于光滑表面,发散沟槽表面对数区速度型存在明显的上移($\Delta u^+=3.03$),而汇聚表面则出现显著的下移($\Delta u^+=-5.67$),进一步说明前者具有减阻特性而后者使壁面阻力增加.

    图 3给出了3种表面边界层内流向、法向湍流度分布及雷诺应力的变化曲线.需要指出的是此处采用了光滑平板表面湍流边界层厚度$\delta_{\rm s}$对法向高度进行无量纲化,这一无量纲方式可以展示同一物理高度处沟槽对边界层湍流脉动特性的影响[13].图 3(a)表明除在非常靠近壁面的位置($y/\delta_{\rm s} < 0.04$)外,汇聚表面边界层内流向湍流度较之于光滑表面明显增大,这是因为边界层内流体沿法向向上的运动将近壁面湍流脉动较大的流体带离壁面,从而使得在相同高度处汇聚表面边界层内流体的湍流强度偏大,且湍流强度的峰值位置更加远离壁面,而发散表面则恰好相反,这与Nugroho等[13]采用热线技术测得的结果相吻合.图 3(b)图 3(c)则分别给出了三种沟槽表面边界层内法向湍流度和雷诺应力的变化曲线,与流向湍流度分布类似,除非常靠近壁面处之外,汇聚表面的法向湍流脉动和雷诺应力均较光滑表面偏大,而发散表面则恰好相反,这种脉动特性的差异与沟槽的排列方向密切相关.

    图  3  3种表面边界层内湍流度及雷诺应力变化曲线
    Figure  3.  Turbulent intensities and Reynolds stress over smooth, convergent and divergent riblets surfaces

    前面从流场统计特性的角度分析了汇聚和发散沟槽表面与光滑表面边界层的差异,这种差异必定是边界层内流动结构变化的体现[7].在湍流边界层中,可以利用伽利略分解提取流场中存在的涡结构,这种方法可以识别具有特定对流速度$U_{\rm c}$的涡结构,但是由于在整个边界层内各法向高度处的对流速度存在较大差异,所以通过单一的伽利略分解难以得到流场中全部的涡结构.通过计算速度梯度张量得到流场中不同位置的旋涡强度($\lambda_{\rm ci}$), 可以有效地区别旋转和剪切的作用,因而被广泛应用于对旋涡的识别[19].同时,利用当地涡量的正负确定涡的旋转方向,相应表达式如下

    $$ \varLambda_{\rm ci} =\lambda_{\rm ci} \dfrac{\omega_z}{|\omega_z|} $$ (3)

    式中$\omega_z$代表当地的涡量,由此可将顺向涡($\varLambda_{\rm ci} < 0$)和逆向涡($\varLambda_{\rm ci}>0$)清晰地予以区别和显示.

    图 4(a)给出了伽利略分解得到的瞬时速度场,对流速度$U_{\rm c}=0.85U_\infty$.对流速度选取过程如下:通过改变$ U_{\rm c}$可以得到不同对流速度下流场涡结构分布及伽利略分解的瞬时速度场,发现在图 4(a)所示瞬时,对流速度为$0.85U_\infty$时所识别的涡结构最清晰,故取对流速度$U_{\rm c}=0.85U_\infty$.相应地,这些涡结构均可由$\varLambda_{\rm ci}$场获得,不同于伽利略分解只能识别单一对流速度的涡结构的限制,不同对流速度的展向涡均可由$\varLambda_{\rm ci}$场识别,它们吻合很好,由此体现了$\varLambda_{\rm ci}$场的伽利略不变性. Natrajan等[20]将得到的展向顺向涡和逆向涡的空间特性与Zhou等[19]提出的发卡涡模型进行比较,证明了二者的对应关系.他们指出流场中大多数顺向涡和逆向涡涡对(图 4(a)所示A/F和E/G)对应于同一个发卡涡结构,该涡对是由测量平面在发卡涡颈部位置所截得.单一的顺向涡则是由测量平面恰好截于发卡涡的头部所得,因此大多数顺向涡即为发卡涡的涡头,由此可以解释顺向涡和逆向涡为何大多成对出现并且前者的数量要远多于后者.在靠近顺向涡结构的上游偏低位置处易诱导产生较强的喷射现象(Q2运动),而在其下游则会诱导产生较强的扫掠现象(Q4运动),在图 4(a)中顺向涡E的上、下游处可清晰地观察到较强的Q2和Q4运动.

    图  4  伽利略分解与$\varLambda_{\rm ci}$确定的展向涡结构分布
    Figure  4.  Distributions of the spanwise vortices identified with Galileo decomposition and $\varLambda_{\rm ci}$ field

    定量化提取流场中广泛分布的展向涡结构需要考虑两方面的因素,即涡强度的阈值和单个涡的尺寸.在本文中采用$\varLambda_{\rm ci}(x, y) \geq 1.5\varLambda_{\rm ci}^{\rm rms}$这一阈值来确定顺向涡结构的边界[18-19],考虑到展向涡结构的尺寸以及本实验的空间分辨率,假定当单个涡的流向和法向尺寸均超过5个网格节点时才将其视为一个独立的展向涡.单个顺向涡的空间范围在流向和法向均超过$20 y^+$ (其中$y^+=v/u_\tau$),同时在该空间范围内各个点均满足$\varLambda_{\rm ci}(x, y)\leq -1.5 \varLambda_{\rm ci}^{\rm rms}(y)$.

    为研究顺向涡结构在不同高度的分布规律,定义$\varPi_{\rm p}(y)$为整个流向测量范围内以$y$为中心法向四个网格间距范围内的顺向涡的数量.在统计过程中确保所有的顺向涡只被统计一次,由此可以得到顺向涡在不同高度的分布规律.鉴于顺向涡与发卡涡涡头的对应关系,这一分布规律实际上反映了发卡涡在不同法向位置的分布特征.

    光滑表面、汇聚沟槽表面、发散沟槽表面边界层内顺向涡数量$\varPi_{\rm p}$随无量纲高度$y/\delta_{\rm s}$的变化如图 5所示.可见不论实验平板的表面如何改变,顺向涡的数量$\varPi_{\rm p}$在湍流边界层的对数区均存在一个极大值.在$y^+ >100$范围,本文与Wu和Christensen[21]所得光滑平板湍流边界层结果的变化趋势相一致,其分布密度都是随高度增加呈现明显的单调递减趋势. Wu和Christensen[21]并未给出$y^+ < 100$范围内顺向涡结构的统计特性,因而无法得到该范围内顺向涡的数量存在峰值,说明本实验结果更全面.事实上,由于大多数顺向涡被认为是发卡涡的涡头,由发卡涡的形成和演化过程可知,大多数发卡涡的涡头位于较高的法向位置而并非是紧靠近壁面处,这就导致顺向涡的分布不可能随着法向高度的增加呈现单调递减.另外,只有具有足够大的空间范围(大于$20y^+$)并且满足一定旋涡强度才会被视作是一个独立的顺向涡,在非常靠近壁面处难以满足这一判据.

    图  5  顺向涡数量随法向高度的变化
    Figure  5.  Population trend of prograde vortices in boundary layer

    图 5(a)可知,3种模型表面边界层内顺向涡的数量均在对数区达到极大值,其中汇聚沟槽表面极大值位置更加远离壁面,而发散沟槽表面$\varPi_{\rm p}$极大值更靠近壁面.另外,三种表面形状对顺向涡数量的影响与阻力的改变类似,减阻时顺向涡的数量减少,发散表面边界层内顺向涡数量较光滑表面减少约9%;而阻力增加对应于顺向涡数量的增加,汇聚表面边界层内顺向涡数量比光滑表面增加约20%.对于汇聚沟槽模型,中心对称面两侧的沟槽均向中间收敛,导致其沟槽内部和紧靠壁面处的流体向对称面汇聚,从而在汇聚沟槽模型的对称面内诱导产生一个由壁面起沿法向向上的流动(common-flow-up运动);与之相反,在发散沟槽模型中沟槽均由对称面向两侧延伸,沟槽内部和靠近壁面的流体有向左右两侧运动的趋势,从而在对称面内诱导产生由较高位置流向壁面的运动(common-flow-down运动).图 6给出了瞬时脉动速度场图像,图 6(a)表明该瞬时汇聚表面边界层内流动大致沿法向向上,图 6(b)则给出了发散表面边界层内流体沿法向向下的运动.通过对所有时刻的速度场数据统计可得,汇聚表面边界层内的流动大多如图 6(a)所示沿法向向上,而在发散表面边界层内则大多如图 6(b)所示沿法向向下.通过瞬时速度场可体现汇聚、发散沟槽表面诱导产生的common-flow运动.

    图  6  瞬时脉动速度场中的common-flow运动
    Figure  6.  Common-flow motion in instantaneous fluctuation velocity field

    图 3指出这种运动使得汇聚和发散表面边界层内的湍流脉动特性发生变化,苏健等[22]指出顺向涡的形成和分布与湍流脉动密切相关,湍流脉动的增强有助于顺向涡的形成.图 4所示的瞬时速度场也表明顺向涡附近具有较强的脉动特性,可诱导形成较强的喷射和扫掠运动,由汇聚沟槽表面诱导产生的沿法向向上的运动使湍流脉动增强,促进了顺向涡的形成和发展,从而使汇聚表面顺向涡的数量得以增加,相反由于湍流脉动的减小使得在发散表面内顺向涡的数量较光滑表面减少.

    另一方面,3种表面在相同来流条件下的雷诺数$Re_\tau$存在较大的差异(见表 1),Wu和Christensen[21]指出光滑表面湍流边界层中顺向涡的分布存在一个雷诺数相关特性,无论是湍流平板边界层还是槽道湍流,$\varPi_{\rm p}(y)Re^{-3/2}_\tau$在内尺度法向高度$100 < y^+ \leq 0.5\delta^+$的范围内几乎完全重合,且此重合特性与雷诺数无关.图 7给出了本实验3种表面顺向涡数量与$Re_\tau$的相关曲线,可见3种表面$\varPi_{\rm p} (y)Re_\tau^{-3/2}$随$y^+$的变化呈现较大的差异,这种差异体现了汇聚和发散沟槽表面对边界层内涡结构分布的显著影响. Wu和Christensen[21]关注光滑表面湍流边界层在不同雷诺数下顺向涡的分布规律,其雷诺数的变化是通过改变自由来流条件实现的,而本实验中3种模型的来流条件完全相同,雷诺数的差异是由不同的沟槽表面所致.由此可见,汇聚和发散沟槽表面对边界层内涡结构数量和分布的改变要远大于$Re_\tau$的影响.

    图  7  $\varPi_{\rm p} Re_\tau^{-3/2}$随$y^+$的变化
    Figure  7.  Variation of $\varPi_{\rm p} Re_\tau^{-3/2}$ with $y^+$

    Christensen和Adrian[23]通过条件平均的方法得到在湍流边界层中存在一系列不同尺度的发卡涡包结构,并且这些结构在边界层相干结构中占据主导地位.根据2.2节的介绍可知沟槽的排列会改变顺向涡结构的分布,而顺向涡又与发卡涡涡头密切相关,由此推测沟槽排列会影响发卡涡包结构的分布特性. Marusic[24]指出脉动速度相关系数$\rho_{uu}$的流向延伸长度同该位置处相干结构沿流向的分布特征相对应,通过两点相关计算可以体现汇聚、发散沟槽表面对边界层内相干结构的影响.脉动速度两点相关系数计算如下式

    $$ \rho_{ij} (x_{\rm r}, y; y_{\rm ref})= \dfrac{\langle u'_i (x, y_{\rm ref}, t)u'_j(x+x_{\rm r}, y, t)\rangle } {\| u'_i(x, y_{\rm ref}, t) \| \| u'_j (x+x_{\rm r}, y, t) \| } $$ (4)

    式中$y_{\rm ref}$代表相关计算选取的参考高度,$x_{\rm r}$代表两点流向的空间间距,$\langle \cdot \rangle$和$\| \cdot \|$分别表示内积和向量的二范数.为了消除相邻两帧速度场数据相关性的影响,每隔十帧取样一帧进行互相关计算,因此互相关计算的实际采样频率为50 Hz,采样间隔是He等[25]进行互相关计算采样间隔的2倍,能够保证相关系数计算的可靠性,同时样本总量足够大,可保证相关系数的收敛.

    图 8 (a)$\sim$图 8(c)给出了光滑和汇聚、发散沟槽模型在流向$\!$-$\!$-$\!$法向平面内相同参考高度$y_{\rm ref}=0.2 \delta_{\rm s}$处流向速度脉动相关系数$\rho_{uu}$的分布云图.相关系数沿流向的倾斜角度同Marusic[24]和Adrian等[26]对发卡涡包倾角的描述相符. Wu和Christensen [27]对比了光滑表面和非规则粗糙表面湍流边界层内流向速度脉动两点相关系数云图,指出二者的差异很小; 而图 8指出在本实验中汇聚和发散沟槽表面的相关系数云图与光滑表面存在明显的差异,这表明沟槽与非规则粗糙表面对边界层内相关特性影响的差异.不同于汇聚和发散沟槽表面这类规则排布的表面粗糙形式,非规则粗糙表面并未在近壁面诱导形成沿法向的运动,因此边界层内大尺度相干结构并未发生显著改变,相关系数云图自然也就不存在明显差异.二者的不同从另一方面体现了表面粗糙形式对边界层内相干结构分布的重要影响,由此可知相比于非规则粗糙表面,汇聚和发散沟槽表面这种规则排布的粗糙形式对边界层内流动结构的影响更具规律性.

    图  8  $y_{\rm ref}=0.2\delta_{\rm s}$处流向速度脉动相关系数$\rho_{uu}$云图
    Figure  8.  Contour of $\rho_{uu}$ at $y_{\rm ref}=0.2\delta_{\rm s}$

    为定量比较光滑表面和汇聚、发散沟槽表面的相关系数$\rho_{uu}$, 图 9(a)给出了在参考高度$y_{\rm ref}=0.2 \delta_{\rm s}$处3种表面相关系数沿流向的变化曲线,可知在相关系数相同时汇聚沟槽表面$\rho_{uu}$的流向延伸长度较光滑表面具有明显的增加,而在发散沟槽表面内该长度则显著减小,图 9(b)给出了3种表面在$y_{\rm ref}=0.2\delta_{\rm s}$处$\rho_{uu}$沿法向的变化曲线,通过对比图 9(a)图 9(b)可知3种表面内相关系数分布沿流向的差异相比于法向差异要大得多.

    图  9  $y_{\rm ref}=0.2\delta_{\rm s} $处$\rho_{uu}$的变化
    Figure  9.  Variation of $\rho_{uu}$ at $y_{\rm ref}=0.2\delta_{\rm s}$

    为研究外区相干结构的分布规律,Christensen和Wu[28]取相关系数$\rho_{uu}=0.5$时所选参考高度上与参考点的流向距离作为流向脉动速度相关系数的延伸长度,即定义$L_x=2r_x |_{\rho_{uu}=0.5}$,并以此作为边界层内沿流向延伸的大尺度相干结构的印记.在选定的参考位置处汇聚沟槽表面内的延伸长度$L_x$较光滑表面要大20%,而在发散沟槽表面内$L_x$则要比光滑表面小25%.考虑到相关系数流向延伸长度与相干结构的相关性,可将$L_x$视作发卡涡包结构流向尺度的重要标志,汇聚和发散沟槽表面内$L_x$的差别恰好与图 5所示的顺向涡结构数量的差异相对应,即汇聚沟槽表面内$L_x$较大同时其顺向涡的数量较光滑表面增加,而发散沟槽表面内$L_x$较小,同时其顺向涡的数量也偏少.这种一致性可解释为随着顺向涡数量的增多,边界层内大尺度相干结构(发卡涡包结构)更容易形成并且其沿流向的范围更大,通过两点相关计算得到的流向尺度$L_x$自然就会增大.而发散沟槽表面正好相反,即顺向涡的数量和流向尺度$L_x$均较光滑表面减小. 3种表面顺向涡数量$\varPi_{\rm p}$与流向延伸长度$L_x$的差异体现了沟槽排布对边界层内发卡涡包等大尺度相干结构的重要影响.图 10给出了光滑和汇聚、发散沟槽表面$L_x$在边界层不同高度的变化曲线,可知在近壁区$L_x$随法向高度的增加迅速增大,约在$y/\delta_{\rm s}=0.4$时$L_x$趋于常数.这一变化规律与Christensen和Wu[28]给出的槽道湍流的结果一致,他们指出在近壁面$L_x$随高度的增加逐渐增大,当$0.2 h < y < 0.6h$ ($h$为槽道高度)时$L_x$基本保持不变.通过比较图 10的三条曲线可知,与近壁面较小的差异不同,3种表面$L_x$的差异在远离壁面处非常明显,尤其当$y/\delta_{\rm s}>0.15$时汇聚表面的$L_x$较光滑表面偏大,而发散表面则偏小.

    图  10  $L_x$随法向高度的变化
    Figure  10.  Variation of $L_x$ in boundary layer

    上述差异与Nugroho等[13]基于泰勒冻结假设得到的预乘谱图像结果一致,需要指出的是,Nugroho等[13]采用热线测速方法,从能谱分析的角度就沟槽排布对大尺度相干结构的影响进行了研究.而本实验采用时间解析的二维粒子图像测速技术,能够直接得到涡结构的变化和空间相关系数等特征量的差异,因而可为汇聚和发散沟槽表面对边界层相干结构的影响提供更加直观的证据.

    前面重点讨论了汇聚和发散沟槽表面对边界层内相干结构的影响,而流动结构的改变必然会导致剪切应力的变化,表 1中3种表面$u_\tau$的差异和图 2中对数区速度型的上下平移均可以体现这种改变.由于顺向涡可诱导产生较强的剪切应力[22],由图 5可知不同表面边界层内顺向涡的数量存在较大差异,因此本文对3种实验情况下顺向涡结构对剪切应力的贡献进行了比较.在湍流边界层中总的剪切应力表达式如下

    $$ \tau(y)=\mu \dfrac{\partial U}{\partial y}-\rho \langle u'v'\rangle $$ (5)

    式中等号右侧两项分别代表黏性切应力和湍流雷诺应力.相似地,顺向涡包含的剪切应力由$\tau_{\rm p}(y)$来表示,特定高度$y$处的顺向涡结构对总体剪切应力的贡献由下式给出

    $$ s_{\rm p}(y)=\dfrac{\tau_{\rm p}(y)}{\tau(y)} $$ (6)

    本文中顺向涡包含剪切应力的计算参照Wu和Christensen[21]提出的方法,首先根据前面所得顺向涡的分布对各个时刻流场中的全部网格节点进行二值化取值,得到与速度场时间序列数目相同的二值化矩阵序列,该二值化矩阵将顺向涡所在位置处的网格节点置为1,其余位置全部置为0,这样就用该矩阵序列标记了每一时刻顺向涡的分布位置.然后将所有时刻的速度场矩阵与该瞬时对应的二值化矩阵进行点乘,得到新的速度场.在顺向涡处的速度场与原速度场相同,在其余位置速度均为0.然后将得到的新的速度场进行与原速度场相同的计算,即可得到顺向涡所包含的剪切应力$ \tau_{\rm p}(y)$,进而可分析顺向涡所包含的剪切应力对总体剪切应力的贡献[21].

    图 11 (a)为不同雷诺数下光滑表面边界层内顺向涡对总体剪切应力的贡献,其中$Re_\tau=483$为本实验光滑表面边界层数据,$Re_\tau=1 400$, 2 350, 3 450为Wu和Christensen[21]零压力梯度光滑平板湍流边界层的实验结果,他们指出随雷诺数的增大顺向涡对剪切应力的贡献逐渐减小,对比本实验$Re_\tau$的数据可知,这一变化规律在低雷诺数下仍然成立.由图 11(a)可知,$S_{\rm p}$在非常靠近壁面处存在一个极大值,而后缓慢减小,$y/\delta$在0.2$\sim$0.6的范围内时$S_{\rm p}$基本保持不变,而后逐渐增加.顺向涡中心的剪切应力在总体剪切应力中所占比例较小.以$Re_\tau=483$为例,$y/\delta$在0.2$\sim$0.6范围内时顺向涡中心区域的剪切应力约占总体剪切应力的9%,其贡献值并不大.但是由瞬时速度场可知,顺向涡周围存在较强的剪切作用,这些喷射和扫掠运动可产生较大的剪切应力,因此我们将统计的顺向涡区域的范围扩大,即包括由其诱导产生的在其周围的运动,分析扩大后的区域对剪切应力的贡献.考虑到顺向涡及由其诱导产生的运动的空间范围以及本实验的空间分辨率,结合涡的边界阈值$|\varLambda_{\rm ci}|\geq 1.5 \varLambda^{\rm rms}_{\rm ci}$将统计范围以顺向涡边缘为界向外延伸3个网格节点范围,并且确保每个顺向涡区域仅被统计一次,所得扩大后的顺向涡区域对剪切应力的贡献由图 11(b)给出.扩大后的面积为原顺向涡中心区域面积的2倍,而其对于总体剪切应力的贡献则为其中心区域的4$\sim$5倍,由此可知尽管顺向涡中心区域对边界层内剪切应力的贡献较小,但是由其诱导产生的位于其周围的喷射和扫掠运动具有较强的剪切作用,亦即顺向涡对边界层内的剪切应力具有重要影响.

    图  11  不同雷诺数下光滑表面边界层内顺向涡对总体剪切应力贡献
    Figure  11.  Fraction of mean shear contributed by prograde vortices over smooth surface

    不同表面边界层内同一高度处顺向涡所含剪切应力与光滑表面边界层总剪切应力之比定义为

    $$ S_{\rm p}(y)=\dfrac{\tau_{\rm p}(y)}{\tau_{\rm s}(y)} $$ (7)

    式中$\tau_{\rm s}(y)$为光滑表面$y$高度处总的剪切应力,结果如图 12所示.可见汇聚表面的顺向涡所含的剪切应力在整个边界层内都大于光滑表面,尤其是在远离壁面处($y/\delta_{\rm s} >0.6$)这种差异尤为明显,这与汇聚表面在该处顺向涡数目的增多和雷诺应力的显著增加相对应.对于发散表面,由于顺向涡数目较光滑表面偏小,且其雷诺应力也较小,故对剪切应力的贡献减小.在远离壁面处由于发散表面和光滑表面的雷诺应力及顺向涡数目的差异较小,故顺向涡所含的剪切应力与光滑表面无明显差异.汇聚沟槽表面阻力的增加与顺向涡数目的增多密切相关,而在发散表面顺向涡数目较光滑表面偏小,故剪切应力及阻力也随之减小.

    图  12  顺向涡对剪切应力的贡献
    Figure  12.  Mean shear contained within prograde vortices

    顺向涡属于大尺度含能结构,对应于湍动能的产生项.由汇聚、发散沟槽表面诱导产生的沿法向的流动会使顺向涡的分布发生改变,顺向涡及其诱导产生的较强的剪切作用会影响湍动能的产生,湍动能的耗散随之改变,进而可使阻力发生改变.因此通过改变表面形状或采用其他控制方法,改变湍流边界层中顺向涡的分布并减少其数量,可达到减阻的目的.

    本实验采用时间解析的二维粒子图像测速技术对光滑表面及汇聚、发散沟槽表面湍流边界层中统计特性及涡结构进行了研究,结果表明汇聚和发散沟槽表面会对湍流边界层统计特性和顺向涡等结构产生显著影响.

    (1) 与光滑平板湍流边界层统计特性相比,发散沟槽表面和汇聚沟槽表面对湍流边界层的影响正好相反.发散沟槽表面使对数区速度分布曲线明显上移、边界层厚度及摩擦速度减小、边界层内流向和法向湍流度及雷诺应力降低.可见,沟槽表面影响整个边界层流动,而发散沟槽表面具有较好的减阻能力和应用前景.

    (2) 不论实验平板的表面如何改变,顺向涡数量$\varPi_{\rm p}$在湍流边界层的对数区均存在一个极大值,且相对于光滑表面,发散沟槽表面所对应的极大值位置更靠近沟槽壁面,而在汇聚沟槽表面则有远离壁面的趋势.此外,在汇聚沟槽表面情况下流向大尺度相干结构更容易形成,且流向延伸其范围更大.

    (3) 尽管顺向涡中心部分对剪切应力的贡献较小,但由其诱导产生的喷射和扫掠运动会在边界层内产生很强的剪切作用.顺向涡的增加将导致阻力的增大.因此,通过改变表面形状,减少湍流边界层中顺向涡的数目,可以达到减阻的效果.

  • 图  1   实验模型

    Figure  1.   Experimental model

    图  2   三种模型平均速度型及速度亏损曲线对比

    Figure  2.   Comparison of mean velocity and velocity defect profiles over three surfaces

    图  3   3种表面边界层内湍流度及雷诺应力变化曲线

    Figure  3.   Turbulent intensities and Reynolds stress over smooth, convergent and divergent riblets surfaces

    图  4   伽利略分解与$\varLambda_{\rm ci}$确定的展向涡结构分布

    Figure  4.   Distributions of the spanwise vortices identified with Galileo decomposition and $\varLambda_{\rm ci}$ field

    图  5   顺向涡数量随法向高度的变化

    Figure  5.   Population trend of prograde vortices in boundary layer

    图  6   瞬时脉动速度场中的common-flow运动

    Figure  6.   Common-flow motion in instantaneous fluctuation velocity field

    图  7   $\varPi_{\rm p} Re_\tau^{-3/2}$随$y^+$的变化

    Figure  7.   Variation of $\varPi_{\rm p} Re_\tau^{-3/2}$ with $y^+$

    图  8   $y_{\rm ref}=0.2\delta_{\rm s}$处流向速度脉动相关系数$\rho_{uu}$云图

    Figure  8.   Contour of $\rho_{uu}$ at $y_{\rm ref}=0.2\delta_{\rm s}$

    图  9   $y_{\rm ref}=0.2\delta_{\rm s} $处$\rho_{uu}$的变化

    Figure  9.   Variation of $\rho_{uu}$ at $y_{\rm ref}=0.2\delta_{\rm s}$

    图  10   $L_x$随法向高度的变化

    Figure  10.   Variation of $L_x$ in boundary layer

    图  11   不同雷诺数下光滑表面边界层内顺向涡对总体剪切应力贡献

    Figure  11.   Fraction of mean shear contributed by prograde vortices over smooth surface

    图  12   顺向涡对剪切应力的贡献

    Figure  12.   Mean shear contained within prograde vortices

    表  1   3种平板模型湍流边界层参数列表

    Table  1   Parameter lists in turbulent boundary layers over three surfaces

    下载: 导出CSV
  • [1]

    Walsh MJ. Riblets as a viscous drag reduction technique. AIAA Paper, 1983, 21(4):485-486 doi: 10.2514/3.60126

    [2] 王晋军.沟槽面湍流减阻研究综述.北京航空航天大学学报, 1998(1):31-34 http://d.wanfangdata.com.cn/Periodical/bjhkhtdxxb199801009

    Wang Jinjun. Review and prospects in turbulent drag reduction over riblets surface. Journal of Beijing University of Aeronautics and Astronautics, 1998(1):31-34 (in Chinese) http://d.wanfangdata.com.cn/Periodical/bjhkhtdxxb199801009

    [3] 王晋军, 兰世隆, 苗福友.沟槽面湍流边界层减阻特性研究.中国造船, 2001, 42(4):1-5 http://d.wanfangdata.com.cn/Periodical/zgzc200104001

    Wang Jinjun, Lan Shilong, Miao Fuyou. Drag-reduction characteristics of turbulent boundary layer flow over riblets surfaces. Shipbuilding of China, 2001, 42(4):1-5 (in Chinese) http://d.wanfangdata.com.cn/Periodical/zgzc200104001

    [4]

    Martin S, Bhushan B. Fluid flow analysis of a shark-inspired microstructure. Journal of Fluid Mechanics, 2014, 756:5-29 doi: 10.1017/jfm.2014.447

    [5]

    Bechert DW, Bruse M, Hage W. Experiments on drag-reducing surfaces and their optimization with an adjustable geometry. Journal of Fluid Mechanics, 1997, 338:59-87 doi: 10.1017/S0022112096004673

    [6]

    Bacher E, Smith C. A combined visualization-anemometry study of the turbulent drag reducing mechanisms of triangular micro-groove surface modifications. AIAA Paper, 1985:85-0546

    [7] 常跃峰, 姜楠.沟槽壁湍流多尺度相干结构实验研究.航空动力学报, 2008, 23(5):788-795 http://d.wanfangdata.com.cn/Periodical/hkdlxb200805003

    Chang Yuefeng, Jiang Nan. Experimental study on coherent structure passive control and drag reduction in turbulent boundary layer with grooved surface. Journal of Aerospace Power, 2008, 23(5):788-795 (in Chinese) http://d.wanfangdata.com.cn/Periodical/hkdlxb200805003

    [8] 胡海豹, 宋保维, 黄桥高等.水下湍流减阻途径分析.摩擦学学报, 2010, 30(6):620-629 http://d.wanfangdata.com.cn/Periodical/mcxxb201006018

    Hu Haibao, Song Baowei, Huang Qiaogao, et al. Analysis about the approaches of underwater turbulence drag reduction. Tribology, 2010, 30(6):620-629 (in Chinese) http://d.wanfangdata.com.cn/Periodical/mcxxb201006018

    [9]

    Klumpp S, Meinke M, Schröder W. Numerical simulation of riblet controlled spatial transition in a zero-pressure-gradient boundary layer. Flow, Turbulence and Combustion, 2010, 85(1):57-71 doi: 10.1007/s10494-010-9251-x

    [10]

    Wang JJ, Lan SL, Lian QX. Effect of the riblets surface on the boundary layer development. Chinese Journal of Aeronautics, 1996, 9(4):257-260

    [11]

    Zhang J, Tian H, Yao Z, et al. Mechanisms of drag reduction of superhydrophobic surfaces in a turbulent boundary layer flow. Experiments in Fluids, 2015, 56(9):1-13

    [12]

    Koeltzsch K, Dinkelacker A, Grundmann R. Flow over convergent and divergent wall riblets. Experiments in Fluids, 2002, 33(2):346-350 doi: 10.1007/s00348-002-0446-3

    [13]

    Nugroho B, Hutchins N, Monty, JP. Large-scale spanwise periodicity in a turbulent boundary layer induced by highly ordered and directional surface roughness. International Journal of Heat & Fluid Flow, 2013, 41(41):90-102

    [14]

    Choi KS. Near-wall structure of a turbulent boundary layer with riblets. Journal of Fluid Mechanics, 1989, 208(208):417-458

    [15]

    Mehta RD, Bradshaw P. Longitudinal vortices imbedded in turbulent boundary layers part 2. Vortex pair with 'common flow' upwards. Journal of Fluid Mechanics, 2006, 188:529-546

    [16]

    Christensen KT. The influence of peak-locking errors on turbulence statistics computed from PIV ensembles. Experiments in Fluids, 2004, 36(3):484-497 doi: 10.1007/s00348-003-0754-2

    [17]

    Champagnat F, Plyer A, Besnerais G. Fast and accurate PIV computation using highly parallel iterative correlation maximization. Experiments in Fluids, 2011, 50 (4):1169-1182 doi: 10.1007/s00348-011-1054-x

    [18]

    Clauser FH. Turbulent boundary layers in adverse pressure gradients. Journal of Aeronautic Science, 1954, 21:91-108 doi: 10.2514/8.2938

    [19]

    Zhou J, Adrian RJ. Balachandar S, et al. Mechanisms for generating coherent packets of hairpin vortices in channel flow. Journal of Fluid Mechanics, 1999, 387:353-396 doi: 10.1017/S002211209900467X

    [20]

    Natrajan VK, Wu Y, Christensen KT. Spatial signatures of retrograde spanwise vortices in wall turbulence. Journal of Fluid Mechanics, 2007, 574:155-167 doi: 10.1017/S0022112006003788

    [21]

    Wu Y, Christensen KT. Population trends of spanwise vortices in wall turbulence. Journal of Fluid Mechanics, 2006, 568:55-76 doi: 10.1017/S002211200600259X

    [22] 苏健, 田海平, 姜楠.逆向涡对超疏水壁面减阻影响的TRPIV实验研究.力学学报, 2016, 48(5):1033-1039 http://lxxb.cstam.org.cn/CN/abstract/abstract146001.shtml

    Su Jian, Tian Haiping, Jiang Nan. TRPIV experimental investigation of the effect of retrograde vortex on drag-reduction mechanism over superhydrophobic surfaces. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(5):1033-1039 (in Chinese) http://lxxb.cstam.org.cn/CN/abstract/abstract146001.shtml

    [23]

    Christensen KT, Adrian RJ. Statistical evidence of hairpin vortex packets in wall turbulence. Journal of Fluid Mechanics, 2001, 431:433-443 doi: 10.1017/S0022112001003512

    [24]

    Marusic I. On the role of large-scale structures in wall turbulence. Physics of Fluids, 2001, 13:735-743 doi: 10.1063/1.1343480

    [25]

    He GS, Pan C, Feng LH, et al. Evolution of lagrangian coherent structures in a cylinder-wake disturbed flat plate boundary layer. Journal of Fluid Mechanics, 2016, 792:274-306 doi: 10.1017/jfm.2016.81

    [26]

    Adrian RJ, Meinhart CD, Tomkins CD. Vortex organization in the outer region of the turbulent boundary layer. Journal of Fluid Mechanics, 2000, 422(13):1-54

    [27]

    Wu YH, Christensen KT. Spatial structure of a turbulent boundary layer with irregular surface roughness. Journal of Fluid Mechanics, 2010, 655(1):380-418

    [28]

    Christensen KT, Wu YH. Characteristics of Vortex Organization in the Outer Layer of Wall Turbulence. Begel House Inc, 2007

图(12)  /  表(1)
计量
  • 文章访问数:  1852
  • HTML全文浏览量:  257
  • PDF下载量:  474
  • 被引次数: 0
出版历程
  • 收稿日期:  2017-07-10
  • 网络出版日期:  2017-10-18
  • 刊出日期:  2017-11-17

目录

/

返回文章
返回