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

基于有限元仿真的不同胸型胸圍截面壓力位移變化研究

2023-12-15 11:05:32劉宿慧顧明月潘怡婷鄒奉元
絲綢 2023年12期
關(guān)鍵詞:有限元模型

劉宿慧, 顧明月, 潘怡婷, 鄒奉元,b,c

(浙江理工大學(xué) a.服裝學(xué)院; b.絲綢文化傳承與產(chǎn)品設(shè)計(jì)數(shù)字化技術(shù)文化和旅游部重點(diǎn)實(shí)驗(yàn)室; c.浙江省服裝工程技術(shù)研究中心,杭州 310018)

人體在受到外界壓力后會(huì)發(fā)生形變,形變程度是反映壓力舒適性的重要指標(biāo)[1]。因此,探究壓力與體表形變的規(guī)律對(duì)壓力舒適性影響具有重要意義。近年來(lái),許多學(xué)者利用有限元數(shù)值模擬對(duì)服裝壓力與相應(yīng)位移之間的關(guān)系進(jìn)行研究。Sun等[2-3]利用逆問(wèn)題的迭代求解方法預(yù)測(cè)了乳房超彈材料系數(shù),基于該系數(shù)建立了人體生物力學(xué)有限元(FE)模型,分析文胸對(duì)乳房皮膚的壓力分布。Wang等[4]以跑步緊身褲為研究對(duì)象,分析了腿部骨骼、肌肉、皮膚、服裝面料與體表位移之間的關(guān)系,預(yù)測(cè)了穿著過(guò)程中緊身褲對(duì)小腿肌肉位移的影響。王燕珍等[5]通過(guò)對(duì)人體皮膚的動(dòng)態(tài)測(cè)量,研究穿著緊身褲時(shí)服裝織物彈性與皮膚變形之間的相關(guān)性。Liu等[6]建立標(biāo)準(zhǔn)女體有限元模型,研究了服裝壓力對(duì)胸部的影響及服裝胸圍應(yīng)變與織物楊氏模量間的關(guān)系。以上研究主要針對(duì)單一體型或歸一化人體進(jìn)行研究,沒(méi)有考慮體型差異造成的壓力與位移變化。因此,本文結(jié)合胸型形態(tài)差異,利用有限元方法對(duì)不同胸型的壓力位移變化規(guī)律進(jìn)行研究。

由于胸型差異的影響,女體胸圍截面各部受壓狀態(tài)及形態(tài)變化有所不同,不同體型的人體受到外力時(shí)會(huì)呈現(xiàn)不同變化規(guī)律[7]。因此,用有限元仿真得出體表形態(tài)變化規(guī)律時(shí),需將胸型差異考慮在內(nèi)。常麗霞等[8]基于乳深、乳平圍兩個(gè)特征變量對(duì)女性乳房的立體程度進(jìn)行表征,將乳房形態(tài)細(xì)分為9類(lèi)。Liu等[9]通過(guò)對(duì)人體108個(gè)測(cè)量項(xiàng)目進(jìn)行主成分分析,篩選出乳房深和胸寬厚比兩個(gè)指標(biāo),對(duì)乳房形態(tài)進(jìn)行分類(lèi)。Shi等[10]分析年齡和BMI對(duì)胸罩尺碼系統(tǒng)的影響,提出將年齡和BMI作為胸型分類(lèi)的輔助標(biāo)準(zhǔn)。探究不同胸型下女體受法向壓力與體表內(nèi)向位移之間的關(guān)系,對(duì)于壓力舒適性研究具有實(shí)際和應(yīng)用價(jià)值。

綜上,本文為獲取不同胸型的人體數(shù)據(jù),采用三維掃描儀測(cè)量女體軀干數(shù)據(jù),提取重建坐標(biāo)系下的乳房參數(shù),根據(jù)乳房深與乳房寬比值對(duì)胸型進(jìn)行分類(lèi),構(gòu)建不同胸型的有限元接觸模型,模擬不同胸型在不同受壓狀態(tài)下的形態(tài)變化,探究不同胸型的壓力—位移變化規(guī)律,建立胸圍截面處的壓力、位移和角度間函數(shù)關(guān)系,以期為緊身服裝壓力舒適性研究提供參考。

