CN105760585A - 一种基于高阶摄动理论的模糊声压快速求解方法 - Google Patents

Abstract

本发明公开了一种基于高阶摄动理论的模糊声压快速求解方法，步骤如下：结构振动声辐射问题的有限元建模；考虑模型的模糊不确定变量，建立声辐射问题的模糊有限元方程；选取截集水平，利用截集运算，将模糊变量转化为区间变量；利用子区间分解技术，将不确定度超过5％的区间变量分解为多个子区间，进而得到一组子区间有限元方程；基于高阶摄动理论对所有的子区间有限元方程进行求解，得到声压的区间变化范围；重组所有截集水平下的声压区间，得到模糊声压的隶属度函数。本发明可系统化解决含有模糊不确定变量的构振动声辐射问题，在保证计算效率的同时，进一步提高了摄动法的计算精度，这是一般商用软件所不能实现的。

Description

一种基于高阶摄动理论的模糊声压快速求解方法

技术领域

本发明属于机械工程领域，具体涉及一种基于高阶摄动理论的模糊声压快速求解方法。

背景技术

结构振动声辐射是噪声的主要来源之一，是工程设计时需要考虑的重要方面，通常希望降低有害振动和噪声的幅值以提升系统的安全性和舒适性。对于实际工程当中的许多薄壁结构，除了考虑内部噪声由结构振动引起外，还通常考虑内部声压对结构的反作用。这类考虑结构和声场耦合作用的振动系统，在汽车、船舶等工业领域中普遍存在。随着计算机软硬件的发展，以有限元技术为代表的数值计算方法在结构振动声辐射研究中的作用和地位不断提高，并在工程领域内发挥着越来越大的作用。

然而，实际工程结构的材料属性、外部载荷和边界条件等模型输入参数不可避免的受到多种不确定因素的影响。对于复杂系统而言，即使很小的不确定输入，对最终的系统响应也可能会产生明显的扰动。以概率理论为基础的随机分析方法在不确定声学研究中已经取得了许多成果，但概率模型的建立需要大量的样本信息来事先确定参数的概率密度函数。而在数值分析的初始阶段，获得足够的样本数据往往花费较大或代价过高。人们发现，工程问题除了存在随机现象外，还存在大量的模糊现象。某些事物的概念或参数的数值是难以确定的，只能根据实验数据或主观经验确定一个大致的范围。如此一来，模糊模型在不确定建模方面表现出了很强的方便性和经济性。模糊理论在结构静力学分析中已取得了一些研究成果，但在声学领域中的应用才刚刚起步。另外，随着工程领域中对计算精度要求的不断提高，传统的模糊摄动理论所导致的计算精度低下问题越来越严重。因此，如何改进摄动方法以快速准确预测系统响应的模糊不确定特征，是目前学术领域的一个研究热点，对于弥补现有结构振动声辐射分析方法的不足，具有重要的工程应用价值。

发明内容

本发明所要解决的技术问题为：克服现有技术在结构振动声辐射问题求解中存在的不足，充分考虑声辐射问题中的模糊不确定因素，基于子区间分解技术和高阶摄动理论，提出了一种快速有效预测模糊声压隶属度函数的数值计算方法，可系统化解决含有模糊不确定变量的构振动声辐射问题，在保证计算效率的同时，进一步提高了摄动法的计算精度。

：本发明为解决上述技术问题采用的技术方案为：一种基于高阶摄动理论的模糊声压快速求解方法，包括以下步骤：

步骤一：针对结构振动声辐射问题，选取合适的单元类型对结构和声场进行离散，计及结构和声场的耦合作用，建立声辐射有限元离散方程；

步骤二：充分考虑声辐射问题中的不确定因素，并用模糊变量对其进行定量化表示，根据步骤一中的离散方程可以建立声辐射问题的模糊有限元方程；

步骤三：选取截集水平，利用截集运算，可以将步骤二中的模糊变量转化为区间变量；

步骤四：利用子区间分解技术，将步骤三中得到的不确定度超过5％的区间变量分解为多个子区间，使每个子区间的不确定度小于5％，进而将步骤二中的模糊有限元方程转化为一组子区间有限元方程；

步骤五：基于高阶摄动理论对步骤四中的所有子区间有限元方程进行求解，利用区间组合运算得到声压的区间变化范围，进而可得到所有截集水平对应的声压区间；

步骤六：将步骤五中得到的所有截集水平下的声压区间进行重组，最终得到模糊声压的隶属度函数。

其中，所述步骤三中截集水平的选取并不是固定不变的；根据模糊不确定变量隶属度函数的分布类型来确定所需截集水平的数量规模和数值大小。

其中，所述步骤四中子区间分解策略并不是固定不变的；根据每个区间变量不确定度的大小以及计算精度的要求灵活选取分解形式，只要保证每个子区间的不确定度小于5％即可。

上述各步骤具体包括以下过程：

步骤一：针对结构振动声辐射问题，选取合适的单元类型对结构和声场进行离散；除了考虑结构对声场的作用外，还充分计及声场对薄壁结构的反作用。在频域分析中，可以建立此声辐射问题的耦合有限元离散方程：

KU＝F

其中 $K=\left(\begin{array}{cc}{K}_s-{\mathrm{\ω}}^2{M}_s& -C\\ {\mathrm{\ρ}}_a{\mathrm{\ω}}^2{C}^T& K_a-{\mathrm{\ω}}^2{M}_a\end{array}\right)$ 为耦合系数矩阵，F＝(F_bF_q)^T为耦合载荷向量，U＝(up)^T为耦合系统响应向量；M_s和K_s表示结构质量矩阵和刚度矩阵，M_a和K_a表示声场质量矩阵和刚度矩阵；F_b为结构体载荷向量，F_q为声场附加载荷向量；u为结构的位移向量，p为声压向量；ω表示频率，C为此声辐射系统的耦合矩阵，ρ_a表示声场介质的密度。

