边界层微分方程

边界层微分方程（精选8篇）

边界层微分方程 第1篇

在计算流体力学(computational fluid dynamics,CFD)和计算气动声学(computational aeroacoustics,CAA)中,计算区域不可能无限大,因此需要人为地将计算区域和远场分隔开来。为了保持计算稳定,需要在分界处添加数值吸收边界(有时也称为辐射边界、非反射边界等)条件。这种数值边界起两个作用:一是让计算区域内的各种波(声波、涡波、熵波)能无反射地通过边界;二是让外界的任何扰动都不能传入计算区域。吸收边界一直是计算气动声学领域的一个研究热点和难点。现在主要的吸收边界有Thompson提出的特征边界、Giles提出的基于Fourier分析的边界、Tam和Webb提出的辐射边界[1]、理想匹配层(perfectly matched layer,PML)边界[2]等。在所有的这些边界中, PML边界的效果最佳。但Tam等[3]指出在有平均流的情况下PML边界也不能完全吸收所有反射波[4]。Hu通过坐标变换方法改进了原有的PML边界,提出了在有非均匀平均流的情况下稳定的PML边界条件,并证明了PML边界的分离形式与非分离形式是等价的。此外,他还提出了非线性欧拉方程和纳维-斯托克斯方程的PML边界条件[5,6,7,8]。 但Hu的PML边界条件还存在以下两个问题:一是变换系数β是个经验参数,缺乏严格证明;二是在采用非线性欧拉方程的PML边界条件时,需要事先假定一个平均流场,这样不利于编程实现。本文主要解决以上两个问题,并给出了计算实例予以验证。

1 线性欧拉方程变换系数推导

假定平均流只位于x方向,速度为um,密度为ρm,压力为pm,则线性欧拉方程可表达为

式中,ρ′为扰动(相对于平均流而言)的密度;u′为扰动x方向的速度;v′为扰动y方向的速度;p′为扰动压力;γ为空气绝热指数。

联立式(1)中的后三式,消去u′、v′可得

$(\frac{\partial }{\partial t}+u{}_{\text{m}}\frac{\partial }{\partial x}){}^2{p}^{\prime }-c{}^2(\frac{\partial {}^2{p}^{\prime }}{\partial y{}^2}+\frac{\partial {}^2{p}^{\prime }}{\partial x{}^2})=0$ (2)

$c=\sqrt{\gamma p{}_{\text{m}}/\rho {}_{\text{m}}}$

式中,c为声速。

式(2)即为对流的波动方程。由于空气对流的作用,波的群速度和相速度不一定沿着同一方向。如果此时在PML中引入阻尼,则不能完全消除反射波,从而影响数值稳定性。为了解决这个问题,需要对式(2)进行适当的线性坐标变换。坐标变换有两种方式,一种是保持时间不变,进行空间变换,即

其中,$\overline{t}$、$\overline{x}$分别为新坐标系下的时间和位移,但在计算中,所有的计算网格都已经划分完毕,如果再采用这种变换,需要重新映射网格,势必引起更多的麻烦。所以只能够对时间实行变换,即

其中,β为一常数,其他的所有量都保持不变。这种变换的目的是为了使式(2)具有与平均流为零情况下的波动方程一样的特点,然后再引入PML边界就不会引起数值不稳定。

经式(4)变换后,有

其中,下标t、x分别表示对原坐标下的t、x求偏导。将式(5)代入式(2)可得

$\begin{array}{l}\frac{\partial {}^2{p}^{\prime }}{\partial \overline{t}{}^2}+2u{}_{\text{m}}(\frac{\partial {}^2{p}^{\prime }}{\partial \overline{t}\partial \overline{x}}+\beta \frac{\partial {}^2{p}^{\prime }}{\partial \overline{t}{}^2})+\\ u{}_{\text{m}}^2(\frac{\partial {}^2{p}^{\prime }}{\partial \overline{x}{}^2}+2\beta \frac{\partial {}^2{p}^{\prime }}{\partial \overline{t}\partial \overline{x}}+\beta {}^2\frac{\partial {}^2{p}^{\prime }}{\partial \overline{t}{}^2})-\\ c{}^2(\frac{\partial {}^2{p}^{\prime }}{\partial \overline{x}{}^2}+2\beta \frac{\partial {}^2{p}^{\prime }}{\partial \overline{t}\partial \overline{x}}+\beta {}^2\frac{\partial {}^2{p}^{\prime }}{\partial \overline{t}{}^2}+\frac{\partial {}^2{p}^{\prime }}{\partial y{}^2})=0(6)\end{array}$

如果能将式(6)中的混合导数项消去,那么其在本质上与静止时的声音控制方程一样。所以应该保证:

um+βu2m-βc2=0 (7)

也即

$\beta =\frac{u{}_{\text{m}}}{c{}^2-u{}_{\text{m}}^2}$ (8)

因此坐标变换为

即通过式(9)的坐标变换,由式(1)决定的声波在传播的过程中与无平均流时的声波在本质上具有相同的特点。所以在采用非线性欧拉方程的PML边界条件时,需要先通过式(9)的坐标变换。

如果速度量纲以声速c为单位,平均流的马赫数Ma为um/c, 则β可以改写为

$\beta =\frac{{\rm M}a}{1-{\rm M}a{}^2}$ (10)

这个系数正是Hu[6]在推导线性欧拉方程的PML方程时用的经验系数。此处从声学方程的角度对这个系数进行了严格的推导。

2 笛卡儿坐标系下全欧拉方程的PML方程

笛卡儿坐标系下的全欧拉方程为

现假定所有的平均流只限于x方向,并将所有的流动变量分离成平均量与扰动量之和,那么

式中,ρm、um、pm分别为平均流的密度、速度和压力;ρ′、u′、p′分别为扰动量的密度、速度和压力。

由于假定平均流只限于x方向,vm=0,那么将式(12)代入式(11)可得

如果写成矩阵的形式,式(13)可表示为

$\frac{\partial {\mathit{U}}^{\prime }}{\partial t}+\mathit{A}\frac{\partial {\mathit{U}}^{\prime }}{\partial x}+\mathit{B}\frac{\partial {\mathit{U}}^{\prime }}{\partial y}=0$ (14)

