收藏 分销(赏)

水下振荡涡流场声散射与调制特性.pdf

上传人:自信****多点 文档编号:585572 上传时间:2024-01-04 格式:PDF 页数:7 大小:5.90MB
下载 相关 举报
水下振荡涡流场声散射与调制特性.pdf_第1页
第1页 / 共7页
水下振荡涡流场声散射与调制特性.pdf_第2页
第2页 / 共7页
水下振荡涡流场声散射与调制特性.pdf_第3页
第3页 / 共7页
亲,该文档总共7页,到这儿已超出免费预览范围,如果喜欢就下载吧!
资源描述

1、本文针对水下振荡涡流场的声散射与调制问题建立了理论模型。采用间断伽辽金数值方法模拟了水下振荡涡流场对声波的散射与调制作用,研究了不同马赫数、不同涡流场尺度与声波波长的尺度波长比和不同涡流场振荡频率与声波频率的频率比条件下水下振荡涡流场对声场散射与调制特性。研究结果表明:声波经过振荡涡流场会产生散射和调制现象;振荡涡流场对声场的散射和调制作用均随着尺度波长比、马赫数增大而增强;散射和调制作用随着频率比降低而增强,随着频率比的降低,散射和调制作用增强速度变小,直至频率比的改变对声场的散射和调制作用几乎没有影响。水下振荡涡流场的声散射与调制能够使散射声场含有涡流特征信息,对水下目标探测与识别具有重要

2、意义。关键词:流声耦合;振荡涡流场;散射特性;调制特性;包络;马赫数;尺度波长比;频率比DOI:10.11990/jheu.202205001网络出版地址:https:/ 文献标志码:A 文章编号:1006-7043(2023)08-1249-07Acoustic scattering and modulation characteristics of underwater oscillating flow fieldsJING Chenxuan3,SHI Shengguo1,2,3,YANG Desen1,2,3(1.Acoustic Science and Technology Labor

3、atory,Harbin Engineering University,Harbin 150001,China;2.Key Laboratory of Marine Information Acquisition and Security(Harbin Engineering University),Ministry of Industry and Information Technology,Harbin 150001,China;3.College of Underwater Acoustic Engineering,Harbin Engineering University,Harbin

4、 150001,China)Abstract:A theoretical model is established to solve the problem of acoustic scattering and modulation of underwa-ter oscillating eddy current fields.The discontinuous Galerkin method is used to numerically simulate the scattering and modulation of the underwater oscillating eddy curre

5、nt field on the acoustic waves.The numerical simulation in-cluded the scattering and modulation characteristics of underwater oscillating eddy current fields on sound fields un-der different Mach numbers,scale-wavelength ratios of eddy current fields to sound wave wavelengths,and frequen-cy ratios o

6、f eddy current fields to sound wave frequencies.The results show that sound waves will generate scatter-ing and modulation phenomena when passing through the oscillating eddy current field.The scattering and modula-tion effects of the oscillating eddy current field on the sound field are enhanced wi

7、th an increase in the scale-wave-length ratio and Mach number.Meanwhile,the scattering and modulation effects increase with the decrease of the frequency ratio,and the enhancement speed of the scattering and modulation effects reduces with the decrease in the frequency ratio until the frequency rati

8、o variation has minimal effect on the scattering and modulation of the sound field.The acoustic scattering and modulation of underwater oscillating eddy current fields can help the scat-tered acoustic field obtain eddy current characteristic information,which is of considerable importance to underwa

9、-ter target detection and recognition.Keywords:flow-acoustic coupling;oscillating eddy current field;scattering characteristics;modulation characteris-tics;envelope;Mach number;scale-wavelength ratio;frequency ratio收稿日期:2022-05-01.网络出版日期:2023-05-23.基金项目:中央高校基本科研专项资金(3072022TS0501).作者简介:荆晨轩,男,博士研究生;时

10、胜国,男,教授,博士生导师;杨德森,男,教授,博士生导师.通信作者:时胜国,E-mail:shishengguo .声波在涡流场传播时,涡流场与声波会发生流声耦合作用,产生声散射现象1-2。回转体在绕流场中会产生交替脱落的卡门涡街,可演化出具有周期性结构的振荡尾涡,在自然界和实际工程中广泛存在3-5。研究水下振荡涡流与声波的耦合作用规律,对结构物尾流识别及探测、认识复杂流场与声波相互作用有着重要意义,具有工程应用与科研价值。哈 尔 滨 工 程 大 学 学 报第 44 卷二维涡旋声散射问题是声波在复杂流场中传播的基本问题,在水声学、气动声学中都有广泛应用6-8。Colonius 等9基于线性纳维

