CN105426594A - 一种基于时空场和条件插值的平稳均质风场快速模拟方法 - Google Patents

Abstract

本发明公开了一种基于时空场和条件插值的平稳均质风场快速模拟方法，基于时空场和条件插值的平稳均质对风场进行高效高精度较高快速模拟，其特征在于，建立多点均质平稳随机过程与多点均质平稳时空随机场之间的联系，确定两者之间转化的条件以及基于时空场模拟出的多点风场样本所具有的统计与各态历经特；采用POD技术对非等间距点的风速时程进行条件插值，采用本发明的方法，不需要传统方法中经常使用的Cholesky分解以及改进模拟方法中相干函数需服从指数分布的假定。此外，还可以使用二维FFT技术极大的提高模拟效率。数值算例分析表明，本发明具有易于使用、精度较高以及模拟效率很高的特点，可有效解决随机风场模拟点数很大时，谱表示方法存在的模拟效率低下的问题。

Description

一种基于时空场和条件插值的平稳均质风场快速模拟方法

技术领域

本发明属于随机信号模拟领域，具体涉及一种基于时空场和条件插值的平稳均质风场快速模拟方法。

背景技术

虽然频域分析方法被广泛应用于工程实际，但是MonteCarlo模拟方法则被广泛应用于非线性、系统随机性以及其他相关的问题。例如，在大跨桥梁和高层建筑的风致动力响应分析中，MonteCarlo模拟方法更容易考虑结构和气动非线性。此外，该方法也经常用于评价其它诸如频域分析方法的准确性。该方法的一个重要环节则是模拟诸如地震动、风场以及脉动波等随机激励的样本。在众多模拟方法中，谱表示方法由于具有严密的理论推导以及容易使用等特点而被广泛应用^[2-5]。随着FFT技术的引入，平稳高斯过程和随机场的计算效率得到显著提高^[6]。随后，Deodatis和Shinozuka对平稳随机过程和随机场的模拟进行了更细致的探讨与研究^[7,1]。

尽管谱分解方法在平稳高斯随机过程的模拟当中得到了成功的应用，但是仍然值得进一步提高其模拟效率，尤其是在诸如高层建筑和大跨桥梁等模拟点数较多的情况。例如，在大跨桥梁的风致抖振响应分析中，为更好的重构主梁的风场，可能需要成百上千的风场模拟点。此外，在高速列车经过桥梁时，也可能需要上千个模拟点进行风车桥耦合系统的响应分析。对于建立在诸如沿海和平坦等均匀地貌的桥梁而言，主梁上的风速场可以视为多点均质的随机过程。即每个风速点的功率谱密度函数相等，互谱密度函数只为距离的函数，而与位置无关。在实际应用中，可能需要一些非等间距点的模拟时程进行结构动力响应分析，这会在一定程度上增加模拟的计算复杂度。

采用传统谱分解方法模拟上述多点随机过程可能需要很长的计算时间和内存。例如，Deodatis提出的双指标频率方法可以产生各态历经的样本时程。但是，当模拟点数较多时，该方法需要进行很多次功率谱矩阵的分解。此外，该方法模拟出的样本周期也很长，从而进一步降低了模拟效率以及增加了响应分析的计算量^[7]。为了提高模拟效率，目前比较普遍的方法则是减少功率谱矩阵Cholesky分解所需要的时间。这是因为谱矩阵的分解需要耗费很多的时间和计算内存。Yang等人和Cao等人提出了无需进行Cholesky分解的风场模拟闭合公式，但该方法只有在相干函数服从指数分布时才有效^[8-9]。研究表明采用指数模型描述风场相干函数特性可能并不合适，为此许多学者也提出了更有意义的非指数模型^[10-11]。此外，基于闭合公式的风场模拟方法只有在模拟点为等间距分布时才有效。Ding等人推导了单指标频率的模拟公式，该方法可以减小Cholesky矩阵分解的计算量，从而提高模拟效率^[12]。Huang等人通过将谱矩阵的相位分离出来以提高谱矩阵的分解效率^[13]。除了Cholesky分解以外，POD或者特征向量分解也被广泛应用于多点平稳随机过程、随机场以及复杂风场的模拟^[14-17]。综上所述，尽管研究者对多点随机风场的模拟进行了很广泛的研究，但当模拟点数很大时，仍然需要进一步提高其模拟效率。

