CN103559366A - 基于光滑有限元边界元法的车身结构声场耦合预测方法 - Google Patents

Abstract

本发明提出了一种基于光滑有限元边界元法的车身结构声场耦合预测方法，包括如下步骤：生成车身结构-声腔网格模型；利用光滑有限元法对板壳单元中应力场进行光滑处理；采用混合积分方法对板壳中剪切应变的剪切项在自然坐标系下单独积分；构造车身板壳结构动力学方程的光滑Galerkin弱形式并将其离散后得到其动力学方程；利用边界元法构造声场仿真模型；根据耦合界面上位移和压力连续，得到结构-声场耦合的光滑有限元边界元法模型；利用耦合模型进行仿真预测。本发明在车身结构-声场耦合问题中，能够得到更好的计算效果，更高有效的分析频带宽，而且该方法对模型的质量要求更低，这样能降低更多的前处理时间，在工程应用中前景广阔。

Description

基于光滑有限元边界元法的车身结构声场耦合预测方法

技术领域

本发明涉及车身结构-声场耦合预测方法，具体涉及一种基于光滑有限元-边界元法的车身结构-声场耦合预测方法。

背景技术

目前，车身的结构-声场耦合预测分析是车身NVH（Noise、Vibration、Harshness，噪声、振动与声振粗糙度）性能CAE（computerAidedEngineering，计算机辅助工程）分析中的一项重要工作，它对车身NVH性能预测以及指导NVH性能开发都有十分重要的意义。

目前，车身结构-声场耦合预测方法有解析法和数值法，其中解析法一般只针对简单问题和模型；数值法是目前主流的预测方法，包括耦合的有限元法和有限元/边界元法，并在许多商业软件如Nastran、Sysnosie等中广泛应用。许多工程师对这些方法使用比较娴熟，在车身NVH性能开发与分析中广泛采用。但这类方法也存在一些问题，不管是耦合的有限元法和有限元/边界元法，结构-声场耦合分析中的结构域分析均采用有限元法，由于有限元法模型过于“刚硬”，存在数值色散效应，其导致预测结果受模型网格尺寸大小和计算频率的高低的影响较大。

为满足车身NVH问题预测结果的可靠性，工程师需要在建模过程中对模型质量进行检查，这需要耗费较多的前处理时间和人力，因此，有必要对提出一些新的预测方法以改进预测结果。另外，由于工程师们已经习惯了有限元、边界元网格模型，这就要求新的预测方法最好是基于有限元、边界元网格模型，这样能更好地保证模型的通用性。

发明内容

为了克服上述现有技术中存在的缺陷，本发明的目的是提供一种基于光滑有限元边界元法的车身结构声场耦合预测方法,该方法在结构域中采用光滑有限元法降低结构模型硬度，减少数值色散效应，提高了车身结构频率响应分析的精度和分析频率范围。

为了实现本发明的上述目的，本发明提供了一种基于光滑有限元边界元法的车身结构声场耦合预测方法，包括如下步骤：

S1，生成车身结构-声腔网格模型；

S2，利用光滑有限元法中分区光滑技术对板壳单元中应力场ε进行光滑处理，

$\mathrm{\ϵ}={\{{\mathrm{\ϵ}}_{\mathrm{xx}},{\mathrm{\ϵ}}_{\mathrm{yy}},{\mathrm{\γ}}_{\mathrm{xy}},{\mathrm{\γ}}_{\mathrm{xz}},{\mathrm{\γ}}_{\mathrm{yz}}\}}^T=\left\{\begin{array}{c}{\mathrm{\ϵ}}_m\\ 0\end{array}\right\}+\left\{\begin{array}{c}{\mathrm{z\ϵ}}_b\\ 0\end{array}\right\}+\left\{\begin{array}{c}0\\ {\mathrm{\ϵ}}_s\end{array}\right\}$

其中，ε_xx为x向弯曲应力，ε_yy为y向弯曲应力，γ_xy为x-y平面剪切应力，γ_xz为x-z平面剪切应力，γ_yz为y-z平面剪切应力，ε_m为膜应变，zε_b为弯曲应力，ε_s为剪切应力；

将相应单元划分成不重叠的四个光滑域，对车身板壳中膜应力和弯曲应力进行光滑处理，得到光滑应力

；光滑膜应变向量

和光滑弯曲应变向量表示如下，

${\stackrel{\‾}{\mathrm{\ϵ}}}_{\mathrm{bc}}\left(x\right)={\∫}_{{\mathrm{\Ω}}_c}{\mathrm{\ϵ}}_b\left(x\right)H\left(x\right)\mathrm{d\Ω}$

${\stackrel{\‾}{\mathrm{\ϵ}}}_{\mathrm{mc}}\left(x\right)={\∫}_{{\mathrm{\Ω}}_C}{\mathrm{\ϵ}}_m\left(x\right)H\left(x\right)\mathrm{d\Ω}$