11、斯托克斯方程(linear navier Stokes,LNS),选取了多个基本涡流结构,采用直接数值模拟(direct numerical simulation,DNS)开展了声散射研究,并为验证模型算例提供了数值参考。Cheinet10基于线性欧拉方程(linear Euler e-quation,LEE),采用二阶中心时域有限差分(finite difference time domain,FDTD)格式数值计算了声波穿过涡流场的散射情况,为研究大气湍流对声传播的影响提供了一种数值方法。Ke 等11采用高阶加权本质无振荡(weight edessentially non-oscillat

12、ory,WENO)格式数值计算了声波穿过多种涡结构的散射情况,获得的结果与 Colonius 采用直接数值模拟的计算结果相吻合。Clair12 通过建立平面波涡流散射的半解析模型及时域有限差分格式数值求解了线性欧拉方程,研究了运动涡对散射声场的影响特性。文献13-14通过线性紧致格式直接数值求解二维欧拉方程,获得了声波穿过均熵涡流的散射特性,并分析了无量纲尺度量涡强度和长度尺度比对散射声场的影响。流声耦合研究主要集中在气动声学领域,在水下流场声散射研究方面,张咏鸥等15基于摄动声学理论建立流声耦合作用下的声散射模型,采用时域有限差分格式进行数值求解,对水下涡流场前向声散射特性进行了分析,探究了

13、散射声场形态函数和指向性。本文针对经典 Gaussian 涡流模型,采用间断伽辽金(discontinuous Galerkin,DG)方法数值求解了运动流场声散射模型。在验证数值方法有效性的基础上,研究了水下振荡涡流场对声场的散射和调制作用,分析了不同马赫数、尺度波长比和频率比条件下振荡涡流场对声场的散射与调制特性。1 流声耦合理论模型1.1 控制方程 为了建立水下振荡涡流场的声散射调制模型,从流体力学基本方程组出发16:ut+(u)u=0t+(u)=0(1)式中:、u 分别为流体密度、速度。假设声学量为扰动量,声学量级远远小于对应的流体变量:u=u0(x)+u(x,t)p=p0(x)+p(

14、x,t)=0(x)+(x,t)(2)式中:0、u0、p0分别为不考虑声扰动的流场密度、速度矢量和压力,、u、p 分别为声场对应的密度变化量、速度矢量和声压,代入流体力学基本方程组并略去高阶小量,得到16:ut+u0u+uu0+10p-20p0=0t+u0+u0+u0+0u=0(3)方程(3)为水下流声耦合作用的声场控制方程,通过对该方程的数值求解,能够数值模拟水下振荡涡流场对声波的散射和调制作用。为保证数值计算精度,本文采用 DG 方法17对控制方程进行求解,其方法结合了传统有限元方法和有限体积方法的优点,具有精度高、耗散和色散小等特点,且良好的稳定性和收敛性,容易实现高阶格式。时间步进采用四