发明内容

鉴于此，本发明提出一种可有效解决模拟点数较大时的风场模拟方法。具体涉及一种基于时空场和POD条件插值的平稳均质风场快速模拟方法，所述的模拟方法包含以下几部分内容：首先基于时空场的概念对符合转换条件的均质平稳随机风场进行了模拟，并对该模拟方法的各态历经性、模拟效率以及模拟精度进行了详细的探讨。然后在模拟出的等间距样本时程的基础上，采用POD插值方法模拟非等间距点的样本时程。采用本发明模拟均质平稳随机风场主要有以下两个优点：不需要传统方法中经常使用的Cholesky分解技术以及改进模拟方法中的相干函数需服从指数分布的假定；可以使用二维FFT技术极大的提高模拟效率。

为实现上述目的，本发明采用的技术方案为：一种基于时空场和条件插值的平稳均质风场快速模拟方法，基于时空场和条件插值的平稳均质对风场进行高效高精度较高快速模拟，其特征在于，建立多点均质平稳随机过程与多点均质平稳时空随机场之间的联系，确定两者之间转化的条件以及基于时空场模拟出的多点风场样本所具有的统计与各态历经特；采用POD技术对非等间距点的风速时程进行条件插值，以解决实际应用中可能存在的非等间距点风场模拟的问题，包括如下具体步骤：

1)获取平稳均质风场随机过程的功率谱密度矩阵

获取零均值l维均质平稳随机风场的功率谱密度矩阵S⁰(ω)，如下式所示：

$S^0\left(\mathrm{\ω}\right)=\[S_{jk}^0\left(\mathrm{\ω}\right)\],j,k=1,2,...,l---\left(1\right)$

式中：为互功率谱密度；由于该风场为均质风场，因此成立；此时，可表示为：

$S_{jk}^0\left(\mathrm{\ω}\right)=S_{jj}^0\left(\mathrm{\ω}\right){\mathrm{\γ}}_{jk}\left(\mathrm{\ω}\right)e^{{\mathrm{i\θ}}_{jk}\left(\mathrm{\ω}\right)},j\≠k---\left(2\right)$

式中：为和之间的相干函数；θ_jk(ω)为相干函数的相位角且满足如相干函数为实数，且为距离的偶函数，则上述多点平稳均质风场的模拟可转化为二维均质时空场的模拟；相应的，互功率谱密度可进一步表示为：

$S_{jk}^0\left(\mathrm{\ω}\right)=S_{jj}^0\left(\mathrm{\ω}\right){\mathrm{\γ}}_{jk}({\mathrm{\ξ}}_0,\mathrm{\ω}),j\≠k---\left(3\right)$

式中：ξ₀为和之间的距离；

2)获取转化后的均质时空场的功率谱密度

设f⁰(x,t)为转化后的二维、连续均质时空场，则f⁰(x,t)和f⁰(x+ξ,t)之间的互谱密度函数满足下式：

$S_f^0(\mathrm{\ξ},\mathrm{\ω})=S_{jj}^0\left(\mathrm{\ω}\right)\mathrm{\γ}(\mathrm{\ξ},\mathrm{\ω})---\left(4\right)$

式中：ξ为距离变量；γ(ξ,ω)为多点随机过程的相干函数表达式，基于此，连续均质时空场的二维功率谱密度函数可由下式获得：

$S_f^0(\mathrm{\κ},\mathrm{\ω})=\frac{1}{2\mathrm{\π}}{\∫}_{-\mathrm{\∞}}^{\mathrm{\∞}}{S}_f^0(\mathrm{\ξ},\mathrm{\ω})e^{-i\mathrm{\κ}\mathrm{\ξ}}d\mathrm{\ξ}---\left(5\right)$

其中：κ为相对于ξ的波数；

3)多点均质平稳风场的模拟

设[f(x₁,t),f(x₂,t),…,f(x_l,t)]^T为对应于离散点x₁,x₂,…,x_l的模拟时空场，这些离散点的模拟时程即可看成为多点均质平稳风场在离散点x₁,x₂,…,x_l的模拟时程样本[p₁(t),p₂(t),…,p_l(t)]^T；

由于相干函数为距离的实偶函数，转化后的二维均质平稳风场为象限对称的均质风场模拟公式如下式所示：