步骤二：由于结构材料制造工艺的限制以及测量的误差，结构振动声辐射问题中的不确定因素不可避免。用n个模糊变量β₁,β₂,...,β_n来表示所有的不确定输入参数，并记为向量的形式β＝(β_i)_n＝(β₁,β₂,...,β_n)，于是根据步骤一中的离散方程可以进一步建立声辐射问题的模糊有限元方程：

K(β)U(β)＝F(β)

步骤三：在0到1范围内选取截集水平λ，利用截集运算可以将步骤二中的模糊变量β_i转化为区间变量

$\begin{array}{cc}{\mathrm{\β}}_{i,\mathrm{\λ}}^I=\[{\underset{\‾}{\mathrm{\β}}}_{i,\mathrm{\λ}},{\stackrel{\‾}{\mathrm{\β}}}_{i,\mathrm{\λ}}\]={\mathrm{\β}}_{i,\mathrm{\λ}}^c+{\mathrm{\Δ\β}}_{i,\mathrm{\λ}}^I={\mathrm{\β}}_{i,\mathrm{\λ}}^c+{\mathrm{\Δ\β}}_{i,\mathrm{\λ}}{\mathrm{\δ}}_i^I& i=1,2,...,n\end{array}$

其中和为区间变量的下界和上界，和表示区间中点和半径，为标准区间变量 表示区间变量在其中点处的波动变化。进一步，可以将λ截集水平下的所有区间变量记为向量的形式：

${\mathrm{\β}}_{\mathrm{\λ}}^I={\left({\mathrm{\β}}_{i,\mathrm{\λ}}^I\right)}_n=({\mathrm{\β}}_{1,\mathrm{\λ}}^I,{\mathrm{\β}}_{2,\mathrm{\λ}}^I,...,{\mathrm{\β}}_{n,\mathrm{\λ}}^I)$

步骤四：引入符号表示步骤三中所得区间变量的不确定度。当第k个区间变量的不确定度γ_k,λ>5％时，利用子区间分解技术，可以将此区间变量分解为N_k个子区间，使每个子区间的不确定度小于5％：

$\begin{array}{cc}{\left({\mathrm{\β}}_{k,\mathrm{\λ}}^I\right)}_j=\[{\underset{\‾}{\mathrm{\β}}}_{k,\mathrm{\λ}}+\frac{2(j-1){\mathrm{\Δ\β}}_{k,\mathrm{\λ}}}{N_k},{\underset{\‾}{\mathrm{\β}}}_{k,\mathrm{\λ}}+\frac{2{\mathrm{j\Δ\β}}_{k,\mathrm{\λ}}}{N_k}\]& j=1,2,...,N_k\end{array}$

其中表示区间变量的第j个子区间。

对所有区间变量的子区间进行任意组合，当不确定参数相互独立时，总共可以得到N＝N₁N₂…N_n个子区间组合。为了表述方便，引入符号统一表示所有的子区间向量，如此一来，可将步骤二中的模糊有限元方程转化为一组子区间有限元方程：

$K\left({\mathrm{\α}}_{\mathrm{\λ}}^I\right)U_{\mathrm{\λ}}^I\left({\mathrm{\α}}_{\mathrm{\λ}}^I\right)=F\left({\mathrm{\α}}_{\mathrm{\λ}}^I\right)$

其中为系统响应向量U在子区间向量影响下的区间表示。

步骤五：基于高阶摄动理论对步骤四中的所有子区间有限元方程进行求解，利用区间组合运算得到声压的区间变化范围，进而可得到所有截集水平对应的声压区间。首先，基于一阶泰勒展式，将子区间有限元方程中的系数矩阵和载荷向量表示为：

$K\left({\mathrm{\α}}_{\mathrm{\λ}}^I\right)=K\left({\mathrm{\α}}_{\mathrm{\λ}}^c\right)+\underset{i=1}{\overset{n}{\Σ}}\frac{\∂K\left({\mathrm{\α}}_{\mathrm{\λ}}^c\right)}{\∂{\mathrm{\α}}_{i,\mathrm{\λ}}}{\mathrm{\Δ\α}}_{i,\mathrm{\λ}}{\mathrm{\δ}}_i^I=K_{\mathrm{\λ}}^c+{\mathrm{\ΔK}}_{\mathrm{\λ}}^I$

$F\left({\mathrm{\α}}_{\mathrm{\λ}}^I\right)=F\left({\mathrm{\α}}_{\mathrm{\λ}}^c\right)+\underset{i=1}{\overset{n}{\Σ}}\frac{\∂F\left({\mathrm{\α}}_{\mathrm{\λ}}^c\right)}{\∂{\mathrm{\α}}_{i,\mathrm{\λ}}}{\mathrm{\α}}_{i,\mathrm{\λ}}{\mathrm{\δ}}_i^I=F_{\mathrm{\λ}}^c+{\mathrm{\ΔF}}_{\mathrm{\λ}}^I$

其中：

$\begin{array}{cccc}{K}_{\mathrm{\λ}}^c=K\left({\mathrm{\α}}_{\mathrm{\λ}}^c\right)& {\mathrm{\ΔK}}_{\mathrm{\λ}}^I=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{\∂K\left({\mathrm{\α}}_{\mathrm{\λ}}^c\right)}{\∂{\mathrm{\α}}_{i,\mathrm{\λ}}}{\mathrm{\Δ\α}}_{i,\mathrm{\λ}}{\mathrm{\δ}}_i^I& F_{\mathrm{\λ}}^c=F\left({\mathrm{\α}}_{\mathrm{\λ}}^c\right)& {\mathrm{\ΔF}}_{\mathrm{\λ}}^I=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{\∂F\left({\mathrm{\α}}_{\mathrm{\λ}}^c\right)}{\∂{\mathrm{\α}}_{i,\mathrm{\λ}}}{\mathrm{\Δ\α}}_{i,\mathrm{\λ}}{\mathrm{\δ}}_i^I\end{array}$

利用区间矩阵求逆运算，系统子区间响应可以表示为：