15、阶龙格-库塔格式。为了减少边界反射对计算的影响,在边界处设置声吸收层。1.2 涡旋模型 二维涡流声散射问题是声波在复杂流场中传播的基本问题,下面针对经典 Gaussian 涡流模型9进行计算和分析:vr=0,v=2r(1-exp(-(r/L)2)=1.256 431,=2.8Lv(4)式中:vr为径向速度;v为切向速度;r 为观测点到涡核中心的距离;L 为涡核尺度,切向速度在在 r/L=1处达到最大;为参数常量。物理模型如图 1 所示,一个波长为 的平面声波在传播过程中穿越一Gaussian 涡,左右设置为吸收层,避免反射的影响,声场计算域是边长为 10 10 的正方形区域。图 1 涡旋声散射

16、模型示意Fig.1 Schematic diagram of the physical model1.3 模型验证 为了考察本文方法的可行性和计算结果的准确性,选取 Colonius9算例结果并进行对比,如图 2 所示,图 2 中 prms为时域散射声压的均方根值:prms=p2sdt/T=(p-pi)2dt/T(5)0521第 8 期荆晨轩,等:水下振荡涡流场声散射与调制特性式中:ps为散射声压;p 为总声压;pi为入射声压。由图 2 可知,本文方法得到的结果与文献中的计算结果吻合度很高,证明了本文方法的有效性。图 2 涡旋声散射模型验证Fig.2 Verification of vorte

17、x acoustic scattering mode假定 Gaussian 涡流速度结构产生周期性振荡:vr=0v=2r(1-exp(-(r/L)2)cos(2f1t)(6)水声学问题不同于空气声学,一般马赫数较小,且自然流场多为长周期波,因此选取流场振荡频率与入射声波频率之比较小18。马赫数 Ma 为 0.01,入射声波频率 f0为 375 Hz,波长为,周期为 T0,流场振荡频率 f1为 15 Hz,周期 T1为 25T0,流场振荡周期 T1远大于声场周期 T0,涡流尺度 L 为/4,观测点选取在距涡流中心 2.5 处。在此基础上,本文通过声散射指向性曲线进行时间无关性和网格无关性验证。为

18、了避免截取信号时间长短对分析结果的影响,进行时间无关性验证,分别选取时间长短为 4T1、6T1、8T1进行求解,如图 3所示。可以看到计算时间长短达到 4T1及以上时,计算得到的散射指向性具有良好的一致性,分析选取的时间长度为 4T1。对应网格在一个波长范围内分别设置 12、18、24 个节点求解,进行网格无关性验证,如图 4 所示。可以看到每波长网格数达到 12 个及以上时,计算得到的声散射指向性曲线具有良好的一致性,为了兼顾计算量和散射声压云图分辨率,下面进行分析选取的网格数为每波长 18 个网格。2 数值计算结果与分析2.1 振荡涡流对声场的散射和调制 根据图 3、4 可知声场最大散射方

19、向大约处于30处,因此选择方向为 30处的时域声信号,分析涡流有无振荡时对声场的影响。观测点声压随时间变化曲线如图 5、6 所示。通过 Hilbert 变换提取声压信号包络,进行包络谱分析,如图 7 所示,包络线谱与涡流振荡频率 f1一致。对观测点声压信号进行频谱分析,如图 8 所示,涡流场振荡时声压信号频率变为 f0、f0f1、f02f1,涡流场无振荡时声压信号频率依然为入射声波频率 f0。可以看出声波经过振荡涡流场后,涡流场的振荡频率调制到声波频率上,本文将这种现象称为流场与声场发生调制作用,而涡流无振荡时,不会发生调制作用。图 3 时间无关性验证Fig.3 Time independen

20、t verification图 4 网格无关性验证Fig.4 Grid independence verification图 5 振荡涡流场观测点声压时域图Fig.5 Time domain diagram of sound pressure at obser-vation point of oscillating eddy current field图 6 无振荡涡流场观测点声压时域图Fig.6 Time domain diagram of sound pressure at obser-vation point of non-oscillating eddy current field15

21、21哈 尔 滨 工 程 大 学 学 报第 44 卷图 7 观测点声压包络谱图Fig.7Sound pressure envelope spectrum of observation point涡流场振荡时散射声压不同时刻云图如图 9 所示,声波在运动流体中传播受到流场速度和速度梯度的影响,从而产生散射声压1,19,产生调制作用的主要原因是涡流场发生振荡后,由于涡流场速度发 生周期性变化,从而引起散射声压产生对应频率的周期性变化,使得入射声波频率调制上涡流场的振荡频率。图 8 观测点声压频谱Fig.8 Spectrogram of sound pressure at observation po

22、int图 9 不同时刻振荡涡流场散射声压云图Fig.9 Cloud map of scattered sound pressure of oscillating eddy current field at different times 为了更为直观地评价流声耦合作用,分别采用散射截面和调制深度描述流场对声场的影响程度。采用散射截面 描述声场指向性及散射强度特性13:=120p2rms(r,)p2i(r,)rd(7)式中:为观测角度;r 为观测半径。以涡流中心为圆心,在圆周上对散射有效声压与入射声压之比的平方进行积分,为用来衡量流场对声波散射作用的强弱。采用调制深度 m20描述声场某点处的包络

23、特征,用来衡量涡流场对声场的调制作用强弱,涡流场无振荡时,声场不存在调制深度。调制深度 m 为:m=rmax-rminpmax+pmin 100%(8)式中:rmax为已调制波的最大振幅;rmin为最小振幅;pmax为载波最大振幅;pmin为最小振幅。通过改变涡流场特征分别改变马赫数 Ma、涡流尺度与声波波长比(下面称尺度波长比)L/、涡流振 荡 频 率 与 入 射 声 波 频 率 比(下 面 称 频 率比)f1/f0,研究上述涡流场参数对声场散射和调制作用的影响。2.2 马赫数影响分析 在低马赫数条件下,对 Ma 分别为 0.001 25、0.002 5、0.005、0.01、0.02 时的

24、流声耦合问题进行计算,不同马赫数流场声散射指向性随方位变化关系如图 10 所示,对于所计算范围内,随着马赫数增大,不同马赫数下的声散射指向性曲线形状基本一致。散射截面和调制深度随马赫数变化规律曲线分别如图 11、12 所示,可以看出在低马赫数条件下(Ma0.05),对于所计算的马赫数范围,散射截面和调制深度随着马赫数的增加而增加,散射截面与马赫数具有接近平方的关系,调制深度与马赫数具有接近正比的关系。2.3 尺度波长比影响分析 研究不同尺度波长比 L/条件下振荡涡流场对声场散射和调制作用的影响。对尺度波长比 L/分别为 0.062 5、0.125、0.25、0.5、1 时的流声耦合2521第

25、8 期荆晨轩,等:水下振荡涡流场声散射与调制特性问题进行计算,不同尺度波长比声散射指向性随方位变化关系如图 13 所示,对于所计算范围内,当尺度波长比较小时,指向性图存在明显的主瓣和旁瓣,随着尺度波长比的增大,旁瓣逐渐消失,主瓣逐渐平滑,该规律与刚性球体声散射类似15,21。在尺度波长比较小时,流场速度梯度较大,因而产生较强的散射指向性,而尺度波长比较大时,流场速度梯度较小,散射指向性减弱,不同角度声压变化更为平缓。图 10 不同马赫数声散射指向性Fig.10Directivity diagrams of acoustic scattering at dif-ferent Mach numbe

26、rs图 11 散射截面随马赫数变化规律Fig.11 The variations of scattering cross section with Mach number图 12 调制深度随马赫数变化规律Fig.12 Variations of modulation depth with Mach number散射截面和调制深度随尺度波长比变化规律曲线分别如图 14、15 所示。在所计算尺度波长比范围内,散射截面和调制深度会随着尺度波长比的增加而增加。散射截面与尺度波长比具有接近平方的关系,调制深度与尺度波长比具有接近正比的关系。图 13 不同尺度波长比声散射指向性图Fig.13Sound s

27、cattering directivity maps of wavelength ratio at different scales图 14 散射截面随尺度波长比变化规律Fig.14Variation of scattering cross section with scale-wavelength ratio图 15 调制深度随尺度波长比变化规律Fig.15Variation of modulation depth with scale-wave-length ratio2.4 频率比影响分析 研究不同频率比 f1/f0条件下涡流场对声场散射和调制作用的影响。对频率比 f1/f0分别为0.0

28、08、0.016、0.04、0.08、0.2 时的流声耦合问题进行计算,不同频率比声散射指向性随方位变化关系如图 16 所示,对于所计算范围内,可以看出散射指向性随着频率比的改变的规律,同尺度波长比具有一定相似性,即随着频率比的降低,旁瓣逐渐减少,主瓣逐渐平滑,不同之处在于频率比降低到一定程度之后,散射指向性几乎不再随频率比变化而变化。3521哈 尔 滨 工 程 大 学 学 报第 44 卷图 16 不同频率比声散射指向性Fig.16Directivity diagram of sound scattering at differ-ent frequency ratios对于所计算范围内,不同频

29、率比条件下的流场对声场的影响存在一极限值,在此将对应的频率比称为极限频率比(f1/f0)lim,高于极限频率比时,散射截面和调制深度随频率比减小而增大,流场对声场的影响随频率比减小而增强,降低到极限频率比后,散射截面和调制深度几乎不随频率比变化而变化。为了进一步研究不同频率比时的声场变化规律,分析了频率比 f1/f0在0.0040.2 散射截面和调制深度随频率比的变化情况,分别如图 17、18 所示。散射截面随频率比的变化情况呈现 3 个阶段的变化,如图中、,极限频率比(f1/f0)lim约为0.01,在第阶段,随着频率比降低,散射截面迅速增大。逐渐接近极限频率比时,进入阶段,此时散射截面随着

30、频率比降低依然增大,但增大幅度变缓。超过极限频率比后,在阶段,频率比的改变对散射声场几乎没有影响。图 17 散射截面随频率比变化规律Fig.17 The variation law of scattering cross section with frequency ratio2.5 极限频率比时声散射指向性分析 达到极限频率比后,对涡流场声散射指向性进行分析。涡流场有无振荡时声散射指向性随方位变化关系如图 19 所示,其中频率比为 0 时代表涡流场无振荡时的声场计算结果,圆圈表示无振荡时散射声压乘以系数 2/2 的计算结果。当达到极限频率比时,将流场无振荡时散射声压指向性乘以系数2/2 后,

31、与达到极限频率比后的散射声压指向性图吻合,上述现象可以解释为:若将无振荡时的散射声压视为直流信号,当涡流场发生振荡时,散射声压随流场周期性变化而产生周期性变化,从而成为正弦信号,上述现象符合正弦信号的有效值为直流信号幅值的 2/2 规律,说明此时流场对声场的调制作用主要是由于流场的周期性变化引起散射声压周期性变化导致的。而频率比较高时的变化规律,说明振荡涡流场与声场还会发生其他耦合作用,使散射声场结构产生细微改变,实际中流场振荡频率与入射声波频率之比远小于极限频率比,因此在此不再进一步探究频率比较高时的现象。考虑海洋中声波衰减问题,可以选择合适的频率比声信号进行流场探测。图 18 调制深度随频

32、率比变化规律Fig.18The variation law of modulation depth with fre-quency ratio图 19 涡流场有无振荡时声散射指向性比较Fig.19 Comparison of sound scattering directivity with or without oscillation of the eddy current field3 结论 1)随着马赫数增大,散射指向性曲线形状基本一致,随着尺度波长比的增大,散射指向性图旁瓣逐渐消失,主瓣逐渐平滑。2)在所计算范围内,散射截面和调制深度会随着马赫数和尺度波长比的增加而增加,散射截面与马赫