$\begin{array}{c}f(x,t)=\sqrt{2}\underset{n_1=0}{\overset{N_1-1}{\Σ}}\underset{n_2=0}{\overset{N_2-1}{\Σ}}\{\sqrt{2S_{ff}^0({\mathrm{\κ}}_{n_1},{\mathrm{\ω}}_{n_2})\mathrm{\Δ}\mathrm{\κ}\mathrm{\Δ}\mathrm{\ω}}\[\cos ({\mathrm{\κ}}_{n_1}x+{\mathrm{\ω}}_{n_2}t+{\mathrm{\Φ}}_{n_1{n}_2}^{\left(1\right)})\\ +\cos ({\mathrm{\κ}}_{n_1}x-{\mathrm{\ω}}_{n_2}t+{\mathrm{\Φ}}_{n_1{n}_2}^{\left(2\right)})\]\}\end{array}---\left(6\right)$

式中：Δκ＝κ_u/N₁，n₁＝0,1,…,N₁-1；Δω＝ω_u/N₂，n₂＝0,1,…,N₂-1；ε＝1/m，m是大于1的正整数；κ_u和ω_u分别为截止波数和截止频率；和分别为独立和均匀分布在[0,2π]的随机相位角变量；当具体的相位角和代入式(6)时，则可以得到单个模拟样本时程f⁽ⁱ⁾(x,t)；为了避免采样定理中的模态混叠，距离增量Δx和时间增量Δt必须分别满足Δx≤2π/(2κ_u)和Δt≤2π/(2ω_u)的条件；

当ε＝0时，式(6)则退化为随机场的经典公式；模拟出的样本以波长L₀＝2π/Δκ和时间T₀＝2π/Δω为一个周期；此外，为了保证模拟出的样本时程具有各态历经性，必须加入如下约束：

$\begin{array}{cc}{S}_{ff}^0(0,\mathrm{\&omega;})=S_{ff}^0(\mathrm{\&kappa;},0)=0& 0\end{array}\&le;\mathrm{\&kappa;}<\mathrm{\&infin;};0\&le;\mathrm{\&omega;}<\mathrm{\&infin;}---\left(7\right)$

4)等间距点的二维FFT高效模拟

设多点风场中前n个点的风速时程[p₁(t),p₂(t),…,p_n(t)]^T为等间距分布点的风速时程，此时风速时程的模拟可以采用二维FFT极大的提高模拟效率，，公式(6)表达为

其中：Re表示实部；对于ε＝0，p₁＝0,1,...，M₁-1，p₂＝0,1,...,M₂-1；对于ε≠0，p₁＝0,1,...,m×M₁-1，p₂＝0,1,...,m×M₂-1；q₁和q₂分别为p₁/M₁和p₂/M₂的余数；和分别表示为

$C_{q_1{q}_2}=\underset{n_1=0}{\overset{M_1-1}{\&Sigma;}}\underset{n_2=0}{\overset{M_2-1}{\&Sigma;}}\{B_{n_1{n}_2}\exp \&lsqb;i\frac{2{\mathrm{\&pi;n}}_1{q}_1}{M_1}+i\frac{2{\mathrm{\&pi;n}}_2{q}_2}{M_2}\&rsqb;\}---\left(9\right)$

${\tilde{C}}_{q_1{q}_1}\underset{n_1=0}{\overset{M_1-1}{\&Sigma;}}\underset{n_2=0}{\overset{M_2-1}{\&Sigma;}}\{{\tilde{B}}_{n_1{n}_2}\exp \&lsqb;i\frac{2{\mathrm{\&pi;n}}_1{q}_1}{M_1}-i\frac{2{\mathrm{\&pi;n}}_2{q}_2}{M_2}\&rsqb;\}---\left(10\right)$

其中和分别可由下式获得

$B_{n_1{n}_2}=2\sqrt{S_{ff}^0({\mathrm{\&kappa;}}_{n_1},{\mathrm{\&omega;}}_{n_2})\mathrm{\&Delta;}\mathrm{\&kappa;}\mathrm{\&Delta;}\mathrm{\&omega;}}\exp \&lsqb;{\mathrm{i\&phi;}}_{n_1{n}_2}^{\left(1\right)\left(i\right)}\&rsqb;---\left(11\right)$