$U_{\mathrm{\λ}}^I\left({\mathrm{\α}}_{\mathrm{\λ}}^I\right)=U_{\mathrm{\λ}}^c+{\mathrm{\ΔU}}_{\mathrm{\λ}}^I={(K_{\mathrm{\λ}}^c+{\mathrm{\ΔK}}_{\mathrm{\λ}}^I)}^{-1}(F_{\mathrm{\λ}}^c+{\mathrm{\ΔF}}_{\mathrm{\λ}}^I)$

其中为子区间响应的中值，表示系统响应在中值处的波动变化。

利用纽曼级数，区间矩阵的逆可以表示为：

${(K_{\mathrm{\λ}}^c+{\mathrm{\ΔK}}_{\mathrm{\λ}}^I)}^{-1}={\left(K_{\mathrm{\λ}}^c\right)}^{-1}+{\left(K_{\mathrm{\λ}}^c\right)}^{-1}\underset{r=1}{\overset{\mathrm{\∞}}{\Σ}}{(-{\mathrm{\ΔK}}_{\mathrm{\λ}}^I{\left(K_{\mathrm{\λ}}^c\right)}^{-1})}^r$

基于泰勒展式，上式中的项可以进一步改写为：

$\underset{r=1}{\overset{\mathrm{\∞}}{\mathrm{\Σ}}}{\left(-{\mathrm{\ΔK}}_{\mathrm{\λ}}^I{\left(K_{\mathrm{\λ}}^c\right)}^{-1}\right)}^r=\underset{r=1}{\overset{\mathrm{\∞}}{\mathrm{\Σ}}}{\left(-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{\∂K\left({\mathrm{\α}}_{\mathrm{\λ}}^c\right)}{\∂{\mathrm{\α}}_{i,\mathrm{\λ}}}{\mathrm{\Δ\α}}_{i,\mathrm{\λ}}{\mathrm{\δ}}_i^I{\left(K_{\mathrm{\λ}}^c\right)}^{-1}\right)}^r=\underset{r=1}{\overset{\mathrm{\∞}}{\mathrm{\Σ}}}{\left(-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\δ}}_i^I{K}_{i,\mathrm{\λ}}\right)}^r$

其中 $K_{i,\mathrm{\λ}}=\frac{\∂K\left({\mathrm{\α}}_{\mathrm{\λ}}^c\right)}{\∂{\mathrm{\α}}_{i,\mathrm{\λ}}}{\mathrm{\Δ\α}}_{i,\mathrm{\λ}}{\left(K_{\mathrm{\λ}}^c\right)}^{-1}$

当r取不同值时，上式中的求和项可具体表示为：

${\left(-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\δ}}_i^I{K}_{i,\mathrm{\λ}}\right)}^{r=1}=-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\δ}}_i^I{K}_{i,\mathrm{\λ}}$

${\left(-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\δ}}_i^I{K}_{i,\mathrm{\λ}}\right)}^{r=2}=\underset{i,j=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\δ}}_i^I{\mathrm{\δ}}_j^I{K}_{i,\mathrm{\λ}}{K}_{j,\mathrm{\λ}}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\left({\mathrm{\δ}}_i^I{K}_{i,\mathrm{\λ}}\right)}^2+\underset{i\≠j}{\overset{n}{\underset{i,j=1}{\mathrm{\Σ}}}}{\mathrm{\δ}}_i^I{\mathrm{\δ}}_j^I{K}_{i,\mathrm{\λ}}{K}_{j,\mathrm{\λ}}$

......

其中记号表示除i＝j＝k外的所有组合情况。

如此一来，区间矩阵逆的纽曼级数就可改写为：

当K_i,λ的谱范数小于1时，利用矩阵级数的收敛性，同时忽略交叉项，可得：

$\begin{array}{c}.{\left(K_{\mathrm{\λ}}^c+{\mathrm{\ΔK}}_{\mathrm{\λ}}^I\right)}^{-1}\≈{\left(K_{\mathrm{\λ}}^c\right)}^{-1}+{\left(K_{\mathrm{\λ}}^c\right)}^{-1}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{r=1}{\overset{\mathrm{\∞}}{\mathrm{\Σ}}}{\left(-{\mathrm{\δ}}_i^I{K}_{i,\mathrm{\λ}}\right)}^r\\ ={\left(K_{\mathrm{\λ}}^c\right)}^{-1}+{\left(K_{\mathrm{\λ}}^c\right)}^{-1}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{-{\mathrm{\δ}}_i^I{K}_{i,\mathrm{\λ}}}{I+{\mathrm{\δ}}_i^I{K}_{i,\mathrm{\λ}}}={\left(K_{\mathrm{\λ}}^c\right)}^{-1}+{\left(K_{\mathrm{\λ}}^c\right)}^{-1}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{E}_{i,\mathrm{\λ}}^I\end{array}$

其中：

$E_{i,\mathrm{\λ}}^I=\frac{-{\mathrm{\δ}}_i^I{K}_{i,\mathrm{\λ}}}{I+{\mathrm{\δ}}_i^I{K}_{i,\mathrm{\λ}}}=-I+\frac{I}{I+{\mathrm{\δ}}_i^I{K}_{i,\mathrm{\λ}}}=E_{i,\mathrm{\λ}}^c+{\mathrm{\ΔE}}_{i,\mathrm{\λ}}^I=E_{i,\mathrm{\λ}}^c+{\mathrm{\ΔE}}_{i,\mathrm{\λ}}{\mathrm{\δ}}_{e_i}^I$

$\begin{array}{cc}{E}_{i,\mathrm{\λ}}^c=\frac{1}{2}(\frac{K_{i,\mathrm{\λ}}}{I-K_{i,\mathrm{\λ}}}-\frac{K_{i,\mathrm{\λ}}}{I+K_{i,\mathrm{\λ}}})& {\mathrm{\ΔE}}_{i,\mathrm{\λ}}=\frac{1}{2}|\frac{K_{i,\mathrm{\λ}}}{I-K_{i,\mathrm{\λ}}}+\frac{K_{i,\mathrm{\λ}}}{I+K_{i,\mathrm{\λ}}}|\end{array}$