$\begin{array}{l}{\mathit{U}}^{\prime }=\left[\begin{array}{c}{\rho }^{\prime }\\ u^{\prime }\\ v^{\prime }\\ p^{\prime }\end{array}\right]\\ \mathit{A}=\left[\begin{array}{cccc}u{}_{\text{m}}+u^{\prime }& \rho {}_{\text{m}}+{\rho }^{\prime }& 0& 0\\ 0& u{}_{\text{m}}+u^{\prime }& 0& \frac{1}{\rho {}_{\text{m}}+{\rho }^{\prime }}\\ 0& 0& u{}_{\text{m}}+u^{\prime }& 0\\ 0& \gamma (p{}_{\text{m}}+p^{\prime })& 0& u{}_{\text{m}}+u^{\prime }\end{array}\right]\\ \mathit{B}=\left[\begin{array}{cccc}{v}^{\prime }& 0& \rho {}_{\text{m}}+{\rho }^{\prime }& 0\\ 0& v^{\prime }& 0& 0\\ 0& 0& v^{\prime }& \frac{1}{\rho {}_{\text{m}}+{\rho }^{\prime }}\\ 0& 0& \gamma (p{}_{\text{m}}+p^{\prime })& v^{\prime }\end{array}\right]\end{array}$

根据Hu[5,6,7]的方法,可以按照以下三步推导PML的方程。首先将式(14)进行如式(9)所示的坐标变换,可得

$(\mathit{{\rm I}}+\beta \mathit{A})\frac{\partial {\mathit{U}}^{\prime }}{\partial \overline{t}}+\mathit{A}\frac{\partial {\mathit{U}}^{\prime }}{\partial \overline{x}}+\mathit{B}\frac{\partial {\mathit{U}}^{\prime }}{\partial y}=0$ (15)

式中,I为单位矩阵。

假定时间上U′为单个Fourier成分,即${\mathit{U}}^{\prime }={\stackrel{⌢}{\mathit{U}}}^{\prime }\text{e}{}^{-\text{i}\omega t}$,则在频域中,式(15)变为

$-\text{i}\omega (\mathit{{\rm I}}+\beta \mathit{A}){\stackrel{⌢}{\mathit{U}}}^{\prime }+\mathit{A}\frac{\partial {\stackrel{⌢}{\mathit{U}}}^{\prime }}{\partial \overline{x}}+\mathit{B}\frac{\partial {\stackrel{⌢}{\mathit{U}}}^{\prime }}{\partial y}=0$ (16)

然后,通过复数变易法,引入阻尼[6],即

$\frac{\partial }{\partial \overline{x}}\to \frac{1}{1+\frac{\text{i}\sigma }{\omega }}\frac{\partial }{\partial \overline{x}}$ (17)

其中,σ为一个吸收系数,它可以为一个常数,也可以为在空间上变化的函数。则式(16)变为

式(18)便是频域中的PML方程,但所有的计算都是在时域中进行的,所以最后一步就是将式(18)再变回时域并变回原坐标。将式(18)变回时域为

最后,将式(19)变回原坐标为

如果略去式(20)中的高阶小量,则可直接得到线性欧拉方程的PML方程,即

$\begin{array}{l}{\mathit{U}}^{\prime }=\left[\begin{array}{c}{\rho }^{\prime }\\ u^{\prime }\\ v^{\prime }\\ p^{\prime }\end{array}\right]\\ {\mathit{A}}^{\prime }=\left[\begin{array}{cccc}u{}_{\text{m}}& \rho {}_{\text{m}}& 0& 0\\ 0& u{}_{\text{m}}& 0& \frac{1}{\rho {}_{\text{m}}}\\ 0& 0& u{}_{\text{m}}& 0\\ 0& \gamma p{}_{\text{m}}& 0& u{}_{\text{m}}\end{array}\right]\\ {\mathit{B}}^{\prime }=\left[\begin{array}{cccc}0& 0& \rho {}_{\text{m}}& 0\\ 0& 0& 0& 0\\ 0& 0& 0& \frac{1}{\rho {}_{\text{m}}}\\ 0& 0& \gamma p{}_{\text{m}}& 0\end{array}\right]\end{array}$

式(21)与文献[8]中的方程相比,形式大大简化,更易于编程实现。

3 柱坐标系下全欧拉方程的PML方程

如果不考虑流体在周向的流动,则柱坐标系下的全欧拉方程为

$\begin{array}{l}\frac{\partial \rho }{\partial t}+v\frac{\partial \rho }{\partial r}+u\frac{\partial \rho }{\partial x}+\rho (\frac{\partial v}{\partial r}+\frac{v}{r}+\frac{\partial u}{\partial x})=0\\ \frac{\partial u}{\partial t}+v\frac{\partial u}{\partial r}+u\frac{\partial u}{\partial x}+\frac{1}{\rho }\frac{\partial p}{\partial x}=0\\ \frac{\partial v}{\partial t}+v\frac{\partial v}{\partial r}+u\frac{\partial v}{\partial x}+\frac{1}{\rho }\frac{\partial p}{\partial r}=0\\ \frac{\partial p}{\partial t}+v\frac{\partial p}{\partial r}+u\frac{\partial p}{\partial x}+\gamma p(\frac{\partial v}{\partial r}+\frac{v}{r}+\frac{\partial u}{\partial x})=0\end{array}\}$

(22)

现假定平均流沿着x方向,与笛卡儿坐标系下全欧拉方程类似,将所有流体变量分解为平均量与扰动量之和,则式(22)可以改写成

$\frac{\partial {\mathit{U}}^{\prime }}{\partial t}+\mathit{A}\frac{\partial {\mathit{U}}^{\prime }}{\partial x}+\mathit{B}\frac{\partial {\mathit{U}}^{\prime }}{\partial r}+\frac{1}{r}\mathit{C}{\mathit{U}}^{\prime }=0$ (23)

$\mathit{C}=\left[\begin{array}{cccc}0& 0& \rho {}_{\text{m}}+{\rho }^{\prime }& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& \gamma (p{}_{\text{m}}+p^{\prime })& 0\end{array}\right]$

其中,矩阵U′、A、B与笛卡儿坐标系下相同。

采用第一节推导的变换系数β,然后按照第二节的三步变换法,可以得到柱坐标系下全欧拉方程的PML方程如下:

其中,矩阵A、B、C与式(23)一致。如果略去高阶小量,也可以得到柱坐标系下线性欧拉方程的PML方程。

4 应用算例