33、数和尺度波长比具有接近平方的关系,调制深4521第 8 期荆晨轩,等:水下振荡涡流场声散射与调制特性度与马赫数和尺度波长比具有接近正比的关系。3)不同频率比条件下的涡流场对声场的影响存在一极限频率比,高于极限频率比时,流场对声场的散射和调制作用随频率比降低而增强,散射指向性图旁瓣逐渐消失,随着频率比的降低,散射和调制作用增强速度变小,降低到极限频率比后,频率比的改变对声场的散射和调制作用几乎没有影响。参考文献:1 AUREGAN Y,MAUREL A,PAGNEUX V,et al.Sound-Flow Interactions M.Softcover reprint of hardcover

34、 1st ed.2002 edition:Springer,2010.2 程建春.声学原理M.北京:科学出版社,2012:384.3 ZDRAVKOVICH M M.Flow around circular cylinders:Volume 2:Applications M.Oxford university press,1997:27.4 董双岭,吴颂平.圆柱绕流尾迹流态特征和涡的演化过程分析J.北京航空航天大学学报,2009,35(8):933-937.DONG Shuangling,WU Songping.Analysis of wake pat-tern and vortex evol

35、ution in flow past circular cylinderJ.Journal of Beijing University of Aeronautics and Astronau-tics,2009,35(8):933-937.5 CATALANO P,WANG Meng,IACCARINO G,et al.Nu-merical simulation of the flow around a circular cylinder at high Reynolds numbersJ.International journal of heat and fluid flow,2003,24