I为单位矩阵，为区间矩阵对应的标准区间变量。

将上述含有部分高阶项的纽曼级数代回到子区间响应表达式中，利用摄动理论可得：

$U_{\mathrm{\λ}}^c={\left(K_{\mathrm{\λ}}^c\right)}^{-1}{F}_{\mathrm{\λ}}^c+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\left(K_{\mathrm{\λ}}^c\right)}^{-1}{E}_{i,\mathrm{\λ}}^c{F}_{\mathrm{\λ}}^c={\left(K_{\mathrm{\λ}}^c\right)}^{-1}\left(I+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{E}_{i,\mathrm{\λ}}^c\right)F_{\mathrm{\λ}}^c$

$\begin{array}{c}{\mathrm{\ΔU}}_{\mathrm{\λ}}^I={\left(K_{\mathrm{\λ}}^c\right)}^{-1}\left(I+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{E}_{i,\mathrm{\λ}}^c\right){\mathrm{\ΔF}}_{\mathrm{\λ}}^I+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\left(K_{\mathrm{\λ}}^c\right)}^{-1}{\mathrm{\ΔE}}_{i,\mathrm{\λ}}^I{F}_{\mathrm{\λ}}^c\\ =\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\left(K_{\mathrm{\λ}}^c\right)}^{-1}\left(I+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{E}_{i,\mathrm{\λ}}^c\right)\frac{\∂F\left({\mathrm{\α}}_{\mathrm{\λ}}^c\right)}{\∂{\mathrm{\α}}_{j,\mathrm{\λ}}}{\mathrm{\Δ\α}}_{j,\mathrm{\λ}}{\mathrm{\δ}}_j^I+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\left(K_{\mathrm{\λ}}^c\right)}^{-1}{\mathrm{\ΔE}}_{i,\mathrm{\λ}}{F}_{\mathrm{\λ}}^c{\mathrm{\δ}}_{e_i}^I\end{array}$

利用关于标准区间变量和的单调性，可以快速得到子区间响应的半径ΔU_λ：

${\mathrm{\ΔU}}_{\mathrm{\λ}}=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\left|{\left(K_{\mathrm{\λ}}^c\right)}^{-1}\left(I+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{E}_{i,\mathrm{\λ}}^c\right)\frac{\∂F\left({\mathrm{\α}}_{\mathrm{\λ}}^c\right)}{\∂{\mathrm{\α}}_{j,\mathrm{\λ}}}{\mathrm{\Δ\α}}_{j,\mathrm{\λ}}\right|+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left|{\left(K_{\mathrm{\λ}}^c\right)}^{-1}{\mathrm{\ΔE}}_{i,\mathrm{\λ}}{F}_{\mathrm{\λ}}^c\right|$

利用区间表示方法和区间组合运算法则，可以先后得到子区间向量和区间向量影响下的系统响应变化范围：

$U_{\mathrm{\λ}}^I\left({\mathrm{\α}}_{\mathrm{\λ}}^I\right)=\[{\underset{\‾}{U}}_{\mathrm{\λ}}\left({\mathrm{\α}}_{\mathrm{\λ}}^I\right),{\stackrel{\‾}{U}}_{\mathrm{\λ}}\left({\mathrm{\α}}_{\mathrm{\λ}}^I\right)\]=\[U_{\mathrm{\λ}}^c-{\mathrm{\ΔU}}_{\mathrm{\λ}},U_{\mathrm{\λ}}^c+{\mathrm{\ΔU}}_{\mathrm{\λ}}\]$

$U_{\mathrm{\λ}}^I\left({\mathrm{\β}}_{\mathrm{\λ}}^I\right)=\underset{i=1,2,...,N}{\∪}{U}_{\mathrm{\λ}}^I\left({\mathrm{\α}}_{i,\mathrm{\λ}}^I\right)=\[\underset{i=1,2,\mathrm{...},N}{min}{\underset{\‾}{U}}_{\mathrm{\λ}}\left({\mathrm{\α}}_{i,\mathrm{\λ}}^I\right),\underset{i=1,2,...,N}{max}{\stackrel{\‾}{U}}_{\mathrm{\λ}}\left({\mathrm{\α}}_{i,\mathrm{\λ}}^I\right)\]$

其中表示由子区间任意组合所得到的N个子区间参数向量。

根据表达式U＝(up)^T从向量中提取声压项，就可以得到λ截集水平下的声压区间变化范围

对选定的所有截集水平重复上述操作，进而可得到所有截集水平对应的声压区间。

步骤六：将步骤五中得到的所有截集水平下的声压区间进行重组，利用模糊分解定理可以得到模糊参数向量β影响下模糊声压的隶属度函数p(β)：

$p\left(\mathrm{\β}\right)=\underset{j=1,...,m}{\∪}\left\{{\mathrm{\λ}}_j{p}_{{\mathrm{\λ}}_j}^I\right({\mathrm{\β}}_{{\mathrm{\λ}}_j}^I\left)\right\}$

其中m为所选取的截集水平的数量。

本发明与现有技术相比的优点在于：

(1)与传统的结构振动声辐射分析方法相比，所提出的计算方法充分考虑实际工程中材料属性和外部载荷的模糊不确定性，计算结果对声场分析具有更重要的指导意义。

(2)在参数的预处理过程中，采用了子区间分解技术，可有效降低区间变量的不确定度，从而保证了泰勒展式和摄动法的计算精度。

(3)通过矩阵级数的收敛条件，保留了纽曼级数中的部分高阶项，有效提高了区间矩阵求逆运算的精度。

(4)本发明操作简单，实施方便，在保证计算效率的基础上提高了计算精度。

附图说明

图1为本发明的一种基于高阶摄动理论的模糊声压快速求解方法流程图；