为了验证全欧拉方程的PML方程的正确性,现利用其计算喷气发动机进气道(压气机前)的流场。假定气体在周向的速度为零,则该三维问题就转化为一个二维问题,另外忽略气体的黏性(低马赫数下这种假设是可行的)。整个计算区域如图1所示。其中,B4B5之间的粗黑线表示机匣(casing),B8B9之间的粗线表示转子(spinner),B1B9为对称轴,B1B2B3B4表示计算区域的外边界,B5B6B7B8为外加的PML区域,真实的流场并不需要这一区域。

计算时以机匣的内径D为单位长度,以空气中的声速c为单位速度,以D/c为单位时间,以空气密度ρ0为单位密度,以ρ0c2为单位压力。控制方程采用式(22),空间离散采用色散关系保持格式(DRP),时间积分采用优化的多步时间格式[9]。所用边界条件如下:在远离进气道的B1-B2-B3-B4边界上,采用辐射边界条件[9];在机匣和转子表面采用外伸点固壁滑移边界条件(ghost point methods)[10];在B5、B6、B7、B8所围成的区域中,采用PML作为出流边界条件,即先确定PML中的平均流在x方向的速度,然后通过等熵关系式计算出压力、密度、等流体变量,将这些量作为式(24)中平均流的流体变量。通过PML与内部区域不断交换数据,可使与PML交界的内部流场的各变量接近或者达到PML中设定的平均流场变量值。

由于DRP格式要求计算网格不能扭曲,所以不能直接在图1 所示的计算区域中进行计算。需要先通过保形映射将机匣和转子映射到计算平面,然后在计算空间中划分矩形网格,并将式(22)转化成计算空间中的形式求解,最后将计算的结果再通过逆变换转换到图1所示的物理空间中[11]。

马赫数为0.3时的流线图如图2和图3所示。从流线图PML边界条件可以看出,周围的气体都吸入到了进气道中,在气流出口处没有反射现象发生,这说明了全欧拉方程的正确性,同时也说明了利用非线性PML作为出流边界条件是可行的。

当然也可以利用有限体积法来计算流场,这里主要是用此例来证明全欧拉方程PML边界条件的正确性。

参考文献

[1]Hixon R,Shih S H,Mankbadi R R.Evaluation ofBoundary Conditions for Computational Aeroaocus-tics[J].AIAA Paper 95-0160.

[2]Hu Q.On Absorbing Boundary Conditions for Line-arized Euler Equations by a Perfectly Matched Lay-er[J].Journal of Computational Physics,1996,129(1):201-219.

[3]Nizampatnam L S,Hoffmann K A,Papadakis M,etal.Investigation of Boundary Conditions for Compu-tational Aeroacoustics[J].AIAA Paper 99-0357.

[4]Tam C K W,Auriault L,Cambuli F.PerfectlyMatched Layer as an Absorbing Boundary Conditionfor the Linearized Euler Equations in Open andDucted Domains[J].Journal of ComputationalPhysics,1998,144(1):213-234.

[5]Hu F Q.A Stable,Perfectly Matched Layer for Lin-earized Euler Equations in Unsplit Physical Varia-bles[J].Journal of Computational Physics,2001,173(2):455-480.

[6]Hu F Q.A Perfectly Matched Layer AbsorbingBoundary Condition for Linearized Euler Equationswith a Non-uniform Mean-flow[J].Journal ofComputational Physics,2005,208(2):469-492.

[7]Hu F Q.On the Construction of PML AbsorbingBoundary Condition for the Non-linear EulerEquations[J].AIAA Paper 2006-0798.

[8]Hu F Q,Li X D,Lin D K,et al.Absorbing BoundaryConditions for Nonlinear Euler and Navier-stokesEquations Based on the Perfectly Matched LayerTechnique[J].Journal of Computational Physics,2008,227(9):4398-4424.

[9]Tam C K W,Webb J C.Dispersion-relation-pre-serving Finite Difference Schemes for Computa-tional Acoustics[J].Journal of ComputationalPhysics,1992,107(2):262-281.

[10]Tam C K W,Dong Z.Wall Boundary Conditionsfor High-order Finite-difference Schemes inComputational Aeroacoustics[J].Theoretical andComputational Fluid Dynamics,1994,6(5/6):303-322.

边界层微分方程 第2篇

北京及华北平原边界层大气中污染物的汇聚系统--边界层输送汇

由地面、3个大气压风带(1 000,925,850 hPa)、气压场、温湿廓线、能见度、暖湿平流以及激光雷达,航测分析,提出太行山山前、燕山山前输送汇.输送汇及其摆动常造成华北平原及北京地区区域大气污染物汇聚,是形成重污染区的主要形式.输送汇配置各类风向风带区,冬秋季节多为区域性短而浅的边界层输送风带,夏季多为天气尺度型的.、远而厚的边界层输送风带.弱气压场或均压场背景下地方性山风及山前串状城市热岛群形成的热力性、动力性低压环流,是输送汇形成的主要原因.

作 者：苏福庆 任阵海 高庆先 张志刚  作者单位：国家环保总局,气候变化影响研究中心,北京,100012 刊 名：环境科学研究  ISTIC PKU英文刊名：RESEARCH OF ENVIRONMENTAL SCIENCES 年，卷(期)： 17(1) 分类号：X513 关键词：输送汇   边界层   边界层输送风带  

基于分数阶微分的边界扫描 第3篇

关键词：分数阶微分,边缘检测,微分阶数,掩膜模板

图像的边缘检测是图像分割的一种重要手段,是图像的一种紧描述,为目标识别和图像解释提供了一种有价值的特征参数[1]。在图像处理中,基于一阶和二阶微分的边缘提取算子是边缘检测常用的方法,而如何减少噪声的影响是该方法需要解决的问题[2,3,4,5]。理论研究表明,对信号进行阶微分运算,当微分阶数较小,时,信号高频分量被大幅提升,中低频有所加强,甚低频进行了非线性的保留。本文从经典分数阶微分定义中的分数阶差分方程出发[6],给出分数阶微分的掩膜模板,实验表明,基于分数阶微分的边缘检测可以有效地提取图像的边缘信息。

1 分数阶微分运算用于图像处理

1.1 微分运算对信号的作用

现有任一能量型函数f(t)∈L2(R),设其傅里叶变换为。假设f(t)的整数k(k∈Z+)阶微分存在,则其傅里叶变换为,其中,称为k阶微分乘子函数。的指数形式为:

