模拟风沙运动的离散颗粒动力学模型(孙其诚,王光谦)

热度218票 浏览110次 【共0条评论】【我要评论 时间:2001年7月01日 09:50

模拟风沙运动的离散颗粒动力学模型水利论文*r{K$i:}\*\T

孙其诚王光谦水利论文U7f$PRy0F_?
(清华大学 水利水电工程系)

yFE,e/zK+` N`0

摘要:根据风沙运动的特点,建立了基于单沙粒动力学的离散颗粒动力学模型。依据该模型,模拟了9 000个沙粒在风力作用下由静止不动到充分发展的稳定态的全过程,细致研究了沙床处于稳定态时床面形态、跃移轨迹以及跃移质/蠕移质间交换规律,所得结果有助于理解风沙运动的内部机理,这是目前常用的模型所做不到的。研究表明:离散颗粒动力学模型能够真实的模拟风沙运动现象,研究范围涉及风沙运动规律和床面形态的发生与发展等关键的问题,因而比其它模型更具优势。水利论文{ qP.LX%j

关键词:风沙运动;风沙流结构;模型;计算机模拟

{N2M:ku^e0

基金项目国家杰出青年科学基金(59525914);中国博士后科学基金资助项目。水利论文1ne5l v*[U
作者简介:孙其诚(1970-),男,清华大学水利系博士后,主要从事气固两相流的实验研究和数值计算。
1w$N;l%kw#q;J0收稿日期:2000-09-06

xr0YFR0

  近几年春季,我国华北地区沙尘天气频频发生, 并涉及部分南方省市, 对民航交通及日常生活带来严重影响, 直接经济损失数以亿计, 因此加强风沙运动规律的研究已经迫在眉睫。国内外的众多学者对此进行了大量的室内外实验研究和理论上的探索[1,2],但是由于风沙运动过程的内在复杂性和影响因素的多样性,至今仍有许多方面没有搞清楚,比如沙粒的起动机制等。我国则更是存在应用技术多、基础研究少、过程研究少的突出问题。水利论文 a`W8hC F3W7Lo

  20世纪80年代以来,越来越多的研究者开始对风沙运动过程进行数值模拟,提出的模型形式各异。研究发现,几乎所有的模型首先假设风沙运动处于稳定态,所涉及到的空气含沙量、跃移沙粒的速度分布函数和轨迹、跃移沙粒与蠕移沙粒的能量、动量交换、蠕移沙粒的溅射系数等反映风沙运动内在机理的及其关键的因素无不是通过随机采样产生[3,4]。这些模型只注重最终的统计结果,忽略了中间过程中每一时刻沙粒的位置、速度、能量耗散等非常重要的信息。事实上,风沙运动过程受外界因素的影响总是处于非稳定态,那么模型中采用的参数就是一种不精确的平均,模拟结果相对粗糙。目前,风沙运动的研究方法逐渐倾向于微观分析,并与宏观综合相结合;研究深度也不仅限于对现象的描述,而是更着重于揭示内在原因及其演化历程,并力图尽可能准确的预测其未来的发展。因而,现有的模型对于全面反映风沙运动的非均匀特征和动态变化规律具有一定的局限性,不适合风沙运动过程基础研究发展的需要。

-d L!W` qY*y0

  随着CPU速度的提高和高效率计算方法的开发,从20世纪90年代中期,从单颗粒的受力分析入手,并跟踪颗粒轨迹的离散颗粒模型[5,6](discrete element model; granular dynamics model)首先在气固流态化、颗粒技术等领域发展起来,应用范围越来越广,为此国际著名期刊《Powder Technology》出版了介绍离散颗粒模型最新研究进展的专辑(2000年第4期)。在风沙物理学领域,在对沙粒的基本规律有了深刻认识的基础上,我们就可以采用越来越精细的离散颗粒模型,越来越准确的模拟风沙运动的全过程。本文建立了适用于风沙运动本身特点的离散颗粒动力学模型,模拟了9 000个沙粒在风力作用下从静止到充分发展的全过程,对推移质的含沙量分布和输沙量、跃移质与蠕移质的交换以及床面形态的发生与发展等关键问题做了详尽的研究,所得结果与已知观测结果符合较好。水利论文 s&`,w V$a u E9a!X

1 离散颗粒动力学模型水利论文4czFf v#G |;b k}

  沙粒在运动过程中不断地与其它沙粒发生碰撞,由于碰撞时动量的交换瞬时完成,其冲力远大于重力、气流曳力和沙粒间的摩擦力,因此在这一瞬间仅考虑碰撞引起的速度的变化;在两次碰撞之间,每个沙粒在曳力和重力的作用下运动。根据沙粒在上述不同过程的受力特点,其运动分解为受冲力支配的瞬时碰撞过程和受曳力和重力控制的非瞬时漂移过程,从而建立了对沙粒-沙粒以及沙粒-气流两种作用分别处理的离散颗粒动力学模型(细节见文献7,8)。

p2vn5~nw1H0

1.1 沙粒区域分布与近地层流场水利论文pESS | J-Q4J;Lt

  任意堆积的沙粒的空间位置都可以分为图1所示的3个区域:下部的沙粒密实排布,处于两虚线之间的沙粒则排布疏松(厚度为一个沙粒粒径),而上方的沙粒漂浮在空中。空间位置的不同对应着本质不同的受力情况:上部的沙粒主要受到重力和气流曳力的作用,称之为近地气流层;中部的沙粒则涉及到起动问题,需考虑沙粒间的摩擦力、重力和曳力,称之为地表层;下部的沙粒由于地表层的掩盖,不受气流的曳力,称之为屏蔽层。水利论文F/^CkMy2P X

  野外观测和风洞实验结果表明,在发生风沙运动以前,近地层风速遵循对数分布[1,9]水利论文2W d ]K7z/v)|;Nt

|Image748.gif (855 bytes)|=u0lny/y0; y0=D/30

J9J u3t z6{ K m0

(1)

f7z#D(I d+{ y#Z0

  其中y, y0分别为床面高程和粗糙度,D为沙粒粒径,均以mm计。在本模型中,按照式(1)给定近地层流场,并假定流线与床面平行,如图1实线所示。水利论文$D/G,QppV!Dp

1.2 沙粒的起动水利论文*{,a"@ [Xv(`

地表层中沙粒的受力情况及作用力表达式已经清楚,但是不同的研究者采用的起动条件不尽相同[10]。本文采用的起动条件为

:D(XpsP-xMXzLmz0

FD=μ·mg水利论文/a;W]3Gp8wIY @6P

(2)水利论文 yb]A!Wp/c

  其中FD为沙粒受到的曳力:FD=1/8πD2CDρg|Image748.gif (855 bytes)|,μ是沙粒的摩擦系数。

!~V\2jQx#l9z0

cM_*|2d|S+BU:O@0

图1 沙粒分层与风速分布示意图水利论文1[/up)V Tlve
Fig.1 Layer configuration and assumption水利论文J,P U:m cd p/J7n
of local wind velocity

:M^$`'z XHU0

2k%`&UN/f0

图2 沙粒碰撞示意图
QA1B.CB y }q(h8U0Fig.2 Collision diagram of sand particles

3s^'h%U*n-n0

1.3 沙粒运动分解

3Te-m T1^g0

1.3.1 碰撞过程水利论文KG n0JWC9zZ

  假设沙粒是球形的准刚体,不考虑沙粒的变形;碰撞为两体碰撞,且发生在质心平面上。颗粒的碰撞首先假设为两体碰撞,碰撞点仅为两颗粒的接触点。如图2所示, 在发生碰撞的两沙粒质心平面建立碰撞坐标系, 原点选在碰撞点上, y轴和x轴分别沿法向和切向。碰撞过程中产生瞬时法向冲量和切向冲量, 两沙粒动量和角动量守恒。设碰撞结束后平均法向冲量和切向冲量分别为Py, Px, 则速度和角速度改变可表示为下列方程

;Ww[5Cz.x0

沙粒水利论文lD+~O!M%f/O_S(P W-_

q0j\KXG0

  (3)水利论文W T a+GC!}"un

沙粒水利论文2q)q9U6~9i\$K

(4)

其中 ma,mb是沙粒a,b的质量;(υx,a0y,a0a0,(υx,b0y,b0b0)分别为沙粒a,b碰前的速度和角速度;(υx,ay,aa),(υx,by,bb)分别为沙粒a,b碰后的速度和角速度;Ia,Ib分别为沙粒a,b的转动惯量,Ia=1/2maR2a,Ib=1/2mbR2b。 水利论文3i5F~P Ghh&a9I r

  Wang和Mason[11]把颗粒的碰撞划分为5种模式,得到了相应的Py和Px的分析解。本文采用的二维圆片形颗粒仅涉及其中的两种: 滑动碰撞(sliding collision,物理意义为接触点的切向相对速度始终不为零), 粘滞碰撞(sticking collision, 物理意义为摩擦力相对很大, 以至于刚刚碰撞, 接触点的切向相对速度就变为零)。 具体结果见表1, e和μ分别为颗粒的恢复系数和摩擦系数。水利论文.n!iwu~p(a*JY

表1 法向和切向冲量平均值的分析解水利论文p g3Zut"Nz
Table 1 Analytical expression of impulse momentum水利论文k$^0J.b*DepN


碰撞模式条件解析解备注

Sliding collisionPd>(1+e)PqPx=-sμPy
9T6_ u[ J R3j.F'hGr#E0Py=-(1+e)C0/B2
S0=[υx,a0a0Ra]-[υx,b0b0Rb]
'Ec/Yj+wE[0C0y,a0y,b0
7| E+[iAx&H(Bm I0B1=1/ma+1/mb+R2a/Ia+R2b/Ib
3K/[%J&d+RA0B2=1/ma+1/mb
J!_qB5x"FAm*I0Pd=sB2S0Pq=-sμB1C0
(?%G0\Lht0Image751.gif (1342 bytes)

Sticking collisionPd≤(1+e)PqPx=-S0/B1水利论文$Nc2_ RZ2@erz7g
Pv=-(1+e)C0/B2

1.3.2 漂移过程水利论文2Aq5m9deE

在两次碰撞之间沙粒在曳力和重力作用下运动,方程为

3MF] `eH2M0E h0

dImage752.gif (851 bytes)/dt=Image753.gif (860 bytes)+Image754.gif (861 bytes)D/m水利论文~9Bg0@R

(5)水利论文,[iK,Kcn

其中

K:e;l(Vh j&Ja0

Image754.gif (861 bytes)D=1/8πD2CDρg|Image748.gif (855 bytes)-Image752.gif (851 bytes)|(Image748.gif (855 bytes)-Image752.gif (851 bytes))水利论文_1B,@#~`

(6)水利论文 ['a0X(f6Uq

Image755.gif (1607 bytes)水利论文Q:u |@O"I^^[

(7)水利论文"O,at8x5T F

Re=ερg|Image748.gif (855 bytes)-Image752.gif (851 bytes)|D/μg

J%a h:G4B0
(8)

1.4 气体运动模型

%@3e%Tw?b0

本文中气体的运动规律采用气固相耦合的Navier Stokes方程描述气相连续方程为水利论文&K7s6nUm8C

Image756.gif (1211 bytes)

a;r Kz'Kwe0

(9)

)[,V!gkby3y2m0

气相动量方程为

V Bd%C `,Z;r0

Image756.gif (1211 bytes)水利论文u0?)I%L2}3[1@ j W+Mt

(10)水利论文psXj9l fi {Y8J

其中气相应力张量水利论文e ^DM ^K

τ=[(-λg-2/3μg)(Image552.gif (855 bytes)·Image748.gif (855 bytes))Image758.gif (852 bytes)+μg(Image552.gif (855 bytes)Image748.gif (855 bytes))+(Image552.gif (855 bytes)Image748.gif (855 bytes))T]

7gl+O$WN0

(11)水利论文9b lBa-y!@

这里μg为气体剪切动力粘性系数,由于流化床中气体的马赫数小于1,所以设为不可压缩的气体,气体密度ρg为常数,并设气体的λg=0。相间动量交换系数β根据Ergun方程及Wen和Yu(1966)提出的公式确定,即水利论文 F ZY8U_!l

ε<0.8时,β=150(1-ε)2/εμg/d2p+1.75(1-ε)ρg/dp|Image748.gif (855 bytes)-Image752.gif (851 bytes)|

f:[5T&x6~$^0
(12)

ε≥0.8时,β=3/4CDε(1-ε)/dpρg|Image748.gif (855 bytes)-Image752.gif (851 bytes)-2.65水利论文#pc"[S7Tc

(13)

 在等温条件下,当颗粒的速度Image752.gif (851 bytes)已知时,方程(9), (10)的未知量为气体压强P,气相在x,y方向的速度分量ux, uy,空隙率ε及固相在x, y方向的速度分量vx, vy将根据离散颗粒的速度运动来计算。方程(9),(10)的标量形式如下水利论文a4l9pU4GM;C8[Q

Image759.gif (1425 bytes)水利论文'H qq-e1O?$M^t

(14)

Image760.gif (1698 bytes)水利论文L]qsK'gM

(15)

Image761.gif (3235 bytes)

C0TCR!L9a[0
(16)

在上式中固相速度由所在控制单元体内的颗粒平均速度给定

;TCik*V+mN"{ Y(OaXZ0

Image762.gif (1106 bytes)水利论文'ud:H+Y-i~DdO

(17)

其中N为所在控制单元体内的颗粒数。

5e}@FE0

  方程(14), (15), (16)的边界条件如下:水利论文!X#u6K ][1F

  (a) 速度边界条件: 入口处按照式(1)给定气速; 出口处切向速度分量和法向速度分量在法向上的梯度为零;

2nf^p;o4a|0

  (b) 压强条件: 入口处和出口处均为标准大气压。水利论文c+h7N8jTRVI

1.5 计算步骤水利论文[&QiK"n

  (1) 在固定时间步长内, 数值求解两相耦合的N-S方程,确定气体的流动;

$N Qwf2Hkg4_0

  (2) 计算床面形态,并依据当地风速判断地表层中的沙粒是否起动;水利论文QF"A rbY!f

  (3) 判断沙粒是否发生碰撞,如发生,则根据式(3)式(4)确定碰后的速度;计算地表层及在气流中沙粒受到的曳力和重力,根据式(5)计算时间步长Δt后的速度和位移;

)V*n2b)H4P2^8@0

  (4) 根据变化的颗粒速度,修正N S方程中的相间作用系数;

zz*pc(y0

  (5) 重复上述步骤, 就可以动态模拟沙粒运动的全过程。

{ h1kx4uU0

2 风沙运动现象的模拟结果与讨论水利论文$mh%S1yi,QLW

  为了检验离散颗粒动力学模型的正确性,我们对风沙运动现象做了模拟,各类参数见表2;气流从左边界吹入;左右边界采用周期性边界条件,来近似自然界中风沙运动的实际情况;沙粒的初始排布见图3(a)。

A9q|0Z'Qf5[0

2.1 沙粒运动趋稳过程水利论文?jbm7o7T9Rk m

  图3是模拟得到的风沙运动现象。可以看到风力逐渐吹起并带走部分沙粒,这些沙粒受重力作用而逐渐下落,并在空中和床面与其它沙粒发生碰撞,引发其它沙粒发生跃移或蠕移。这样众多沙粒相互作用、共同运动,一方面不断受到风力作用而加速,另一方面又不断消耗能量,最终形成充分发展的稳定态。水利论文0dh1U R{"N%]:}

表2 沙粒、气体物性及计算参数水利论文/lSQg6Q |
Table 2 Parameters used in the simulation水利论文MU2lK%ir/r,F-t


沙  粒水利论文3T8ish${JfI.]9ye\

气流水利论文3[%a~Ur

计算参数水利论文!F0H%Rd%Ba6tB


直径D=1mm
q/Y%N6?I$b;_Q0密度ρg=1.29kgm-3水利论文Ct8x9a6WR;[V
恢复系数e=0.80水利论文.~+?'o*?F X|
摩擦系数μ=0.12水利论文1W^o3J;i9|_&B
数目N=9000 (300×30)
6Bv| [9_2^%F,j |0初始堆积300mm×30mm(见图2a)
密度ρp=2 650kgm-3水利论文A]+l pH L!j%_
动力粘度μg=1.8×105Nsm-2水利论文*?$fVY t-d
u0=0.664ms-1
*H ao1Eo&|0(对应2m处风速为7.1ms-1)
时间步长 Δt=10-6s水利论文P*H8r _!Bw A/vB
计算区域 300mm×100mm水利论文d2E _A I:\o5UX_8Q
流场网格节点 150×50水利论文7S9N H;B$VU
流场网格步长 2mm×2mm
p|D JRY0沙粒网格步长
1mm×1mm

水利论文 E'u1o(R P9H_

图3 风沙运动现象的模拟
h2D2z!C`*d5v0Fig.3 Snapshots of simulated sand movement

G7WLmbI!y0

W3f K+S&c%l;rAf0

图4 沙粒运动趋稳过程
V&^CT0b'i$w9N%}$st0Fig.4 Averaged velocity per particle水利论文D#@Hl |QLa AQ4g*F
against time
水利论文(MAa(f?)D^2AbLqI

  图4是所有沙粒的均方根速度()随时间的变化关系,它量化了风沙系统中风力加速和能量耗散两个因素相互对抗、相互协调的程度。可以看出,整个过程可以划分为两个阶段:加速段和稳定段。在前12秒的加速段内,沙粒主要是在风力作用下加速运动,但是随着时间的推移,风力加速逐渐减小,而沙粒速度不断增大,耗散的能量也因频繁的非弹性碰撞而增大。在12秒左右这两个因素相互协调,沙粒得到充分发展,达到稳定段。图3 (b)是加速段时沙粒的运动情况,图3 (c)是稳定后的运动情况。
  注意到即使风沙运动达到稳定后,沙粒的均方根速度也只有7.05mm/s,这主要是因为绝大部分沙粒处于床面和屏蔽层,速度几近于零。同时,模拟中达到稳定的时间比在沙风洞实验中所需要的时间少的多,这主要因为实验中外界因素难以很好的控制,比如风洞壁面的影响、来流稳定性和速度均匀性的影响等,而在计算机模拟中不存在这些问题,所以较快达到稳定。

2.2 床面形态描述水利论文 gCg3e\-Zo HC-I

水利论文9r/\N"ZQDD I

图5 模拟得到的床面形态水利论文.^1JC#`U'O] V4?
Fig.5 Simulated ripples at dynamic equilibrium

AN3R q!^8sfa0

  由于本模型中沙粒为离散的,计算量巨大,目前只能模拟几千个沙粒,所以如何利用有限的沙粒来描述床面形态就显的非常关键。本文中流场网格步长2mm×2mm,当沙粒密实排布时可以包含4个沙粒,我们定义包含4个沙粒的网格构成整个沙床,最顶部的包含4个沙粒的网格构成了床面顶端形态。由此得到的数据往往表现为离散的点,直接依序连结的折线来反映床面形态常常过于粗糙,因而必须采用光滑方法对原始数据进行修匀。本文采用Savitzky-Golay 滤波方法来光滑床面数据。图5是图3(c) 对应的沙床形态,细线为原始数据,粗线为光滑后的最终床面形态,平均床高为12mm。可以看出:床面在风力作用下由平整演变为高低起伏的波状,具有一定的周期性,波长与波高之比近似为6,这一结果处于合理的观测结果范围内[9]水利论文 l6k)mNc;}I

2.3 跃移质运动轨迹

3a HS%|'i*^I8y0

  图6是t=30.36s~30.52s间隔内部分沙粒的跃移运动轨迹。可以看出跃移运动具有特殊的抛物线形状,一些沙粒在空中与其它沙粒碰撞后弹起,因而轨迹在发生转折。还可以推算出图中所示跃移沙粒的速度约为1.5m/s,约为计算中采用的风速的1/5,这是由于气流流经沙面时,相当一部分能量消耗于地表的摩擦和沙粒间的非弹性碰撞之中,从而使得沙粒速度大大减小。但是,它比在沙风洞中的实验数据稍大(实验中沙粒速度一般比风速小一个数量级),这是因为该模型暂时没有考虑沙粒的旋转,它也消耗很大比例的气流动量,使得沙粒速度进一步减小。需要说明的是,为了便于演示,图6仅选用了几个容易识别的沙粒,并不是所有的跃移沙粒都具有图中所示轨迹,模拟中发现更多的跃移沙粒的轨迹极其复杂,其速度也要比图中所选的几个沙粒的速度小。这个结果是其它模型所无法得到的。水利论文|| geP*}g9k R/j

2.4 跃移质与蠕移质的交换

o4w z-?Zj0

  跃移沙子受重力作用而下落,虽然不断与其它沙粒发生碰撞,之后又被弹起,但是仍然有不少沙粒在能量消耗殆尽时,落到床面成为蠕移沙粒。与此同时,床面的部分蠕移沙粒受到较大的冲击而迅速跃起,成为跃移沙粒。两者相互矛盾,构成了风沙运动内部的基本过程。当跃起量大于下落量时,更多的沙粒被风带走,床面逐渐被吹蚀;当下落量大于跃起量时,沙粒逐渐堆积到床面;当两者总量达到平衡时,形成稳定的床面形态。所以,跃移质/蠕移质交换规律的研究对于了解风沙运动内在机理很重要。但是直接的实验观测极其困难,在离散颗粒动力学模型中却能够很方便地“测量”。图7是23.48~23.60 s内床面不同位置跳离床面成为跃移质的沙粒数目(a线),落到床面成为蠕移质的沙粒数目(b线,定义为负值),两者几乎相等,平均值为6(-6),这样在本文的计算条件下两者交换大约50个/秒。

TMcd)UP"bI)F0

/`'^#ppeh n/@cE$]0

图6 部分沙粒跃移轨迹水利论文 f2K"m&X s"j Z4r.@
Fig.6 Selected saltating paths in the air
水利论文L@F5Q:}TX

M~0X N5D-o:A~0

图7 跃移质与蠕移质交换的数目水利论文+\&egP N!\Ubc
Fig.7 Exchange between saltation and surface
3W~ M+GY0creep along surface3

O gESQ0

结论

KqI6C et0S2s!l0

  离散颗粒动力学直接给定沙粒的真实物理特性,能够模拟风沙运动从非稳态到稳态的发展过程,可以反映风沙运动的动态变化特性;与其它模型相比,该模型的研究范围较广,涉及风沙运动规律和床面形态的发生与发展等关键问题,因而本模型更适用于风沙运动的研究。水利论文D D q6]V6i

参 考 文 献

J7`{oI9Po RJg tw0

[1] 钱宁, 万兆惠。泥沙运动力学[M]。北京: 科学出版社,1983. 493~516.水利论文 h'u:uo Dx7fkMWh

[2] 吴正。中国沙漠与海岸沙丘研究[M]。北京: 科学出版社, 1997. 1~8.[3] Haff P K, Anderson R S. Grain scale simulation of loose sedimentary beds: the example of grain bed impacts in aeolian saltation[J]。Sedimentology, 1993, 40: 175~198.水利论文)\O X0wGU*c;O

[4] Landry W, Wener B T. Computer simulation of self organized wind ripple patterns[J]。Physica D, 1994, 77: 238~360.水利论文Oi|Zs4| a^4m

[5] Tsuji Y, Kawaguchi T, Tanaka T. Discrete particle simulation of two dimensional fluidized beds[J]。Chemical Engineering Science, 1993, 77(1): 79~87.水利论文6\q oqh0}4\*} q

[6] Hoomans B P B, Kuipers J A M, Briels W J, et al. Discrete particle simulation of bubble and slug formation in a two-dimensional gas fluidized bed: A hard sphere approach[J]。Chemical Engineering Science, 1996, 51(1): 99~108.水利论文-^(~Kzb}f

[7] Ouyang Jie, Li Jinghai. Particle motion resolved discrete model for simulating gas solid fluidization[J]。Chemical Engineering Science, 1999, 54(13): 2077~2083.

i't;K7T#U:gYG0

[8] 孙其诚, 李静海。鼓泡流化床中气泡行为的模拟[J]。化工冶金。2000, 21(2):216~218.水利论文7UE,{Mf!bJ

[9] 朱朝云, 丁国栋, 杨明远。风沙物理学[M]。北京: 中国林业出版社,1992. 1~25.水利论文c~?y%y#a1r*JIhh

[10] 韩其为, 何明民。泥沙起动规律及起动流速[M]。北京: 科学出版社, 1999. 3~9.水利论文a,F)t*y*k I-BB3U

[11] Wang Y, Mason M T. Two dimensional rigid body collisions with friction[J]。ASME Journal of Applied Mechanics, 1992, 59:635~642.

6T5Mc]mlyW`"g,S0
TAG: 动力学 风沙 颗粒 王光谦 孙其诚
顶:23 踩:32
【已经有166人表态】
19票
极差
17票
很差
18票
较差
21票
稍差
26票
稍好
19票
较好
23票
很好
23票
极好
下一篇:黄河下游河槽萎缩与防洪(黄金池)
上一篇:风洞中挟沙气流水平集沙实验研究(倪晋仁, 李振山)
查看全部回复【已有0位网友发表了看法】

广告投放

广告投放