图2为本发明的汽车结构振动声辐射问题的有限元模型示意图；

图3为300Hz频率下观测点处模糊声压的隶属度函数示意图；

图4为700Hz频率下观测点处模糊声压的隶属度函数示意图。

具体实施方式

下面结合附图和实施例对本发明做进一步说明。

本发明适用于含有模糊不确定参数的结构振动声辐射问题。本发明实施方式以某汽车结构内部舱室的声场分析为例，具体说明所述的一种基于高阶摄动理论的模糊声压快速求解方法。另外，此汽车结构内部舱室的声场数值计算方法可以推广到其他含有模糊参数的结构振动声辐射问题的声场预测中。

一种基于高阶摄动理论的模糊声压快速求解方法的计算过程如图1所示，考虑结构和声场的耦合作用，对结构振动声辐射问题进行有限元建模，引入模糊不确定变量，建立模糊有限元方程，利用子区间技术，对某些截集水平下的区间变量做精细化处理，基于高阶摄动理论对所有的子区间有限元方程进行求解，得到声压的区间变化范围，并利用模糊分解定理对其进行重组，最终得到模糊声压的隶属度函数。可分为如下几个步骤进行：

步骤一：建立汽车结构和内部声场的有限元模型，如图2所示：汽车结构如车顶2、前窗玻璃3、仪器板4、车身7、后窗玻璃9用二维四边形壳单元来离散，座椅6用三维六面体固体单元来离散，内部声场8用三维六面体流体单元来离散。车顶中心处1施加有幅值为F₀的简谐激励。在驾驶员耳朵所在位置提取一个节点5，作为舱室声压的观测点。除了考虑结构对声场的作用外，还充分计及声场对薄壁结构的反作用。在频域分析中，建立汽车结构振动声辐射问题的耦合有限元离散方程：

KU＝F

其中 $K=\left(\begin{array}{cc}{K}_s-{\mathrm{\ω}}^2{M}_s& -C\\ {\mathrm{\ρ}}_a{\mathrm{\ω}}^2{C}^T& K_a-{\mathrm{\ω}}^2{M}_a\end{array}\right)$ 为耦合系数矩阵，F＝(F_bF_q)^T为耦合载荷向量，U＝(up)^T为耦合系统响应向量；M_s和K_s表示结构质量矩阵和刚度矩阵，M_a和K_a表示声场质量矩阵和刚度矩阵；F_b为结构体载荷向量，F_q为声场附加载荷向量；u为结构的位移向量，p为声压向量；ω表示频率，C为此声辐射系统的耦合矩阵，ρ_a表示声场介质的密度。

步骤二：汽车前后窗采用密度为ρ₁，弹性模量为E₁的玻璃材料；车顶和车身结构采用密度为ρ₂，弹性模量为E₂的金属材料；仪器板和座椅采用密度为ρ₃，弹性模量为E₃的泡沫材料；舱室内空气密度记为ρ_a，声音在空气中的传播速度记为c；汽车车顶承受幅值为F₀的简谐激励。由于材料制造工艺的限制以及测量的误差，所有模型输入参数均含有一定的模糊不确定性，且隶属度函数满足三角分布规律，即ρ₁＝(2500,3000,3500)kg/m³，E₁＝(55,70,85)GPa，ρ₂＝(7000,8500,10000)kg/m³，E₂＝(160,200,240)GPa，ρ₃＝(0.9,1.1,1.3)kg/m³，E₃＝(0.0024,0.003,0.0036)GPa，ρ_a＝(1.00,1.225,1.45)kg/m³，c＝(270,340,410)m/s，F₀＝(0.8,1.0,1.2)N。将所有模糊变量用向量符号统一表示β＝(β_i)₉＝(ρ₁,ρ₂,ρ₃,E₁,E₂,E₃,ρ_a,c,F₀)，根据步骤一中的离散方程可以进一步建立声辐射问题的模糊有限元方程：

K(β)U(β)＝F(β)

步骤三：在0到1范围内选取11个截集水平λ_j＝(j-1)×0.1j＝1,...,11，为了表示方便，将其统一记为λ。利用截集运算可以将步骤二中的模糊变量β_i转化为区间变量

$\begin{array}{cc}{\mathrm{\β}}_{i,\mathrm{\λ}}^I=\[{\underset{\‾}{\mathrm{\β}}}_{i,\mathrm{\λ}},{\stackrel{\‾}{\mathrm{\β}}}_{i,\mathrm{\λ}}\]={\mathrm{\β}}_{i,\mathrm{\λ}}^c+{\mathrm{\Δ\β}}_{i,\mathrm{\λ}}^I={\mathrm{\β}}_{i,\mathrm{\λ}}^c+{\mathrm{\Δ\β}}_{i,\mathrm{\λ}}{\mathrm{\δ}}_i^I& i=1,2,...,9\end{array}$

其中和为区间变量的下界和上界，和表示区间中点和半径，为标准区间变量 表示区间变量在其中点处的波动变化。进一步，可以将λ截集水平下的所有区间变量记为向量的形式：

${\mathrm{\β}}_{\mathrm{\λ}}^I={\left({\mathrm{\β}}_{i,\mathrm{\λ}}^I\right)}_9=({\left({\mathrm{\ρ}}_1\right)}_{\mathrm{\λ}}^I,{\left({\mathrm{\ρ}}_2\right)}_{\mathrm{\λ}}^I,{\left({\mathrm{\ρ}}_3\right)}_{\mathrm{\λ}}^I,{\left(E_1\right)}_{\mathrm{\λ}}^I,{\left(E_2\right)}_{\mathrm{\λ}}^I,{\left(E_3\right)}_{\mathrm{\λ}}^I,{\left({\mathrm{\ρ}}_a\right)}_{\mathrm{\λ}}^I,c_{\mathrm{\λ}}^I,{\left(F_0\right)}_{\mathrm{\λ}}^1)$