将上式阶数k推广至任意阶v(v∈R2),可得算子Dv,相应的分数阶微分f(v)(t)在频域的形式,其中乘子d赞v(ω)在频域的指数形式为:

信号的分数阶微分从信号处理的角度可以理解为广义的调幅调相,其中振幅随频率的变化呈分数阶幂指数变化,而相位为频率的广义Hilbert变换。由(1)(2)两式可画出一阶、二阶和分数阶微分的幅频特性曲线(略)。分析可得,微分运算有提升信号高频、消弱信号甚低频的作用,且二阶微分明显强于一阶微分。可见,分数阶微分一方面可加强信号高中频分量,另一方面可对信号甚低频分量进行非线性保留。

1.2 分数阶微分用于图像边缘处理

对于图像,低频对应图像的平滑区,即图像的非边缘区;中频对应图像的纹理信息;高频对应图像的边缘区和噪声区。用整数阶微分进行处理,图像中灰度变化不大的纹理细节信息将会遭到大幅的线性衰减,导致其结果几乎为零,不能很好的检测出平滑区域的纹理细节信息;而用分数阶微分进行处理,这些信息则可在一定程度上得到非线性保留。可见,分数阶微分比整数阶微分更有利于提升图像平滑区中的中、高频信息。

2 分数阶微分算子的实现

目前为止,分数阶微分并没有一个统一的标准定义,比较著名的是G-L定义。该定义将微积分的阶数由经典定义的整数扩展到分数阶推导而来。

其中,函数,其值有限,图像信号灰度的最大变化量也是有限的。由于两相邻像素是灰度变化发生的最短距离,对于二维数字图像,x和y轴方向上持续时间的度量就要以像素为单位。若一元信号f(t)的持续期t∈[a,t],将其按单位h=1进行等分,得,则信号f(t)分数阶微分的差分表达式为:

一般来说,在M×N大小的图像f上,用m×n大小的滤波器掩膜进行线性滤波:

其中a=(m-1)/2且b=(n-1)/2。为了得到一幅完整的经过滤波处理的图像,需要对图像中的所有像素点进行处理,也就是说要依次对x=0,1,2,…,M-1和y=0,1,2,…,N-1应用上式。

若二维数字图像信号中x和y的持续期间为x∈[x1,x2]和y∈[y1,y2],对图像的微分运算可以使用3×3模板、5×5模板或者7×7模板中当前像素点的4个方向或者8个方向进行,根据(6)式可得模板一至模板六。如图2所示。

模板四是3×3模板当前像素点的8个不同方向,各个方向上的分量依次为:1,-v;

模板五是5×5模板当前像素点的8个不同方向,各个方向上的分量依次为:1,-v,(v×v-v)/2;

模板六是7×7模板的8个不同方向,各个方向上的分量依次为:

3 图像边缘提取实验仿真

对于二维灰度图像,边缘和噪声对应图像中局部特性的不连续点,即对应邻域像素的灰度发生了较大变化,为高频信号。两者也存在区别,即边缘具有邻域一致性、结构性和方向性,且相对于噪声,边缘具有较大的能量和范围。我们可以利用临近像素点信息减少或抵消噪声的作用,适当扩大信号处理的范围,从而加强边缘信号。

3.1 不同模板实验效果分析

通过模板一、模板二和模板三实验效果分析表明,当分数阶为0.5时边缘既清晰噪声又比较少。

通过模板四、模板五和模板六实验效果分析表明,阶数的大小对图像边缘的检测有直接的影响,随着阶数的增加,边缘逐渐清晰。

3.2 分数阶固定效果分析

通过实验效果可以看出,当阶数固定时,随着迭代次数的增加,图像逐渐不清晰,出现双边缘现象,但图像边缘亮度变化基本不大。当分数阶阶数小于0.5时,第一次得到的效果图亮度最好,以后的每次迭代都会削弱图像的亮度信息;当阶数大于0.5时,迭代次数越多,图像越亮,但会伴有噪声的出现。

4 结论

仿真实验表明:基于分数阶微分的边缘检测算法在不同的分数阶阶数时具有不同的提取图像边缘信息的能力,总体来说随着阶数的增大,边缘提取的能力逐渐增强;但是边缘的有效提取虽然随着模板大小的增大而有所提升,但不总是成正比例提升。但基于分数阶微分的边缘检测相对于传统的一阶和二阶微分算子的边缘检测效果有很好的抗噪性。

参考文献

[1]才辉,张广新,张浩,等.一种新的基于多信息测度融合的边缘检测方法[J].浙江大学学报:工学版,2008,42(10):1671-1675.

[2]Jeong H,Kim C.Adaptive determination of filter scales for edge detection[J].IEEE Trans on Pattern Analysis and Machine Intelli-gence,1992,14(5):579-585.

[3]Turgut A Y,Yemez E A,Bulem S.Multidirectional and multiscale edge detection via M-band vavelet transform[J].IEEE Trans ImageProcessing,1996,5(9):1370-1377.

[4]Elder J H,Zucher S W.Local scale control for edge detection and blur estimation[J].IEEE Trans on Pattern Analysis and Machine In-telligence,1998,20(7):699-716.

[5]Andrew P P.Directional filtering in edge detection[J].IEEE Trans Imaging Processing,1998,7(4):611-615.

边界层微分方程 第4篇

考虑边界波浪方向的缓坡方程自适应求解模型

缓坡方程是描述近岸波浪运动较好的数学模型之一.在发展的自适应有限元求解缓坡方程的基础上,采用迭代求解的方法,确定波浪相对于边界的入射方向,从而对边界条件进行改进,建立了求解缓坡方程的数值计算模型.典型算例表明,考虑波浪相对于边界的.入射角度后,模型可以更好地模拟吸收波浪边界,同时对多向波对双突堤的绕射进行了模拟研究,与试验结果比较表明,所建立的数值计算模型能够适用于多向不规则波传播过程的模拟研究.

作 者：柳淑学 孙冰 LIU Shu-xue SUN Bing  作者单位：大连理工大学,海岸和近海工程国家重点实验室,辽宁,大连,116024 刊 名：海洋工程  ISTIC PKU英文刊名：THE OCEAN ENGINEERING 年，卷(期)： 25(1) 分类号：P731.2 关键词：波向角   自适应有限元法   缓坡方程   边界条件  

边界层通用气象数据采集系统设计 第5篇

1 系统的拓扑结构