1 胸圍截面幾何模型構(gòu)建

本文通過(guò)三維掃描儀獲取女體表面的空間信息,利用逆向工程軟件對(duì)三維數(shù)據(jù)進(jìn)行預(yù)處理并與內(nèi)層結(jié)構(gòu)數(shù)據(jù)裝配,分胸型生成人體胸部的幾何模型,為后續(xù)分類(lèi)建立有限元模型提供數(shù)據(jù)支撐。

1.1 數(shù)據(jù)獲取

為獲取各種胸型的體表形態(tài)數(shù)據(jù),本文利用美國(guó)[TC]2NX-16三維掃描儀獲取了136名年齡在18~25歲女性受試者的點(diǎn)云數(shù)據(jù)。實(shí)驗(yàn)環(huán)境溫度控制在20 ℃左右,相對(duì)濕度控制在60%左右,實(shí)驗(yàn)對(duì)象保持上半身全裸,下半身穿著貼體淺色內(nèi)褲的狀態(tài)[11]。

1.2 胸型分類(lèi)

由于采集到的人體數(shù)據(jù)差異性大且存在誤差數(shù)據(jù),需要對(duì)數(shù)據(jù)進(jìn)行篩選、精簡(jiǎn)、降噪、對(duì)稱(chēng)等前處理操作來(lái)確保數(shù)據(jù)的有效性。在保證人體特征點(diǎn)清晰的前提下,消除噪點(diǎn)影響,從而獲得較為光滑的人體數(shù)據(jù)。剔除誤差數(shù)據(jù)后,保留114名人體數(shù)據(jù)。

通過(guò)提取乳房輪廓特征點(diǎn)來(lái)構(gòu)建胸部的局部坐標(biāo)系。本文提取乳底點(diǎn)(BBP)、乳房外緣點(diǎn)(OBP)、乳房?jī)?nèi)緣點(diǎn)(IBP)和乳點(diǎn)(BP)作為乳房特征點(diǎn)[12],將BBP點(diǎn)、OBP點(diǎn)、IBP點(diǎn)連接,建立YZ坐標(biāo)平面,經(jīng)過(guò)BP點(diǎn)且法線為OBP點(diǎn)和IBP點(diǎn)連線的面作為XZ坐標(biāo)平面,經(jīng)過(guò)BP點(diǎn)且與YZ平面和XZ平面垂直的面作為XY坐標(biāo)平面,三個(gè)坐標(biāo)平面的交點(diǎn)為原點(diǎn)O,使用對(duì)齊坐標(biāo)功能構(gòu)建局部坐標(biāo),如圖1(a)所示。將BP點(diǎn)到O點(diǎn)的距離作為乳房深a,將OBP點(diǎn)和IBP點(diǎn)的距離作為胸寬b,如圖1(b)所示。

圖1 局部坐標(biāo)系

為建立有限元模型,探究不同胸型下的受壓情況,需對(duì)受試者進(jìn)行胸型分類(lèi)。本文選取能夠較好體現(xiàn)胸型形狀差異的胸深寬比(Depth Width Ratio,DWR)來(lái)對(duì)胸型進(jìn)行分類(lèi)[9]。此方法將乳房深a與乳房寬b的比值作為分類(lèi)指標(biāo),分類(lèi)時(shí)胸型覆蓋率較高[13]。對(duì)樣本進(jìn)行分類(lèi),能夠體現(xiàn)胸部形態(tài)差異及胸型特點(diǎn)。根據(jù)DWR值的大小范圍將114名人體數(shù)據(jù)分為6類(lèi)胸型,DWR值越大,乳房形狀越突出且起伏越大。胸型及對(duì)應(yīng)的樣本數(shù)如表1所示。

表1 胸型分類(lèi)及樣本數(shù)

1.3 幾何曲線裝配

截取114個(gè)人體模型過(guò)BP點(diǎn)的水平截面作為胸圍截面,在截面曲線上以胸寬為橫軸,胸厚為縱軸構(gòu)建坐標(biāo)系,利用Matlab編程對(duì)截面曲線進(jìn)行平滑和提取數(shù)據(jù)坐標(biāo)點(diǎn)處理,每5°定一個(gè)坐標(biāo)點(diǎn),共劃分為72個(gè)點(diǎn),將每類(lèi)胸型下的坐標(biāo)數(shù)據(jù)進(jìn)行平均化處理[14]。本文將女體胸圍截面曲線作為皮膚層,得到6種胸型的皮膚層曲線。