${\tilde{B}}_{n_1{n}_2}=2\sqrt{S_{ff}^0({\mathrm{\&kappa;}}_{n_1},{\mathrm{\&omega;}}_{n_2})\mathrm{\&Delta;}\mathrm{\&kappa;}\mathrm{\&Delta;}\mathrm{\&omega;}}\exp \&lsqb;{\mathrm{i\&phi;}}_{n_1{n}_2}^{\left(2\right)\left(i\right)}\&rsqb;---\left(12\right)$

在上式中，Δx和Δt需满足ΔxΔκ＝2π/M₁和ΔtΔω＝2π/M₂的条件；为了避免模态混叠，M₁≥2N₁；M₂≥2N₂必须成立；值得一提的是，当ε＝0时，需在n₁＝0,1,...M₁-1和n₂＝0,1,...M₂-1时成立；此外，在N₁≤n₁≤M₁-1或N₂≤n₂≤M₂-1时，必须成立；

5)非等间距点的POD条件插值；

对于非等间距点的风速时程模拟，可以采用公式(6)直接进行三角函数的叠加计算，这种方法适用于非等间距模拟点数较少的情况；当模拟点数较多时，采用POD方法对非等间距点的风速时程进行条件插值，从而近似模拟非等间距点的风速时程；

设为模拟出的一部分等间距分布的多点风速时程，其中然后，可将P⁽ⁱ⁾(t)投影在一系列优化正交基上，如下式所示：

$P^{\left(i\right)}\left(t\right)=\underset{j=1}{\overset{\stackrel{\&OverBar;}{n}}{\&Sigma;}}{\mathrm{\&Phi;}}_j{a}_j\left(t\right)---\left(13\right)$

式中：a_j(t)为j阶主坐标且可由计算；优化正交基可通过求解以下特征值方程获得：

R(0)Φ_j＝λ_jΦ_j(14)

式中：R(0)为样本时程P⁽ⁱ⁾(t)的互协方差矩阵；λ_j为j阶特征值，一旦特征值按照降序排列，低阶模态将包含大部分的能量；相应的，重构的样本时程可由下式计算：

${\hat{P}}^{\left(i\right)}\left(t\right)=\underset{j=1}{\overset{q}{\&Sigma;}}{\mathrm{\&Phi;}}_j{a}_j\left(t\right),q\&le;\stackrel{\&OverBar;}{n}---\left(15\right)$

POD将会被应用于非等间距点的风速时程插值；在获取每个优化正交基以后，可对所有非等间距点的正交基进行插值，并最终得到插值后的所有正交基注意到，此时正交基的维数已包括非等间距点的个数；包含非等间距点的风速时程样本可由下式获得：

$\tilde{P}\left(t\right)=\underset{j=1}{\overset{q}{\&Sigma;}}{\widetilde{\mathrm{\&Phi;}}}_j{a}_j\left(t\right),q\&le;\stackrel{\&OverBar;}{n}---\left(16\right)$

通过上述插值，可以得出所有非等间距点的风速时程，进而可模拟出整个多点均质风场的样本时程。

从上述模拟步骤可以看出，基于时空场的模拟方法不需要传统方法中所需的Cholesky分解。因此，相干函数服从指数分布的假定可以忽略。此外，二维FFT的使用可以极大的提高风场的模拟效率。当模拟点数很多时，模拟效率提高的更加明显。

本发明的有益效果为：

采用本发明来模拟均质平稳随机风场时，不需要传统方法中经常使用的Cholesky分解以及改进模拟方法中相干函数需服从指数分布的假定。此外，还可以使用二维FFT技术极大的提高模拟效率。数值算例分析表明，该方法具有易于使用、精度较高以及模拟效率很高的特点。因此，该方法可有效解决随机风场模拟点数很大时，谱表示方法存在的模拟效率低下的问题。

附图说明

图1为实例1主梁上的风场模拟点分布图；

图2为实例1模拟出的等间距点样本时程；其中，图(a)为点1的模拟时程；(b)为点76的模拟时程；

图3为实例1估计自相关函数和互相关函数与理论值的对比图(等间距点)；其中，(a)为模拟点1的自相关函数；(b)模拟点1和76之间的互相关函数；

图4为实例1模拟出的非等间距点样本时程；其中：(a)点A的模拟时程；(b)点B的模拟时程。

图5为实例1估计自相关函数和互相关函数与理论值的对比图(非等间距点)；其中，(a)模拟点A的自相关函数；(b)模拟点A和B之间的互相关函数；