36、(4):463-469.6 ORR M H,HAURY L R,WIEBE P H,et al.Backscatter of high-frequency(200 kHz)acoustic wavefields from o-cean turbulenceJ.The journal of the acoustical society of America,2000,108(4):1595-1601.7 PROKHOROV V E,CHASHECHKIN Y D.Experimental study of the scattering of ultrasound from wake flows

37、in a linearly stratified fluidJ.Acoustical physics,2005,51(1):95-104.8 FERZIGER J H.Low-frequency acoustic scattering from a trailing vortexJ.The journal of the acoustical society of America,1974,56(6):1705-1707.9 COLONIUS T,LELE S K,MOIN P.The scattering of sound waves by a vortex:numerical simulat

38、ions and analyt-ical solutionsJ.Journal of fluid mechanics,1994,260:271-298.10 CHEINET S,EHRHARDT L,JUV D,et al.Unified modeling of turbulence effects on sound propagationJ.The journal of the acoustical society of America,2012,132(4):2198-2209.11 KE Guoyi,LI Wen,ZHENG Z C.Vortex scattering effects o

39、n acoustic wave propagationC/Proceedings of the 21st AIAA/CEAS Aeroacoustics Conference.Reston,Virginia:AIAA,2015:AIAA2015-3267.12 CLAIR V,GABARD G.Spectral broadening of acoustic waves by convected vorticesJ.Journal of fluid mechan-ics,2018,841:50-80.13 马瑞轩,王益民,张树海,等.旋涡声散射特性的尺度效应 数 值 研 究 J.物 理 学 报,

