999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

混凝土中柱形裝藥的爆炸破壞分區 及應力波衰減規律

2025-09-15 00:00:00周鑫馮彬陳力
爆炸與沖擊 2025年9期

中圖分類號:0382 國標學科代碼:13035 文獻標志碼:A

Study on failure zones and attenuation law of stress waves in concrete induced by cylindrical charge explosion

ZHOU Xin,FENG Bin,CHENLi (Engineering Research Center ofSafety and Protection of Explosion amp; Impact of Ministry of Education, Southeast University,Nanjing211189,Jiangsu,China)

Abstract:Inblast-resistantstructraldesignforconventional weapos,previoussudiesonblas-iducedstresswaves nsolid media have predominantly focusedonsoilandrockmedia(i.e.,groundshock isues),whereas researchonthe propagationand atenuation lawsof stresswaves inconcrete remains relativelylimited.Based onthe Karagozian and Case concrete(KCC) constitutive model in conjunction with the multi-material ALE(MMALE)algorithm, the propagation laws of stress waves in concrete induced bycylindrical chargeexplosion were numericallyinvestigated.Firstly,the applicabilityofthecostitutive modelparameters andnumerical algorithm were validatedbycomparing theresults with theexisting experiments. Subsequently,thepeak stresswas employed asacriterion todelineate theexplosive damage zones in theconcrete surrounding thecharge.Aditionally,theatenuationlawsofexplosionstress waves ineachdamagezone were disussed.Finallyteefect ofburialdepthwastakenintofurtherconsidered,andaformulaforcalculatingthepeak stressinconcreteinduced by cylindrical chargeexplosion was established.It was foundthattheattenuation patesof blast-induced stresswaves differ significantlyineach explosion failure zone.Thestress wavesin the near-field zone(quasi-fluid and crushing zones) demonstrates amorerapid atenuationrate compared tothat inthe mid-field zone (transition and fracture zones).Furthermore, an increase intheaspectratioof thecylindricalchargeleads toanaccelerationintheatenuationof thenormal peak stre. Moreover,theestablished formula forcalculating the peak stressof blast-induced stresswaves enablesaccurateandrapid determinationofthe normal peak stressgenerated bycylindricalcharges withvaryinggeometries andburialdepths,whichcan be served as a valuable reference for blast-resistant design of concrete structures.

Keywords: concrete; cylindrical charges; failure zones; explosion stress waves; peak stress

鋼筋混凝土結構是抗爆結構中最主要的結構類型,混凝土材料廣泛應用于各類民用建筑、防護工事及遮彈層等防護屏障。以往針對常規武器的抗爆結構設計中,對固體介質中爆炸應力波的研究多針對土壤和巖石介質(即地沖擊問題),對混凝土中爆炸應力波的傳播與衰減規律研究較少。一方面是因為以往常規武器的精度較低,很難直接命中結構目標,地面結構的防護設計研究多基于結構外部爆炸場景;另一方面,對地下結構來說,以往常規武器的侵徹能力一般難以突破上方的覆土層和遮彈層而直接命中結構,因此地下結構抗爆通常僅考慮武器在結構外圍介質中爆炸的情況[1]。然而,隨著精確制導技術、深鉆地戰斗部、超高音速武器等先進技術的發展,戰斗部直接命中結構的概率大大增加。此外,當前對結構類目標毀傷評估的需求日增,相對于僅需考慮一般偏保守情況的結構設計計算,毀傷評估計算更需要明確各類極端加載情況下結構的破壞機理和毀傷程度。因此,對混凝土介質中常規武器爆炸應力波效應的研究在當今顯得尤為重要。