步骤四：引入符号表示步骤三中所得区间变量的不确定度。当第k个区间变量的不确定度γ_k,λ>5％时，利用子区间分解技术，可以将此区间变量分解为N_k个子区间，使每个子区间的不确定度小于5％。这里以截集水平λ₈＝0.7和λ₆＝0.5为例，分别用两个子区间和四个子区间对所有区间变量进行划分，如表1所示。

表1不同截集水平下区间变量的子区间划分

引入符号统一表示所有的子区间向量，如此一来，可将步骤二中的模糊有限元方程转化为一组子区间有限元方程：

$K\left({\mathrm{\α}}_{\mathrm{\λ}}^I\right)U_{\mathrm{\λ}}^I\left({\mathrm{\α}}_{\mathrm{\λ}}^I\right)=F\left({\mathrm{\α}}_{\mathrm{\λ}}^I\right)$

其中为系统响应向量U在子区间向量影响下的区间表示。

步骤五：基于高阶摄动理论对步骤四中的所有子区间有限元方程进行求解，利用区间组合运算得到声压的区间变化范围，进而可得到所有截集水平对应的声压区间。首先，基于一阶泰勒展式，将子区间有限元方程中的系数矩阵和载荷向量表示为：

$K\left({\mathrm{\α}}_{\mathrm{\λ}}^I\right)=K\left({\mathrm{\α}}_{\mathrm{\λ}}^c\right)+\underset{i=1}{\overset{n}{\Σ}}\frac{\∂K\left({\mathrm{\α}}_{\mathrm{\λ}}^c\right)}{\∂{\mathrm{\α}}_{i,\mathrm{\λ}}}{\mathrm{\Δ\α}}_{i,\mathrm{\λ}}{\mathrm{\δ}}_i^I=K_{\mathrm{\λ}}^c+{\mathrm{\ΔK}}_{\mathrm{\λ}}^I$

$F\left({\mathrm{\α}}_{\mathrm{\λ}}^I\right)=F\left({\mathrm{\α}}_{\mathrm{\λ}}^c\right)+\underset{i=1}{\overset{n}{\Σ}}\frac{\∂F\left({\mathrm{\α}}_{\mathrm{\λ}}^c\right)}{\∂{\mathrm{\α}}_{i,\mathrm{\λ}}}{\mathrm{\α}}_{i,\mathrm{\λ}}{\mathrm{\δ}}_i^I=F_{\mathrm{\λ}}^c+{\mathrm{\ΔF}}_{\mathrm{\λ}}^I$

其中：

$\begin{array}{cccc}{K}_{\mathrm{\λ}}^c=K\left({\mathrm{\α}}_{\mathrm{\λ}}^c\right)& {\mathrm{\ΔK}}_{\mathrm{\λ}}^I=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{\∂K\left({\mathrm{\α}}_{\mathrm{\λ}}^c\right)}{\∂{\mathrm{\α}}_{i,\mathrm{\λ}}}{\mathrm{\Δ\α}}_{i,\mathrm{\λ}}{\mathrm{\δ}}_i^I& F_{\mathrm{\λ}}^c=F\left({\mathrm{\α}}_{\mathrm{\λ}}^c\right)& {\mathrm{\ΔF}}_{\mathrm{\λ}}^I=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{\∂F\left({\mathrm{\α}}_{\mathrm{\λ}}^c\right)}{\∂{\mathrm{\α}}_{i,\mathrm{\λ}}}{\mathrm{\Δ\α}}_{i,\mathrm{\λ}}{\mathrm{\δ}}_i^I\end{array}$

利用区间矩阵求逆运算，系统子区间响应可以表示为：

$U_{\mathrm{\λ}}^I\left({\mathrm{\α}}_{\mathrm{\λ}}^I\right)=U_{\mathrm{\λ}}^c+{\mathrm{\ΔU}}_{\mathrm{\λ}}^I={(K_{\mathrm{\λ}}^c+{\mathrm{\ΔK}}_{\mathrm{\λ}}^I)}^{-1}(F_{\mathrm{\λ}}^c+{\mathrm{\ΔF}}_{\mathrm{\λ}}^I)$

其中为子区间响应的中值，表示系统响应在中值处的波动变化。

利用纽曼级数，区间矩阵的逆可以表示为：

${(K_{\mathrm{\λ}}^c+{\mathrm{\ΔK}}_{\mathrm{\λ}}^I)}^{-1}={\left(K_{\mathrm{\λ}}^c\right)}^{-1}+{\left(K_{\mathrm{\λ}}^c\right)}^{-1}\underset{r=1}{\overset{\mathrm{\∞}}{\Σ}}{(-{\mathrm{\ΔK}}_{\mathrm{\λ}}^I{\left(K_{\mathrm{\λ}}^c\right)}^{-1})}^r$

基于泰勒展式，上式中的项可以进一步改写为：

$\underset{r=1}{\overset{\mathrm{\∞}}{\mathrm{\Σ}}}{\left(-{\mathrm{\ΔK}}_{\mathrm{\λ}}^I{\left(K_{\mathrm{\λ}}^c\right)}^{-1}\right)}^r=\underset{r=1}{\overset{\mathrm{\∞}}{\mathrm{\Σ}}}{\left(-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{\∂K\left({\mathrm{\α}}_{\mathrm{\λ}}^c\right)}{\∂{\mathrm{\α}}_{i,\mathrm{\λ}}}{\mathrm{\Δ\α}}_{i,\mathrm{\λ}}{\mathrm{\δ}}_i^I{\left(K_{\mathrm{\λ}}^c\right)}^{-1}\right)}^r=\underset{r=1}{\overset{\mathrm{\∞}}{\mathrm{\Σ}}}{\left(-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\δ}}_i^I{K}_{i,\mathrm{\λ}}\right)}^r$

其中 $K_{i,\mathrm{\λ}}=\frac{\∂K\left({\mathrm{\α}}_{\mathrm{\λ}}^c\right)}{\∂{\mathrm{\α}}_{i,\mathrm{\λ}}}{\mathrm{\Δ\α}}_{i,\mathrm{\λ}}{\left(K_{\mathrm{\λ}}^c\right)}^{-1}$