该系统由边界层通用气象要素采集模块、GPRS通信网络、防火墙、数据处理中心、应用服务器、客户端组成。其拓扑结构如图1 1所示。

2 系统结构框架图设计

2.1 系统结构

以AT89S52单片机为基础,采用两级模拟开关为核心的传感器,配合上高速的模/数转换器,实现系统自动识别传感器输出信号,利用通过RS-232总线|I2C串行总线交换数据、接收控制指令,由中控系统解析处理后,分别将接收到的气象检测数据以及定位模块信息发送至通信模块传输到应用服务器,通过应用程序中实现信息的可视化。

3 系统硬件设计

3.1 主控模块

主控模块的主要功能是提供控制总线、输入输出总线、数据总线端口,同时具备数据的存储与分析功能。边界层远离地面,能量来源主要依靠太阳能,所以低耗就显得尤为重要,所以在主控模块中,本系统选用的单片机是AT89S52,是一款低功耗、高性能CMOS8位微控制器,并8K容量的可编程Flash存储器,相比于与传统的51单片机更有超强的抗干扰能力。该系统中,单片机主要负责控制总线使用权的分配,对采集到的气象数据、定位数据进行处理,并通过I2C串行总线|RS232总线传至通信模块,以便发送回地面,完成气象数据采集工作。

3.2 以传感器为核心的采集模块

边界层通用气象要素采集系统设计中,既要考虑到经典电路中所要求的准确性,由于受高空条件的限制,在设计时还要尽量做到低功耗,所以在本系统中采用的通用式采集电路能够一路多用,大大降低采集模块的功耗。为满足大阳能供电的低耗的要求,本系统采用的通用式气象要素采集方法,将自动判别传感器输出信号的类型,继而调整电路对应信号的测量,为保证精度,需要软件来支持判定信号类型(温度、湿度、气压、风向等)。

在本系统中,温度、湿度、气压等气象数据采用IC2总线来读入,温度传感器采用DS18B20,其能够适配各种单片机或系统机,测温范围为-55℃到+125℃,测量分辨率为0.0625℃。采用单总线专用技术,既可通过串行口线过其它I/O口线与,也可通微机接口。用户可分别设定各路温度的上、下限。湿度传感器采用HS1101,互换性好,高可靠性,长期稳定性好,响应快。精确度±1.5%,工作温度范围 -40~140℃,量程1~99%RH。压力传感器采用STP-TW,产品精度高,且功耗很低,有过电压、反极性等多重保护,绝压0—110pa,供电电源24v DC,工作温度-40—80℃。系统中所选用的传感器均能长时间工作在大气边界层中,且保持精度,能有效、准确测量温度、湿度、压强等气象要素。

3.4 通信模块

本系统在设计时的主要目的及时能够向比较多的用户能够及时、准确的了解到由边界层的气象状况所推测出的天气情况,所以我们选择通过GPRS与使用者进行联系。GPRS作为新型的传输方式,以封包的方式,无需使用整个通道,这就使使用者所需要的流量最少。并且GPRS的接入时间短,传输速率高,平均速率可达70Kbps,这为数据通信的适时性提供了保障。GPRS可以说是GSM的延续,所以其覆盖范围很大,可以满足全球各地的数据使用要求。通过这种传输方式,边界层处的气象数据能稳定、准确的传到上位机中,对于用户获取信息也比较直接,快捷。

3.5 电源模块

由于边界层处的条件限制,电源要给系统各个模块供电,电源的设计就尤为重要,所以作为单元的模块来设计。考虑到高空太阳能充足,可采用太阳能电池供电方式。这种方式节能环保,且能长时间使用,可降低产品回收的成本,但是也使低耗成为系统设计需要考虑的重要因素。

4 边界层通用气象数据采集模块软件设计

边界层通用气象数据采集模块软件系统的主要任务:

1) 初始化系统各个模块,识别各输出输入信号,判断信号类型;

2) 控制采集模块,实现对指令要求的气象要素信息进行采集;

3) 完成对各通道的数据进行分析,配合硬件对采集到的气象信息进行计算、分析、存储;

4) 实现系统各模块的工作状态的检测并显示;

5) 发送指令,实现与应用服务器的通信。

5 结束语

郑州市春季边界层风气候变化研究 第6篇

郑州在纬向地处丘陵与平原交界处, 径向地处南北气候过渡带。空中水汽输送及气象自然灾害产生形成的微过程与郑州高空气象资料信息系统关系密切。

边界层气候为气候变化主要组成部分。随着气候资源的开发利用日益深入, 对边界层风的科研需求日益明显。

通过对郑州市边界层风的气候研究, 为低空高空风气候业务开展打下基础, 为河南省风能资源开发提供科学依据, 为防御气象减灾提供科学依据, 为高空业务发展奠定基础。以郑州高空站资料为主, 结合中国高空气候特征, 研究本地边界层气候的变化情况。

1.1 资料来源

利用郑州1955-2010年55年的高空气象序列化资料, 并进行均一化处理。主要数据来源于高空气象记录表, 经月报表数据、电子数字化数据核对后纳入采用。

地面气象资料来源于本站档案馆保存的数据, 具有一定的数据质量可信度。

1.2 研究方法

在郑州1955~2012年高空气象资料系统化处理的基础上, 利用聚类分析技术, 结合高空气候资料数据库, 对高空边界层风数据进行处理, 研究其气候变化特征。

类聚分析主要是结合高空风垂直切变特征、季节变化特征、天气过程特征、大气辐射因子、雾霾及灾害性天气过程按群的观念进行类别分析。

在统计方面结合高空风资料的实际情况, 采用基础的统计方式, 以便于交流表达。统计分析过程中的主要表达式有以下几个:

(1) 径向风

风向为风的来向与正北方向一致, 北风为正。风速表达式为:

式中, V径向———径向风风速, V———风速, θ———风向。

(2) 纬向风

风向为风的来向与正东方向一致, 东风为正。风速表达式为:

式中, V纬向———径向风风速, V———风速, θ———风向。

(3) 均值

(4) 样本方差

(5) 相关系数统计分布

其中为Pearson相关系数, 简化后为

(6) 类聚分析定距型变量采用切比雪夫距离

(7) 类聚分析定序型变量采用卡方距离

2 资料系统分析

2.1 资料系统