图6为实例2各种方法的定义；

图7为实例2各种方法模拟效率对比。

具体实施方式

为了使本发明的目的、技术方案和有益效果更加清楚，下面将结合附图和实施例，对本发明的实施例进行详细的说明，以方便技术人员理解。

实施例1：实际风场的模拟

1)目标风场描述

本发明采用一个千米级悬索桥的主梁风场模拟作为数值案例。该桥梁的跨度布置为390、1080、1080和390m，其桥跨布置图如图1所示。为了简便，本案例只对两个中间主跨的顺风向风场进行模拟。主梁上总共有301个等间距模拟点，因此两点之间的距离为7.2m。其中三个非等间距模拟点A，B和C分别位于模拟点75和76以及77和78之间，如图1所示。

桥梁主梁上的顺风向脉动风谱采用Kaimal风谱模型，如下式所示：

$S\left(\mathrm{\&omega;}\right)=\frac{200}{4\mathrm{\&pi;}}\frac{z}{U}\frac{u^2}{{\&lsqb;1+50\stackrel{\&OverBar;}{\mathrm{\&omega;}}\&rsqb;}^{5/3}}---\left(17\right)$

式中：z＝65.5为桥梁主梁的高度；U＝38.9m/s为主梁高度处的平均风速；u_*为摩阻速度，且定义为u_*＝kU/ln(z/z₀)。其中k≈0.4，粗糙长度z₀＝0.01。此外，相干函数模型选用Davenport指数函数模型，如下式所示：

$\mathrm{\&gamma;}(\mathrm{\&xi;},\mathrm{\&omega;})=\exp (-\frac{\mathrm{\&lambda;}\mathrm{\&xi;}\mathrm{\&omega;}}{2\mathrm{\&pi;}U})---\left(18\right)$

其中：λ＝7为衰减因子；ξ为两点之间的距离。

相应的，转化后时空场的互谱密度函数可表示为

$S(\mathrm{\&xi;},\mathrm{\&omega;})=S\left(\mathrm{\&omega;}\right)\exp (-\frac{\mathrm{\&lambda;}\mathrm{\&xi;}\mathrm{\&omega;}}{2\mathrm{\&pi;}U})---\left(19\right)$

当相位角为0时，相干函数为偶函数，因此多点均质随机过程的模拟可转化为时空随机场的模拟。根据公式(5)，可得出转化时空随机场的空间谱密度函数，如下式所示：

$S(\mathrm{\&kappa;},\mathrm{\&omega;})=S\left(\mathrm{\&omega;}\right)\frac{\mathrm{\&alpha;}\mathrm{\&omega;}}{\mathrm{\&pi;}\&lsqb;{\left(\mathrm{\&alpha;}\mathrm{\&omega;}\right)}^2+{\mathrm{\&kappa;}}^2\&rsqb;}---\left(20\right)$

式中：α＝λ/(2πU)。当ω和κ趋于0时，S(κ,ω)的值趋近于无穷大。即在原点时，S(κ,ω)的值为一个冲击函数。但是，功率谱函数在频率和波数二个轴上的二重积分为有限值18.8m²/s²，即该二维功率谱下的能量为有限的。

2)等间距点的模拟结果

在模拟多点随机过程时，截止频率选为4πrad/s，相应的时间增量Δt为0.25s。频率离散数N₂设为2048，而波数离散数N₁分别选择256，512和1024三个值进行对比。在模拟中，选用频率偏移ε＝1/2以提高模拟精度。采用50个样本去估计模拟随机过程的统计互相关函数。严格来说，模拟出的样本不服从随机过程概念下的各态历经性。但是，该样本具有随机场意义下的各态历经性。因此，基于时间平均计算方法可首先得到每个样本的估计互相关函数。然后，采用50个样本进行统计平均，以使统计互相关函数看起来比较平滑。

图2为模拟出的点1和点76的风速时程，可以看出，两者之间的相关性较弱，这是由于这两模拟点之间的距离较远。图3则为点1风速时程的估计自相关函数以及点1和点76风速时程的估计互相关函数与理论值的对比。从图中可以看出，当N₁＝512时，估计值与理论值仍然存在一定误差。但是，当N₁＝1024，模拟出的随机过程具有满意的精度。由此可见，如要获取更高的精度，可以采取频率偏移以及较大的波数离散数。当然，这会使得模拟出的样本周期变大。但是，由于本发明的模拟效率很高，由样本周期变大而造成的计算量增大并不重要。

