2. 广州汽车集团股份有限公司 汽车工程研究院, 广州 511434
2. Automotive Research & Deve-lopment Center, Guangzhou Automobile Group Co. Ltd., Guangzhou 511434, China
随着人们消费水平的提升和消费观念的升级, 汽车NVH (Noise, Vibration, Harshness)性能中的车内噪声水平越来越受到消费者的关注。汽车在高速行驶时, 由外部气流引起的气动噪声是车内噪声的主要来源之一[1]。
在汽车的开发过程中, 由于风洞实验成本高, 数值模拟已成为研究车辆气动噪声的主要手段, 而准确的噪声分析是建立在精确的流场数据基础之上的。因此, 提高流场计算精度十分必要。
大涡模拟(LES)是计算非定常流场的常用方法。该方法利用滤波函数对流场中的涡进行尺度划分, 对大尺度的涡结构进行直接求解, 对小尺度的涡结构则采用亚格子模型进行建模, 而汽车气动噪声的产生与流场中小尺度涡旋的形成与运动有着重要的关系。滤波处理后的N-S方程将产生表征大、小尺度相互作用的亚格子应力项, 对亚格子应力项的求解通过构建亚格子模型进行。不同亚格子模型的选取会获得不同的流场结果, 进而导致不同的声场结果。
1 数值计算理论 1.1 流场数值计算方法汽车运动的马赫数小于0.3, 其绕流可看作不可压缩流, 不可压N-S方程经过滤波后得到LES控制方程:
(1)
(2)
式中, τij为亚格子应力项, 表示可求解的大尺度涡和不可求解的小尺度涡的动量交换:
(3)
为使方程封闭, 需构建亚格子应力模型, 即τij的表达式。采用Smagorinsky-Lilly、Dynamic Smagorinsky、WALE、WMLES、WMLES S-Omega以及Dynamic Kinetic Energy Subgrid-Scale Model 6种亚格子模型进行数值模拟。亚格子模型的τij表达式可写为:
(4)
式中,μt为湍流涡粘度,
(5)
在Smagorinsky-Lilly模型[2]中, 湍流涡粘度:
(6)
其中,
(7)
式中, k为冯卡尔曼常数; d为近壁面的距离; Cs为Smagorinsky常数, 通常取0.1;Δ为亚网格尺寸, 由网格尺寸V确定, Δ=V1/3; Ls为亚格子尺度的混合长度。
动态Smagorinsky(Dynamic Smagorinsky)亚格子尺度模型[3]中, Cs基于运动尺度求解得出。该模型添加了网格滤波和检验滤波函数, 进一步得出湍流应力张量:
(8)
式中, Tij为滤波后的亚格子应力项。
通过Lilly最小二乘分析, 最终得到Cs表达式为:
(9)
式中:
(10)
局部涡粘度壁面自适应模型[4](Wall-Adapting Local Eddy-Viscosity Model, WALE)中湍流涡粘度为:
(11)
式中, 亚格子尺度的混合长度Ls和应变率张量Sijd的表达式为:
(12)
(13)
式中, 常数Cw=0.325,
(14)
代数形式的壁面函数模型[5](Algebraic Wall-Modeled LES Model, WMLES)在边界层内部使用雷诺平均方法进行模拟, 在边界层外部使用LES进行模拟。WMLES可应用于更高雷诺数的流场模拟, 其湍流涡粘度:
(15)
式中, dw为计算点与壁面的距离; S为应变率; 常数k=0.41, CSmag=0.2;Δ为亚网格尺度:
(16)
式中, hmax是六面体单元的最大边长, hwn是壁面法线方向网格间距, 常数Cw=0.15。
由于无法为恒定的剪切流提供零涡流粘度, WMLES模型在分离剪切层时会产生过大的涡粘性。在式(15)中, 若将S替换为|S-Ω|的差分形式, 则可以改进WMLES模型, 这种改进模型称为WMLES S-Omega模型, 其中Ω为旋涡强度。
动态动能亚格子模型[6](Dynamic Kinetic Energy Subgrid-Scale Model, KET)在传统Smagorinsky-Lilly模型上, 进一步考虑亚格子尺度湍动能的传输, 能更好地模拟亚格子尺度的湍流。亚格子湍流动能ksgs为:
(17)
湍流涡粘度μt通过亚格子湍流动能计算得出:
(18)
式中滤波尺寸Δf=V1/3, Ck根据流场状态获得。该模型的亚格子应力表达式为:
(19)
采用莱特希尔声类比理论(Lighthill’s acoustic analogy)对气动声学进行求解。从流体力学的连续方程和动量方程推导得出莱特希尔方程:
(20)
其中:
(21)
式中, ρ'为流体密度波动量, c0为声音在均质介质中的速度, Tij为莱特希尔张量, ui为速度矢量沿i方向的分量, pi、ρ为流体的压力和密度, p0、ρ0分别为参考压力和密度,
应用格林函数进行积分可对莱特希尔声类比方程求解, 但对于复杂边界, 格林函数往往难以获得。Oberai等人提出一种新的求解方法[7], 通过获得频域下的莱特希尔方程变分形式, 采用有限元法进行求解。
如图 1所示, 气流流经的固体区域为ΩS, 固体区域的边界为ΓS。流体区域为Ωf, 流体是非均匀的或者湍流的, 受N-S方程控制。在Ωa和Ω∞区域, 流体状态是均匀的。另外, 湍流区域(声源区)的外边界为Γf, 声学截断面Γa为无反射边界。n为区域Ωt=Ωf∪Ωa向外的法向向量。对于任意速度矢量u, 可用un= u · n =uini来表示其法向分量。
|
| 图 1 气动声学问题求解计算域示意图 Fig.1 Computational domain for aeroacoustic problem solving |
基于莱特希尔方程, 利用转换等式:
(22)
得到莱特希尔方程的频域形式:
(23)
式中, ψ=-b/iω, b为总熵, i(非下标)为虚数单位。
基于Oberai提出的变分公式, 可得:
(24)
式中, δψ为测试函数。利用格林函数进行变分转换, 同时结合动量方程, 最终得到弱解形式:
(25)
采用有限元法进行内声场求解。在一个由边界面Ω所包围、声源分布为q的流体域V中, 任意空间点(x, y, z)的声压p受二阶赫姆霍兹方程控制:
(26)
式中, k=ω/c=2πf/c为波数, c为声速, ρ为流体密度, ω为圆频率。
将式(26)转换为等效的加权残差积分方程:
(27)
式中,
(28)
式中nV为区域V向外的法向向量。
在有限元方法中, 连续体V离散成一系列有限小单元以及节点ne, 在每个单元内, 利用形函数Nie, 声压p可用一定数量np的
(29)
式中, nf为整个离散域的节点总数, N为全局形函数的(1×nf)向量, 
(30)
式中, B为全局形函数梯度分量的(3×nf)矩阵。将式(29)和(30)代入到式(28)中, 得到:
(31)
式中, K为声学刚度矩阵, M为声学质量矩阵, Q i为声学激励向量, C为声学阻尼矩阵。未知节点的声压值
在密封良好情况下, 泄露噪声、风振噪声的影响可以忽略。此时车内气动噪声可归结为以下两种机理[8]:(1)空气中的旋涡收缩与变形及其相互作用, 以及气流与壁面相互作用会产生气动噪声源, 多种气动噪声源向四周传递噪声, 其中到达玻璃壁面的声压会引起壁面的振动从而向模型内部传递噪声; (2)当气流流经玻璃表面时, 引起了壁面振动, 从而向模型内部传递噪声。
虽然作用在玻璃表面上的湍流压力的压力级数要比声压的压力级数大, 但是在传递效率上, 湍流压力的传递损失却比声压的传递损失要大。所以, 在研究外部流场诱发汽车内部气动噪声时, 需要同时考虑上述2种作用机理。
2 计算模型与设置 2.1 计算模型采用的计算模型为现代汽车公司建立的简易汽车模型Hyundai Simplified Model(HSM), 参考相关文献[9-10]建立HSM模型。如图 2所示, 模型长2 m, 宽1 m, 高1 m, 壁面从外到内分别为12 mm铝、40 mm吸声材料、4 mm隔膜以及20 mm吸声材料, 玻璃均为4 mm厚。材料属性如表 1和2所示。
|
| 图 2 HSM模型 Fig.2 Hyundai Simplified Model |
| 材料 | 密度/(kg·m-3) | 杨氏模量/GPa | 泊松比 |
| 铝 | 2700 | 69 | 0.33 |
| 隔膜 | 2000 | 0.04 | 0.45 |
| 车窗玻璃 | 2500 | 70 | 0.22 |
| 空气 | 1.2 | - | - |
| 密度/(kg·m-3) | 声速/(m·s-1) | 孔隙率 | 孔曲率 |
| 827 | 343 | 0.879 | 3.31 |
实验中, 通过垫片、密封剂、螺栓等将车窗玻璃安装在铝制外壳上。为了验证车窗玻璃边界条件对噪声传播和数值模拟的影响, 文献[11]对车窗玻璃表面内外的频率响应函数(FRF)、玻璃传递函数和车内混响时间进行了测量, 结果表明, 玻璃内外的FRF结果几乎相同。可以认为, 玻璃安装到铝制外壳上是两种材料的简单附着, 仿真模型的简化对计算结果的影响可以忽略。
2.2 流场计算参考现代汽车公司风洞实验段的外形[12]建立的流场计算域如图 3所示, 计算域长20 m, 宽14 m, 高10 m, 模型距离入口4 m。模型表面网格尺寸为5 mm, 计算域外表面网格尺寸为500 mm。利用密度盒对模型附近网格进行加密, 2个加密区的最大网格尺寸为20和100 mm, 其中, 第一个密度盒包含了湍流较为剧烈的声源区。在模型表面划分5层边界层网格, 第一层厚0.2 mm, 增长率为1.2。整个计算域网格数量约为1500万, 边界条件和参数设置见表 3。
|
| 图 3 计算域示意图 Fig.3 Computational domain |
| 边界 | 类型 | 数值 |
| Inlet | Velocity inlet | 110 km/h |
| Outlet | Pressure outlet | 0 Pa |
| Others | Wall | - |
采用Fluent进行流场计算。利用稳态结果作为瞬态的初始值进行计算, 以减少瞬态计算达到收敛的时间。分别使用上述6种亚格子模型进行瞬态计算, 以获得流场的速度和压力结果。考虑到气动噪声集中在中低频段, 本次计算的最高频率初步定为1000 Hz。根据奈奎斯特采样定理, 确定采样时长为0.08 s, 计算时间步长为2×10-4 s。
2.3 声场计算采用Actran进行声学计算, 使用积分插值的方法将流场信息映射至声学网格。外声场模型如图 4所示, 分为声源区、声传播区和无限元边界, 无限元边界为无反射边界。实验风洞为声学风洞, 壁面进行了吸声处理, 采用无反射边界和实验保持一致。根据分析, 频率的波长内至少包含6个网格单位的要求, 确定声场的网格尺寸为40 mm。采用四面体网格对外声场进行划分, 最终网格数约为140万。
|
| 图 4 外声场模型 Fig.4 External sound field model |
图 5为内声场模型局部网格。玻璃和隔膜使用三棱柱网格划分, 其余部分使用四面体网格划分, 网格尺寸和外声场相同。
|
| 图 5 内声场局部网格 Fig.5 Local grid of internal sound field |
为了和实验数据[9, 11]对比, 在模型内部建立2个监测点, 监测点位置如图 6所示。
|
| 图 6 监测点位置 Fig.6 Position of monitoring points |
由流场决定的CFD截止频率fc为:
(32)
式中, ε为湍动能耗散率, Δx为网格尺度。
图 7为CFD截止频率分布图, 由图可知在湍流核心区域, CFD截止频率大部分在1000 Hz左右, 符合计算要求。
|
| 图 7 CFD截止频率分布 Fig.7 CFD cut-off frequency distribution |
图 8和9为1/3倍频程频率中心点的声压级。总体来看, 声压级Lp随频率f的上升而逐渐降低, 不同亚格子模型计算得到的结果的变化趋势和实验一致, 但在数值上有偏差。上述2个监测点的总声压级以及和实验值的误差如表 4和5所示。
|
| 图 8 监测点1声压级对比 Fig.8 Comparison of sound pressure level at monitoring point 1 |
|
| 图 9 监测点2声压级对比 Fig.9 Comparison of sound pressure level at monitoring point 2 |
| SLM | DSLM | WALE | WMLES | WMLES S-Omega | KET | 实验值 | |
| 总声压级/dB | 90.1 | 90.8 | 80.4 | 77.2 | 80.2 | 85.9 | 71.6 |
| 误差/% | 25.8 | 26.8 | 12.3 | 7.8 | 12.0 | 19.9 | - |
| SLM | DSLM | WALE | WMLES | WMLES S-Omega | KET | 实验值 | |
| 总声压级/dB | 90.3 | 90.4 | 83.7 | 77.5 | 76.1 | 86.2 | 73.1 |
| 误差/% | 23.5 | 23.7 | 14.5 | 6.0 | 4.1 | 17.9 | - |
从总声压级看, 只有WMLES模型的结果误差在10%以内, 其余模型误差都较大。所以采用WMLES模型计算总声压级的精度更高。
由频谱曲线可知, 100~160 Hz, 6种亚格子模型得到的计算结果均大于实验值; 160~500 Hz, 计算结果均小于实验值; 500~1000 Hz, 计算结果和实验结果基本一致。在200 Hz时, 实验曲线出现峰值, 原因可能是该频率点为模型的固有频段, 共振使得声压级增大, WMLES和WMLES S-Omega模型准确地捕捉到了这一信息。
声场结果和流场结果有着直接的联系, 高频段结果和流场中小旋涡的变化及相互作用有关。出现低频段计算结果高于实验值、高频段计算结果低于实验值的原因可能是受计算资源的限制, 数值模拟没能精确捕捉到流场中的小旋涡运动, 导致部分高频段能量计入低频段。
在近壁面附近, 湍流尺度随壁面距离的减小而线性减小。所以, 壁面附近小尺度的旋涡占主要部分。由1.3节可知, 壁面附近的流场对车内噪声有着显著的影响。壁面附近气流模拟的精确性将直接影响车内噪声的计算结果。
LES在处理近壁面的流场时需要非常精细的网格尺度和非常小的时间步长(y+=1)。WMLES模型对LES进行了改进, 在处理近壁面流场时使用RANS进行模拟, 其余部分使用LES进行模拟。RANS在求解近壁面流场时, 不需要LES所要求的非常精细的网格尺度, 使其在求解高雷诺数下近壁面湍流更加有效。
在LES中, 亚格子模型表示了大、小旋涡的相互作用。几种亚格子模型的区别在于对湍流涡粘度μt的建模方式不同。湍流模拟中, μt对小尺度旋涡的形成与运动有很大的影响。在对μt建模过程中, 和网格尺寸有关的亚网格尺度Δ是一个重要的参数。其余亚格子模型采用网格体积的1/3次方定义Δ, WMLES模型改变了Δ的定义, 综合考虑了网格边长、网格中心点离壁面距离和网格法向尺寸, 此外建模时还考虑了y+值, 使得WMLES模型对壁面附近的旋涡捕捉更为精确, 而其余模型计算得到的湍流衰减过大, 这一点在文献[5]中已有体现。200 Hz时, WMLES及其改进模型WMLES S-Omega能够捕捉到峰值, 和其能够准确捕捉流场信息有关。
4 结论基于不同亚格子模型, 对HSM模型进行了流场和声场的计算, 并与实验结果进行了对比分析, 结果表明:采用相同网格尺度时, 基于WMLES亚格子模型的LES数值模拟能较为准确地获得车外气动声源。
| [1] |
Buchheim R, Dobrzynski W, Mankau H, et al. Vehicle interior noise related to external aerodynamics[J]. International Journal of Vehicle Design, 1982, 3(4): 398-410. |
| [2] |
Smagorinsky J. General circulation experiments with the primitive equtions: Ⅰ. the basic experiments[J]. Monthly Weather Review, 1963, 91(3): 99-164. DOI:10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2 |
| [3] |
Gallerano F, Napoli E. A dynamic subgrid-scale tensorial eddy viscosity model[J]. Continuum Mechanics and Thermody-namics, 1999, 11(1): 1-14. DOI:10.1007/s001610050101 |
| [4] |
Nicoud F, Ducros F. Subgrid-scale stress modelling based on the square of the velocity gradient tensor[J]. Flow, Turbulence and Combustion, 1999, 62(3): 183-200. DOI:10.1023/A:1009995426001 |
| [5] |
Shur M L, Spalart P R, Strelets M K, et al. A hybrid RANS-LES approach with delayed-DES and wall-modelled LES capabilities[J]. International Journal of Heat & Fluid Flow, 2008, 29(6): 1638-1649. |
| [6] |
Kim W W, Menon S. Application of the localized dynamic subgrid-scale model to turbulent wall-bounded flows[R]. AIAA 1997-0210, 1997.
|
| [7] |
Oberai A A, Roknaldin F, Hughes T J R. Computational procedures for determining structural-acoustic response due to hydrodynamic sources[J]. Computer Methods in Applied Mechanics and Engineering, 2000, 190(3-4): 345-361. DOI:10.1016/S0045-7825(00)00206-1 |
| [8] |
Yao H D, Davidson L. Generation of interior cavity noise due to window vibration excited by turbulent flows past a generic side-view mirror[J]. Physics of Fluids, 2018, 30(3): 36104. DOI:10.1063/1.5008611 |
| [9] |
Mendonca F G, Connelly T, Bonthu S, et al. CAE-based prediction of aero-vibro-acoustic interior noise transmission for a simple test vehicle[R]. SAE Technical Paper 2014-01-0592, 2014.
|
| [10] |
Khondge A, Lee M. Numerical investigation of sunroof buffeting for hyundai simplified model[J]. Korea Noise Vibration Engineering Association Papers, 2014, 24(3): 180-188. DOI:10.5050/KSNVE.2014.24.3.180 |
| [11] |
Cho M, Kim H G, Oh C, et al. Benchmark study of numerical solvers for the prediction of interior noise transmission excited by a-pillar vortex[C]//Proc of the 43th International Congress and Exposition on Noise Control Engineering (Inter-noise). 2014.
|
| [12] |
Kim M S, Lee J H, Kee J D, et al. Hyundai full scale aero-acoustic wind tunnel[R]. SAE Technical Paper 2001-01-0629, 2001.
|