郑州高空风探测系统建于1955年。截止目前, 高空风原始数据采集经历光学经纬仪、马拉黑、910雷达、701雷达、C波段雷达、L波段雷达等探测装备, 高空风原始数据处理经历计算尺、表格法、绘图板法、初级可编程计算器、高级可编程计算器、微型计算机等发展阶段, 采集手段和数据处理能力经历滚动式发展。高空风数据可靠性进一步提高。

郑州的高空风资料处理主要分为柱面坐标系统高空风数据和球面坐标系统高空风数据。以光学经纬仪、马拉黑为代表的柱面坐标系统高空风数据, 仅采集仰角、方位角, 风向风速的其他值通过递推计算求得。自1964年郑州站应用910雷达以来, 球面坐标系统高空风数据开始用于高空风求取, 原始数据有仰角、方位角、斜距。

球面坐标系统高空风数据的完整性优于柱面坐标系统高空风数据。

2.2 资料生成背景

高空风资料生成的关键环节有原始数据采集系统和高空风资料处理系统。

高空风原始数据采集系统为高空风探测基础。高空风原始数据采集依据探测手段的不同分为人工采集和器械采集, 而器械采集分为有源目标和无源目标的采集。现行常规高空风观测属于有源目标观测, 而风廓线雷达属于无源目标采集系统。对于常规高空风观测, 郑州站采用目标轨迹跟踪技术, 分为近设备跟踪和远设备跟踪两种。对于远设备跟踪 (距雷达天线距离一般为180米左右) , 计算高空风过程中进行了有效订正, 有效采集数据一般为放球30秒后。

依据数据处理的自动化程度, 高空风资料处理系统分为原始初级高空风处理系统和信息化处理系统。

3 雾霾沙尘与边界层风相关性分析

春季为郑州雾霾天气出现频率较高的季节, 边界层风具有明显的气候特点。

表1为2008年4月5日至8日边界层风。从表中可以看出, 高空有较强偏西气流, 边界层下部较强的偏南气流为扬沙天气的流场配置要素之一, 而前期的边界层风的稳定度差为判别的参量之一。

该表比较典型地反映出了郑州雾发生的气象场配置。

表1明显反映出, 地面为静风或偏南的微风, 边界层下部为偏南风, 高空虽为偏西气流, 但存在向东风顺向偏转的趋势的情况下, 浓雾产生。对比地面气象观测资料, 2008年4月7日5时至6时产生大雾, 能见度仅有600米。

4 结论及讨论

4.1 结论

综合以上分析可以得出如下结论:

(1) 郑州风速随高度变化与艾克曼风速指数基本吻合, 郑州南偏西气流风向随高度变化以顺时针变化为主征, 变化范围在248~348之间。

(2) 由于受地形等要素影响, 风速风向变化存在一定的波动性。

(3) 在时空上, 郑州春季高空风向变化呈现比较明显的循环性变化特征, 上边界层以西北气流为主。郑州近50年边界层气候资料显示, 除2000年海拔3 000米和2010年距地900米天气系统过境前风向出现切变外, 边界层风向多为西北气流。

(4) 郑州的降水过程产生前, 郑州边界层风变化存在明显特征, 边界层上层西南气流引导, 低空由东北气流转向偏东南气流, 继而旋转为西南气流, 与高空降水带衔接后降水产生, 属于典型的垂直结构的顺时针变化气旋特点。

(5) 春季为郑州雾霾天气出现频率较高的季节。其边界层风具有明显的气候特点。高空有较强偏西气流, 边界层下部较强的偏南气流为扬沙天气的流场配置要素之一。而前期的边界层风的稳定度差为判别的参量之一。

(6) 地面为静风或偏南的微风, 边界层下部为偏南风, 高空虽为偏西气流, 但存在向东风顺向偏转的趋势的情况下, 雾产生的概率较大。

4.2 讨论与思考

通过对郑州高空边界层风的研究, 我们得到如下启示:

高空气象资料经历从无到有, 从比较凌乱到规范性, 资料索取手段从落后、原始到科学、系统、现代化的发展过程, 高空资料的均一性问题是目前资料业务的制约瓶颈。高空探测站点偏少, 探测时次密度稀疏, 不能满足现代气象业务的需求, 应引起气象工作者的重视。

高空气象信息系统的技术资源配置相对滞后, 业务需求与业务现状脱节, 高空业务现状应尽快改进。解决办法是强化新技术资源在基础需求强的高空业务方面的配置。

边界层微分方程 第7篇

关键词：大气边界层,激光雷达,EMD

0 引言

大气边界层 (Atmospheric Boundary Layer, ABL) 是指靠近地球表面, 受地面摩擦阻力等表面强迫力影响并且响应时间尺度约为1小时或更小的大气层[1]。它的高度随天气过程、地形和地表特征而变化, 一般在2公里范围内。大气边界层是与人类关系密切的大气层, 是大气污染物的主要集中层, 边界层高度的大小, 决定了排入大气中污染物的浓度值。所以针对气候及大气环境的研究需要, 开展对大气边界层高度的观测势在必行。同时, 对边界层高度的观测要求也越来越高。

激光雷达作为一种新型大气主动遥感探测工具[2,3,4], 已经广泛应用于大气气溶胶和污染物、云等领域的测量和研究工作。激光雷达利用气溶胶作为示踪粒子, 所测量到的光功率信号正比于大气的气溶胶含量。与此相对应, 大气边界层中含有丰富的气溶胶 (较自由对流层的气溶胶含量更大) , 这是因为气溶胶粒子主要来源于地球表面, 从而导致对激光的更大散射, 因此可以利用激光雷达较容易的探测到这两层之间的边界。

本文利用激光雷达探测系统饿实测数据, 采用了梯度法中的二阶求导法来提取大气边界层高度, 并在分析处理中应用了经验模式分解 (EMD) 方法, 给出了典型天气情况下的提取结果, 并与诊断式求得的高度值进行比对, 吻合度较高。

1 梯度法反演大气边界层高度

因为边界层顶与自由对流层气溶胶浓度差距较大, 所以回波信号的梯度最大值点对应的高度边界层厚度[5,6,7]。常用的梯度计算方法为一阶导数法或二阶导数法。计算激光雷达回波距离校正信号关于高度的二阶导数, 并计算其绝对最小值, 最小值对应的高度定义为过渡层, 即混合层和自由对流层交界的中间位置。一阶导数的最小值对应高度是混合层的顶部。

激光雷达接收到的米散射回波信号可以式 (1) 所示的激光雷达方程表示:

由 (1) 式可得:

激光雷达距离平方校正回波信号P (z) z2在一定程度上反映了大气气溶胶浓度随探测高度变化的情况。由于覆盖逆温的作用, 大量的大气气溶胶粒子富集在大气边界层以内, 这样大气边界层大自由大气层之间的大气气溶胶浓度就会发生变化。P (z) z2的梯度D (z) 的变化即代表着大气气溶胶浓度梯度的变化[8]。

D (z) 定义为D (z) =d[P (z) z2]/dz该梯度变化的最大位置就是大气边界层的高度。

2 EMD方法在信号处理中的应用

经验模式分解 (EMD) 方法的本质是对信号进行平稳化处理, 根据特征时间尺度把一个复杂的数据序列分解成有限项内在模函数 (Instrinsic Mode Functions, IMF) [9,10]。最低频率的IMF分量通常情况下代表原始数据的整体趋势或均值, 确定IMF的标准为: (1) 整个数据集中, 极值点的数目和跨零点的数目必须相等或至多只差一个; (2) 任一点上, 由局部极值定义的上下包络线的均值必须等于零。

作为一种应用, EMD分解方法可以有效地提取一个数据序列的趋势或去掉该数据序列的均值。IMFS是通过筛选过程从数据中提取出的:找出原序列x (t) 的各个局部极大值、极小值, 用三阶样条函数分别进行插值, 相应得到原序列x (t) 的上包络序列值xmax (t) 及下包络序列值xmin (t) ;上下包络线的均值被指定为m10;原序列x (t) 和均值m10的差值h10是第一个组成。

x (t) -m10=h10

理想情况下, h10应该是一个IMF, 但由于非线性数据包络线的平均值可能不同于实际的局部均值, 导致非对称波的存在, 因此重复这一筛选过程就非常必要了。把h10当作原序列, 重复以上处理过程k次, 直至满足内在模函数的定义, 求出第一个内在模函数为止。

h1k=h1 (k-1) -m1k;

c1=h1k;

第一个IMF分量c1应包含原始信号中最短的周期分量。第二个IMF函数c2通过对余数项r1=x (t) -c1进行筛选得到。通常第n个IMF分量通过筛选rn=rn-1-cn得到。筛选过程直到下列任一标准满足: (1) 当分量cn或剩余分量rn比预定值小; (2) 当剩余分量rn变成单调函数时, 从中不能再筛选出基本模式分量。对于有趋势的数据, 其剩余分量就应该是该趋势。如果把分解后的各分量合并起来, 就得到原序列

EMD分解方法的特点是在分解的过程中, 内在模函数保留了原始数据本身的特性, 各项内在模函数的和等于原始数据序列, 这就允许可以只用IMFS来分析和处理数据, 保留了数据本身的特性。

3 EMD方法在信号处理中的应用

本次实验中, 激光雷达以脉冲能量50m J出射激光, 重复频率1Hz, 所接收到的大气回波信号经6分钟时间的平均后存储一次, 距离分辨率为15m。

图1为2014年11月03日晚20:00的大气回波信号 (未校准) , 考虑到大气边界层的高度通常为2km以下, 通常我们只需处理从地面到2km这一段高度的数据。但本文为了使比较结果更形象明显, 故对采样数据全部进行处理。

对该段数据先进行预处理, 接着用上述介绍过的EMD方法分解, 按照频率由高到低得到16个IMFS, 见图2。从各个内在模函数可以看出, 其分解出的前6项IMF代表原始数据中较高频的组成, 对信号的能量贡献最少, 我们可以考虑将其消除, 只处理低频的部分。图3为消除高频项后的信号与原始信号的比较。通过比较可发现, 二者基本吻合, 消除高频项后的信号较原始信号平滑的多, 且仍能显著地描述出原始信号的变化趋势。将经过以上处理的信号采用Menut求取二阶导数的方法, 认为二阶导数最小值所对应的高度为大气边界层高度。从图4可以看出, 由EMD预先处理再求二阶导数, 信号强度最小值对应的高度为1290m。

在此本文把通过以上方法求得的边界层高度值与根据常规气象资料计算出的值比较, 边界层的计算公式采用Benkley and Schulman提供的中性和稳定情况下的关系式:, 式中u*为摩擦速度 (m/s) ;f为地转参数 (1/s) ;a为常数, 由气象资料决定, 不同研究者给出的a值不同。本文中, 我们根据u*值的范围选择a, 文中有两种情况, a=0.31或a=0.56。利用青岛气象台提供的气象资料及以上关系式, 我们得到11月03日18:00时的大气边界层厚度值为1053m。

下面, 我们把2014年10月和11月测量的部分结果在表中作比较, 同时表中亦列出仅根据Menut的距离校准信号的二阶微分最小确定的大气边界层厚度。

从以上实验结果可计算出:激光雷达回波信号用EMD处理结合二阶导数最小值对应的大气边界层高度与气象资料所确定的高度值的相关性高达76.32%, 依据二阶微分最小确定的大气边界层高度与气象资料所确定的高度值的相关性为52.11%。显然, 前者要绝对优于后者。也就是说, 用EMD结合二阶微分最小方法确定大气边界层高度相比较与仅利用二阶微分最小的方法来确定大气边界层高度, 结果得到明显改善, 进一步证明了经验模式分解方法在大气边界层厚度的激光雷达提取应用中的可行性。同时, 由于边界层高度的诊断关系式受探测时大气边界层稳定度的影响, 所以与EMD结合二阶微分最小确定边界层高度的方法相比, 要相对复杂得多

EMD结合二阶微分最小确定边界层高度的方法, 仅根据激光雷达所接收到的回波信号 (与大气中气溶胶含量成正比) , 利用经验模式分解与Menut的二阶微分求导二者的结合, 即可得到大气边界层高度值, 方法简单而有效。

参考文献

[1]Stull R B.Introduction to Boundary Layer Meteorology[M].Norwell:Kluwer Academic Publisher, 1988

[2]Menut, L., Flamant, C., Pelon, J., Flamant, P.H., Urban boundary-layer height determination from lidar measurements over the Paris area, Applied Optics, Vol.38, No.6, 945-954, 1999[Z].

[3]Steyn D G, Bottenheim J W, Thomson R B.Over-view of tropospheric ozone in the lower fraser valley and the Pacific 93 field study[J].J.Atmos.Environ., 1997, 31:2025-2035.