当r取不同值时，上式中的求和项可具体表示为：

${\left(-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\δ}}_i^I{K}_{i,\mathrm{\λ}}\right)}^{r=1}=-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\δ}}_i^I{K}_{i,\mathrm{\λ}}$

${\left(-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\δ}}_i^I{K}_{i,\mathrm{\λ}}\right)}^{r=2}=\underset{i,j=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\δ}}_i^I{\mathrm{\δ}}_j^I{K}_{i,\mathrm{\λ}}{K}_{j,\mathrm{\λ}}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\left({\mathrm{\δ}}_i^I{K}_{i,\mathrm{\λ}}\right)}^2+\underset{i\≠j}{\overset{n}{\underset{i,j=1}{\mathrm{\Σ}}}}{\mathrm{\δ}}_i^I{\mathrm{\δ}}_j^I{K}_{i,\mathrm{\λ}}{K}_{j,\mathrm{\λ}}$

......

其中记号表示除i＝j＝k外的所有组合情况。

如此一来，区间矩阵逆的纽曼级数就可改写为：

当K_i,λ的谱范数小于1时，利用矩阵级数的收敛性，同时忽略交叉项，可得：

$\begin{array}{c}.{\left(K_{\mathrm{\λ}}^c+{\mathrm{\ΔK}}_{\mathrm{\λ}}^I\right)}^{-1}\≈{\left(K_{\mathrm{\λ}}^c\right)}^{-1}+{\left(K_{\mathrm{\λ}}^c\right)}^{-1}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{r=1}{\overset{\mathrm{\∞}}{\mathrm{\Σ}}}{\left(-{\mathrm{\δ}}_i^I{K}_{i,\mathrm{\λ}}\right)}^r\\ ={\left(K_{\mathrm{\λ}}^c\right)}^{-1}+{\left(K_{\mathrm{\λ}}^c\right)}^{-1}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{-{\mathrm{\δ}}_i^I{K}_{i,\mathrm{\λ}}}{I+{\mathrm{\δ}}_i^I{K}_{i,\mathrm{\λ}}}={\left(K_{\mathrm{\λ}}^c\right)}^{-1}+{\left(K_{\mathrm{\λ}}^c\right)}^{-1}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{E}_{i,\mathrm{\λ}}^I\end{array}$

其中：

$E_{i,\mathrm{\λ}}^I=\frac{-{\mathrm{\δ}}_i^I{K}_{i,\mathrm{\λ}}}{I+{\mathrm{\δ}}_i^I{K}_{i,\mathrm{\λ}}}=-I+\frac{I}{I+{\mathrm{\δ}}_i^I{K}_{i,\mathrm{\λ}}}=E_{i,\mathrm{\λ}}^c+{\mathrm{\ΔE}}_{i,\mathrm{\λ}}^I=E_{i,\mathrm{\λ}}^c+{\mathrm{\ΔE}}_{i,\mathrm{\λ}}{\mathrm{\δ}}_{e_i}^I$

$\begin{array}{cc}{E}_{i,\mathrm{\λ}}^c=\frac{1}{2}(\frac{K_{i,\mathrm{\λ}}}{I-K_{i,\mathrm{\λ}}}-\frac{K_{i,\mathrm{\λ}}}{I+K_{i,\mathrm{\λ}}})& {\mathrm{\ΔE}}_{i,\mathrm{\λ}}=\frac{1}{2}|\frac{K_{i,\mathrm{\λ}}}{I-K_{i,\mathrm{\λ}}}+\frac{K_{i,\mathrm{\λ}}}{I+K_{i,\mathrm{\λ}}}|\end{array}$

I为单位矩阵，为区间矩阵对应的标准区间变量。

将上述含有部分高阶项的纽曼级数代回到子区间响应表达式中，利用摄动理论可得：

$U_{\mathrm{\λ}}^c={\left(K_{\mathrm{\λ}}^c\right)}^{-1}{F}_{\mathrm{\λ}}^c+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\left(K_{\mathrm{\λ}}^c\right)}^{-1}{E}_{i,\mathrm{\λ}}^c{F}_{\mathrm{\λ}}^c={\left(K_{\mathrm{\λ}}^c\right)}^{-1}(I+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{E}_{i,\mathrm{\λ}}^c)F_{\mathrm{\λ}}^c$

$\begin{array}{c}{\mathrm{\ΔU}}_{\mathrm{\λ}}^I={\left(K_{\mathrm{\λ}}^c\right)}^{-1}\left(I+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{E}_{i,\mathrm{\λ}}^c\right){\mathrm{\ΔF}}_{\mathrm{\λ}}^I+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\left(K_{\mathrm{\λ}}^c\right)}^{-1}{\mathrm{\ΔE}}_{i,\mathrm{\λ}}^I{F}_{\mathrm{\λ}}^c\\ =\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\left(K_{\mathrm{\λ}}^c\right)}^{-1}\left(I+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{E}_{i,\mathrm{\λ}}^c\right)\frac{\∂F\left({\mathrm{\α}}_{\mathrm{\λ}}^c\right)}{\∂{\mathrm{\α}}_{j,\mathrm{\λ}}}{\mathrm{\Δ\α}}_{j,\mathrm{\λ}}{\mathrm{\δ}}_j^I+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\left(K_{\mathrm{\λ}}^c\right)}^{-1}{\mathrm{\ΔE}}_{i,\mathrm{\λ}}{F}_{\mathrm{\λ}}^c{\mathrm{\δ}}_{e_i}^I\end{array}$

利用关于标准区间变量和的单调性，可以快速得到子区间响应的半径ΔU_λ：