式中，ε_b(x)为膜应变，ε_m(x)为弯曲应变，Ω_c为光滑域，Η(x)为光滑函数，满足 $H\left(x\right)=\left\{\begin{array}{cc}1/A_C& x\∈{\mathrm{\Ω}}_C\\ 0& x\∉{\mathrm{\Ω}}_C\end{array}\right.,$ 其中 $A_C={\∫}_{{\mathrm{\Ω}}_C}\mathrm{d\Ω};$

S3，采用混合积分方法对板壳中剪切应变的剪切项在自然坐标系下单独积分；

S4，构造车身板壳结构动力学方程的光滑Galerkin弱形式；

${\∫}_{\mathrm{\Ω}}\mathrm{\δ}{\stackrel{\‾}{\mathrm{\ϵ}}}_m^T{D}_m{\stackrel{\‾}{\mathrm{\ϵ}}}_m\mathrm{d\Ω}+{\∫}_{\mathrm{\Ω}}{\mathrm{\δ}\stackrel{\‾}{\mathrm{\ϵ}}}_b^T{D}_b{\stackrel{\‾}{\mathrm{\ϵ}}}_b\mathrm{d\Ω}+{\∫}_{\mathrm{\Ω}}{\mathrm{\δ\ϵ}}_s^T{D}_s{\mathrm{\ϵ}}_s\mathrm{d\Ω}+{\∫}_{\mathrm{\Ω}}\mathrm{\δ}{u}^T{\mathrm{\ρ}}_{s}t\stackrel{\·\·}{u}\mathrm{d\Ω}+{\∫}_{\∂\mathrm{\Ω}}\mathrm{\δ}{u}^T{t}_s\mathrm{dS}-{\∫}_{\mathrm{\Ω}}{\mathrm{\δu}}^T{b}_s\mathrm{d\Ω}=0$

其中，

为虚光滑膜应变，D_m为膜本构系数，

为虚光滑弯曲应变，D_b为弯曲本构系数，为虚剪切应变，D_s为剪切本构系数，δu^T为虚位移，dS为微分面积，ρ_s为材料密度，t为单元厚度，t_s为表面载荷牵引力，b_s为体积力，

为加速度，

为光滑膜应变，

为光滑膜弯曲应变。

S5，将结构的光滑Galerkin弱形式离散后得到其动力学方程：

$\stackrel{\‾}{K}u-M\stackrel{\·\·}{u}=F_f+F_b$

式中,

为板单元的光滑刚度矩阵；M为板单元的质量矩阵，F_f为表面载荷列阵，F_b为体积力列阵;

S6，利用边界元法构造声场仿真模型，并将声场边界区分为耦合边界和非耦合边界两部分，

$\left[\begin{array}{cc}{H}_{11}& H_{12}\\ H_{21}& H_{22}\end{array}\right]\left[\begin{array}{c}{p}_a\\ p_b\end{array}\right]=-{\mathrm{\ρ\ω}}^2\left[\begin{array}{cc}{G}_{11}& G_{12}\\ G_{21}& G_{22}\end{array}\right]\left[\begin{array}{c}{u}_n^a\\ 0\end{array}\right]$

式中，ω为圆频率，H_mn和G_mn分别为分块矩阵，m取值为1或2，n取值为1或2，p_a、p_b表示耦合和非耦合部分边界的节点声压向量；

表示耦合边界a处的节点法向位移分量；

S7，根据耦合界面上位移和压力连续，得到结构-声场耦合的光滑有限元边界元法模型，

$\left[\begin{array}{ccc}\stackrel{\‾}{K}-{\mathrm{\ω}}^{2}M& L& 0\\ {\mathrm{\ρ}}_f{\mathrm{\ω}}^2{G}_{11}{n}_a& H_{11}& H_{12}\\ {\mathrm{\ρ}}_f{\mathrm{\ω}}^2{G}_{21}{n}_a& H_{12}& H_{22}\end{array}\right]\left[\begin{array}{c}U\\ p_a\\ p_b\end{array}\right]=\left[\begin{array}{c}F\\ 0\\ 0\end{array}\right]$

其中，

为结构场的光滑刚度矩阵，L为耦合矩阵，ρ_f为流体密度，n_a为耦合边界法向量，U为位移向量，F为载荷向量；

S8，根据光滑有限元-边界元法对车身结构-声场耦合模型进行仿真预测。

在本发明的一种优选实施中，

所述步骤S2中，车身结构单元中光滑域内光滑膜应变

和光滑弯曲应变为：

${\stackrel{\‾}{\mathrm{\ϵ}}}_{\mathrm{mC}}\left(x\right)=\underset{I=1}{\overset{4}{\mathrm{\Σ}}}{\stackrel{\‾}{B}}_{\mathrm{mI}}^C\left(x\right)d_I^{\mathrm{shell}}$ ${\stackrel{\‾}{\mathrm{\ϵ}}}_b^C\left(x\right)=\underset{I=1}{\overset{4}{\mathrm{\Σ}}}{\stackrel{\‾}{B}}_{\mathrm{bI}}^C\left(x\right)d_I^{\mathrm{shell}}$

式中，d^shell为壳单元位移向量，表示为 $d^{\mathrm{shell}}=\left[\begin{array}{cccc}{d}_1^{\mathrm{shell}}& d_2^{\mathrm{shell}}& d_3^{\mathrm{shell}}& d_4^{\mathrm{shell}}\end{array}\right];$

为壳单元节点I在第C光滑域的光滑膜应力矩阵：

${\stackrel{\‾}{B}}_{\mathrm{mI}}^C\left(x\right)=\frac{1}{A_C}{\∫}_{{\mathrm{\Γ}}_C}\left[\begin{array}{ccccc}{N}_I{n}_x& 0& 0& 0& 0\\ 0& N_I{n}_y& 0& 0& 0\\ N_I{n}_y& N_I{n}_x& 0& 0& 0\end{array}\right]\mathrm{d\Γ},$

为四边形壳单元节点I的形函数，I＝1,2,3，4，表示为 $N_I^{\mathrm{shell}}=\mathrm{diag}\left\{\begin{array}{ccccc}{N}_I& N_I& N_I& N_I& N_I\end{array}\right\},$ 其中，n_x、n_y分别为x、y向的法向分量，Γ_C为边界，N_I为等参单元形函数；

为壳单元节点I在第C光滑域的光滑弯曲应变矩阵

${\stackrel{\‾}{B}}_{\mathrm{bI}}^C\left(x\right)=\frac{1}{A_C}{\∫}_{{\mathrm{\Γ}}_C}\left[\begin{array}{ccccc}0& 0& 0& N_I{n}_x& 0\\ 0& 0& 0& 0& N_I{n}_y\\ 0& 0& 0& N_I{n}_y& N_I{n}_x\end{array}\right]\mathrm{d\Γ}.$

在本发明的另一种优选实施中，

步骤S3中车身板壳结构单元中剪切应变为：

其中BsI表示与壳单元节点I对应的剪切应变矩阵：

$B_{\mathrm{sI}}=\left[\begin{array}{ccccc}0& 0& N_{I,x}& 0& N_I\\ 0& 0& N_{I,y}& -N_I& 0\end{array}\right],$

其中，N_I,x为形函数对x的偏导，N_I,y为形函数对y的偏导。

在本发明的再一种优选实施中，车身板壳结构单元中光滑应力

为：

$\stackrel{\‾}{\mathrm{\ϵ}}=\left\{\begin{array}{c}{\stackrel{\‾}{\mathrm{\ϵ}}}_m\\ 0\end{array}\right\}+\left\{\begin{array}{c}z{\stackrel{\‾}{\mathrm{\ϵ}}}_b\\ 0\end{array}\right\}+\left\{\begin{array}{c}0\\ {\mathrm{\ϵ}}_s\end{array}\right\}.$

在本发明的一种优选实施中，车身结构域离散后的动力学方程为：

$({\stackrel{\‾}{K}}_s-{\mathrm{\ω}}^2{M}_s)U=F_f+F_b$

式中，

为车身板壳结构单元的光滑刚度矩阵，表示为：

${\stackrel{\‾}{K}}_s={\stackrel{\‾}{K}}^m+{\stackrel{\‾}{K}}^b+K^s=\underset{C=1}{\overset{\mathrm{SC}}{\mathrm{\Σ}}}{\left({\stackrel{\‾}{B}}_{\mathrm{mI}}^C\right)}^T{D}_m{\stackrel{\‾}{B}}_{\mathrm{mI}}^C{A}_C+\underset{C=1}{\overset{\mathrm{SC}}{\mathrm{\Σ}}}{\left({\stackrel{\‾}{B}}_{\mathrm{bI}}^C\right)}^T{D}_b{\stackrel{\‾}{B}}_{\mathrm{bI}}^C{A}_C+{\∫}_{\mathrm{\Ω}}{\left(B^s\right)}^T{D}^s{B}^s\mathrm{d\Ω}k=\mathrm{1,2},...\mathrm{SC}$

其中，

为光滑膜刚度矩阵，

为光滑弯曲刚度矩阵，K^s为剪切刚度矩阵，B^s为剪切应变矩阵，D^s为剪切本构系数，SC为光滑域个数。

M_s为单元质量矩阵，采用集中质量矩阵形式，

M_s＝diag{m₁ m₂ m₃ m₄}

$m_I=\mathrm{diag}\left\{\begin{array}{cccccc}{\mathrm{\ρ}}_s{A}_{\mathrm{CI}}t& {\mathrm{\ρ}}_s{A}_{\mathrm{CI}}t& {\mathrm{\ρ}}_s{A}_{\mathrm{CI}}t& \frac{{\mathrm{\ρ}}_s{A}_{\mathrm{CI}}{t}^3}{12}& \frac{{\mathrm{\ρ}}_s{A}_{\mathrm{CI}}{t}^3}{12}& \frac{{\mathrm{\ρ}}_s{A}_{\mathrm{CI}}{t}^3}{12}\end{array}\right\},I=\mathrm{1,2,3,4}$

式中，A_CI表示第C个光滑域的面积；m_I表示与节点I相对应的第C个光滑域的质量矩阵；ρ_s表示结构密度，所述C从1至4内取整数。

本发明车身NVH性能预测问题中采用板壳结构的光滑有限元模型，并采用分区应力光滑技术对车身结构域单元进行应力光滑处理，将光滑处理后的应力带至动力学方程的伽辽金弱形式中，获得光滑有限元的离散动力学方程，并根据结构声场耦合边界条件、结合声场边界元模型，预测车身结构-声场耦合问题。

本发明采用分区光滑处理技术改善结构模型硬度，减少数值模型色散效应，能获得比有限元法更好的计算结果，光滑有限元对网格质量要求比有限元要低，有效的分析频率更高。这一优点为本发明解决车身结构-声场耦合问题的预测提供了技术基础。利用本发明分析车身结构-声场耦合系统时，在结构域中采用光滑有限元降低结构模型硬度，减少数值色散效应，提高了车身结构频率响应分析的精度和分析频率范围，进而获得更准确的计算结果。

本发明在车身结构-声场耦合问题中，能够得到更好的计算效果，更高有效的分析频带宽，而且该方法对模型的质量要求更低(如网格扭曲度、尺寸)，这样能降低更多的前处理时间，在工程应用中前景广阔。

本发明的附加方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。

附图说明

本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解，其中：

图1是本发明基于光滑有限元边界元法的车身结构声场耦合预测方法的流程图；

图2是本发明一种优选实施方式中采用的车身壳结构单元示意图；

图3是本发明一种优选实施方式中车身结构声场耦合简化图；

图4是是本发明一种优选实施方式中光滑有限元模型光滑域的划分图；

图5是本发明一种优选实施方式中车身板状结构-声场耦合图；

图6是本发明一种优选实施方式中车身结构-声场耦合模型图；

图7是本发明一种优选实施方式中计算车身结构-声场耦合的结果。

具体实施方式

下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，仅用于解释本发明，而不能理解为对本发明的限制。

本发明提供了一种基于光滑有限元边界元法的车身结构声场耦合预测方法，如图1所示，下面结合图1对本发明的预测方法进行详细说明。

本发明的基于光滑有限元边界元法的车身结构声场耦合预测方法包括如下步骤：

S1，生成车身结构-声腔网格模型。如图2所示，车身结构声场耦合模型由车身结构模型和车内空腔模型组成，并且结构域和空腔域在耦合界面上应满足位移和压力连续的条件。对车身结构而言，其主要由板状结构构成，可采用壳单元进行模拟。根据第一弯曲剪切变形理论，车身板壳结构的位移分量u、v和w分别表示为：

u(x,y,z)＝u₀(x,y)+zθ_x(x,y)

v(x,y,z)＝v₀(x,y)+zθ_y(x,y)

w(x,y,z)＝w₀(x,y)

式中u₀,v₀和w₀分别表示壳中面x,y和z三个方向的位移；θ_y和θ_x分别表示xoz和yoz平面内的转角，如图3所示。

S2，本发明中所用的壳单元为四边形单元，需要将车身结构和空腔划分成四边形单元，本发明对车身结构离散为N_e个四边形壳单元，包含N_d个节点。设壳单元节点位移向量为u_shell＝[u,v,w,θ_x,θ_y]^T，则单元中面上的位移近似为：

$u_{\mathrm{shell}}^h\left(x\right)=\underset{I=1}{\overset{4}{\mathrm{\Σ}}}{N}_I^{\mathrm{shell}}{d}_I^{\mathrm{shell}}=N^{\mathrm{shell}}{d}^{\mathrm{shell}}$

将各单元划分成四个光滑域，如图4所示。按照分区应力光滑方法将各壳结构划分成四个光滑域，并对光滑域内膜应力和弯曲应力进行光滑处理，分别获得其光滑应力梯度矩阵；剪切应力采用MITC4单元中的混合积分方法对剪切项在自然坐标系下单独积分。

本发明利用光滑有限元法中分区光滑技术对板壳单元中应力场ε进行光滑处理，

$\mathrm{\ϵ}={\{{\mathrm{\ϵ}}_{\mathrm{xx}},{\mathrm{\ϵ}}_{\mathrm{yy}},{\mathrm{\γ}}_{\mathrm{xy}},{\mathrm{\γ}}_{\mathrm{xz}},{\mathrm{\γ}}_{\mathrm{yz}}\}}^T=\left\{\begin{array}{c}{\mathrm{\ϵ}}_m\\ 0\end{array}\right\}+\left\{\begin{array}{c}{\mathrm{z\ϵ}}_b\\ 0\end{array}\right\}+\left\{\begin{array}{c}0\\ {\mathrm{\ϵ}}_s\end{array}\right\}$

其中，ε_xx为x向弯曲应力，ε_yy为y向弯曲应力，γ_xy为x-y平面剪切应力，γ_xz为x-z平面剪切应力，γ_yz为y-z平面剪切应力，ε_m为膜应变，zε_b为弯曲应力，ε_s为剪切应力；

将相应单元划分成不重叠的四个光滑域，对车身板壳中膜应力和弯曲应力进行光滑处理，得到光滑应力；光滑膜应变向量

和光滑弯曲应变向量

表示如下，

${\stackrel{\‾}{\mathrm{\ϵ}}}_{\mathrm{bc}}\left(x\right)={\∫}_{{\mathrm{\Ω}}_c}{\mathrm{\ϵ}}_b\left(x\right)H\left(x\right)\mathrm{d\Ω}$

${\stackrel{\‾}{\mathrm{\ϵ}}}_{\mathrm{mc}}\left(x\right)={\∫}_{{\mathrm{\Ω}}_C}{\mathrm{\ϵ}}_m\left(x\right)H\left(x\right)\mathrm{d\Ω}$

式中，ε_b(x)为膜应变，ε_m(x)为弯曲应变，Ω_c为光滑域，Η(x)为光滑函数，满足 $H\left(x\right)=\left\{\begin{array}{cc}1/A_C& x\∈{\mathrm{\Ω}}_C\\ 0& x\∉{\mathrm{\Ω}}_C\end{array}\right.,$ 其中 $A_C={\∫}_{{\mathrm{\Ω}}_C}\mathrm{d\Ω}.$

在本发明的一种更加优选的实施中，车身结构单元中光滑域内光滑膜应变

和光滑弯曲应变

为：

${\stackrel{\‾}{\mathrm{\ϵ}}}_{\mathrm{mC}}\left(x\right)=\underset{I=1}{\overset{4}{\mathrm{\Σ}}}{\stackrel{\‾}{B}}_{\mathrm{mI}}^C\left(x\right)d_I^{\mathrm{shell}}$ ${\stackrel{\‾}{\mathrm{\ϵ}}}_b^C\left(x\right)=\underset{I=1}{\overset{4}{\mathrm{\Σ}}}{\stackrel{\‾}{B}}_{\mathrm{bI}}^C\left(x\right)d_I^{\mathrm{shell}}$

式中，d^shell为壳单元位移向量，表示为 $d^{\mathrm{shell}}=\left[\begin{array}{cccc}{d}_1^{\mathrm{shell}}& d_2^{\mathrm{shell}}& d_3^{\mathrm{shell}}& d_4^{\mathrm{shell}}\end{array}\right];$

为壳单元节点I在第C光滑域的光滑膜应力矩阵：

${\stackrel{\‾}{B}}_{\mathrm{mI}}^C\left(x\right)=\frac{1}{A_C}{\∫}_{{\mathrm{\Γ}}_C}\left[\begin{array}{ccccc}{N}_I{n}_x& 0& 0& 0& 0\\ 0& N_I{n}_y& 0& 0& 0\\ N_I{n}_y& N_I{n}_x& 0& 0& 0\end{array}\right]\mathrm{d\Γ},$

为四边形壳单元节点I(I＝1,2,3，4)的形函数，表示为 $N_I^{\mathrm{shell}}=\mathrm{diag}\left\{\begin{array}{ccccc}{N}_I& N_I& N_I& N_I& N_I\end{array}\right\},$ 其中，n_x、n_y分别为x、y向的法向分量，Γ_C为边界，N_I为等参单元形函数；

为壳单元节点I在第C光滑域的光滑弯曲应变矩阵

${\stackrel{\‾}{B}}_{\mathrm{bI}}^C\left(x\right)=\frac{1}{A_C}{\∫}_{{\mathrm{\Γ}}_C}\left[\begin{array}{ccccc}0& 0& 0& N_I{n}_x& 0\\ 0& 0& 0& 0& N_I{n}_y\\ 0& 0& 0& N_I{n}_y& N_I{n}_x\end{array}\right]\mathrm{d\Γ}.$

在本发明的另一种优选实施中，车身板壳结构单元中剪切应变为：

其中B_sI表示与壳单元节点I对应的剪切应变矩阵：

$B_{\mathrm{sI}}=\left[\begin{array}{ccccc}0& 0& N_{I,x}& 0& N_I\\ 0& 0& N_{I,y}& -N_I& 0\end{array}\right],$

其中，N_I,x为形函数对x的偏导，N_I,y为形函数对y的偏导。

在本发明的再一种优选实施中，车身板壳结构单元中光滑应力

为：

$\stackrel{\‾}{\mathrm{\ϵ}}=\left\{\begin{array}{c}{\stackrel{\‾}{\mathrm{\ϵ}}}_m\\ 0\end{array}\right\}+\left\{\begin{array}{c}z{\stackrel{\‾}{\mathrm{\ϵ}}}_b\\ 0\end{array}\right\}+\left\{\begin{array}{c}0\\ {\mathrm{\ϵ}}_s\end{array}\right\}.$

S3，采用混合积分方法对板壳中剪切应变的剪切项在自然坐标系下单独积分。具体计分方法可按照现有的积分方法进行。

S4，构造车身板壳结构动力学方程的光滑Galerkin弱形式；

${\∫}_{\mathrm{\Ω}}\mathrm{\δ}{\stackrel{\‾}{\mathrm{\ϵ}}}_m^T{D}_m{\stackrel{\‾}{\mathrm{\ϵ}}}_m\mathrm{d\Ω}+{\∫}_{\mathrm{\Ω}}{\mathrm{\δ}\stackrel{\‾}{\mathrm{\ϵ}}}_b^T{D}_b{\stackrel{\‾}{\mathrm{\ϵ}}}_b\mathrm{d\Ω}+{\∫}_{\mathrm{\Ω}}{\mathrm{\δ\ϵ}}_s^T{D}_s{\mathrm{\ϵ}}_s\mathrm{d\Ω}+{\∫}_{\mathrm{\Ω}}\mathrm{\δ}{u}^T{\mathrm{\ρ}}_{s}t\stackrel{\·\·}{u}\mathrm{d\Ω}+{\∫}_{\∂\mathrm{\Ω}}\mathrm{\δ}{u}^T{t}_s\mathrm{dS}-{\∫}_{\mathrm{\Ω}}{\mathrm{\δu}}^T{b}_s\mathrm{d\Ω}=0$

其中，

为虚光滑膜应变，D_m为膜本构系数，

为虚光滑弯曲应变，D_b为弯曲本构系数，

为虚剪切应变，D_s为剪切本构系数，δu^T为虚位移，dS为微分面积，ρ_s为材料密度，t为单元厚度，t_s为表面载荷牵引力，b_s为体积力，

为加速度，

为光滑膜应力，

为光滑弯曲应力。

S5，将结构的光滑Galerkin弱形式离散后得到其动力学方程：

$\stackrel{\‾}{K}u-M\stackrel{\·\·}{u}=F_f+F_b$

式中,

为板单元的光滑刚度矩阵；M为板单元的质量矩阵，F_f为表面载荷列阵，F_b为体积力列阵。

在本发明的一种优选实方式中，车身结构域离散后的动力学方程为：

$({\stackrel{\‾}{K}}_s-{\mathrm{\ω}}^2{M}_s)U=F_f+F_b$

式中，

为车身板壳结构单元的光滑刚度矩阵，表示为：

${\stackrel{\‾}{K}}_s={\stackrel{\‾}{K}}^m+{\stackrel{\‾}{K}}^b+K^s=\underset{C=1}{\overset{\mathrm{SC}}{\mathrm{\Σ}}}{\left({\stackrel{\‾}{B}}_{\mathrm{mI}}^C\right)}^T{D}_m{\stackrel{\‾}{B}}_{\mathrm{mI}}^C{A}_C+\underset{C=1}{\overset{\mathrm{SC}}{\mathrm{\Σ}}}{\left({\stackrel{\‾}{B}}_{\mathrm{bI}}^C\right)}^T{D}_b{\stackrel{\‾}{B}}_{\mathrm{bI}}^C{A}_C+{\∫}_{\mathrm{\Ω}}{\left(B^s\right)}^T{D}^s{B}^s\mathrm{d\Ω},k=\mathrm{1,2},...\mathrm{SC}$

其中，

为光滑膜刚度矩阵，

为光滑弯曲刚度矩阵，K^s为剪切刚度矩阵，B^s为剪切应变矩阵，D^s为剪切本构系数，SC为光滑域个数。

M_s为单元质量矩阵，采用集中质量矩阵形式，

M_s＝diag{m₁ m₂ m₃ m₄}

$m_I=\mathrm{diag}\left\{\begin{array}{cccccc}{\mathrm{\ρ}}_s{A}_{\mathrm{CI}}t& {\mathrm{\ρ}}_s{A}_{\mathrm{CI}}t& {\mathrm{\ρ}}_s{A}_{\mathrm{CI}}t& \frac{{\mathrm{\ρ}}_s{A}_{\mathrm{CI}}{t}^3}{12}& \frac{{\mathrm{\ρ}}_s{A}_{\mathrm{CI}}{t}^3}{12}& \frac{{\mathrm{\ρ}}_s{A}_{\mathrm{CI}}{t}^3}{12}\end{array}\right\},I=\mathrm{1,2,3,4}$

式中，A_CI表示第C个光滑域的面积；m_I表示与节点I相对应的第C个光滑域的质量矩阵；ρ_s表示结构密度，C从1至4内取整数。

S6，利用边界元法构造声场仿真模型，并将声场边界区分为耦合边界和非耦合边界两部分，其离散边界积分方程如下：

$\left[\begin{array}{cc}{H}_{11}& H_{12}\\ H_{21}& H_{22}\end{array}\right]\left[\begin{array}{c}{p}_a\\ p_b\end{array}\right]=-{\mathrm{\ρ\ω}}^2\left[\begin{array}{cc}{G}_{11}& G_{12}\\ G_{21}& G_{22}\end{array}\right]\left[\begin{array}{c}{u}_n^a\\ 0\end{array}\right]$

式中，ω为圆频率，H_mn和G_mn分别为分块矩阵，m取值为1或2，n取值为1或2，p_a、p_b表示耦合和非耦合部分边界的节点声压向量；

表示耦合边界a处的节点法向位移分量；

S7，根据位移连续和压力连续的耦合边界条件，将车身结构域模型和空腔域模型耦合，如图5所示。获得光滑有限元边界元法的预测模型：

$\left[\begin{array}{ccc}\stackrel{\‾}{K}-{\mathrm{\ω}}^{2}M& L& 0\\ {\mathrm{\ρ}}_f{\mathrm{\ω}}^2{G}_{11}{n}_a& H_{11}& H_{12}\\ {\mathrm{\ρ}}_f{\mathrm{\ω}}^2{G}_{21}{n}_a& H_{12}& H_{22}\end{array}\right]\left[\begin{array}{c}U\\ p_a\\ p_b\end{array}\right]=\left[\begin{array}{c}F\\ 0\\ 0\end{array}\right]$

其中，为结构场的光滑刚度矩阵，L为耦合矩阵，ρ_f为流体密度，n_a为耦合边界法向量，U为位移向量，F为载荷向量。

S8，根据光滑有限元边界元法对车身结构-声场耦合模型进行仿真预测。

本发明中所提出的预测流程通过分区应力光滑技术能有效降低模型“硬度”，降低对结构网格模型质量和尺寸要求，进而提高预测精度，提高模态频率和频响分析的精度。

运用本发明中方法对车身结构-声场耦合问题的概念模型进行仿真，预测其声学频率响应，如图6所示。该实施方式主要评价概念阶段车身地板结构对车内声场的贡献影响。为验证光滑有限元边界元法的有效性，地板结构模型划分成168个四边形壳单元，网格尺寸较大；声腔边界四边形单元为764个。参考结果则由精细网格模型(地板900个四边形壳单元、声腔边界6724个四边形单元)通过商业软件的数值实验得到，可以视为近似结果。

该车身结构-声场耦合问题的材料参数条件如下：材料参数为：弹性模量E=3.0E6Mpa，泊松比μ=0.3，密度ρ_s=7800kg/m3，厚度为1mm。车内空气声腔的参数为：密度ρ_f=1.25kg/m3，声速c=343m/s；载荷激励为底盘硬点传递给车身地板的单位谐波力，频率范围50-150Hz；地板边界固支，声腔中非耦合边界设置成刚性，耦合边界设置成强耦合方式。运用光滑有限元边界元法预测驾驶员右耳处的声压频率响应，结果如图7所示。

从图中可以看出：对于粗糙网格质量的计算模型，光滑有限元边界元法在分析车身结构-声场耦合问题时具有比有限元边界元法更高的精度，其结果接近参考结果值，表明光滑有限元边界元法具有很好的预测车身结构声场耦合问题。

基于以上方法介绍和数值实验分析，证明了基于光滑有限元-边界元法非常适合于车身结构-声场耦合问题预测，其预测精度比传统方法要高；对模型的质量要求更低。

在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

尽管已经示出和描述了本发明的实施例，本领域的普通技术人员可以理解：在不脱离本发明的原理和宗旨的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由权利要求及其等同物限定。

Claims (5)

1.一种基于光滑有限元边界元法的车身结构声场耦合预测方法，其特征在于，包括如下步骤：

S1，生成车身结构-声腔网格模型；

S2，利用光滑有限元法中分区光滑技术对板壳单元中应力场ε进行光滑处理，

$\mathrm{\ϵ}={\{{\mathrm{\ϵ}}_{\mathrm{xx}},{\mathrm{\ϵ}}_{\mathrm{yy}},{\mathrm{\γ}}_{\mathrm{xy}},{\mathrm{\γ}}_{\mathrm{xz}},{\mathrm{\γ}}_{\mathrm{yz}}\}}^T=\left\{\begin{array}{c}{\mathrm{\ϵ}}_m\\ 0\end{array}\right\}+\left\{\begin{array}{c}{\mathrm{z\ϵ}}_b\\ 0\end{array}\right\}+\left\{\begin{array}{c}0\\ {\mathrm{\ϵ}}_s\end{array}\right\}$

其中，ε_xx为x向弯曲应力，ε_yy为y向弯曲应力，γ_xy为x-y平面剪切应力，γ_xz为x-z平面剪切应力，γ_yz为y-z平面剪切应力，ε_m为膜应变，zε_b为弯曲应力，ε_s为剪切应力；

将所述板壳单元划分成不重叠的四个光滑域，对车身板壳中膜应力和弯曲应力进行光滑处理，得到光滑应力

；光滑膜应变向量

和光滑弯曲应变向量

表示如下，

${\stackrel{\‾}{\mathrm{\ϵ}}}_{\mathrm{bc}}\left(x\right)={\∫}_{{\mathrm{\Ω}}_c}{\mathrm{\ϵ}}_b\left(x\right)H\left(x\right)\mathrm{d\Ω}$

${\stackrel{\‾}{\mathrm{\ϵ}}}_{\mathrm{mc}}\left(x\right)={\∫}_{{\mathrm{\Ω}}_C}{\mathrm{\ϵ}}_m\left(x\right)H\left(x\right)\mathrm{d\Ω}$

式中，ε_b(x)为弯曲应变，ε_m(x)为膜应变，Ω_c为光滑域，Η(x)为光滑函数，满足 $H\left(x\right)=\left\{\begin{array}{cc}1/A_C& x\∈{\mathrm{\Ω}}_C\\ 0& x\∉{\mathrm{\Ω}}_C\end{array}\right.,$ 其中 $A_C={\∫}_{{\mathrm{\Ω}}_C}\mathrm{d\Ω};$

S3，采用混合积分方法对板壳中剪切应变的剪切项在自然坐标系下单独积分；

S4，构造车身板壳结构动力学方程的光滑Galerkin弱形式；

${\∫}_{\mathrm{\Ω}}\mathrm{\δ}{\stackrel{\‾}{\mathrm{\ϵ}}}_m^T{D}_m{\stackrel{\‾}{\mathrm{\ϵ}}}_m\mathrm{d\Ω}+{\∫}_{\mathrm{\Ω}}{\mathrm{\δ}\stackrel{\‾}{\mathrm{\ϵ}}}_b^T{D}_b{\stackrel{\‾}{\mathrm{\ϵ}}}_b\mathrm{d\Ω}+{\∫}_{\mathrm{\Ω}}{\mathrm{\δ\ϵ}}_s^T{D}_s{\mathrm{\ϵ}}_s\mathrm{d\Ω}+{\∫}_{\mathrm{\Ω}}\mathrm{\δ}{u}^T{\mathrm{\ρ}}_{s}t\stackrel{\·\·}{u}\mathrm{d\Ω}+{\∫}_{\∂\mathrm{\Ω}}\mathrm{\δ}{u}^T{t}_s\mathrm{dS}-{\∫}_{\mathrm{\Ω}}{\mathrm{\δu}}^T{b}_s\mathrm{d\Ω}=0$

其中，

为虚光滑膜应变，D_m为膜本构系数，

为虚光滑弯曲应变，D_b为弯曲本构系数，

为虚剪切应变，D_s为剪切本构系数，δu^T为虚位移，dS为微分面积，ρ_s为材料密度，t为单元厚度，t_s为表面载荷牵引力，b_s为体积力，

为加速度，

为光滑弯曲应力，

为光滑剪切应力。

S5，将结构的光滑Galerkin弱形式离散后得到其动力学方程：

$\stackrel{\‾}{K}u-M\stackrel{\·\·}{u}=F_f+F_b$

式中,

为板单元的光滑刚度矩阵；M为板单元的质量矩阵，F_f为表面载荷列阵，F_b为体积力列阵;

S6，利用边界元法构造声场仿真模型，并将声场边界区分为耦合边界和非耦合边界两部分，

$\left[\begin{array}{cc}{H}_{11}& H_{12}\\ H_{21}& H_{22}\end{array}\right]\left[\begin{array}{c}{p}_a\\ p_b\end{array}\right]=-{\mathrm{\ρ\ω}}^2\left[\begin{array}{cc}{G}_{11}& G_{12}\\ G_{21}& G_{22}\end{array}\right]\left[\begin{array}{c}{u}_n^a\\ 0\end{array}\right]$

式中，ω为圆频率，H_mn和G_mn分别为分块矩阵，m取值为1或2，n取值为1或2，p_a、p_b表示耦合和非耦合部分边界的节点声压向量；

表示耦合边界a处的节点法向位移分量；

S7，根据耦合界面上位移和压力连续，得到结构-声场耦合的光滑有限元边界元法模型，

$\left[\begin{array}{ccc}\stackrel{\‾}{K}-{\mathrm{\ω}}^{2}M& L& 0\\ {\mathrm{\ρ}}_f{\mathrm{\ω}}^2{G}_{11}{n}_a& H_{11}& H_{12}\\ {\mathrm{\ρ}}_f{\mathrm{\ω}}^2{G}_{21}{n}_a& H_{12}& H_{22}\end{array}\right]\left[\begin{array}{c}U\\ p_a\\ p_b\end{array}\right]=\left[\begin{array}{c}F\\ 0\\ 0\end{array}\right]$

其中，为结构场的光滑刚度矩阵，L为耦合矩阵，ρ_f为流体密度，n_a为耦合边界法向量，U为位移向量，F为载荷向量；

S8，根据光滑有限元-边界元法对车身结构-声场耦合模型进行仿真预测。

2.如权利要求1所述的基于光滑有限元边界元法的车身结构声场耦合预测方法，其特征在于：所述步骤S2中，车身结构单元中光滑域内光滑膜应变和光滑弯曲应变

为：

${\stackrel{\‾}{\mathrm{\ϵ}}}_{\mathrm{mC}}\left(x\right)=\underset{I=1}{\overset{4}{\mathrm{\Σ}}}{\stackrel{\‾}{B}}_{\mathrm{mI}}^C\left(x\right)d_I^{\mathrm{shell}}$ ${\stackrel{\‾}{\mathrm{\ϵ}}}_b^C\left(x\right)=\underset{I=1}{\overset{4}{\mathrm{\Σ}}}{\stackrel{\‾}{B}}_{\mathrm{bI}}^C\left(x\right)d_I^{\mathrm{shell}}$

式中，d^shell为壳单元位移向量，表示为 $d^{\mathrm{shell}}=\left[\begin{array}{cccc}{d}_1^{\mathrm{shell}}& d_2^{\mathrm{shell}}& d_3^{\mathrm{shell}}& d_4^{\mathrm{shell}}\end{array}\right];$

为壳单元节点I在第C光滑域的光滑膜应力矩阵：

${\stackrel{\‾}{B}}_{\mathrm{mI}}^C\left(x\right)=\frac{1}{A_C}{\∫}_{{\mathrm{\Γ}}_C}\left[\begin{array}{ccccc}{N}_I{n}_x& 0& 0& 0& 0\\ 0& N_I{n}_y& 0& 0& 0\\ N_I{n}_y& N_I{n}_x& 0& 0& 0\end{array}\right]\mathrm{d\Γ},$

为四边形壳单元节点I的形函数，I＝1,2,3，4，表示为 $N_I^{\mathrm{shell}}=\mathrm{diag}\left\{\begin{array}{ccccc}{N}_I& N_I& N_I& N_I& N_I\end{array}\right\},$ 其中，n_x、n_y分别为x、y向的法向分量，Γ_C为边界，N_I为等参单元形函数；

为壳单元节点I在第C光滑域的光滑弯曲应变矩阵

${\stackrel{\‾}{B}}_{\mathrm{bI}}^C\left(x\right)=\frac{1}{A_C}{\∫}_{{\mathrm{\Γ}}_C}\left[\begin{array}{ccccc}0& 0& 0& N_I{n}_x& 0\\ 0& 0& 0& 0& N_I{n}_y\\ 0& 0& 0& N_I{n}_y& N_I{n}_x\end{array}\right]\mathrm{d\Γ}.$

3.如权利要求1所述的基于光滑有限元边界元法的车身结构声场耦合预测方法，其特征在于：步骤S3中车身板壳结构单元中剪切应变为：

其中B_sI表示与壳单元节点I对应的剪切应变矩阵：

$B_{\mathrm{sI}}=\left[\begin{array}{ccccc}0& 0& N_{I,x}& 0& N_I\\ 0& 0& N_{I,y}& -N_I& 0\end{array}\right],$

其中，N_I,x为形函数对x的偏导，N_I,y为形函数对y的偏导。

4.如权利要求1所述的基于光滑有限元边界元法的车身结构声场耦合预测方法，其特征在于：车身板壳结构单元中光滑应力为：

$\stackrel{\‾}{\mathrm{\ϵ}}=\left\{\begin{array}{c}{\stackrel{\‾}{\mathrm{\ϵ}}}_m\\ 0\end{array}\right\}+\left\{\begin{array}{c}z{\stackrel{\‾}{\mathrm{\ϵ}}}_b\\ 0\end{array}\right\}+\left\{\begin{array}{c}0\\ {\mathrm{\ϵ}}_s\end{array}\right\}.$

5.如权利要求1所述的基于光滑有限元边界元法的车身结构声场耦合预测方法，其特征在于：车身结构域离散后的动力学方程为：

$({\stackrel{\‾}{K}}_s-{\mathrm{\ω}}^2{M}_s)U=F_f+F_b$

式中，

为车身板壳结构单元的光滑刚度矩阵，表示为：

${\stackrel{\‾}{K}}_s={\stackrel{\‾}{K}}^m+{\stackrel{\‾}{K}}^b+K^s=\underset{C=1}{\overset{\mathrm{SC}}{\mathrm{\Σ}}}{\left({\stackrel{\‾}{B}}_{\mathrm{mI}}^C\right)}^T{D}_m{\stackrel{\‾}{B}}_{\mathrm{mI}}^C{A}_C+\underset{C=1}{\overset{\mathrm{SC}}{\mathrm{\Σ}}}{\left({\stackrel{\‾}{B}}_{\mathrm{bI}}^C\right)}^T{D}_b{\stackrel{\‾}{B}}_{\mathrm{bI}}^C{A}_C+{\∫}_{\mathrm{\Ω}}{\left(B^s\right)}^T{D}^s{B}^s\mathrm{d\Ωk}=\mathrm{1,2},...\mathrm{SC}$

其中，

为光滑膜刚度矩阵，

为光滑弯曲刚度矩阵，K^s为剪切刚度矩阵，B^s为剪切应变矩阵，D^s为剪切本构系数，SC为光滑域个数。

M_s为单元质量矩阵，采用集中质量矩阵形式，

M_s＝diag{m₁ m₂ m₃ m₄}

$m_I=\mathrm{diag}\left\{\begin{array}{cccccc}{\mathrm{\ρ}}_s{A}_{\mathrm{CI}}t& {\mathrm{\ρ}}_s{A}_{\mathrm{CI}}t& {\mathrm{\ρ}}_s{A}_{\mathrm{CI}}t& \frac{{\mathrm{\ρ}}_s{A}_{\mathrm{CI}}{t}^3}{12}& \frac{{\mathrm{\ρ}}_s{A}_{\mathrm{CI}}{t}^3}{12}& \frac{{\mathrm{\ρ}}_s{A}_{\mathrm{CI}}{t}^3}{12}\end{array}\right\},I=\mathrm{1,2,3,4}$

式中，A_CI表示第C个光滑域的面积；m_I表示与节点I相对应的第C个光滑域的质量矩阵；ρ_s表示结构密度，所述C从1至4内取整数。

Images