軟組織層包括輸乳管竇、乳腺、肌肉等內(nèi)部結(jié)構(gòu)[6],位于皮膚層內(nèi)部,其坐標(biāo)數(shù)據(jù)根據(jù)皮膚層和皮膚層厚度計(jì)算可得。女體胸部皮膚層厚度為1.5 mm[14],利用Matlab編程,分別將6類(lèi)胸型皮膚層坐標(biāo)值導(dǎo)入,分類(lèi)獲取胸圍截面處軟組織層坐標(biāo)值,計(jì)算如下式所示。

(1)

式中:x0代表幾何模型皮膚層曲線上任意一個(gè)點(diǎn)的x坐標(biāo),y0代表幾何模型皮膚層曲線上任意一個(gè)點(diǎn)的y坐標(biāo),x代表軟組織層曲線上與皮膚層相對(duì)應(yīng)的點(diǎn)的x坐標(biāo),y代表軟組織層曲線上與皮膚層相對(duì)應(yīng)的點(diǎn)的y坐標(biāo);d代表皮膚層厚度,θ代表相應(yīng)點(diǎn)的角度。

骨骼層曲線參照中國(guó)可視化女體斷層圖片[15],重建幾何模型中的骨骼層相對(duì)位置,進(jìn)行順滑、對(duì)稱(chēng)處理,獲得骨骼層曲線坐標(biāo),如圖2所示。

圖2 女體胸圍截面橫向剖面

根據(jù)女體真實(shí)結(jié)構(gòu),重建皮膚層、軟組織層和骨骼層曲線坐標(biāo)系,將每層曲線的原點(diǎn)坐標(biāo)對(duì)齊,對(duì)截面模型進(jìn)行幾何裝配,構(gòu)建的幾何模型作為有限元建模的基礎(chǔ)。本文以B胸型為例,每層曲線及幾何裝配模型如圖3所示。

圖3 幾何模型裝配曲線

2 胸圍截面有限元建模與仿真

本文基于6種胸型的幾何模型,建立有限元模型,施加5種不同大小的壓力載荷,模擬不同受壓狀態(tài)下6類(lèi)胸型所受壓力與體表位移的變化規(guī)律。在確保模型精度的情況下簡(jiǎn)化模型,對(duì)模型進(jìn)行假設(shè)并提高運(yùn)算精度確保模型有效性[16]。表2為本文對(duì)模型做出的假設(shè)。

表2 模型假設(shè)條件

把幾何模型導(dǎo)入Abaqus中,創(chuàng)建胸圍截面各層部件,將胸圍截面看作各向同性的彈性體。皮膚層、軟組織層、骨骼層材料參數(shù)如表3所示[17]。

表3 材料參數(shù)

網(wǎng)格類(lèi)型為CPE3,種子大小為2。設(shè)置部件相互作用,使皮膚層、軟組織層及骨骼層在邊界上連接,但卻各自獨(dú)立,皮膚、軟組織、內(nèi)部骨骼間互相黏接,當(dāng)受到外部應(yīng)力時(shí)無(wú)相對(duì)滑移,此時(shí)將皮膚層和軟組織層接觸定義為變形體與變形體的黏結(jié),將軟組織層和骨骼層接觸定義為變形體與剛體的黏結(jié)[18]。由于骨骼在受壓后不發(fā)生位移,將邊界條件設(shè)定為骨骼固定,即x方向和y方向位移均為0。

本文基于6類(lèi)胸型的不同形狀特點(diǎn),研究胸圍截面上受體表特征差異影響的不同點(diǎn)處受壓位移變化特征。為得出不同胸型受壓力影響易產(chǎn)生皮膚形變的位置,對(duì)胸圍截面施加不同大小的均勻壓力。因此,在緊身胸衣舒適范圍[19]內(nèi)施加2、2.5、3、3.5、4.0 kPa的5種法向壓力載荷,作用于皮膚表面。探究當(dāng)受壓發(fā)生變化時(shí),各胸型受壓位移變化規(guī)律。邊界條件及壓力施加情況如圖4所示。本文以B類(lèi)胸型為例,有限元模型如圖5所示。