40、2021,70(10):104301.MA Ruixuan,WANG Yimin,ZHANG Shuhai,et al.Nu-merical investigation of scale effect on acoustic scattering by vortex J.Acta physica sinica,2021,70(10):104301.14 王益民,马瑞轩,武从海,等.旋涡声散射的空间尺度特性 数 值 研 究 J.物 理 学 报,2021,70(19):194302.WANG Yimin,MA Ruixuan,WU Conghai,et al.Nu-merical study on

41、 spatial scale characteristics of sound scat-tering by a static isentropic vortexJ.Acta physica sini-ca,2021,70(19):194302.15 杜浩,熊鳌魁,张咏鸥.水下涡流场前向声散射特性分析J.声学学报,2020,45(1):55-61.DU Hao,XIONG Aokui,ZHANG Yongou.Forward sound scattering from an underwater vortex J.Acta acustica,2020,45(1):55-61.16 PIERCE

42、 A D.Acoustics:An introduction to its physical principles and applicationsM.Cham:Springer Interna-tional Publishing,2019.17 吕宏强,张涛,孙强,等.间断伽辽金方法在可压缩流数值模拟中的应用研究综述J.空气动力学学报,2017,35(4):455-471.LYU Hongqiang,ZHANG Tao,SUN Qiang,et al.Appli-cations of discontinuous Galerkin method in numerical simulations

43、 of compressible flows:a review J.Acta aerodynamica sinica,2017,35(4):455-471.18 侍茂宗.物理海洋学M.济南:山东教育出版社,2004:255.19 BRILLANT G,CHILL F,PINTON J.Transmission of sound through a single vortexJ.The European physical journal B-condensed matter and complex systems,2004,37(2):229-239.20 任晶.宽带相控参量阵的数字系统实现D.

44、成都:电子科技大学,2021.REN Jing.Digital system realization of broadband phased parametric arrayD.Chengdu:University of Electronic Science and Technology of China,2021.21 马黎黎,王仁乾.海洋波导中刚性球及旋转椭球前向散射场时频特征的畸变J.声学学报,2014,39(4):407-416.MA Lili,WANG Renqian.Distortion of time-frequency characteristics of forward sc

45、attering from a rigid sphere/prolate-spheroid in an oceanic waveguide J.Acta acustica,2014,39(4):407-416.本文引用格式:荆晨轩,时胜国,杨德森.水下振荡涡流场声散射与调制特性J.哈尔滨工程大学学报,2023,44(8):1249-1255.JING Chenxuan,SHI Shengguo,YANG Desen.Acoustic scattering and modulation characteristics of underwater oscillating flow fieldsJ.Journal of Harbin Engineering University,2023,44(8):1249-1255.5521

展开阅读全文
相似文档                                   自信AI助手自信AI助手
猜你喜欢                                   自信AI导航自信AI导航
搜索标签

当前位置:首页 > 学术论文 > 论文指导/设计

移动网页_全站_页脚广告1

关于我们      便捷服务       自信AI       AI导航        获赠5币

©2010-2024 宁波自信网络信息技术有限公司  版权所有

客服电话:4008-655-100  投诉/维权电话:4009-655-100

gongan.png浙公网安备33021202000488号   

icp.png浙ICP备2021020529号-1  |  浙B2-20240490  

关注我们 :gzh.png    weibo.png    LOFTER.png 

客服