3)非等间距点的模拟结果

如上所述，本发明采用三个非等间距点来说明POD条件插值方法。其中，点A位于点75和点76的中点。模拟点B和模拟点C为等间距分布在点77和78之间。采用非等间距点周围的66到85点的模拟时程进行条件插值，其中插值方法为线性插值。

图4为非等间距点A和B的样本时程，可以看出两个样本时程具有较强的相关性。图5则为点A风速时程的估计自相关函数以及点A和点B风速时程的估计互相关函数与理论值的对比，可以看出，两者之间吻合很好，说明采用POD进行非等间距点的风速时程插值具有较高的精度。

实施例2：模拟效率对比

为了显示本发明的计算效率，将本发明与经典的多点随机过程谱表示方法进行对比。其中，谱分解方法分别采用经典的Cholesky分解方法^[7]和闭合公式方法^[8-9]。同时，在具体模拟时，分别采用双指标频率和单指标频率方法进行模拟。综合上述因素，可以定义四种模拟方法，如图6所示。本发明可以定义为方法5，且采用ε＝1/2进行模拟。

在对比时，所有模拟方法的离散频率数都设为1024，而模拟点数则从16变化到1024，其它参数均采用前述实施案例里面的参数。所有的算法都采用MATLAB进行编程，并运行在CPU为Intel(R)Xeon(R)E5-2609v2processor(2.50GHz)，内存为32GB的64位电脑上。各种方法的计算时间如图7所示，星号表示计算时间为估计值，时间的单位为s。此外，图中比值1定义为方法2与方法5计算时间的比值，比值2定义为方法4与方法5计算时间的比值。

由图7可以看出，本发明的计算效率明显高于其它方法的计算效率。此外，随着模拟点数的增大，本发明的计算效率提高的更加明显。例如，当模拟点数为16时，比值1和2分别为62.8和6，本发明的模拟效率提高不是很明显。但是，当模拟点数为128时，比值1和2则分别为16374.8和166.2，计算效率得到明显的提高。事实上，比值2可能更具参考性，因为当模拟点数很大时，双指标频率方法模拟出的样本周期过长。

综上所述，本发明涉及一种基于时空场和条件插值的平稳均质风场快速模拟方法，该模拟方法主要包括以下几个部分：首先基于时空场的概念对符合转换条件的均质平稳随机风场进行了模拟，并对该模拟方法的各态历经性、模拟效率以及模拟精度进行了详细的探讨。然后在模拟出的等间距样本时程的基础上，采用POD插值方法模拟非等间距点的样本时程。采用本发明模拟均质平稳随机风场主要有以下两个优点：第一，不需要传统方法中经常使用的Cholesky分解以及改进模拟方法中相干函数需服从指数分布的假定；第二，可以使用二维FFT技术极大的提高模拟效率。数值案例分析结果表明，与其它方法相比，本发明的模拟效率得到了极大的提高。此外，该方法的模拟精度也比较满意。因此，本发明可有效解决当模拟点数很大时，平稳均质风场的谱表示模拟方法存在的效率低下问题。

需要说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (1)

1.一种基于时空场和条件插值的平稳均质风场快速模拟方法，基于时空场和条件插值的平稳均质对风场进行高效高精度较高快速模拟，其特征在于，建立多点均质平稳随机过程与多点均质平稳时空随机场之间的联系，确定两者之间转化的条件以及基于时空场模拟出的多点风场样本所具有的统计与各态历经特；采用POD技术对非等间距点的风速时程进行条件插值，以解决实际应用中可能存在的非等间距点风场模拟的问题，包括如下具体步骤：

1)获取平稳均质风场随机过程的功率谱密度矩阵

获取零均值l维均质平稳随机风场的功率谱密度矩阵S⁰(ω)，如下式所示：

式中：为互功率谱密度；由于该风场为均质风场，因此成立；此时，可表示为：

式中：为和之间的相干函数；θ_jk(ω)为相干函数的相位角且满足如相干函数为实数，且为距离的偶函数，则上述多点平稳均质风场的模拟可转化为二维均质时空场的模拟；相应的，互功率谱密度可进一步表示为：

式中：ξ₀为和之间的距离；

2)获取转化后的均质时空场的功率谱密度