[4]Davis K J, Len Schow D H, Oncley S P, et al.Role of entrainment in surfaceatmosphere interaction over the boreal forest[J].J.Geophy Res, 1997, 102 (D24) :29219-29230.

[5]Arya S P S.Parameterizing the height of the stable atmospheric boundary layer.JAM., 20, 1981:1192-1202[Z].

[6]Estournel C, Guedalia D.Improving the diagnostic relation for the nocturnal boundary layer height.BLM, 53, 1990:191-198[Z].

[7]Surridge A D, Swanepoel D J.On the evolution of the height and temperature difference across the nocturnal stable Bondary Layer.BLM, 40, 1987:87[Z].

[8]王珍珠, 李炬, 钟志庆, 等.激光雷达探测北京城区夏季大气边界层[J].应用光学, 2008, 29 (1) :96-100.

[9]Norden E.Huang The empirical mode decomposition and Hilbert spectrum for nonlinear and non-stationary time series analysis.Proc.R.Soc.Lond A 1998, 454 (1971) :903-995[Z].

翼型结构边界层吹气控制的数值研究 第8篇

关键词：翼型,边界层,吹气控制,开缝,数值模拟

0 引言

众所周知, 翼型结构在国民生产生活中有着广泛的应用, 一方面, 我们需要利用流体通过翼型结构来获得能量[1];在另一些方面, 我们希望尽量减小流体通过翼型结构时产生的阻力, 以节省能源、稳定翼型结构的工作状态[2-3]。人们发现, 当翼型前来流方向 (即攻角) 大于某一值 (即临界攻角) 时, 稳定压力面的绕流流场将会发生非定常分离和失速, 翼型的升力系数突然降低, 阻力系数迅速增大。就目前而言, 在众多的减阻控制技术中, 层流流动控制是一项非常有效的减小摩擦阻力的控制技术, 它通过采取控制措施抑制边界层内部各种不稳定扰动的发展和放大, 使失稳的边界层变得稳定, 从而延迟边界层转捩的发生。因此, 将翼型开缝吹气技术应用于改善翼型气动特性有重要的意义。

1 计算模型和数值计算方法

对G4-73No.8D翼型进行了等比例的缩小, 翼型弦长C=20cm。设定计算域四周边界距翼型表面的距离均为15倍。为了在满足计算精度的前提下尽可能减少网格数量, 提高计算效率, 翼型附近网格加密处理如图1所示, 网格由翼型表面向四周扩散。

本文以FLUENT软件作为流场求解器, 采用雷诺时均守恒Navier-Stokes方程。湍流模型选取Realizable k-ε两方程模型。入口边界设置为速度进口边界条件, 并给定初始湍流参数。出口边界设置为压力出口, 翼型表面采用无滑移边界条件。

2 数值模拟结果

2.1 开缝前的模拟

空气来流速度为50m/s, 攻角为-30°时模拟得到的翼型附近的速度分布云图如2图所示:

速度分布云图

由图2可以看出, 在-30°这样一个较大的攻角下, 风机翼型腹面形成旋涡从前缘开始几乎已经占据整个翼型, 在翼型表面产生了滞至流体区域, 且随着气流流向, 滞止区域逐渐增大, 从而排挤上游来流, 使整个边界层发生分离;同时可以看出, 在滞止区域内有较大的漩涡产生, 这些现象都极不利于叶片的气动特性, 会对叶片本身的升阻力产生较大影响。

经分析可知, 在攻角较小时, 翼型壁面附近的流场都处于顺压梯度 (dp/dx<0) 区, 此时, 由于顺压差和层外势流的加速, 边界层内的流体始终保持向下游流动, 流体质点沿翼型表面前进不会停滞, 也不会出现边界层分离现象。而在攻角较大时, 翼型壁面附近的流场就会出现逆压梯度 (dp/dx>0) 区, 逆压差和层外势流的减速使得边界层中流动减速, 边界层内滞止区域增大, 使得压强升高, 滞止流体发生回流, 在边界层内形成涡流。

2.2 开缝处理及模拟结果分析

由图2所表现的翼型的大攻角绕流速度分布可以看出, 在此工况下在翼型压力面处的流体微团流速很小, 并且具有很高的压力。故确定开缝方案为:前缘端开缝位置为翼型几何弦长2%处, 在翼型压力面, 为来流进口端;尾缘端开缝位置为翼型几何弦长70%处, 在翼型吸力面, 为流体出口端;开缝宽度为1mm。网格划分方案为, 对开缝内部采用尺寸函数加密;在外部流场与开缝入口出口相接区域内, 由于流体流速流向发生较大变化, 也采用尺寸函数进行加密, 并采用自适应网格进行修正;对于其他边界层部分及远场计算域的网格划分, 均参照未开缝模型进行。

风机开缝翼型在来流速度为50m/s时的模拟结果如图3所示:

对图2与图3做比较分析, 可以发现, 在相同攻角下, 开缝风机翼型的涡流区域明显减小, 而且攻角越大, 减小效果越明显。开缝使得压力面的高压流体经过直接由吸力面吹出, 吹出的气流不仅破坏了物面的连续性, 限制了边界层的发展, 而且通过开缝缝隙的气流的吹除作用, 使得涡流能够迅速脱落, 从而减少了滞止流体的数量, 减小了涡流区的面积。翼型开缝后压力面与吸力面的压力差较开缝前大大的减小。此时高压流体通过缝隙自行进入低压处, 起到了类似于平衡孔的作用, 有效的改善了翼型在恶劣工况下的气动性能。但同时可以看到, 由于开缝位置的相对固定, 单个模型对于翼型升阻力的影响规律可能并不明确。已图3为例, 尽管开缝平衡了压力面和吸力面的压力, 也限值了吸力面后半段滞止区域的过分发展, 但由于攻角较大, 在翼型前段就已经发生较大区域的边界层分离与流体滞止, 因而使得开缝后的吸力面出现了多个不均匀的涡流区, 且可以推断随着时间变化, 这些涡流区将不停的发生变化, 这一现象很不利于翼型升阻力的优化。所以对于开缝方案对升阻力的影响规律还需进一步研究和分析。

参考文献

[1]杨科, 王会社, 徐建中等.开缝式风力机静态失速特性的研究[J].工程热物理学报, 2008, 29 (01) :32-35.

[2]胡丹梅, 李佳, 闫海津.水平轴风力机翼型动态失速的数值模拟[J].中国电机工程学报, 2010, 30 (20) :106-111.