圖4 壓力施加及邊界條件

圖5 B類(lèi)胸型的有限元模型

3 結(jié)果與分析

按照DWR值大小得到了6類(lèi)胸型,不同胸型對(duì)應(yīng)的胸圍截面形狀、軟組織層厚度等差異導(dǎo)致模擬結(jié)果和胸圍截面壓力與位移函數(shù)關(guān)系有所不同。因此,基于胸型差異本文分別對(duì)6類(lèi)胸型仿真的結(jié)果進(jìn)行分析。

3.1 仿真結(jié)果

通過(guò)Abaqus有限元分析軟件計(jì)算,得到在2~4 kPa舒適壓力內(nèi)A、B、C、D、E、F這6類(lèi)胸型受5種不同大小壓力載荷作用下的應(yīng)力及位移變化情況。圖6為受壓3 kPa時(shí),6種胸型所受應(yīng)力情況。由圖6可知,乳房位置應(yīng)力最大,約為2 kPa;其次是后側(cè)位置,約為1.7 kPa。故可知,在外力作用下,曲率較大且位置突出處應(yīng)力較大。

圖6 3 kPa壓力下不同胸型應(yīng)力云圖

圖7為受壓3 kPa時(shí),6類(lèi)胸型的位移云圖。由圖7可知,A胸型的形變較大處集中在靠近胸高點(diǎn)并偏向前中的位置,形變較為集中。從A胸型向F胸型過(guò)渡的過(guò)程中,乳房表面形態(tài)起伏逐漸變大,胸部軟組織層變厚,形變逐漸從胸高點(diǎn)向側(cè)面位置分散增加,但形變最大位置始終在BP點(diǎn)附近,隨胸型變化增大呈現(xiàn)外擴(kuò)趨勢(shì)。

圖7 3 kPa壓力下不同胸型位移云圖

3.2 胸圍截面壓力、位移與角度函數(shù)關(guān)系

體表的受壓變形量是表征人體壓力舒適性及緊身衣設(shè)計(jì)的理論依據(jù)。本文在前期數(shù)據(jù)預(yù)處理階段,通過(guò)相關(guān)性分析得出“壓力/位移”作為整體時(shí)與角度間的相關(guān)性較高,6類(lèi)胸型的擬合度較高,R2均大于0.95。因此,以角度為自變量,“壓力/位移”為因變量,探究不同胸型壓力、位移與角度間的函數(shù)關(guān)系。根據(jù)“壓力/位移”值,可通過(guò)易測(cè)得的壓力數(shù)據(jù)計(jì)算難以常規(guī)測(cè)量的體表位移數(shù)據(jù),建立壓力與體表位移間的關(guān)系。對(duì)A~F這6類(lèi)胸型“(壓力/位移)—角度”的函數(shù)關(guān)系曲線分段進(jìn)行擬合,有助于對(duì)以胸圍為依據(jù)的緊身服裝設(shè)計(jì)進(jìn)行優(yōu)化。由于對(duì)胸型曲線進(jìn)行了平均、對(duì)稱(chēng)等標(biāo)準(zhǔn)化處理,在[180°,360°]的數(shù)據(jù)以x=180°為對(duì)稱(chēng)軸與[0°,180°]對(duì)稱(chēng),因此,本文僅對(duì)[0°,180°]“壓力/位移—角度”的函數(shù)關(guān)系進(jìn)行擬合,并得出曲線函數(shù)及其相應(yīng)擬合度,如圖8所示。擬合函數(shù)及擬合度如表4所示,其中,x為角度,y為“壓力/位移”值。

表4 A~F胸型的“(壓力/位移)—角度”函數(shù)關(guān)系

圖8 A~F胸型“(壓力/位移)—角度”函數(shù)關(guān)系擬合曲線

由表4可知,在[0°,180°]內(nèi),A胸型分兩段進(jìn)行擬合,B~F胸型分三段進(jìn)行擬合,其擬合曲線的擬合度(R2)均在0.99左右。根據(jù)擬合函數(shù),也可以得知不同胸型下特定壓力時(shí)的位移大小。因壓力相對(duì)于體表位移更易獲得,根據(jù)壓力計(jì)算的位移值是評(píng)價(jià)壓力舒適性及緊身服裝塑形美觀效果的重要指標(biāo),利用此擬合函數(shù),可以高效、準(zhǔn)確地得到體表位移值這一關(guān)鍵數(shù)據(jù)。