常規武器戰斗部的裝藥形狀多為圓柱形,因此,待研究問題可轉化為柱形裝藥在混凝土介質中的爆炸應力波問題。目前主要的研究手段為試驗和數值模擬。早期試驗研究方面,以柱形裝藥在空氣自由場中的爆炸試驗居多,重點研究了不同長徑比、不同起爆方式以及不同方位處爆炸波的荷載(超壓峰值和沖量)特征。Plooster等[2]開展了一系列柱形裝藥空氣自由場爆炸試驗,測得了不同裝藥長徑比、方位角和比例距離處的爆炸荷載。Ismail等3]發現柱形裝藥產生的爆炸波非常復雜,超壓時程曲線存在多重峰值。 wu 等[4通過開展空氣自由場爆炸試驗,對比了球型裝藥和柱形裝藥作用于鋼筋混凝土板上的爆炸荷載,試驗結果表明,柱形裝藥正下方處的反射超壓和沖量遠大于相同質量的球型裝藥。Shi等[5]通過爆炸試驗對比分析了裝藥形狀對爆炸荷載的影響,結果表明,裝藥形狀對爆炸近區的爆炸荷載影響較大,當比例距離大于 5.0m/kg1/3 時,其影響可以忽略不計。對混凝土介質中柱形裝藥爆炸試驗的研究相對較少,黃家蓉等測量了混凝土介質中柱形裝藥的爆炸應力波,并對試驗工況進行了數值模擬,模擬結果與試驗結果吻合較好。Gebbeken等[7]通過開展柱形裝藥接觸爆炸試驗,獲取了混凝土狀態方程參數。

隨著數值計算方法的發展與計算效率的提升,數值模擬已成為研究爆炸應力波問題的重要工具。Sherkar等[8通過數值模擬分析了柱形裝藥形狀和起爆點對空氣自由場中爆炸應力波的影響,并認為裝藥形狀對入射超壓峰值和沖量產生影響的臨界比例距離分別為 3.69,2.74m/kg1/3 ,而起爆點對爆炸荷載產生影響的臨界比例距離為 3m/kg1/3 ,超過該距離時,可忽略其對爆炸荷載的影響。Xiao等[同樣也開展了相關的研究工作,結果表明,柱形裝藥一端起爆產生的超壓峰值和沖量最大,分別是等當量中心起爆的球型裝藥的約4.5倍(超壓峰值)和約4.0倍(最大沖量)。Gao等[基于數值模擬研究了長徑比和方位對中心起爆柱形裝藥產生的空氣自由場爆炸荷載的影響,并建立了爆炸超壓峰值和沖量的實用化計算公式。在此基礎上,王明濤等[1]進一步開展了柱形裝藥空中爆炸數值模擬研究,提出了柱形裝藥空中爆炸入射和反射沖擊波荷載的計算方法。Gao等[12]開展了混凝土中柱形裝藥爆炸試驗,并結合數值模擬,在先前提出的混凝土中球形裝藥爆炸應力波的峰值應力計算公式[13]基礎上,建立了混凝土中柱形裝藥爆炸應力波峰值應力的實用化計算公式,但其并未考慮爆炸應力波峰值應力衰減速度隨傳播距離的變化。楊耀宗等[4開展了混凝土中帶殼柱形裝藥爆炸應力波衰減規律的數值模擬研究,建立了帶殼柱形裝藥峰值應力的計算公式,其適用比例爆距為 0.30~1.0m/kg1/3 。

對于柱形裝藥爆炸應力波問題的已有研究以空氣介質中爆炸荷載的分布特征居多,對混凝土介質中的爆炸應力波傳播規律的研究相對較少。與空氣沖擊波相比,混凝土類介質中的應力波傳播與衰減規律更為復雜,其與介質受力特征及介質狀態密切相關[15]。炸藥起爆之后,在混凝土介質中產生應力波并向外傳播。在傳播過程中應力波不斷衰減,由初始的強間斷沖擊波衰減為彈塑性波[10,混凝土介質狀態也由高應力擬流體狀態向低應力固體彈塑性狀態轉變[12-13]。然而,現有的混凝土中應力波衰減規律的研究鮮有考慮介質受力特征和介質狀態對應力波衰減特征的影響;此外,在建立爆炸應力波峰值的實用化計算公式時,已有研究多采用單一衰減指數來統一描述擬流體狀態和彈塑性狀態的混凝土中的應力波衰減規律,其合理性和準確性有待商榷,需進一步探究。

本文中基于已有的柱形裝藥接觸爆炸試驗[,,利用LS-DYNA有限元軟件開展數值模擬研究,依據爆炸應力波特征對柱形裝藥周圍混凝土介質破壞分區進行劃分,分析不同破壞分區中的爆炸應力波衰減規律,并綜合考慮柱形裝藥長徑比、破壞分區(介質狀態與受力特征)以及裝藥埋深對峰值應力的影響,提出柱形裝藥法向峰值應力實用化計算公式。研究結果可為混凝土中爆炸應力波分析及防護工程抗爆設計提供參考。

1數值模型及驗證

基于Gebbeken等[7開展的混凝土靶板接觸爆炸試驗,采用LS-DYNA軟件建立精細化數值模型,并對數值模擬方法和材料模型參數進行驗證。

1.1 有限元模型

試驗工況如圖1所示,方形混凝土靶板邊長為 100cm ,高度為 30cm ,混凝土水灰比為0.45,單軸抗壓強度為 51.2MPa ,密度為 2.35g/cm3 ;靶板內部共布置6個壓力傳感器,傳感器分3層布置,同一層傳感器間隔 80mm ,上下層間隔 20mm ,首層傳感器距靶體上表面 40mm 。圓柱形PETN炸藥直徑與高度均為 75mm ,裝藥質量為 500g. ,起爆點位于裝藥尾部中心正下方 10mm 處。

圖1試驗布置(單位:mm)

Fig.1Schematic diagram of the experiment (unit: mm)

考慮到模型具有較好的對稱性,為提高計算效率,采用二維軸對稱方法建立有限元模型,模型尺寸與試驗一致,如圖2所示?;炷涟械酌婧蛡让娌捎米杂蛇吔纾諝庥蛏媳砻婧蛡让娌捎猛干溥吔纭2捎枚辔镔|ALE算法(MMALE)進行數值模擬,其中炸藥、空氣和混凝土均采用ALE 網格;相較于Lagrange算法和Euler算法,該方法既可以避免網格畸變問題,又可以較好地追蹤物質界面,已廣泛應用于侵徹爆炸等問題的數值模擬[16-19]。

圖2有限元模型

Fig.2 Finite elementmodel

1.2 材料參數

采用Karagozian and Case concrete (KCC)本構模型[20]作為混凝土的材料模型,該模型引入了3個獨立的強度面,即初始強度面、最大強度面和殘余強度面,并綜合考慮了材料損傷、應變率和靜水壓力對屈服應力的影響,可以較好地捕捉復雜應力狀態下的混凝土行為,被廣泛應用于混凝土類材料在爆炸荷載作用下的破壞效應分析[21-23]。Kong等[24]發現KCC模型自動生成的狀態方程曲線與試驗數據較為接近,但是強度面參數僅適用于低靜水壓,對于侵徹爆炸這類高靜水壓問題并不適用,并基于大量混凝土三軸試驗數據重新確定了強度面參數。因此,本文中采用Kong等[24]改進的強度面參數,狀態方程參數由KCC模型自動生成算法獲得。

PETN炸藥材料模型采用*MAT_HIGHEXPLOSIVE BURN,狀態方程采用Jones-Wilkins-Lee(JWL)狀態方程:

式中: p 為爆轟產物壓力, A,B,ω,R1,R2 為狀態方程參數, V 為相對體積, E 為單位體積內能。炸藥材料參數采用參考Xiao等[25]提供的參數,如表 1[25] 所示。

表1炸藥本構模型及狀態方程參數[25]

Table 1 Parameters of constitutive model and EOS for explosive[25]

空氣采用*MATNULL材料模型,狀態方程采用多項式狀態方程:

p=C0+C1μ+C2μ2+C3μ3+(C4+C5μ+C6μ2)E0

式中: p 為空氣壓力; C0?C1?C2?C3?C4?C5?C6 為自定義系數; μ=ρ/ρ0-1,ρ/ρ0 為當前密度與參考密度的比值; E0 為單位參考體積的初始能量;空氣材料參數如表 2[26] 所示。

網格尺寸對數值模擬預測結果有顯著影響,利用上述有限元模型和材料參數對網格收斂性進行分析。圖3(a)和(b)分別為不同網格尺寸時,炸藥正下方 0.2m 處混凝土中的壓力時程曲線和峰值壓力,可以看出當網格尺寸小于 3.0mm 時,壓力時程曲線和峰值壓力均開始收斂。因此,后續數值模擬的網格尺寸均設置為 3.0mm 。

Fig.3Meshconvergenceanalysis

1.3數值模擬結果驗證

圖3網格收斂性分析

Gebbeken等[7]開展了3組相同工況的接觸爆炸試驗,各測點的應力時程曲線如圖4所示,其中\" gauge 1-3′′ 表示第1組試驗、測點3所測的試驗數據??梢钥闯?,試驗數據存在較大的離散性,例如,3組爆炸試驗中的第1層測點處的峰值應力最大值為 16.06GPa (gauge2-4處),而最小值僅為 1.55GPa (gauge 3-4處)。Xiao 等[25]研究表明,在爆炸近區,PETN炸藥的等效TNT當量系數為3.31,因此,試驗中 500gPETN 炸藥與 1655g 相同形狀的 TNT炸藥威力相同?;贖opkinson 定律[27],測點1、4處的峰值應力與相同形狀的TNT炸藥在比例距離為 0.074m/kg1/3 時的峰值應力基本相同。Tu等28開展了類似工況的研究,并提出了TNT接觸爆炸計算模型,確定了柱形裝藥(長徑比為1.25)在C35混凝土中的爆炸應力波峰值應力,其中比例距離為 0.074m/kg1/3 處的峰值應力為 8.22GPa 。對比試驗數據可以看出:第2組試驗中的測點1、4處的峰值應力分別為14.72、16.06GPa,明顯偏大;而第3組試驗數據分別為3.26,1.55GPa ,明顯偏小;第1組試驗數據分別為11.16、9.21GPa,與上述計算值相近,較為合理。

圖4不同測點處的應力時程曲線

Fig.4Stress-time curves at different measurement points

試驗數據和數值模擬結果的對比如圖4所示,可以看出:整體上,數值模擬的波形與試驗數據吻合較好;雖然試驗數據與數值模擬結果存在一些差異,但是數值模擬結果處于試驗數據范圍內。進一步可以發現,數值模擬計算的峰值應力與第1組試驗數據吻合較好。由于同層傳感器的比例距離相同,因此,取同層傳感器的峰值應力平均值代表該比例距離處的峰值應力,進而計算出平均誤差,其值均小于 18% ,如表3所示。此外,需要注意的是,數值模擬得到的應力時程曲線下降段呈現振蕩現象。主要原因是:一方面,接觸爆炸條件下,該區域直接地沖擊與感生地沖擊會發生耦合作用[29;另一方面,爆炸結束前,混凝土中始終會存在壓縮波與稀疏波的相互作用[13.30],從而引起波形的震蕩。與數值模擬結果相比,試驗數據曲線下降段振蕩頻率較弱,這可能是因為傳感器采樣頻率較低,未捕捉到其余峰值應力。由此可見,采用的KCC本構模型和MMALE算法可以較為準確地描述混凝土介質中爆炸應力波的傳播規律。

表3第1組試驗中各測點峰值應力的試驗結果與數值模擬結果的對比

Table 3 Comparison of stress peak of tests with that of numerical simulation in the first group

注:平均誤差=(試驗平均值-數值模擬值)/平均值 ×100% 門

2柱形裝藥作用下混凝土破壞分區劃分及應力波衰減規律

混凝土中柱形裝藥爆炸引起的應力波與裝藥形狀(長徑比)、裝藥類型、埋置深度、比例距離等因素相關,首先對混凝土自由場中的爆炸應力波衰減機理進行分析,之后再進一步考慮裝藥長徑比和埋置深度的影響。常規武器的威力通常以TNT當量來衡量,其參數已被廣泛應用于防護結構抗爆設計,其他裝藥類型當量可通過TNT等效系數進行換算,因此,后續數值模型中裝藥采用TNT,參數見文獻[31]?;谏鲜鲵炞C的本構模型和數值算法,建立了長徑比為1的柱形裝藥封閉爆炸的數值模型,起爆點位于柱形裝藥頂部中心點,裝藥當量和網格信息與上述模型一致,有限元模型如圖5所示。

2.1柱形裝藥作用下混凝土破壞分區劃分

炸藥爆炸之后,裝藥周圍介質在爆炸應力波的劇烈作用下發生汽化、粉碎、破裂等現象[32-33]。依據介質的破壞程度和損傷狀態,裝藥周圍介質可以劃分為不同的區域。李守巨等[34、錢七虎等3和王明洋等[3]將裝藥周圍的巖石劃分為粉碎區、破裂區和彈性區,之后,張志呈等[37]和冷振東等[38]基于巖體破壞特征將破碎區細分為剪切破碎區和裂隙區。近兩年,Mandal等[39]和Gao等[12-13]將封閉空間下裝藥周圍混凝土介質破壞分為5個區域,即近流體區、壓碎區、過渡區、破裂區和彈性區,如圖6所示。

圖5有限元模型 Fig.5Finite elementmodel

圖6裝藥周圍介質破壞分區

Fig.6Failure zones of the surrounding medium around the charge

本文中借鑒Mandal等[39]和Gao等[12-13]提出的方法對混凝土中爆炸破壞分區進行劃分,近流體區混凝土介質中靜水壓力遠大于其剪切強度,呈塑性流動狀態。壓碎區混凝土的強度效應逐漸顯現,但是由于所受作用力遠大于其強度,混凝土被擠壓破碎[13]。Amelsfort 等[40]根據試驗數據建立了混凝土壓碎強度 σcr 與混凝土抗壓強度 fc 關系式:

基于上述關系式確定文中混凝土的壓碎強度范圍為 0.74~1.86GPa ,當爆炸應力波峰值應力處于該范圍內時,混凝土被完全壓碎,形成壓碎區。當峰值應力超過 1.86GPa 時,混凝土呈流體狀態,形成近流體區。利用數值模擬計算結果確定近流體區和壓碎區范圍為 0.08~0.13m/kg1/3 和 0.13~0.17m/kg1/3

過渡區是連接破碎區和破裂區的中間區域,由于泊松效應,該區域的混凝土受到徑向和環向壓應力的同時作用,發生剪切破壞,裂隙呈網狀分布。破裂區的混凝土環向約束消失,環向拉伸應力開始起主導作用,當環向拉應力大于混凝土抗拉強度時,形成徑向裂隙;隨著應力波繼續向外傳播,強度進一步下降,當環向拉伸應力低于混凝土抗拉強度時,裂隙停止擴展,裂隙之外為彈性區。對比過渡區和破裂區混凝土的受力特征,不難發現:過渡區和破碎區混凝土的環向應力存在顯著差異,前者主要承受環向壓應力作用,而后者主要承受環向拉應力,且拉應力大于其抗拉強度,因此,可以基于這一特征來劃分過渡區和破碎區。

圖7為比例距離 Z 為 0.32~1.07m/kg1/3 時,柱形裝藥正下方測點的環向應力 (σx) 時程曲線??梢钥闯觯罕壤嚯x為 0.32~0.55m/kg1/3 時,混凝土環向主導應力由壓應力(正值)轉變為拉應力(負值),且比例距離為 0.46m/kg1/3 時,環向拉應力峰值為 4.9MPa ,大于本文中混凝土抗拉強度( 4.5MPa[41] ),因此,本文中以比例距離 0.46m/kg1/3 作為過渡區與破裂區的分界線。

圖7環向應力時程曲線

Fig.7Circumferential stress-time curves

基于上述分析,長徑比為1的柱形TNT爆炸作用下,混凝土介質近流體區、壓碎區、過渡區和破裂區比例距離范圍分別為 0.08~0.13,0.13~0.17,0.17~0.46,0.46~1.0m/kg1/3 ,Gao等[12]計算的近流體區和壓碎區范圍均稍大于本文計算結果,主要因為柱形裝藥在近區的峰值應力衰減速度比球型裝藥更快[,42],混凝土的損傷云圖及破壞分區劃分如圖8所示。

圖8混凝土的損傷云圖及破壞分區劃分

Fig.8Damage contour and the division of failure zones for concrete targe

2.2柱形裝藥法向各破壞分區應力波衰減規律

圖9為不同比例距離處的柱形裝藥法向 (y 方向)應力時程曲線(以壓應力為正),可以看出,隨著比 例距離的增大,沖擊波迅速衰減為具有升壓時間的塑性波,波形逐漸變緩變長,峰值應力逐漸降低。比 例距離為 0.08~0.12m/kg1/3 時,爆炸應力波為強間斷沖擊波,升壓時間接近 0μs ,峰值應力呈現線性衰 減,由 5.69GPa 衰減為 2.27GPa ;比例距離為 0.18~0.30m/kg1/3 時,爆炸應力波升壓時間由 13.91μs 增至 39.40μs ,峰值應力由 0.64GPa 衰減為 0.22GPa ;比例距離為 0.40~0.60m/kg1/3 時,爆炸應力波升壓時間 由 50.35μs 增至 55.38μs ,峰值應力由 0.14GPa 衰減為 0.07GPa ,此時升壓時間增長和峰值應力衰減速率 均減緩;比例距離為 0.70~1.00m/kg1/3 時,升壓時間由 58.49μs 增至 59.01μs ,峰值應力由 0.06GPa 衰減 為 0.04GPa ,此范圍內塑性波進一步衰減為彈性波,升壓時間基本保持不變,由于柱(球)面波的幾何擴散 效應,峰值應力繼續衰減。

圖9柱形裝藥法向不同測點應力時程曲線

Fig.9 Normal stress-time curves at different measurement points for cylindrical charges

上述分析表明,近流體區和壓碎區的混凝土中應力波為強間斷沖擊波,而到了過渡區和破裂區,沖擊波已經衰減為塑性波,其波形的上升段與下降段均放緩。當比例距離達到 1.0m/kg1/3 時,峰值應力已小于混凝土的抗壓強度,此時爆炸應力波已衰減為彈性波。表4總結了裝藥長徑比為1時,混凝土自由場爆炸破壞分區邊界尺寸及介質受力特征。

表4混凝土自由場破壞分區邊界尺寸及介質受力特征

Table4 Concrete failure zone boundaries in a free field and stressed medium properties

基于Hopkinson相似律[27],固體介質中爆炸應力波峰值應力通常以冪指數形式給出:

式中: σm 為某一測點的峰值應力; k 為應力衰減系數; 為測點至裝藥中心的距離; W 為裝藥質量; Q/W1/3 (20為比例距離,用 Z 表示; n 為衰減指數,其值越大表示峰值應力衰減越快。

在防護工程抗爆設計中,峰值應力的最大值通常是設計者關注的重點,主要出現在常規武器法向[(y方向),因此后續研究中的峰值應力σm 默認指代法向峰值應力。雙對數坐標系下柱形裝藥正下方的峰值應力 σm 與比例距離Z關系如圖10所示,可以看出,近流體區與壓碎區的峰值應力衰減明顯快于過渡區和破裂區,這表明單一衰減指數難以準確描述各個破壞分區的峰值應力衰減規律。王明洋等[43]和吳祥云等[15]基于試驗和數值模擬研究了巖石類介質中爆炸地沖擊傳播規律,也發現裝藥近區 0~0.2m/kg1/3 和中遠區( 0.2~1.0m/kg1/3 的峰值應力衰減規律存在顯著差異。

值得注意的是,現有的混凝土介質中爆炸應力波峰值應力計算公式多數都采用了單一衰減指數,如李重情等[44]和""等[45]提出的球型裝藥峰值應力計算公式的衰減指數 n 為1.739(C50混凝土);高矗等[13]提出的球型裝藥峰值應力計算公式的衰減指數 n 為1.734,并在此基礎上建立了柱形裝藥的峰值應力計算公式的衰減指數為 2.38(0~0.24m/kg1/3"其余區間與球型裝藥一致)[12];楊耀宗等[14]提出的柱形裝藥峰值應力計算公式的衰減指數 n 為1.39(CF120混凝土)。圖10對比了現有計算公式結果與數值模擬結果,需要說明的是,由于衰減系數k 與介質類型相關4,為此本節只選取上述已有計算公式的衰減指數 n ,通過擬合數值模擬數據確定最優的衰減系數k。結果顯示, n 為1.734和1.739時,在 0.15~0.80m/kg1/3"范圍內,計算值與數值模擬結果吻合較好,而范圍外計算值明顯偏小,最大差值達 4.25GPa;n 為1.39時,在 0.24~1.0m/kg1/3"范圍內,計算值與數值模擬結果吻合較好,而在 0~0.24m/kg1/3"范圍內,計算值顯著低估峰值應力;相較于單一衰減指數的計算公式,兩段式[12]( n 為2.38和1.734)的計算公式更優,但是計算值與數值模擬結果的最大差值仍可達 2.97GPa ,這可能是由于未考慮衰減指數隨著傳播距離的變化所致。

為了分析衰減指數隨傳播距離的變化規律,圖11給出了各破壞分區測點的峰值應力與比例距離的關系,并通過最小二乘法擬合確定各破壞分區的峰值應力衰減系數 k 和衰減指數 n 。可以看出,擬合曲線與峰值應力衰減趨勢吻合良好,相關性系數 R2 均大于0.98;整體上,裝藥近區(近流體區和壓碎區)峰值應力衰減速度大于中遠區(過渡區和破裂區),且衰減系數 k 與衰減指數 n 呈負相關。此外,數值結果發現壓碎區衰減指數(3.224)大于近流體區(2.539),這可能是因為壓碎區混凝土介質由高應力擬流體狀態向固體塑性狀態轉變,造成阻力突變,加速了峰值應力的衰減[47]。

圖11各破壞分區豎向峰值應力與比例距離散點圖

Fig.11 Scatter plot of vertical peak stress versus scaled distance for each failure zone

2.3裝藥長徑比對峰值應力衰減規律的影響

已有學者研究表明裝藥長徑比I/d對峰值應力具有顯著影響[5.8.48]。為探究裝藥長徑比對峰值應力的影響規律,分別開展了長徑比為1、2、4、6和8的柱形裝藥在封閉空間下的數值模擬研究?;?.1節中爆炸破壞分區劃分方法以及2.2節中衰減系數 k 和指數 n 計算方法,分別計算出了長徑比2、4、6和8的柱形裝藥法線方向的破壞分區及衰減參數,如表5所示,可以看出:長徑比在 2~8 之間的柱形裝藥工況下的峰值應力衰減系數 k 隨著長徑比增大而遞減,而衰減指數 n 隨著長徑比增大而遞增。

表5不同長徑比的柱形裝藥各破壞分區參數

利用上述表格所列的區間范圍,確定破壞分區隨長徑比變化的規律,如圖12所示,(圖中的虛線為各個破壞分區的分界線,例如 Zca,fl 為空腔與近流體區的分界線)??梢钥闯觯弘S著柱形裝藥長徑比的增大,近流體區、壓碎區和過渡區逐漸變窄。圖13為不同長徑比柱形裝藥底部峰值應力 (p1) 和空腔區邊界處峰值應力 (p2) 的變化趨勢圖,可以發現:一方面,柱形裝藥長徑比增大使得裝藥底部的峰值應力提升,與文獻[28]觀點一致;另一方面,柱形裝藥長徑比增大導致峰值應力衰減系數變大,應力衰減加快[4],導致空腔區邊界處峰值應力 (p2) 隨裝藥長徑比的增大而減小,進一步導致其他破壞分區變窄。

圖12破壞分區范圍隨著長徑比的變化關系 Fig.12 Failure zone boundariesversus ratio of length to diameter

圖13裝藥底部及空腔壁處峰值應力變化趨勢 Fig.13Peak stress in guage1and2 versus ratioof length to diameter

圖14展示了不同長徑比柱形裝藥的峰值應力衰減規律,可以看出,隨著長徑比增大,裝藥正下方峰值應力衰減速度加快,且近流體區和壓碎區的衰減速度明顯快于其他破壞分區,表明裝藥近區峰值應力對裝藥長徑比更為敏感。圖15進一步給出了各破壞分區的衰減系數 k, 衰減指數 n 與長徑比l/d的關系,之后,采用最小二乘法擬合建立 k,n 與 l/d 之間的函數表達式??梢钥闯觯瘮当磉_式與數值模擬結果吻合較好,相關性系數 R2 均大于 0.98 。衰減系數 k 與柱形裝藥長徑比 l/d 呈負相關,其中在近流體區和壓碎區,兩者呈指數關系,而在過渡區和破裂區,兩者呈線性關系。文獻[14]中指出衰減系數 k 的取值與介質類型相關,這也說明了近流體區和壓碎區混凝土介質屬性相近。衰減指數 n 與長徑比 l/d 呈線性正相關,并且近流體區和壓碎區的衰減指數 n 大于過渡區和破裂區,這也說明了裝藥近區的峰值應力衰減更快。

圖14不同長徑比條件下峰值應力與比例距離散點圖 Fig.14Scatter plot of peak stress versus scaled distance under different aspect ratios

3柱形裝藥峰值應力實用化計算公式

3.1不同長徑比的柱形裝藥峰值應力計算公式

基于上述分析,當比例距離大于 1.00mkg1/3 時,柱形裝藥爆炸波已衰減為彈性波。本文中重點關注比例距離小于 1.00m/kg1/3 的混凝土中的峰值應力分布,此范圍內混凝土介質可以分為近流體區、壓碎區、過渡區和破裂區。鑒于近流體區和壓碎區范圍較窄,并隨柱形裝藥長徑比增大而逐漸縮小,且區間內混凝土介質受力特征相似,均以承受靜水壓為主。參考巖石介質的相關研究[34.38]將兩者合并,統稱為粉碎區?;跀抵的M結果,分段擬合得到 k 和 n 的表達式如下:

式中: Zi,j 表示破壞分區 i 與破壞分區 j 分界線,如 Zca,fl 表示空腔與近流體區的分界線,不同長徑比柱形裝藥的 Zi,j 值見表5;通過等式 (4)~(6) ,可快速地計算出柱形裝藥峰值應力。

不同長徑比的柱形裝藥爆炸應力波峰值應力數值模擬結果與計算公式結果對比如圖16所示,可以看出,兩者吻合良好,最大誤差為 10.1% ,說明計算公式可以準確預測不同長徑比的柱形裝藥峰值應力值。

3.2變埋深條件下的柱形裝藥峰值應力計算公式

裝藥埋置深度決定了耦合傳入混凝土中的爆炸能量,通常采用TM5-855-1中的峰值應力耦合系數f[49]來量化分析其對柱形裝藥爆炸應力波峰值應力的影響。Gao等[12]和楊耀宗等[14]認為柱形裝藥長徑比對耦合系數 f 影響較小,對實際工程而言可以忽略。此外,Gao等[12]基于爆轟產物與混凝土之間的耦合機制建立了峰值應力耦合系數 f 與比例距離 Z 之間的簡化模型,如圖17(b)所示,具體表達式如下:

圖16數值模擬與式 (4)~(6) 計算結果對比

Fig.16Comparison of stress peak values between numerical simulations and Eqs. (4)-(6)

中: h/l 表示為相對埋深[12],指柱形裝藥底端至靶體表面的距離 h 與裝藥長度1的比值,如圖17(a)所示。結合式 (4)~(8) ,可以得到變埋深條件下,不同長徑比的柱形裝藥峰值應力計算公式:

圖17相對埋深及耦合系數簡化模型示意圖[12]

Fig.17Schematic diagram of relative burial depth and coupling coefficient f[12]

4結論

基于KCC本構模型和多物質ALE算法,采用LS-DYNA軟件開展了柱形裝藥爆炸應力波在混凝土介質中的衰減規律研究。主要對裝藥周圍介質破壞分區進行了劃分,并探討了各個破壞分區上爆炸應力波衰減規律,以及柱形裝藥長徑比對各破壞分區峰值應力衰減規律的影響,并在此基礎上提出了柱形裝藥峰值應力實用化計算公式,主要結論如下。

(1)采用徑向壓應力和環向拉應力為閾值對裝藥周圍介質進行劃分,可以較好地表征爆炸破壞分區的分布;近流體區和壓碎區爆炸應力波為沖擊波,而過渡區和破裂區為塑性波,并且相較于過渡區和破裂區,近流體區和壓碎區爆炸應力波峰值應力衰減更快。這說明峰值應力的衰減規律無法使用單一衰減指數進行描述,需進行分段描述。

(2)隨著柱形裝藥長徑比增加,爆炸應力波峰值應力衰減加快,衰減指數呈線性遞增,近流體區和壓碎區衰減系數呈指數遞減,過渡區和破裂區衰減系數呈線性遞減;此外,柱形裝藥長徑比增加導致法向近流體區、壓碎區、過渡區和破裂區范圍逐漸減小。

(3)基于混凝土介質中爆炸應力波衰減規律的分析,綜合考慮了各破壞分區的差異性、裝藥長徑比以及埋置深度等因素,提出了柱形裝藥爆炸應力波峰值應力實用化計算公式,可以準確快速地計算出柱形裝藥爆炸應力波的峰值應力。

參考文獻:

[1] KRAUTHAMMER T. Modern protective structures [M].Boca Raton: CRC Press,2008. DOI: 10.1201/9781420015423.

[2] PLOOSTER MN. Blast efects from cylindrical explosive charges: experimental measurements [M].Fort Belvoir: Defense Technical Information Center, 1982:11-18.

[3] ISMAIL MM,MURRAY S G.Studyof the blast waves from theexplosion of nonsphericalcharges[J].Propellants, Explosives, Pyrotechnics,1993, 18(3): 132-138. DO1: 10.1002/prep.19930180304.

[4] WU CQ,FATTORIG,WHITAKERA,etal.Investigationofair-blasteffectsfromspherical-andcylindrical-shaped charges[J]. IntermationalJournal ofProtective Structures,2010,1(3): 345-362.DOI: 10.1260/2041-4196.1.3.345.

[5] SHI YC,WANGN,CUIJ,etal.Experimentalandnumericalinvestigationofchargeshapeefectonblastloadinduced by near-field explosions[J].Process SafetyandEnvironmental Protection2022,165:266-277.DOI:10.1016/jpsep.2022. 07.018.

[6] 黃家蓉,劉光昆,吳飚,等.爆炸沖擊作用下混凝土中動態應力波測試與模擬[J].防護工程,2020,42(4):23-28.DOI: 10.3969/j.issn.1674-1854.2020.04.003. HUAIG J Λ,LIU UA, WUD,et a1. Iesung anu sHuiauon OI uyuav suess wave I conuete uuer cxpIosi anu impact[J]. Protective Enginering,2020,42(4): 23-28.DOI: 10.3969/j.issn.1674-1854.2020.04.003.

[7] GEBBEKEN N,GREULICH S,PIETZSCH A.Hugoniot propertiesfor concrete determined by fullscale detonation experiments andflyer-plate-impact tests[J].Intermational Journal of Impact Enginering,2006,32(12): 2017-2031.DOI: 10.1016/j.ijimpeng.2005.08.003.

[8] SHERKAR P,SHINJ, WHITTAKER A, et al. Influence of charge shape and point of detonation on blast-resistant design[J]. Journal of Structural Engineering, 2016,142(2): 04015109. DOI: 10.1061/(ASCE)ST.1943-541X.0001371.

[9]XIAO W F, ANDRAE M, GEBBEKEN N.Efect of charge shape and initiation configuration of explosivecylinders detonating infreeaironblast-resistantdesign[J].JoualofStructuralEnginering2020,146(8):0402146.DO:101061/ (ASCE)ST.1943-541X.0002694.

[10]GAOC,KONGX Z,FANGQ,etal.Numercal investigationonfreeairblastloads generated fromcenter-initiatedcyldrical charges with variedaspectratio inarbitraryorientation[J]. Defence Technology2022,18(9):1662-1678.DOI:10.1016/dt. 2021.07.013.

[11]王明濤,程月華,吳昊.柱形裝藥空中爆炸沖擊波荷載研究[J].爆炸與沖擊,2024,44(4):043201.DOI:10.11883/bzycj2023-0197. WANG M T, CHENG Y H, WU H.Study on blast loadings of cylindrical charges air explosion[J]. Explosion and Shock Waves, 2024, 44(4): 043201. DOI: 10.11883/bzycj-2023-0197.

[12]GAO C, KONG X Z,FANG Q.Experimentaland numerical investigationonthe attenuation ofblast waves inconcrete induced by cylindricalcharge explosion[J].Interational Joual of ImpactEnginering,2023,174:104491.DO:10.1016/ j.ijimpeng.2023.104491.

[13]高矗,孔祥振,方秦,等.混凝土中爆炸應力波衰減規律的數值模擬研究[J].爆炸與沖擊,2022,42(12):123202.DOI: 10.11883/bzycj-2022-0041. GAO C,KONG X Z,FANG Q, et al. Numerical study on attenuation of stress wave in concrete subjected to explosion[J]. Explosion and Shock Waves,2022, 42(12): 123202. DOI: 10.11883/bzycj-2022-0041.

[14]楊耀宗,孔祥振,方秦,等.混凝土中帶殼柱形裝藥爆炸應力波衰減規律的數值模擬[J].爆炸與沖擊,2024,44(11): 112202. DOI: 10.11883/bzycj-2023-0342. YANG Y Z,KONG X Z,FANG Q,et al. Numerical ivestigationon attenuation of stresswaves inconcrete induced by cylindricalcased charge explosion[J].ExplosioandShock Waves,2024,44:12202.DOI:10.11883/bzycj-20-0342.

[15]吳祥云,曲建波,張光明,等.巖石中不同埋深爆炸自由場直接地沖擊參數的預計方法[C]/崔京浩.第 20屆全國結構工 程學術會議論文集(第I冊).《工程力學》雜志社,2011:262-267. WU X Y,QUJB, ZHANGG M,etal. Prediction methodof the direct ground shock parameters of explosion atdierent buried depths in fre fieldofrock[C/CUIJH.Procedings ofthe Twentieth National Conferenceon Structural Enginering (No.I). Engineering Mechanics Magazine,2011: 262-267.

[16]LIU ZY, ZHAIJZ,SU S. Numerical simulation on conical shaped charge with copper liner in several typical shapes[J]. Materials Research Proceedings,2019,13(3): 7-12. DOI: 10.21741/9781644900338-2.

[17]ABIR M,ARUMUGAMD,DHANA B,etalNumerical simulationofblast wave propagatio inlayered soil featuring soilstructure interaction [Cl//COMPDYN. Proceedings of the 6th Intermational Conference on Computational Methods in Structural Dynamics and Earthquake Enginering Methods in Structural Dynamics and Earthquake Engineering.Rhodes Island,2017: 4752-4765. DO1: 10.7712/120117.5759.16936.

[18]KULAK R F,BOJANOWSKI C.Modeling ofcone penetration test using SPH and MM-ALE aproaches [C//Ansys Company. Proceedings of the 8th European LS-DYNA? Users Conference. Strasbourg, 2011: 1-10.

[19]VAN DORSSELAER N,LAPOUJADE V. A contribution to new ALE 2D method validation [C]// Ansys Company. Proceedings of the 11th International LS-DYNA?Users Conference.Dearborn, 2010: 39-50.

[20]MALVAR L J, CRAWFORD JE,WESEVICHJ W,et al. A plasticity concrete material model for DYNA3D [J]. International Journal of Impact Engineering,1997,19(9/10): 847-873. DOI: 10.1016/S0734-743X(97)00023-7.

[21]TUZG,LUY.Evaluationof typicalconcrete materialmodelsusedihydrocodes forhighdynamicresponsesimulations[J]. International Journal of Impact Engineering, 2009,36(1): 132-146. DOI: 10.1016/j.ijimpeng.2007.12.010.

[22]匡志平,陳少群.混凝土Kamp;C模型材料參數分析與模擬[J].力學季刊,2015,36(3):517-526.DOI:10.15959/j.cnki.0254- v0JS.z01J.05.01y. KUANG Z P, CHEN SQ.Analysis and simulation for the material parameters of Kamp;Cconcrete model[J]. Chinese Quarterly ofMechanics,2015,36(3): 517-526. DO1: 10.15959/j.cnki.0254-0053.2015.03.019.

[23]SU Q, WU H,FANG Q.Calibrationof KCC model for UHPC under impact and blast loadings[J]. Cement and Concrete Composites, 2022, 127: 104401. DOI: 10.1016/j.cemconcomp.2021.104401.

[24]KONG X Z,FANG Q,LI QM,et al.Modified Kamp;Cmodel forcratering andscabbing ofconcrete slabs under projectile impact[J]. International Jourmal ofImpact Engineering, 2017,108: 217-228.DOI: 10.1016/.ijimpeng.2017.02.016.

[25]XIAO WF,ANDRAE M, GEBBEKEN N. Air blast TNT equivalence factors of high explosive material PETN for bare charges[J]. Journal of HazardousMaterials,2019,377: 152-162.DOI: 10.1016/j.jhazmat.2019.05.078.

[26]甘露,陳力,宗周紅,等.近距離爆炸比例爆距的界定標準及荷載模型[J].爆炸與沖擊,2021,41(6):064902.DOI: 10.11883/bzycj-2020-0194. GANL,CHENL, ZONG ZH,etal. Definition ofscaled distanceof close-in explosion and blast load caleulation model[J]. Explosion and Shock Waves,2021, 41(6): 064902. DOI: 10.11883/bzycj-2020-0194.

[27]HOPKINSONB.Britishordnance board minutes[J].Joualof the SocietyforArmyHistorical Research,1915,230(57): 88-107.

[28]TUH,FUNGTC,TANK H,etal.Ananalytical model to predict thecompressive damageofconcreteplates undercontact detonation[J]. Interational Journal ofImpactEngineering,2019,134:103344.DOI: 10.1016/j.ijimpeng.2019.103344.

[29]劉琦,翟超辰,張躍飛,等.地面和埋置爆炸土中地沖擊作用分區數值模擬及試驗研究[J].爆炸與沖擊,2022,42(8): 082201. DOI: 10.11883/bzycj-2021-0326. LIU Q,ZHAI C C, ZHANG Y F,et al. Numerical simulation and test study on ground shock subzones in soil produced by ground and buried explosion[J]. Explosion and Shock Waves,2022,42(8): 082201.DOI1: 10.11883/bzycj-2021-0326.

[30]FORBES J W. Shock wave compression of condensed matter: a primer [M]. Berlin: Springer, 2012.

[31]DOBRATZ B M.LLNL explosives handbook: properties of chemical explosives and explosives and explosive simulants: UCRL-52997 [R]. Lawrence: Livermore National Laboratory,1981.

[32]鄭哲敏,解伯民,劉育魁,等.地下核爆炸流體彈塑性計算方案和若干結果[M]//鄭哲敏.鄭哲敏文集.北京:科學出版社, 2004.

[33]鄭哲敏.爆炸成形模型律[M].北京:科學出版社,2004. ZHENG ZM. Explosion forming model law [M]. Beijing: Science Press,2004.

[34]李守巨,何慶志,費鴻祿.巖石爆破破壞分區的研究[J].爆破,1991(1):16-19. LI S J,HE Q Z,FEI HL. Research on the division of rock blasting damage zones [J]. Blasting,1991(1):16-19.

[35]錢七虎,王明洋.巖土中的沖擊爆炸效應[M].北京:國防工業出版社,2010. QIANQH, WANGMY.Impactand explosion effcts inrock andsoil[M].Beijing:NationalDefense IndustryPres,2010.

[36]王明洋,鄧宏見,錢七虎.巖石中侵徹與爆炸作用的近區問題研究[J].巖石力學與工程學報,2005,24(16):2859-2863. DOI: 10.3321/j.issn:1000-6915.2005.16.008. WANG M Y, DENG H J, QIAN Q H. Study on problems of near cavityof penetration and explosion in rock[J]. Chinese Journal ofRockMechanicsandEnginering,2005,24(16): 2859-2863.DOI:10.3321/j.issn:1000-6915.2005.16.008.

[37]張志呈.定向斷裂控制爆破機理綜述[J].礦業研究與開發,200,20(5):40-42.DOI:10.3969/j.ssn.1005-2763.2000. 05.015. ZHANG Z C.Summaryof the mechanism of directional fracture controlled blasting[J]. Mining Research and Development, 2000,20(5): 40-42. DOI: 10.3969/j.issn.1005-2763.2000.05.015.

[38]冷振東.巖石爆破中爆炸能量的釋放與傳輸機制[D].武漢:武漢大學,2017. LENG ZD. Explosion energy release and transmision mechanism in rock blasting[D]. Wuhan: Wuhan University2017.

[39]MANDAL J, GOEL MD, AGARWAL A K. Surface and buried explosions: an explorative review with recent advances [J]. Archives ofComputational Methods in Engineering,2021,28(7): 4815-4835. DOI: 10.1007/s11831-021-09553-2.

[40]AMELSFORTR, WEERHEIJMJT.The failure modeof concrete slabs due to contact charges[R].John Wiley amp; Sons,1994.

[41]SALAMMR. Analytical expresions foruniaxial tensile strengthofconrete interms ofuniaxialcompresivestrength[J]. Transportation Research Record, 1992(1335): 52-54.

[42]宋守志.條形藥包爆炸時的高速沖擊效應[C]/第四屆全國巖石破碎學術討論會論文集.成都:中國巖石力學與工程學 會,中國金屬學會采礦學會,中國土木工程學會隧道及地下工程學會,1989:4. SONG S Z.High-speed impact effcts oflinear charge explosion[Cl//Proceedings of the 4th National Symposium on Rock Fragmentation. Chengdu: Chinese Society forRock Mechanicsand Enginering,Chinese Society of Metals Mining Society, Chinese Society of Civil Engineering Tunnel and Underground Engineering Society,1989: 4.

[43]王明洋,李杰,鄧國強.超高速動能武器鉆地毀傷效應與工程防護[M].北京:科學出版社,2021. WANGMY,LIJ,DENGG Q.Penetrationanddestruction effectsof hypervelocity kineticenergy weaponsand engieering protection [M].Beijing:Science Press,2021.

[44]李重情,穆朝民,石必明.變埋深條件下混凝土中爆炸應力傳播規律的研究[J].振動與沖擊,2017,36(6):140-145.DOI: 10.13465/j.cnki.jvs.2017.07.021. LI Z Q,MUCM,SHIBM. Investigate onshock stress propagation inconcrete atdiferent depths under blasting[J]. Joual of Vibration and Shock, 2017, 36(6): 140-145. DO1: 10.13465/j.cnki.jvs.2017.07.021.

[45]MU CM,ZHOU H,MA HF.Prediction method for ground shock parameters of explosion in concrete [J]. Constructionand BuildingMaterials,2021,291:123372.DOI:10.1016/j.conbuildmat.2021.123372.

[46]LEONG EC,ANAND S,CHEONG HK,etal.Re-examinationof peak stressand scaled distance due to groundshock[J]. International Journal of Impact Engineering, 2007,34(9): 1487-1499. DOI: 10.1016/j.jimpeng.2006.10.009.

[47]YANKELEVSKY D Z, KARINSKI Y S,FELDGUN V R. Re-examination ofthe shock wave's peak pressure attenuation in soils[J]. IntermationalJournalofImpactEngineering,2011,38(11):864881.DOI:10.1016/j.ijimpeng.2011.05.011.

[48]FAN Y,CHENL,LI Z,etal.Modeling theblast load induced byaclose-inexplosionconsideringcylindricalcharge parameters [J].Defence Technology,2023,24: 83-108.DOI: 10.1016/j.dt.2022.02.005.

[49]US Army Enginer Waterways Experiment Station.Fundamentals of protective design for conventional weapons [M]. Washington:USDepartment of theArmy,1986.

(責任編輯 曾月蓉)

主站蜘蛛池模板: 综合五月天网| 亚洲精品大秀视频| 九九热精品免费视频| 精品无码国产自产野外拍在线| 福利在线不卡| 高清色本在线www| 高清久久精品亚洲日韩Av| 在线观看亚洲国产| 欧美午夜精品| 国产成人乱码一区二区三区在线| 久久久久国色AV免费观看性色| 在线精品视频成人网| 六月婷婷综合| 久久久久九九精品影院| 国产日韩丝袜一二三区| 农村乱人伦一区二区| 欧美中文字幕在线视频| 无码中文字幕乱码免费2| 国产精品区视频中文字幕| 国产精品无码作爱| 亚洲高清中文字幕在线看不卡| 中文字幕亚洲另类天堂| 狠狠五月天中文字幕| 99偷拍视频精品一区二区| 婷婷综合色| 亚洲欧美综合另类图片小说区| 中文字幕在线免费看| 欧美一级在线| 中国国产A一级毛片| 欧美在线一二区| 亚洲成人www| 国产欧美日韩另类精彩视频| 亚洲天堂自拍| 久草中文网| 精品人妻无码中字系列| 91 九色视频丝袜| 97狠狠操| 日韩视频福利| 第一页亚洲| 国产美女一级毛片| 亚洲天堂久久久| 无码内射在线| 欧美色综合网站| 国产久操视频| 99久久99视频| 欧美一区国产| 午夜日本永久乱码免费播放片| 97视频免费看| 亚洲高清资源| 国产丝袜无码一区二区视频| 找国产毛片看| 国产网友愉拍精品视频| 国产网站免费看| 蝌蚪国产精品视频第一页| 亚洲午夜片| 精品少妇三级亚洲| 日韩中文欧美| 亚洲精品国产精品乱码不卞| 国产一线在线| 99er这里只有精品| 午夜一级做a爰片久久毛片| 国产精品制服| 思思99热精品在线| 亚洲天堂视频网| 在线人成精品免费视频| 区国产精品搜索视频| 亚洲欧美日韩成人在线| 亚洲色婷婷一区二区| 欧美综合中文字幕久久| 尤物亚洲最大AV无码网站| 国产经典免费播放视频| 热久久这里是精品6免费观看| 伊人婷婷色香五月综合缴缴情| 日韩av高清无码一区二区三区| 大陆国产精品视频| 人妻少妇久久久久久97人妻| 暴力调教一区二区三区| 亚洲一区二区约美女探花| 99re视频在线| 日日噜噜夜夜狠狠视频| 国产成人在线无码免费视频| 免费啪啪网址|