设f⁰(x,t)为转化后的二维、连续均质时空场，则f⁰(x,t)和f⁰(x+ξ,t)之间的互谱密度函数满足下式：

式中：ξ为距离变量；γ(ξ,ω)为多点随机过程的相干函数表达式，基于此，连续均质时空场的二维功率谱密度函数可由下式获得：

其中：κ为相对于ξ的波数；

3)多点均质平稳风场的模拟

设[f(x₁,t),f(x₂,t),…,f(x_l,t)]^T为对应于离散点x₁,x₂,…,x_l的模拟时空场，这些离散点的模拟时程即可看成为多点均质平稳风场在离散点x₁,x₂,…,x_l的模拟时程样本[p₁(t),p₂(t),…,p_l(t)]^T；

由于相干函数为距离的实偶函数，转化后的二维均质平稳风场为象限对称的均质风场模拟公式如下式所示：

式中：Δκ＝κ_u/N₁，n₁＝0,1,…,N₁-1；Δω＝ω_u/N₂，n₂＝0,1,…,N₂-1；ε＝1/m，m是大于1的正整数；κ_u和ω_u分别为截止波数和截止频率；和分别为独立和均匀分布在[0,2π]的随机相位角变量；当具体的相位角和代入式(6)时，则可以得到单个模拟样本时程f⁽ⁱ⁾(x,t)；为了避免采样定理中的模态混叠，距离增量Δx和时间增量Δt必须分别满足Δx≤2π/(2κ_u和Δt≤2π/(2ω_u)的条件；

当ε＝0时，式(6)则退化为随机场的经典公式；模拟出的样本以波长L₀＝2π/Δκ和时间T₀＝2π/Δω为一个周期；此外，为了保证模拟出的样本时程具有各态历经性，必须加入如下约束：

4)等间距点的二维FFT高效模拟

设多点风场中前n个点的风速时程[p₁(t),p₂(t),…,p_n(t)]^T为等间距分布点的风速时程，此时风速时程的模拟可以采用二维FFT极大的提高模拟效率，，公式(6)表达为

其中：Re表示实部；对于ε＝0，p₁＝0,1,...,M₁-1，p₂＝0,1,...,M₂-1；对于ε≠0，p₁＝0,1,...,m×M₁-1，p₂＝0,1,...,m×M₂-1；q₁和q₂分别为p₁/M₁和p₂/M₂的余数；和分别表示为

其中和分别可由下式获得

在上式中，Δx和Δt需满足ΔxΔκ＝2π/M₁和ΔtΔω＝2π/M₂的条件；为了避免模态混叠，M₁≥2N₁；M₂≥2N₂必须成立；值得一提的是，当ε＝0时，需在n₁＝0,1,...M₁-1和n₂＝0,1,...M₂-1时成立；此外，在N₁≤n₁≤M₁-1或N₂≤n₂≤M₂-1时，必须成立；

5)非等间距点的POD条件插值；

对于非等间距点的风速时程模拟，可以采用公式(6)直接进行三角函数的叠加计算，这种方法适用于非等间距模拟点数较少的情况；当模拟点数较多时，采用POD方法对非等间距点的风速时程进行条件插值，从而近似模拟非等间距点的风速时程；

设为模拟出的一部分等间距分布的多点风速时程，其中然后，可将P⁽ⁱ⁾(t)投影在一系列优化正交基Φ＝[Φ₁,Φ₂,…,Φ_n]上，如下式所示：

式中：a_j(t)为j阶主坐标且可由计算；优化正交基可通过求解以下特征值方程获得：

R(0)Φ_j＝λ_jΦ_j(14)

式中：R(0)为样本时程P⁽ⁱ⁾(t)的互协方差矩阵；λ_j为j阶特征值，一旦特征值按照降序排列，低阶模态将包含大部分的能量；相应的，重构的样本时程可由下式计算：

POD将会被应用于非等间距点的风速时程插值；在获取每个优化正交基Φ_j＝[φ_1j,φ_2j,…,φ_nj]^T以后，可对所有非等间距点的正交基进行插值，并最终得到插值后的所有正交基注意到，此时正交基的维数已包括非等间距点的个数；包含非等间距点的风速时程样本可由下式获得：

通过上述插值，可以得出所有非等间距点的风速时程，进而可模拟出整个多点均质风场的样本时程。