3.3 基于不同胸型的壓力位移變化規(guī)律分析

3.3.1 不同胸型的體表位移變化規(guī)律

圖9為5種壓力作用下人體的位移曲線。由圖9可知,在不同壓力作用下,B、C、D、E、F胸型的曲線整體變化趨勢(shì)基本一致,隨著壓力增大,位移也隨之增大,曲線整體呈現(xiàn)上移趨勢(shì),胸型之間的位移差也是遞增的。在不同壓力作用下,A胸型在0°和180°位置的形變始終大于其他胸型。

B、C、D、E、F胸型形變最大點(diǎn)均在50°和310°左右。而受BP點(diǎn)位置影響,A胸型的形變最大點(diǎn)在35°和325°左右。B、C、D、E、F胸型形狀相對(duì)A胸型較外擴(kuò),每類(lèi)胸型形變最大點(diǎn)和其對(duì)應(yīng)的BP點(diǎn)位置一致,越接近BP點(diǎn),形變就會(huì)越大。

B、C、D、E、F胸型形變最小點(diǎn)在0°左右,即前中位置。A胸型的形變最小點(diǎn)在90°~95°,對(duì)應(yīng)人體側(cè)縫位置。這6類(lèi)胸型在0°附近位移變化按順序遞減,此位置處A胸型位移變化最大,形變值為E胸型的2倍左右。

對(duì)于B、C、D、E、F胸型,僅次于最低點(diǎn)的極小值點(diǎn)出現(xiàn)在160°和200°附近,即180°對(duì)應(yīng)的胸圍截面的后中位置附近。這是由于A胸型女體偏瘦,其90°和270°附近的人體側(cè)縫位置軟組織層厚度較薄且基本小于A胸型胸圍截面的其他位置,導(dǎo)致最小形變量出現(xiàn)在該位置附近;而B(niǎo)~F胸型軟組織層最薄處位于0°附近,同時(shí)后中180°附近位置也相對(duì)較薄的,僅次于0°位置,0°和180°±20°分別與最小位移變化和極小值點(diǎn)的位置相對(duì)應(yīng)。因此,胸型影響下的軟組織層厚度和受壓位移密切相關(guān),受不同胸型形態(tài)起伏及軟組織層厚度規(guī)律影響,相同壓力下,軟組織層較薄的位置,其位移變化較小。

3.3.2 不同壓力下的體表位移變化規(guī)律

圖10為6種胸型在不同壓力下的位移變化規(guī)律。由圖10可知,6種胸型體表的位移變化規(guī)律基本一致,和每類(lèi)胸型的體表形態(tài)特征變化一致。隨著壓力增大,體表各處形變呈遞增趨勢(shì)。從前中到胸高點(diǎn)位置,位移隨體表形態(tài)變化逐漸增大,再過(guò)渡到側(cè)縫位置,隨軟組織厚度減小,形變相應(yīng)減少。從側(cè)縫至后側(cè)120°左右位置,位移值隨軟組織變厚相應(yīng)增加,最后向后中方向總體呈下降趨勢(shì),180°~360°和0°~180°對(duì)稱(chēng)變化。

圖10 A~F胸型位移隨壓力變化曲線

4 結(jié) 論

受不同胸型的截面形態(tài)差異影響,受壓后胸圍截面發(fā)生的形變規(guī)律隨胸型變化呈現(xiàn)不同特點(diǎn)。本文根據(jù)乳房輪廓特征參數(shù)差異將114名人體數(shù)據(jù)細(xì)分為6類(lèi),針對(duì)每類(lèi)胸型特點(diǎn)分別建立有限元數(shù)值模型,對(duì)其施加5種均勻大小的法向壓力,研究胸型差異下胸圍截面壓力與位移變化規(guī)律,并擬合得到了“(壓力/位移)—角度”曲線及其函數(shù)關(guān)系,主要結(jié)論如下:

1) 6類(lèi)胸型受5種不同大小壓力載荷作用下,胸圍截面BP點(diǎn)應(yīng)力最大,其次是后側(cè)位置,曲率較大且位置突出處應(yīng)力較大。

2) 在壓力作用下,A胸型的形變主要集中在靠近胸高點(diǎn)并偏向前中的位置,產(chǎn)生的位移變化比較集中。而從A胸型向F胸型過(guò)渡的過(guò)程中,由于乳房表面起伏形態(tài)逐漸變大,胸部及人體側(cè)面軟組織層變厚,形變隨胸型變化增大且呈現(xiàn)外擴(kuò)趨勢(shì),但形變最大的位置始終在BP點(diǎn)附近。

3) 基于胸型特點(diǎn)構(gòu)建了6類(lèi)胸型壓力、位移和角度間三角函數(shù)關(guān)系,擬合度(R2)均在0.99左右,擬合效果較好。由于實(shí)際皮膚表面形變不易通過(guò)測(cè)量直接獲得,可以在已知壓力數(shù)據(jù)情況下即可求得體表形變數(shù)據(jù),替代了受壓位移繁瑣的測(cè)量,從而簡(jiǎn)化測(cè)量步驟,可為緊身服裝壓力舒適性研究提供參考。

《絲綢》官網(wǎng)下載

中國(guó)知網(wǎng)下載

猜你喜歡
有限元模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
新型有機(jī)玻璃在站臺(tái)門(mén)的應(yīng)用及有限元分析
基于有限元的深孔鏜削仿真及分析
基于有限元模型對(duì)踝模擬扭傷機(jī)制的探討
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
磨削淬硬殘余應(yīng)力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
主站蜘蛛池模板: 无码专区在线观看| 欧美激情综合| 一级全免费视频播放| 日本a∨在线观看| 国产精品嫩草影院av| 国产丰满成熟女性性满足视频| 9999在线视频| 91视频免费观看网站| 国内精自视频品线一二区| 免费视频在线2021入口| 日本黄网在线观看| 亚洲天堂精品视频| 亚洲日本中文字幕乱码中文 | 一本大道香蕉中文日本不卡高清二区| 亚洲女同一区二区| 亚洲永久精品ww47国产| 亚洲欧美日韩综合二区三区| 亚洲欧美日韩动漫| 婷婷成人综合| 熟妇丰满人妻| 99国产精品免费观看视频| 午夜人性色福利无码视频在线观看| 国产另类视频| 国产sm重味一区二区三区| 亚洲毛片一级带毛片基地| 啪啪永久免费av| 久久精品国产精品国产一区| 人人看人人鲁狠狠高清| 久久无码av一区二区三区| 中文字幕在线播放不卡| 欧美综合区自拍亚洲综合绿色| 国产精品久久久久久久伊一| 亚洲日韩在线满18点击进入| 成人免费视频一区二区三区| 91视频免费观看网站| 2022国产无码在线| 国产污视频在线观看| 欧美激情视频一区| 亚洲无码精彩视频在线观看| 成人欧美在线观看| 亚洲最大综合网| 国产精品亚洲五月天高清| 亚洲综合精品第一页| 亚洲中文无码h在线观看 | 国产黑丝一区| 伊人久久青草青青综合| 露脸国产精品自产在线播| 亚洲精品成人福利在线电影| 四虎国产精品永久一区| 青青青国产视频| 免费一级毛片在线播放傲雪网| 国产在线专区| 四虎国产精品永久一区| 精品三级网站| 91色国产在线| 91在线激情在线观看| 中文字幕在线日韩91| 亚洲成人免费在线| 超级碰免费视频91| 精品国产欧美精品v| 55夜色66夜色国产精品视频| 亚洲av无码人妻| 国产精品999在线| 国产在线视频二区| 亚洲日本中文字幕天堂网| 99久久精品无码专区免费| 亚洲成肉网| 在线观看热码亚洲av每日更新| 老司机精品久久| 无码aⅴ精品一区二区三区| 在线观看国产精美视频| 免费在线色| 久久这里只有精品23| 综1合AV在线播放| 热热久久狠狠偷偷色男同| 一级在线毛片| 久久香蕉国产线看观看亚洲片| 午夜国产精品视频| 婷五月综合| 国产人免费人成免费视频| 91精品免费高清在线| 中文字幕免费播放|