${\mathrm{\ΔU}}_{\mathrm{\λ}}=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\left|{\left(K_{\mathrm{\λ}}^c\right)}^{-1}\left(I+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{E}_{i,\mathrm{\λ}}^c\right)\frac{\∂F\left({\mathrm{\α}}_{\mathrm{\λ}}^c\right)}{\∂{\mathrm{\α}}_{j,\mathrm{\λ}}}{\mathrm{\Δ\α}}_{j,\mathrm{\λ}}\right|+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left|{\left(K_{\mathrm{\λ}}^c\right)}^{-1}{\mathrm{\ΔE}}_{i,\mathrm{\λ}}{F}_{\mathrm{\λ}}^c\right|$

利用区间表示方法和区间组合运算法则，可以先后得到子区间向量和区间向量影响下的系统响应变化范围：

$U_{\mathrm{\λ}}^I\left({\mathrm{\α}}_{\mathrm{\λ}}^I\right)=\[{\underset{\‾}{U}}_{\mathrm{\λ}}\left({\mathrm{\α}}_{\mathrm{\λ}}^I\right),{\stackrel{\‾}{U}}_{\mathrm{\λ}}\left({\mathrm{\α}}_{\mathrm{\λ}}^I\right)\]=\[U_{\mathrm{\λ}}^c-{\mathrm{\ΔU}}_{\mathrm{\λ}},U_{\mathrm{\λ}}^c+{\mathrm{\ΔU}}_{\mathrm{\λ}}\]$

$U_{\mathrm{\λ}}^I\left({\mathrm{\β}}_{\mathrm{\λ}}^I\right)=\underset{i=1,2,...,N}{\∪}{U}_{\mathrm{\λ}}^I\left({\mathrm{\α}}_{i,\mathrm{\λ}}^I\right)=\[\underset{i=1,2,\mathrm{...},N}{min}{\underset{\‾}{U}}_{\mathrm{\λ}}\left({\mathrm{\α}}_{i,\mathrm{\λ}}^I\right),\underset{i=1,2,...,N}{max}{\stackrel{\‾}{U}}_{\mathrm{\λ}}\left({\mathrm{\α}}_{i,\mathrm{\λ}}^I\right)\]$

其中表示由子区间任意组合所得到的N个子区间参数向量。

根据表达式U＝(up)^T从向量中提取声压项，就可以得到λ截集水平下的声压区间变化范围

对选定的所有截集水平重复上述操作，进而可得到所有截集水平对应的声压区间。

选定300Hz，400Hz，500Hz，600Hz，700Hz五个频率，上述λ₈＝0.7和λ₆＝0.5两个截集水平下观测点处区间声压的计算结果如表2和表3所示。通过与样本数为10⁷的传统蒙特卡洛抽样方法对比可以看出，本方法的计算误差小于2％，计算精度可以完全满足工程需求。另外，本方法的计算时间要远远少于蒙特卡洛抽样方法，更适用于实际复杂工程问题。

表2截集水平λ₈＝0.7下观测点处区间声压上下界

表3截集水平λ₆＝0.5下观测点处区间声压上下界

步骤六：将步骤五中得到的所有截集水平下的声压区间进行重组，利用模糊分解定理可以得到模糊参数向量β影响下模糊声压的隶属度函数p(β)：

$p\left(\mathrm{\β}\right)=\underset{j=1,...,11}{\∪}\left\{{\mathrm{\λ}}_j{p}_{{\mathrm{\λ}}_j}^I\right({\mathrm{\β}}_{{\mathrm{\λ}}_j}^I\left)\right\}$

以300Hz和700Hz两个频率为例，观测点处模糊声压的隶属度函数如图3和图4所示，横坐标表示声压值，纵坐标表示声压对应的隶属度大小，实线和虚线分别表示蒙特卡洛抽样方法和本发明方法计算得到的结果。可以看出，模糊声压的隶属度函数与模型输入参数一样，仍近似满足三角分布规律。另外，本发明计算得到的隶属度函数曲线与传统蒙特卡洛抽样得到的参考值吻合程度很好，计算结果真实可信。用本发明可以有效解决含有模糊输入参数的汽车结构振动声辐射问题，计算精度高，计算时间短，此功能是一般商用软件所不能实现的。

以上所述的仅为本发明的较佳实施例而已，本发明不仅仅局限于上述实施例，凡在本发明的精神和原则之内所作的局部改动、等同替换、改进等均应包含在本发明的保护范围之内。

Claims (3)

1.一种基于高阶摄动理论的模糊声压快速求解方法，其特征在于包括以下步骤：

步骤一：针对结构振动声辐射问题，选取合适的单元类型对结构和声场进行离散，计及结构和声场的耦合作用，建立声辐射有限元离散方程；

步骤二：充分考虑声辐射问题中的不确定因素，并用模糊变量对其进行定量化表示，根据步骤一中的有限元离散方程建立声辐射问题的模糊有限元方程；

步骤三：选取截集水平，利用截集运算，将步骤二中的模糊变量转化为区间变量；

步骤四：利用子区间分解技术，将步骤三中得到的区间变量中不确定度超过5％的区间变量分解为多个子区间，使每个子区间的不确定度小于5％，进而将步骤二中的模糊有限元方程转化为一组子区间有限元方程；

步骤五：基于高阶摄动理论对步骤四中的所有子区间有限元方程进行求解，利用区间组合运算得到声压的区间变化范围，进而得到所有截集水平对应的声压区间；

步骤六：将步骤五中得到的所有截集水平下的声压区间进行重组，最终得到模糊声压的隶属度函数。

2.根据权利要求1所述的一种基于高阶摄动理论的模糊声压快速求解方法，其特征在于：所述步骤三中截集水平的选取并不是固定不变的；根据模糊不确定变量隶属度函数的分布类型来确定所需截集水平的数量规模和数值大小。

3.根据权利要求1所述的一种基于高阶摄动理论的模糊声压快速求解方法，其特征在于：所述步骤四中子区间分解策略并不是固定不变的；根据每个区间变量不确定度的大小以及计算精度的要求灵活选取分解形式，只要保证每个子区间的不确定度小于5％即可。