亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        電阻接地狀態(tài)下星用電路板深層充電仿真方法

        2017-07-05 15:33:49王建昭陳鴻飛蔡震波
        航天器環(huán)境工程 2017年3期
        關(guān)鍵詞:有限元

        王建昭,陳鴻飛,蔡震波

        (1. 北京空間飛行器總體設(shè)計部,北京 100094;2. 北京大學(xué) 空間物理與應(yīng)用技術(shù)研究所,北京 100871)

        電阻接地狀態(tài)下星用電路板深層充電仿真方法

        王建昭1,陳鴻飛2,蔡震波1

        (1. 北京空間飛行器總體設(shè)計部,北京 100094;2. 北京大學(xué) 空間物理與應(yīng)用技術(shù)研究所,北京 100871)

        文章基于Monte Carlo方法和有限元方法,對雙層和4層電路板覆銅層通過電阻接地時的深層充電進(jìn)行仿真分析,詳細(xì)討論了空間各向同性電子通量模型、電路板背面和中間覆銅層分別通過電阻接地時的計算方法和邊界條件,以及不同接地條件下有限元矩陣方程的建立;最終定量計算了電阻阻值對電路板充電結(jié)果的影響。仿真結(jié)果表明,較之完全接地情況,通過電阻接地會增加充電電場和電勢,最大電勢深度也相應(yīng)變化;電阻接地層電勢和阻值呈線性關(guān)系;當(dāng)接地電阻為109?量級及以下時,其對深層充電的影響可以忽略,驗證了NASA-HDBK-4002A手冊中設(shè)計指南的正確性。

        深層充電;接地電阻;Monte Carlo方法;有限元方法;Geant4

        0 引言

        地球空間輻射環(huán)境對衛(wèi)星在軌安全運(yùn)行有重要影響。其中高能電子(0.1~10 MeV)穿透衛(wèi)星蒙皮和儀器外殼,沉積在介質(zhì)內(nèi)并逐步建立電場的過程,稱為深層充電。一旦電場強(qiáng)度超過介質(zhì)放電閾值,將發(fā)生深層放電,對介質(zhì)和儀器設(shè)備造成破壞。

        國內(nèi)外深層充電研究主要集中于地面試驗[1]、數(shù)值模擬和在軌監(jiān)測[2-3]等方面。其中,數(shù)值模擬具有成本低廉且可以詳細(xì)探究充電過程的優(yōu)勢。對于粒子輸運(yùn)計算,主要有解析方法[4-5]和 Monte Carlo方法[6]:前者方便快捷,后者在模擬粒子數(shù)足夠多時更加精確。對于充電過程計算,主要有時域差分法[7]和有限元法[8]:前者計算方法簡單,后者在邊界條件復(fù)雜時更方便并且可以計算三維情況。

        對于星用電路板,為了防止由于電路板孤立使得其充電電場強(qiáng)度過高,一般采用覆銅層接入衛(wèi)星結(jié)構(gòu)地的方式增加電荷泄放速率,降低充電電場強(qiáng)度。接地條件和入射電子能量的不同對深層充電有很大影響[9],此外雙層接地相對單面接地、多層接地相對雙層接地,充電電場強(qiáng)度進(jìn)一步降低[10-11]。相對電路板中心,電路板邊緣區(qū)域的充電電場更大[8]。

        在深層充放電的地面或星載測量設(shè)備中[12-13],需要測量接地處電流以評估深層放電風(fēng)險。因此電路板不能采用完全接地的方式,而是應(yīng)在覆銅層和地之間接入一個電阻,通過電阻上的壓降來測量覆銅層與結(jié)構(gòu)地之間的放電電流和脈沖,進(jìn)而評估放電風(fēng)險。另外,NASA-HDBK-4002A手冊[14]的設(shè)計指南提出,電路板應(yīng)通過0~10 M?的電阻接地。因此需要定量分析電路板充電特性和接地電阻的關(guān)系,為接地電阻的選擇提供依據(jù)和參考。當(dāng)覆銅層通過電阻接地時,邊界條件更復(fù)雜。不同于此前的時域差分法,本文建立了一種基于Monte Carlo方法和有限元方法的仿真模型,可以計算通過電阻接地時的電路板深層充電電場和電位。

        1 仿真方法

        本文選擇深層充電風(fēng)險較高的典型地球同步軌道(GEO),使用基于Monte Carlo方法的軟件包Geant4編寫的程序,模擬空間高能電子在介質(zhì)中的輸運(yùn),再根據(jù)有限電阻接地的物理模型和邊界條件,用一維有限元方法計算得到電路板不同深度的深層充電情況。

        1.1 電子通量模型

        深層充電多發(fā)生在空間環(huán)境較惡劣時,因此選用ESA的FLUMIC3.0電子通量模型,可以描述內(nèi)外輻射帶電子通量[15]。對于外輻射帶(L>2.5),任意L值處、能量>E的高能電子積分通量為

        其中,F(xiàn)(foy, fsc)是與季節(jié)相位因子foy和太陽活動周相位因子fsc有關(guān)的GEO上能量>2 MeV的電子積分通量。

        對 GEO,由于內(nèi)源場偏離偶極場以及外源場的存在,衛(wèi)星位置處的L值隨時間變化,電子輻射環(huán)境也隨之變化。因此,把軌道分為某一固定時間間隔的點陣,利用ONERA-DESP的IRBEM函數(shù)庫計算這些點的軌道參數(shù),并轉(zhuǎn)化為(B, L)坐標(biāo),再利用FLUMIC模型計算這些點的電子通量能譜,將各點能譜進(jìn)行平均可得到該軌道的平均電子積分通量能譜。

        電路板是放置在機(jī)箱內(nèi)的,因此各向同性的高能電子須入射并穿透屏蔽層后才能沉積到電路板中。設(shè)置衛(wèi)星蒙皮和機(jī)箱外殼的屏蔽厚度均為1 mm等效鋁。計算得到GEO入射電子經(jīng)2 mm鋁屏蔽后的電子微分通量能譜,如圖1所示。

        1.2 電路板模型

        選用文獻(xiàn)[11]中所用的雙層和4層電路板模型,如圖2所示。電路板以FR4(環(huán)氧樹脂玻璃布層壓板)為基材,厚度為1.41 mm,密度為1.78×103kg/m3,其上下表面分別有厚度為 35 μm 的覆銅層;對于4層電路板,距上下表面287 μm處還各有一層17 μm厚的覆銅層。這些覆銅層作完全接地或有限電阻接地處理。

        本文用Geant4程序模擬粒子輸運(yùn),為模擬真實空間的各向同性電子分布,如圖 3所示,選用平面電子發(fā)射源,平面 A為可能的電路板位置,電子通量沿平面源法向x的極角θ呈cos型分布。對于平面源上的一點,沿θ角度發(fā)射的粒子通量為沿x軸通量的cos θ倍,這樣任意放置的平面A所接收的通量不隨θ而變,從而形成空間各向同性粒子分布。

        假設(shè)真實空間電子各向同性電子積分通量為F,m-2·s-1·sr-1;且電路板放置在靠近衛(wèi)星蒙皮一側(cè),可忽略其背面所受輻射,即只有一半電子入射到屏蔽層。對于屏蔽層上一點,θ0是與電路板法線之間的夾角,對半空間積分可得:

        而對于平面源上一點,電子通量呈cos型分布,對半空間積分可得:

        因此,cos型分布導(dǎo)致的各向同性分布通量為真實空間的1/2,即歸一化因子fac=1/2。

        假設(shè)Geant4程序模擬電子數(shù)為N,則相應(yīng)的仿真時間(真實空間輻射同樣個數(shù)電子所需時間)為

        其中,S為屏蔽層的面積。

        1.3 物理模型

        本文計算深層充電所用的物理模型為

        式中:ε為介電常數(shù);φ為電勢;ρ為電荷密度;j為電子束電流,由Geant4程序計算得到;σ為總電導(dǎo)率;E為電場。方程分別為Possion方程、電流連續(xù)性方程和電勢定義方程。

        高分子聚合物在輻射下電導(dǎo)率將明顯增加,

        式中:σd為材料未受輻射時的暗電導(dǎo)率;σr為輻射誘導(dǎo)電導(dǎo)率;D˙為輻射劑量率,rad·s-1,由Geant4程序計算得到;k和?為輻射誘導(dǎo)電導(dǎo)率的系數(shù)和指數(shù),通過實驗測得。另外,充電電場也會導(dǎo)致電導(dǎo)率的增加[16],

        式中:σ0為電場為 0 時的電導(dǎo)率;βF=(e3/πε)1/2,e為電子電量;kB為Boltzmann常數(shù);T為溫度。環(huán)氧樹脂的物理特性見表1。

        表1 環(huán)氧樹脂物理特性Table 1 Physical properties of FR4

        1.4 計算方法和邊界條件

        仿真計算中,用有限元方法求解Possion方程,用時域差分法求解電流連續(xù)性方程,流程如圖4所示。先求出第一時刻 t1=Δt的電荷密度分布 ρ;再根據(jù)Possion方程得到電勢分布φ;轉(zhuǎn)化為電場分布 E后,再求解電流連續(xù)性方程,得到下一時刻t2=t1+Δt的電荷密度分布 ρ。如此反復(fù),直到求解的電場和電勢達(dá)到飽和值。

        1)背面覆銅層有限電阻接地

        考慮電路板覆銅層通過有限電阻接地情況的邊界條件,以背面接電阻為例,如圖5所示。下標(biāo)b代表電路板背面。

        式中:Jb+和Jb-分別為覆銅層上表面兩側(cè)的電子束電流,Jb+= Jb-;σb和Eb分別為覆銅層上表面的總電導(dǎo)率和電場;jb為流經(jīng)電阻的電流。則邊界條件為

        式中:φb為覆銅層上表面的電勢;Rb為覆銅層連接的接地電阻。此邊界條件為第三類邊界條件。

        2)中間覆銅層有限電阻接地

        圖6為中間覆銅層有限電阻接地情況。

        根據(jù)電流守恒,

        式中:Jb1和 Jb2、σb1和 σb2、Eb1和 Eb2分別為覆銅層上下表面的電子束電流、總電導(dǎo)率和電場;jb為流經(jīng)電阻的電流。則邊界條件為

        式中,φb1和φb2分別為覆銅層上下表面的電勢。

        1.5 有限元矩陣方程的建立和邊界條件的強(qiáng)加

        1)有限元矩陣方程的建立

        考察深層充電情況隨電路板深度的變化,只需考慮一維有限元情況。將電路板分為M層,共N個(N=M+1)節(jié)點,如圖7所示。用上標(biāo)表示單元編碼,下標(biāo)表示一個單元內(nèi)的結(jié)點編碼,即對于第i個單元,有:x2(i)=x1(i+1)。

        選擇線性插值函數(shù),電勢φ(i)(x)=a(i)+b(i)x,單元長度為l(i)。用Galerkin方法構(gòu)建有限元矩陣方程(詳細(xì)推導(dǎo)過程參見文獻(xiàn)[17])

        式中:K為 N×N矩陣;φ、b、g為 N×1矩陣;φn表示第 n個節(jié)點的電勢值。對于一維 Possion方程,K為

        其余元素為0;b為

        若電勢連續(xù),則gi=0。

        2)邊界條件

        分別考慮不同接地方式下的邊界條件,具體計算時,不同覆銅層采取不同的邊界條件。

        a)完全接地或覆銅層加偏置電壓

        以正面完全接地或加偏壓 U為例,此時為第一類邊界條件φ|x=0=0或φ|x=0=U,即φ1=b1+g1=0或φ1=0,對應(yīng)的矩陣元素變化為:

        b)背面覆銅層有限電阻接地

        根據(jù) 1.4節(jié)的推導(dǎo),此時為第三類邊界條件

        即 bN+gN=bN–φNε/(σbRb),對應(yīng)的矩陣元素變化為

        c)中間覆銅層有限電阻接地

        假設(shè)第n層為覆銅層,通過有限電阻接地,根據(jù)1.4節(jié)的推導(dǎo),邊界條件為

        假設(shè) σb1=σb2=σb,即 bn+gn=bn+ε(Jn–Jn+1)/σb–φnε/(σbRb),bn+1+gn+1=bn+1+ε(Jn–Jn+1)/σb–φn+1ε/(σbRb),對應(yīng)于矩陣變化,為

        2 結(jié)果分析

        2.1 背面覆銅層有限電阻接地

        考慮雙層電路板正面完全接地處理,背面分為完全接地和有限電阻(1012?)接地2種情況,模擬得到經(jīng)過100 h輻射后的深層充電情況(如圖8所示)。

        相對于完全接地情況,背面有限電阻接地時正面電場增大20.0%,背面電場減小26.9%,最大電勢增大39.5%,最大電勢深度位置也相應(yīng)變化,背面電勢是非0值。這是因為電阻的存在阻礙了電荷的泄放,更因而容易產(chǎn)生高電場和高電勢,從而增加了深層放電的危險。具體數(shù)值見表2。

        表2 雙層電路板飽和電場和電勢Table 2 Saturated electric field and potential of 2-layer PCB

        當(dāng)雙層電路板背面通過電阻不完全接地時,背面電勢和最大電勢隨時間變化的趨勢如圖9所示。充電過程是一個電場和電勢快速增加、逐漸達(dá)到飽和的過程,可以用指數(shù)上升函數(shù) V=Vs[1–exp(-t/τ)]表示,其中,Vs為飽和電勢,τ為充電時間常數(shù)。可見,經(jīng)過100 h的充電,電勢達(dá)到飽和。圖9中還給出了背面電勢和最大電勢的比值隨時間的變化,發(fā)現(xiàn)這一比值固定為0.58。

        為了使有限電阻接地產(chǎn)生的電勢和最大電勢可比擬,前面計算中選用的接地電阻為1012??,F(xiàn)考察背面電阻阻值對深層充電的影響,如圖10所示。背面電勢和電阻呈線性關(guān)系,而當(dāng)電阻大于1011?時最大電勢才有明顯變化,此時背面電勢為-45.7 V,最大電勢比完全接地時增加4.6%。因此,當(dāng)接地電阻較小時(如<1011?),接地電阻對深層充電的影響可以忽略。

        2.2 中間覆銅層有限電阻接地

        考慮4層電路板正面和背面完全接地,中間覆銅層分為完全接地和電阻接地(2個中間金屬層都接1010?的電阻)2種情況,模擬得到經(jīng)過100 h輻射后的深層充電情況(如圖11所示)。相對于完全接地,中間層不完全接地時正面電場增大24.4%,背面電場增大47.3%,最大電勢增大14.3%,最大電勢深度位置也相應(yīng)地變化,具體數(shù)值見表3。4層電路板正面或背面有限電阻接地的情況與 2.1節(jié)中類似。

        表3 4層電路板飽和電場和電勢Table 3 Saturated electric field and potential of 4-layer PCB

        如前所述,在充電過程中,有限電阻端電勢和最大電勢比值不隨時間改變,針對雙層和4層電路板分別考察此比值隨接地電阻大小的變化,如圖12所示。當(dāng)電阻不是很大時(<1011?),電阻端電勢和最大電勢的比值與電阻值呈線性關(guān)系。而相對于雙層電路板,4層電路板的這一比值更大,說明有限電阻接地對4層電路板充電特性的影響更大。

        3 結(jié)束語

        本文建立了一種計算電路板覆銅層有限電阻接地時深層充電情況的方法,在對雙層和4層電路板建模的基礎(chǔ)上,以FLUMIC3.0電子通量模型為輸入,用Monte Carlo軟件包Geant4模擬電子在電路板中的輸運(yùn),先用有限元方法求解泊松方程,再用時域差分法求解電流連續(xù)性方程以得到下一時刻的電荷分布,如此反復(fù)計算直到電場和電勢達(dá)到飽和。

        為了確保仿真結(jié)果的全面性,考察了接地電阻大范圍變化下的深層充電情況。結(jié)果表明,雙層電路板背面通過1012?電阻接地時,經(jīng)過惡劣空間環(huán)境充電100 h,背面電勢為-3.6×102V,最大電勢增大39.5%,相應(yīng)的最大電勢深度較完全接地情況也有變化。相同條件下4層電路板中間層接1010?電阻,中間層電勢分別為-30 V和-14 V,最大電勢增大14.3%。表明有限電阻接地總體上會增加深層充電的電勢,從而加大放電風(fēng)險,這是因為電阻的存在阻礙了沉積電荷的泄放。另外,有限接地層電勢和最大電勢的比值不隨充電時間變化,而且接地電阻不同時,接地電勢和電阻阻值呈線性關(guān)系。以接地層電勢和最大電勢比值為判據(jù),有限電阻接地對4層電路板充電特性影響更大。綜合考慮雙層和4層電路板,當(dāng)接地電阻<109?時,其對深層充電幾乎沒有影響,這也驗證了 NASA-HDBK-4002A手冊中的設(shè)計要求,即電路板應(yīng)通過0~10 M?的電阻接地,因為該范圍大小的電阻接地不會對充電結(jié)果產(chǎn)生顯著影響;當(dāng)接地電阻較大時,其對充電的影響增大。這一數(shù)值結(jié)果可為實驗及探測儀器的設(shè)計提供參考。

        本文構(gòu)建了量化分析覆銅層電路板通過不同電阻接地的深層充電風(fēng)險評估方法,利用該方法可以針對其他材料參數(shù)、結(jié)構(gòu)和接地條件進(jìn)行分析。另外,本文所用的有限元模型只考慮厚度方向的充電特性,而覆銅層寬度對深層充電的影響可以作為下一階段研究方向。

        (References)

        [1] 張振龍, 全榮輝, 韓建偉, 等. 衛(wèi)星部件內(nèi)部充放電實驗與仿真[J]. 原子能科學(xué)技術(shù), 2010, 44(增刊 1):538-544 ZHANG Z L, QUAN R H, HAN J W, et al. Internal charging-discharging test and simulation for satellite components[J]. Atomic Energy Science and Technology,2010, 44(sup 1): 538-544

        [2] FREDERICKSON A R, BRAUTIGAM D H. Mining CRRES IDM pulse data and CRRES environmental data to improve spacecraft charging/discharging models and guidelines: NASA/CR-2004-213228[R], 2004

        [3] RYDEN K A, MORRIS P A, FORD K A, et al.Observations of internal charging currents in medium earth orbit[J]. IEEE Trans on Plasma Science, 2008,36(5): 2473-2481

        [4] TABATA T, ANDREO P, SHINODA K. An algorithm for depth-dose curves of electrons fitted to Monte Carlo data[J]. Radiation Physics and Chemistry, 1998, 53(3):205-215

        [5] FREDERICKSON A R, BELL J T, BEIDL E A. Analytic approximation for charge current and deposition by 0.1 to 100 MeV electrons in thick slabs[J]. IEEE Trans on Nuclear Science, 1995, 42(6): 1910-1921

        [6] HALBLEIB J A, KENSEK R P, VALDEZ G, et al. ITS:the integrated TIGER series of electron/photon Monte Carlo transport codes-version 3.0[J]. IEEE Trans on Nuclear Science, 1992, 39(4): 1025-1029

        [7] JUN I, GARRETT H B, KIM W, et al. Review of an internal charging code, NUMIT[J]. IEEE Trans on Plasma Science, 2008, 36(5): 2467-2472

        [8] TANG X J, YI Z, MENG L F, et al. 3-D internal charging simulation on typical printed circuit board[J].IEEE Trans on Plasma Science, 2013, 41(12): 3448-3452

        [9] 秦曉剛, 賀德衍, 王驥. 基于Geant4 的介質(zhì)深層充電電場計算[J]. 物理學(xué)報, 2009, 58(1): 684-689 QIN X G, HE D Y, WANG J. Geant4-based calculation of electric field in deep dielectric charging[J]. Acta Physica Sinica, 2009, 58(1): 684-689

        [10] 黃建國, 陳東. 不同接地方式的衛(wèi)星介質(zhì)深層充電研究[J]. 物理學(xué)報, 2004, 53(5): 1611-1616 HUANG J G, CHEN D. A study of deep dielectric charging on satellites for different grounding patterns[J].Acta Physica Sinica, 2004, 53(5): 1611-1616

        [11] 王建昭, 陳鴻飛, 于向前, 等. 多層電路板的深層充電研究[J]. 中國科學(xué): 技術(shù)科學(xué), 2015, 45(3): 330-337 WANG J Z, CHEN H F, YU X Q, et al. Study on internally dielectric charging of multilayer circuit board[J]. Science China: Technological Sciences, 2015,45(3): 330-337

        [12] 李存惠, 柳青, 史亮, 等. 一種介質(zhì)材料充放電測試設(shè)備: 中國, CN102162825A[P], 2011-08-24

        [13] KIM W, JUN I, KOKOROWSKI M. Internal electrostatic discharge monitor (IESDM)[J]. IEEE Trans on Nuclear Science, 2010, 57(6): 3143-3146

        [14] NASA. Mitigating in-space charging effects: a guideline:NASA-HDBK-4002A[S], 2011

        [15] RODGERS D J, HUNTER K A, WRENN G L. The FLUMIC electron environment model[C]//8thSpacecraft Charging Technology Conference. Huntsville, AL,2004: NASA/CP-2004-13091

        [16] ADAMEC V, CALDERWOOD J H. Electrical conduction in dielectrics at high fields[J]. J Phys D Appl Phys, 1975, 8(5): 551-560

        [17] 金建銘. 電磁場有限元方法[M]. 西安: 電子科技大學(xué)出版社, 1998: 24-34

        (編輯:張艷艷)

        A method for simulating internal charging of spacecraft circuit board grounded by electric resistance

        WANG Jianzhao1, CHEN Hongfei2, CAI Zhenbo1
        (1. Beijing Institute of Spacecraft System Engineering, Beijing 100094, China;2. Institute of Space Physics and Applied Technology, Peking University, Beijing 100871, China)

        Based on the Monte Carlo method and the finite element method (FEM), a method is proposed to calculate the internal charging process of the two-layer and four-layer printed circuit boards (PCB) with the copper layer grounded by the electric resistance. We discuss the isotropic source of electrons in the Geant4 program, the grounding boundary condition by the electric resistance, and how to build the finite element matrix function. We calculate the impact of the electric resistance on the charging process, and it is shown that the resistance can increase the electric field and the potential as compared with the fully-grounded situation.Accordingly, the maximum electric field value changes. The calculation also shows that the electric potential of the layer which is not fully-grounded is in a linear relation with the resistance. When the resistance is about 109? or smaller, its impact on the internal charging can be ignored. This result validates the applicability of guidelines in NASA-HDBK-4002A.

        internal charging; electric resistance; Monte Carlo method; finite element method; Geant4

        V542

        :A

        :1673-1379(2017)03-0258-07

        10.3969/j.issn.1673-1379.2017.03.006

        王建昭(1989—),男,碩士學(xué)位,研究方向為航天器空間環(huán)境效應(yīng)。E-mail: wangjz1989@126.com。

        2016-11-21;

        2017-05-20

        北京空間飛行器總體設(shè)計部自主研發(fā)基金

        王建昭, 陳鴻飛, 蔡震波. 電阻接地狀態(tài)下星用電路板深層充電仿真方法[J]. 航天器環(huán)境工程, 2017, 34(3):258-264

        WANG J Z, CHEN H F, CAI Z B. A method for simulating internal charging of spacecraft circuit board grounded by electric resistance[J]. Spacecraft Environment Engineering, 2017, 34(3): 258-264

        猜你喜歡
        有限元
        基于擴(kuò)展有限元的疲勞裂紋擴(kuò)展分析
        非線性感應(yīng)加熱問題的全離散有限元方法
        TDDH型停車器制動過程有限元分析
        新型有機(jī)玻璃在站臺門的應(yīng)用及有限元分析
        基于I-DEAS的履帶起重機(jī)主機(jī)有限元計算
        基于有限元模型對踝模擬扭傷機(jī)制的探討
        10MN快鍛液壓機(jī)有限元分析
        磨削淬硬殘余應(yīng)力的有限元分析
        基于SolidWorks的吸嘴支撐臂有限元分析
        箱形孔軋制的有限元模擬
        上海金屬(2013年4期)2013-12-20 07:57:18
        av人妻在线一区二区三区| 色综合自拍| 亚洲色图视频在线观看网站| av天堂一区二区三区| 亚洲综合中文字幕综合| 国模无码一区二区三区| 精品推荐国产精品店| 麻豆久久久国内精品| 在线视频一区二区国产| 少妇被粗大的猛进出69影院| 国产在线无码制服丝袜无码| 成人在线免费视频亚洲| 欧美另类在线视频| 精品熟妇av一区二区三区四区| 国产自产在线视频一区| 久久精品亚洲熟女av蜜謦| 欧美寡妇xxxx黑人猛交| 日韩AVAV天堂AV在线| av网址不卡免费在线观看| 风韵犹存丰满熟妇大屁股啪啪| 久久精品国产成人| 精品国产高清自在线99| 青青草免费在线手机视频| 亚洲综合中文字幕综合| 2021久久精品国产99国产精品| 国产成人av综合亚洲色欲| 国产91大片在线观看| 全免费a级毛片免费看无码| 无码国产午夜福利片在线观看| 国产av一区二区三区丝袜| 亚洲三级香港三级久久| 日本少妇春药特殊按摩3| а中文在线天堂| 性感人妻一区二区三区| 国产精品主播在线一区二区| 天堂aⅴ无码一区二区三区| 毛片无遮挡高清免费久久| 久久色悠悠综合网亚洲| 亚洲av国产av综合av卡| 国产人妖视频一区二区| 毛片av中文字幕一区二区 |