• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    太陽高能粒子(SEP)傳播數值模擬中的太陽風背景場研究

    2016-07-29 02:02:32魏穩(wěn)穩(wěn)沈芳左平兵秦剛楊子才
    地球物理學報 2016年3期
    關鍵詞:太陽風磁場

    魏穩(wěn)穩(wěn),沈芳,左平兵,秦剛,楊子才

    1 中國科學院國家空間科學中心 空間天氣學國家重點實驗室,北京 100190 2 中國科學院大學,北京 100049

    ?

    太陽高能粒子(SEP)傳播數值模擬中的太陽風背景場研究

    魏穩(wěn)穩(wěn)1,2,沈芳1*,左平兵1,秦剛1,楊子才1,2

    1 中國科學院國家空間科學中心 空間天氣學國家重點實驗室,北京1001902中國科學院大學,北京100049

    摘要太陽高能粒子(SEP)事件是一類重要的空間天氣災害性事件,如能準確預報SEP 事件,人們便可以采取必要的防護措施,保障衛(wèi)星、星載設備以及航天員的安全,盡可能地降低經濟損失.因此,其數值預報研究在空間天氣預報研究中占有很重要的地位.SEP 事件中的高能粒子在不同的時間尺度內被耀斑過程或者CME 驅動的激波加速,并且在被擾動后的行星際太陽風中傳輸,這些過程都緊緊依賴于太陽風背景場.因此獲取更加接近物理真實的太陽風背景場是模擬SEP 事件的重要部分,也是提高 SEP 物理模式的關鍵因素之一.我們目前的工作基于張明等發(fā)展的 SEP 在行星際空間傳播的模型,嘗試將Parker太陽風速度解及WIND飛船觀測的磁場實時數據融入模型中,研究不同的太陽風速度以及真實磁場分布對SEP在行星際空間中傳播的影響.通過求解聚焦傳輸方程,我們的模擬結果表明:(1)快太陽風條件下,絕熱冷卻效應項發(fā)揮了更大的作用,使粒子能量衰減的更快,而慢太陽風對粒子的通量變化沒有顯著影響;(2)加入觀測的磁場數據時,粒子的全向通量剖面發(fā)生了比較明顯的變化,具體表現在:通量峰值推遲到達、出現多峰結構、各向異性也發(fā)生一些改變.分析表明真實磁場的極性對粒子在行星際空間中傳播有著重要的影響.

    關鍵詞太陽高能粒子; 磁場; 太陽風; 行星際輸運

    1引言

    太陽高能粒子(Solar Energetic Particles,簡稱SEP) 事件是能量粒子的通量突然增強的事件,美國國家海洋和大氣管理局(NOAA)將SEP 事件定義為能量大于10 MeV 的粒子在連續(xù)15 min以上的時間內數目超過10 pfu(cm-2·s-1·sr-1)的爆發(fā)性事件(Rodríguez-Gasén,2011).SEP 事件主要包括三種類型(Kallenrode,2003):與太陽耀斑爆發(fā)相關聯的脈沖型事件,與日冕物質拋射驅動的激波相關聯的緩變型事件,以及同時具有緩變型和脈沖型事件特征的混合型事件.SEP 事件會對衛(wèi)星和宇航員構成嚴重的威脅,如能準確預報SEP 事件,人們便可以對衛(wèi)星和星載設備采取必要的防護措施,同時也可以為航天器的故障分析和航天員的安全防護提供科學依據,盡可能地降低損失.然而SEP 事件產生和發(fā)展的物理機制非常復雜,目前還沒有被完全解釋清楚,預報SEP 事件的能力也不盡如人意(Lario,2005).因此,正確理解和恰當地動力學描述這些機制,對推進和提高空間天氣預報能力是非常重要的(魏穩(wěn)穩(wěn)等,2015).

    SEP 事件中的高能粒子在不同的時間尺度內被耀斑過程或者CME 驅動的激波加速,并且在被擾動后的行星際磁場和太陽風中傳輸.近幾十年里,很多研究者探究了太陽風背景場對高能粒子在行星際中傳播的影響.Lario等(2003,2008)基于Ulysses以及ACE衛(wèi)星的觀測,對太陽活動高年里較低能段的高能粒子事件進行了分析,認為行星際磁力線結構會影響高能粒子的經緯度分布.并重點分析了大尺度的行星際結構對2004年9月SEP事件的影響,認為其強度和各向異性時間剖面是由一個共轉的低密度、低速、低質子β值的太陽風流決定.大尺度太陽風結構對高能粒子傳播的影響決定不同位置的宇宙飛船觀測到的SEP事件的性質,不同于低β值區(qū)域有利于高能粒子的散射,被壓縮的磁場區(qū)域會阻礙高能粒子的自由流動.磁鏡效應與散射過程能夠對高能粒子進行有效的約束,這些結構相對于觀測者和粒子源項的位置決定SEP事件的特征(Lario et al.,2008).另外,Shen 等(2008)的研究也表明磁云這種復雜結構也能對高能粒子進行有效的約束并顯著影響高能粒子的通量.

    近年來一些研究者采用聚焦傳輸方程模擬粒子在行星際空間中的傳播過程,聚焦傳輸方程包括了許多重要的粒子傳輸機制:沿磁力線的流動、太陽風對流、投擲角擴散、磁場聚焦效應、垂直擴散、絕熱冷卻效應等.模擬結果(Ruffolo,1995;Qin et al.,2006; Zhang et al.,2009)發(fā)現,在聚焦傳播方程中引入絕熱冷卻效應項,將增加粒子強度的衰減率.Lario等(Lario et al.,1997)研究了太陽風對能量為0.1~20 MeV的質子在行星際空間傳播的影響,發(fā)現太陽風對流能使粒子更早地到達觀測者,絕熱冷卻效應則導致粒子通量減小.

    SEP 事件中的高能粒子在被擾動后的行星際磁場和太陽風中傳輸,這一過程緊緊依賴于太陽風和背景場(Barouch and Burlaga,1976;Lario et al.,1997;Malandraki et al.,2007).以往粒子傳輸模型中大多采用理想的Parker螺旋磁場和恒定的太陽風速度,背景場的設置過于簡單(Heras et al.,1995; Li et al.,2003; Zhang et al.,2009; Wang et al.,2012).因此,建立更加接近物理真實的背景場來模擬粒子在三維行星際空間中的傳播是提高SEP 事件模擬結果不可或缺的部分,本文將圍繞背景場模擬展開討論,主要研究太陽風速度以及真實磁場分布對SEP在行星際空間中傳播的影響.

    2模型簡介

    SEP在日球層中的傳播是SEP研究中的一個核心問題,SEP的空間分布等都與之密切相關.研究高能粒子在行星際空間中的傳輸時需要求解聚焦傳輸方程.聚焦傳輸方程的本質是Fokker-Plank方程,用于描述粒子分布函數f在相空間的演化.有些學者(Kallenrode and Wibberenz,1997; Ruffolo,1995; Lario et al.,1997)通過有限差分格式對聚焦傳輸方程進行離散求解.Qin等(2006)和Zhang等(2009)開發(fā)并發(fā)展了一種計算SEP 在三維日球層磁場中傳播的模式,該模式采用后向隨機微分方法解聚焦傳輸方程.用于模擬研究觀測者處于行星際空間不同的經度、緯度以及徑向距離時,SEP 通量和各向異性隨時間的演化過程.

    模型采用三維的聚焦傳輸方程來研究粒子的傳播過程,其形式如下:

    (1)

    (2)

    (3)

    我們根據Qin等和Zhang等(Zhang,1999;Qinetal.,2006;Zhangetal.,2009)采用的后向隨機微分方法來求解聚焦傳輸方程.方程(1)可以改寫為五個一維的后向隨機微分方程:

    求解隨機微分方程采用后向時間的技巧,即讓虛擬粒子從指定的地點r0出發(fā),以一定的初始動量p,從時刻t=t0開始沿后向時間運動,直到粒子退出邊界.我們假設所有的高能粒子都在太陽附近注入,因此,粒子分布函數在源區(qū)滿足邊界條件:

    (4)

    其中,a(φ,θ)是一個關于日面經度φ和日面緯度θ的函數,用來描述SEP源區(qū)強度的空間分布,p是源區(qū)粒子動量,Tc和Tl分別表示源區(qū)粒子注入剖面的上升和衰減時間尺度.我們以下的模擬事例中均假設源區(qū)里面的粒子分布均勻一致,即a(φ,θ)=1.

    方程的解即為所有數值模擬過程在退出點上的平均值:

    (5)

    其中,xe,μe,pe表示虛擬粒子退出邊界時的參量,te是粒子退出邊界的時刻,N是虛擬粒子的總量.

    表1 模型中參數設定

    3模擬結果及分析

    3.1太陽風速度對粒子傳播過程的影響

    根據物理性質和起源區(qū)域的不同,太陽風一般分為快太陽風(v≥550 km·s-1)和慢太陽風(v≤450 km·s-1)(Wang et al.,2000).我們分別計算了快、慢太陽風對粒子通量的影響,快太陽風對不同能量質子的傳播過程以及對不同觀測位置的通量剖面的影響.并且統計分析了快、慢太陽風條件下統計到的粒子的動量分布特征.模型中磁場采用Parker螺旋磁場.

    為了研究發(fā)展的太陽風速度對粒子在行星際空間中傳播的影響,我們將模型中沿徑向恒定不變的太陽風速度Vsw=Vswer(Vsw=400km·s-1)變?yōu)镻arker太陽風速度解(由(6)式解出),Parker太陽風速度解的具體數值實現方法詳見文獻Parker(1958).從(6)式中可以看出,Parker太陽風速度解不再沿徑向不變,而是與日心距離r相關.

    其中Vsw滿足如下方程:

    (6)

    圖1給出了太陽風速度接近慢太陽風的模擬結果.圖1a給出了該事例中太陽風速度隨日心距離r的變化曲線,太陽風速度在1AU處達到383km·s-1.以下圖形中標注為慢太陽風的都是基于圖1a中的Parker慢太陽風速度解.圖1b給出了不同太陽風速度條件下100MeV能量的質子在1AU赤道處的全向通量對比圖.黑色實線代表固定速度400km·s-1的模擬結果,紅色虛線代表固定速度383km·s-1的模擬結果,綠色虛線代表Parker太陽風速度解接近慢太陽風速度下的模擬結果.從圖1b中可以看出這三條通量曲線非常接近,說明慢太陽風相對固定的太陽風速度而言,對粒子的通量變化影響不大.

    圖1 太陽風速度接近慢太陽風的模擬結果(a) 太陽風速度隨日心距離的變化; (b) 不同太陽風條件下能量為100 MeV的質子的全向通量對比.

    圖2給出了太陽風速度接近快太陽風的模擬結果.圖2a中Parker太陽風速度解在1AU處達到644km·s-1.以下圖形中標注為快太陽風的都是基于圖2a中的Parker快太陽風速度解.圖2b給出了不同太陽風速度條件下100MeV能量的質子在1AU赤道處的全向通量對比圖.黑色實線代表固定速度400km·s-1的模擬結果,紅色虛線代表固定速度644km·s-1的模擬結果,綠色虛線代表Parker太陽風速度解接近快太陽風速度下的模擬結果.從圖2b中可以看出,固定太陽風速度變大時,將增加粒子通量曲線的衰減率,而Parker快太陽風條件下,粒子通量曲線的衰減率進一步增加.通過分析我們認為這是由于快太陽風條件下,絕熱冷卻效應項發(fā)揮了更大的作用.在Parker太陽風速度解的模型中,絕熱冷卻項dp/dt的表達式如下:

    圖2 太陽風速度接近快太陽風的模擬結果(a)太陽風速度隨日心距離的變化; (b) 不同太陽風條件下能量為100 MeV的質子的全向通量對比.

    (7)

    由于在SEP的初始相時期粒子密度的梯度大,擴散效應比絕熱冷卻效應更重要,所以在SEP的初始相可以忽略絕熱效應,在SEP的衰減時期,粒子密度空間梯度小,絕熱冷卻效應將起重要的作用(涂傳詒等,1988).在Parker太陽風速度解的條件下,對于某一固定的徑向距離r,太陽風速度增大時.相應的Vsw/r和dVsw/dr也隨之增大,第一項和第二項與太陽風速度的大小呈線性關系,因此,太陽風速度越高,會造成SEP向外對流的越多,絕熱能量損失也隨之增大,絕熱效應也就越顯著.

    我們又對快、慢太陽風速度條件下統計到的粒子的動量進行了統計分析.圖3給出了觀測到的粒子的動量分布圖.黑色實線,紅色虛線分別代表Parker太陽風速度解接近快太陽風速度、慢太陽風速度時的統計結果.我們將觀測到的粒子動量的最大值和最小值(pmax&pmin)等分為10個區(qū)間Δp= 0.1(pmax-pmin),然后統計每個動量區(qū)間內的粒子數占總的粒子數的百分比.表2為不同太陽風條件下的統計結果,我們設置的目標粒子的動量為p=0.445.其中目標粒子的動量定義為觀測者處最終觀測到的粒子動量,通過能量計算可以得到,該能量是人為設定的,是為了便于統計分析某一個具體能量段內的粒子的傳播特征、分布特征等.

    圖3 快慢太陽風條件下統計到的粒子動量分布

    從圖3和表2的結果可以看出:

    表2 不同太陽風條件下的統計結果

    (1) 從所統計到粒子的最低動量來看,它們都比目標粒子對應的動量高,說明粒子在運動過程中都有能量衰減;

    (2) 從所統計到粒子的最高動量來看,快太陽風條件下粒子能量的最大值更大,說明快太陽風條件下粒子能量的衰減可以更加顯著;

    (3) 投放相同的粒子數(投放1×106個粒子)時,快太陽風條件下統計到的粒子數小于慢太陽風條件下的粒子數,也驗證了上述結論;

    (4) 從粒子分布來看,越靠近目標粒子能量段的粒子數越多,且呈現出單調遞減的趨勢,說明大部分統計到的粒子是能量衰減較少的粒子.

    此外,我們又模擬了快太陽風對不同能量的粒子和不同觀測位置處的通量剖面的影響.圖4a、4b分別給出了不同太陽風速度條件下100 MeV能量的質子在觀測者位于0.5AU赤道處以及2AU、緯度80°處的全向通量對比.圖4c給出了不同太陽風速度條件下10 MeV能量的質子在1AU赤道處的全向通量對比圖.黑色實線代表固定速度為400 km·s-1的模擬結果,紅色虛線代表Parker太陽風速度解接近快太陽風速度下的模擬結果.從結果中可以得出,快太陽風對不同能量段的粒子在不同觀測位置處的通量剖面有著相似的影響.

    圖4 快太陽風條件下觀測者處于不同位置時能量為100 MeV的質子的全向通量對比((a)和(b))以及能量為10 MeV的質子的全向通量對比(c)

    3.2用實時的磁場觀測數據

    由1AU處的磁場估算得到.

    (8)

    假設測試粒子從太陽源表面(Rs,φs)離開,到達空間任意位置(r,φ),并且太陽風在每一個卡林頓周期內是不變的,則可以根據宇宙飛船的觀測數據估算Rs處太陽風的參量.參數r,Rs,φ,φs,觀測時間t,區(qū)間[t0,t0+T]的開始時間t0由以下公式計算得到:

    (9)

    模型中使用的是WIND飛船一分鐘精度的磁場數據,我們分別計算了處于太陽活動低年以及太陽活動高年的不同卡林頓自轉周期(簡稱CR)內SEP的全向通量和各向異性的時間剖面.

    太陽活動低年的模擬結果如圖5所示.其中,圖5b是根據WIND飛船在2007年5月3日至2007年5月29日處于CR2056和CR2057內的觀測數據得到的黃道面內的磁場圖.為了更好的可視化效果,我們將任意日心距離r處的磁場進行了標準化處理,具體標準化處理方法為:將同一條Parker螺線上,任意日心距離r處的磁場,均賦為該條螺線在1AU處的磁場數值,然后進行可視化作圖.太陽坐標為(X,Y)=(rcosφ,rsinφ)=(0,0),地球處在黑色圓圈處,并按照順時針方向進行背景磁場的設置.兩條黑線的夾角為統計到的粒子所在源區(qū)的大致分布范圍.圖5a中紅色虛線為在此CR期間內根據類Parker磁場模擬得到的粒子全向通量剖面,黑色實線為原Parker磁場下的模擬結果,圖5c為相應的各向異性剖面.相比之下,我們發(fā)現在類Parker磁場下模擬得到的粒子通量峰值會到達的晚一些,全向通量剖面會出現多個峰值并且在衰減相下降的更快,而且一階各向異性也有明顯的不同.

    圖5 (a) 不同行星際磁場下,對比100 MeV質子的通量; (c) 各向異性; (b) 根據WIND飛船數據得到的CR2056-CR2057內磁場Br分量圖.太陽坐標為(X,Y)=(rcosφ,rsinφ)=(0,0),黑色圓圈為地球

    由式子(10)—(12)可知,當加入磁場的觀測數據后,cosβ和sinβ的大小不會改變,但是它們的符號會隨著磁場的極性而發(fā)生相應的改變,從而會影響投擲角μ值的變化,并最終對動量p的值產生影響.因此,我們認為粒子通量峰值推遲到達,多峰結構以及各向異性的差異主要是由于加入的真實磁場

    具有不同的極性所致.

    針對通量峰值推遲到達,我們比較了Parker螺旋磁場和類Parker螺旋磁場下最大峰值附近的粒子初始能量的大小,發(fā)現Parker磁場下粒子初始能量的平均值大約為0.463,類Parker磁場下粒子初始能量的平均值大約為0.628,而我們統計的目標粒子能量為0.445.由于能量越高,衰減到目標能量所需要的時間越長,所以在類Parker磁場中由于磁場極性對粒子能量的調制,它們到達觀測者所需的時間會更久,從而導致峰值推遲到達.

    針對出現的多峰結構,我們在圖6中給出了第二個較大峰值到達前后不同時間段內統計到的粒子初始動量p的分布范圍.黑色實線、紅色虛線、綠色虛線分別代表第二個峰值到達前兩天,第二個峰值期間以及第二個峰值過后兩天內的粒子動量分布.從圖6可以看出,不同時間段內動量的分布存在明顯差異.在第二個峰出現的時間段內(紅色虛線),p=1附近出現了另一個峰值.通過對通量(為p-γ的函數)的計算,我們發(fā)現這個峰值對該期間的粒子通量的貢獻起決定性作用,并且使得該通量超過了前后兩天的通量大小.通常情況下隨著時間的推移,在衰減相不同時間段內統計到的大部分粒子的動量會不斷增大,即類似圖6中黑色實線的峰值會不斷后移,從而使得全向通量不斷減小.然而,同樣,由于磁場極性對粒子能量的調制,使得粒子動量出現了額外的峰值并發(fā)揮了主要作用,從而導致出現了多峰結構.

    圖6 CR2056-CR2057模擬得到的不同時間段內到達觀測者的粒子初始動量p的分布范圍

    (10)

    (11)

    (12)

    太陽活動高年的模擬結果,我們以CR1969為例展示在了圖7中.類似圖5的格式,圖7a給出了模擬得到的粒子全向通量剖面,圖7c為各向異性剖面,圖7b給出了WIND飛船在該周期內所觀測到的磁場.其中兩條黑色線同樣代表統計到的粒子所在源區(qū)的大致范圍.類似太陽活動小年的情況,從圖7a和圖7c中可以看出,粒子通量峰值也比原Parker磁場到達的晚一些,粒子的全向通量剖面出現了多峰結構,而各向異性在誤差范圍內并沒有太大的差異.

    圖7 (a) 不同行星際磁場下,對比100 MeV質子的通量; (c) 各向異性;(b) 根據WIND飛船數據得到的CR1969的磁場Br分量

    同樣地,我們計算得到類Parker磁場下最大峰值附近所統計到的粒子初始能量的平均值大約為0.562,要大于Parker磁場下粒子能量的平均值0.463,因此粒子在類Parker螺線磁場里運動到達觀測者所需的時間更久,導致峰值推遲到達.但是,相比于太陽活動小年粒子的平均初始能量為0.628,該周期內統計到的粒子初始能量要更小一些,所以相比之下粒子的全向通量剖面峰值到達的也更早一些.至于更小的粒子初始能量,可能與太陽活動高年更復雜的磁場極性變化對粒子能量的調制有關.

    此外,與太陽活動小年相比,在太陽活動高年里模擬得到的粒子的全向通量剖面會出現更多的峰值.圖8給出了第二、三個峰值到達時以及它們到達之前一段內觀測到的粒子初始動量p的分布范圍.從圖8可以看出,在這兩個峰值到達之前,如黑線所示,粒子初始動量大約以0.8為峰值,而且粒子的動量分布范圍更加廣泛.在紅線所示的第二個峰值到 達期間,粒子的初始動量更加集中地分布在以0.6左右為峰值的一個小范圍內,因此根據通量計算公式可知該段時間內粒子的全向通量更大.而在綠線所示的第三個峰值到達期間,粒子的初始動量仍以0.6左右為峰值,但是它們占該時間段內所統計到總粒子數的比例更小,而且來源更加廣泛,這就是為什么第三個峰值要比第二個峰值低的原因.類似于太陽活動小年的情況,在太陽活動高年里磁場對粒子能量的調制會導致多峰結構,而這種更多峰值結構的出現,是由于磁場的極性出現了更多的變化,對粒子能量的調制更加復雜所致.

    圖8 CR1969模擬得到的不同時間段內到達觀測者的粒子初始動量p的分布范圍

    從圖7c中可以看出,類Parker磁場下的各向異性在上升相位時大于零,說明朝磁場相同方向運動(μ=1)的粒子占主導.同樣地,在衰減相期間,由于磁場極性的影響以及粒子投擲角方向的改變,會使得各向異性在誤差范圍內出現一些波動.

    從模擬結果可以看出,太陽活動高年和低年的結果雖然都包含通量峰值推遲到達、出現多峰結構、各向異性發(fā)生改變等特點,但也存在差異.由于太陽活動高年磁場比較復雜,導致通量剖面出現更多的峰,各向異性剖面也存在差異.但二者的模擬結果都充分說明了真實磁場極性對粒子在行星際空間中的傳播有著重要的影響.

    4結論

    本文通過求解聚焦傳輸方程研究了太陽風速度和背景磁場對粒子在行星際空間中傳播的影響,模擬結果表明:

    (1) 將固定的太陽風速度改為變化的Parker太陽風速度解,從模擬結果可以看出慢太陽風對粒子的通量變化沒有顯著影響,而當太陽風速度加快時,粒子通量剖面在后期下降的更快.這是因為快太陽風條件下,絕熱冷卻效應更加顯著,從而使粒子的能量衰減的更快.

    (2) 相對于原Parker磁場模型,加入觀測的磁場數據時粒子的全向通量剖面發(fā)生了比較明顯的變化:通量峰值推遲到達、出現多峰結構,各向異性也發(fā)生一些改變.通過分析,我們認為加入觀測的磁場數據時,磁場有了極性,粒子的投擲角會隨著磁場的極性而發(fā)生相應改變,進而也會對粒子的能量進行調制,從而導致模擬結果出現了上述的變化.無論在太陽活動低年還是高年,粒子全向通量剖面和各向異性剖面都伴隨著這些改變,而太陽活動高年更加復雜的變化可能是由于磁場的極性更加復雜所致.

    我們的模擬結果說明太陽風速度以及背景磁場對粒子在行星際空間中的傳播有著重要的影響.以往粒子傳輸模型中大多采用恒定的太陽風速度和Parker螺旋磁場,背景場的設置太過理想,這可能導致無法準確模擬太陽高能粒子在行星際空間中的演化過程.因此,設置更加真實的背景場對提高我們的模擬結果有著重要的意義.未來的研究工作中我們將嘗試加入由MHD模擬(Fengetal.,2010,2012;Shenetal.,2011,2014)得到的實時變化的太陽風速度和三維磁場分量,建立更加完整、更加接近物理真實的太陽風背景場,以期得到更真實的粒子在三維行星際空間的傳播物理過程,同時可以研究日冕物質拋射相互作用等對高能粒子事件形成和傳播的影響(Gopalswamyetal.,2002;Shenetal.,2008;Lietal.,2012;Shenetal.,2013).

    References

    Barouch E,Burlaga L F.1976.Three-dimensional interplanetary stream magnetism and energetic particle motion.J.Geophys.Res.,81(13): 2103-2110,doi: 10.1029/JA081i013p02103.

    Feng X S,Yang L P,Xiang C Q,et al.2010.Three-dimensional solar WIND modeling from the sun to earth by a SIP-CESE MHD model with a six-component grid.Astrophys.J.,723(1): 300-319,doi: 10.1088/0004-637X/723/1/300.

    Feng X S,Jiang C W,Xiang C Q,et al.2012.A data-driven model for the global coronal evolution.Astrophys.J.,758(1),doi: 10.1088/0004-637X/758/1/62.

    Florens M S L,Cairns I H,Knock S A,et al.2007.Data-driven solar wind model and prediction of type II bursts.Geophys.Res.Lett.,34(4): L04104,doi: 10.1029/2006GL028522.Gopalswamy N S,Yashiro S,Michaek G,et al.2002.Interacting coronal mass ejections and solar energetic particles.Astrophys.J.,572(1): L103-L107,doi: 10.1086/341601.

    Heras A M,Sanahuja B,Lario D,et al.1995.Three low-energy particle events: modeling the influence of the parent interplanetary shock.Astrophys.J.,445(1): 497-508,doi: 10.1086/175714.

    Kallenrode M B.2003.Current views on impulsive and gradual solar energetic particle events.Journal of Physics G: Nuclear and Particle Physics,29(5): 965-981,doi: 10.1088/0954-3899/29/5/316.

    Kallenrode M B,Wibberenz G.1997.Propagation of particles injected from interplanetary shocks: A black box model and its consequences for acceleration theory and data interpretation.J.Geophys.Res.,102(A10): 22311-22334,doi: 10.1029/97JA01677.

    Lario D,Sanahuja B,Heras A M.1997.Modeling the interplanetary propagation of 0.1~20 MeV shock-accelerated protons.I: Effects of the adiabatic deceleration and solar wind convection.Adv.Space Res.,20(1): 115-120,doi: 10.1016/S0273-1177(97)00492-4.

    Lario D,Roelof E C,Decker R B,et al.2003.Solar maximum low-energy particle observations at heliographic latitudes above 75 degrees.Advances in Space Research,32(4): 579-584,doi: 10.1016/S0273-1177(03)00339-9.

    Lario D.2005.Advances in modeling gradual solar energetic particle events.Advances in Space Research,36(12): 2279-2288,doi: 10.1016/j.asr.2005.07.081.

    Lario D,Decker R B,Malandraki O E,et al.2008.Influence of large-scale interplanetary structures on energetic particle propagation: September 2004 event at Ulysses and ACE.J.Geophys.Res.,113(A3): A03105,doi: 10.1029/2007JA012721.

    Li G,Zank G P,Rice W K M.2003.Energetic particle acceleration and transport at coronal mass ejection-driven shocks.J.Geophys.Res.,108(A2): 1082,doi: 10.1029/2002JA009666.Li G,Moore R,Mewaldt R A,et al.2012.A Twin-CME scenario for ground level enhancement events.Space Science Reviews,171(1-4): 141-160,doi: 10.1007/s11214-011-9823-7.

    Malandraki O E,Marsden R G,Tranquille C,et al.2007.Energetic particle observations by Ulysses during the declining phase of solar cycle 23.J.Geophys.Res.,112(A6): A06111,doi: 10.1029/2006JA011876.

    Parker E N.1958.Dynamics of the interplanetary gas and magnetic fields.Astrophys.J.,128: 664-676,doi: 10.1086/146579.

    Qin G,Zhang M,Dwyer J R.2006.Effect of adiabatic cooling on the fitted parallel mean free path of solar energetic particles.J.

    Ruffolo D.1995.Effect of adiabatic deceleration on the focused transport of solar cosmic rays.Astrophys.J.,442(2): 861-874,doi: 10.1086/175489.

    Shen C L,Wang Y M,Ye P Z,et al.2008.Enhancement of solar energetic particles during a shock-magnetic cloud interacting complex structure.Solar Physics,252(2): 409-418,doi: 10.1007/s11207-008-9268-7.

    Shen C L,Li G,Kong X,et al.2013.Compound twin coronal mass ejections in the 2012 May 17 GLE event.Astrophys.J.,763(2): 114-121,doi: 10.1088/0004-637X/763/2/114.

    Shen F,Feng X S,Wu S T,et al.2011.Three-dimensional MHD simulation of the evolution of the April 2000 CME event and its induced shocks using a magnetized plasma blob model.J.Geophys.Res.,116(A4): A04102,doi: 10.1029/2010JA015809.Shen F,Shen C L,Zhang J,et al.2014.Evolution of the 12 July 2012 CME from the Sun to the Earth: Data-constrained three-dimensional MHD simulations.J.Geophys.Res.,119(9): 7128-7141,doi: 10.1002/2014JA020365.

    Wang Y,Qin G,Zhang M.2012.Effects of perpendicular diffusion on energetic particles accelerated by the interplanetary coronal mass ejection shock.Astrophys.J.,752(1): 37,doi: 10.1088/0004-637X/752/1/37.

    Wang Y M,Sheeley N R Jr,Socker D G,et al.2000.The dynamical nature of coronal streamers.J.Geophys.Res.,105(A11): 25133-25142,doi: 10.1029/2000JA000149.

    Wei W W,Shen F,Zuo P B.2015.Research progress on the solar energetic particle model based on magnetohydrodynamic simulation.Progress in Astronomy (in Chinese),33(1): 1-26,doi: 10.3969/j.issn.1000-8349.2015.01.01.

    Zhang M.1999.A Markov stochastic process theory of cosmic-ray modulation.Astrophys.J.,513(1): 409-420,doi: 10.1086/306857.

    Zhang M,Qin G,Rassoul H.2009.Propagation of solar energetic particles in three-dimensional interplanetary magnetic fields.Astrophys.J.,692(1): 109-132,doi: 10.1088/0004-637X/692/1/109.

    附中文參考文獻

    涂傳詒等.1988.日地空間物理學.北京: 科學出版社.

    魏穩(wěn)穩(wěn),沈芳,左平兵.2015.基于磁流體力學模擬的太陽高能粒子物理模式研究進展.天文學進展,33(1): 1-26,doi: 10.3969/j.issn.1000-8349.2015.01.01.

    (本文編輯胡素芳)

    基金項目國家重點基礎研究專項經費資助項目(2012CB825601),中國科學院知識創(chuàng)新工程重大項目(KZZD-EW-01-4),國家自然科學基金(41174150,41174152,41374188,41474152)聯合資助.

    作者簡介魏穩(wěn)穩(wěn),女,1989年生,博士研究生,方向為太陽高能粒子傳播背景場的研究.E-mail: wwwei@spaceweather.ac.cn E-mail: fshen@spaceweather.ac.cn

    *通訊作者沈芳,研究員,主要從事日地空間背景太陽風結構以及行星際擾動傳播過程的三維MHD數值模擬.

    doi:10.6038/cjg20160301 中圖分類號P354 Geophys.Res.,111(A8): A08101,10.1029/2005JA011512.Rodríguez-Gasén R.2011.Modelling SEP events: latitudinal and longitudinal dependence of the injection rate of shock-accelerated protons and their flux profiles [Ph.D.thesis].Barcelona: Universitatde Barcelona.

    收稿日期2015-06-12,2015-08-31收修定稿

    Effects of the solar wind background field on the numerical simulationof the Solar Energetic Particle (SEP) transportation

    WEI Wen-Wen1,2,SHEN Fang1*,ZUO Ping-Bing1,QIN Gang1,YANG Zi-Cai1,2

    1StateKeyLaboratoryofSpaceWeather,NationalSpaceScienceCenter,ChineseAcademyofSciences,Beijing100190,China2UniversityofChineseAcademyofSciences,Beijing100049,China

    AbstractSolar energetic particles (SEPs) pose one of the most serious hazards to spacecraft systems and constrain human activities in space.Thus,it is of importance to forecast SEP events.Several theories and numerical models are applied to simulate SEP events.Each model makes some assumptions to simplify the complex acceleration and transportation processes within such events.In general,SEP will interact with ambient solar wind and background magnetic field during transportation.It is recognized that interplanetary transport effects must be taken into account at any analysis of SEP propagation.In the previous simulation,it always assumed Parker magnetic field and fixed solar wind speed as the input parameters.However,these assumptions are too simple when compared with the real conditions.In order to get better results,it is necessary to use more accurate background conditions.Recently,we change the fixed solar wind speed into spatial-dependent speed profile based on Parker′s theory,and replace the Parker magnetic field with another Parker-like magnetic field based on in situ data at 1 AU.By solving the focused transport equation with simulation of time-backward stochastic processes method,our results show that:(1) Under fast solar wind speed assumption,it is clear that the omnidirectional flux decreases faster than that for the situation with slow solar wind speed in the decay phase.We suggest that it is due to the adiabatic cooling effect.Fast solar wind speed has a significant effect on the adiabatic cooling,which leads the SEPs to lose energy more quickly during transportation.However,slow solar wind speed has less impact on the time profiles of SEP flux and anisotropy.We also compare the time profiles of SEP event observed at different observatories and energies,the results remain the same as previous; (2) When applying in situ data of magnetic field observed by WIND during different Carrington Rotations,the omnidirectional flux time profiles vary greatly,and the main results are as followings: the peak flux appears to be delayed,multi-peak occur,anisotropy also has some differences.We think it results from the magnetic field polarity,which affects the pitch angle,and,furthermore,modulates the momentum.The characteristics are similar in solar minimum and solar maximum,while the peaks seem to be more when solar activity is active.We conclude that the real magnetic field polarity may exert a significant influence during the propagation of SEP.In the future,we will try to use the real-time background conditions which obtain from MHD models in our simulations,in order to make a thorough study of the SEP propagation.

    KeywordsSolar energetic particle; Magnetic field; Solar wind; Interplanetary transport

    魏穩(wěn)穩(wěn),沈芳,左平兵等.2016.太陽高能粒子(SEP)傳播數值模擬中的太陽風背景場研究.地球物理學報,59(3):767-777,doi:10.6038/cjg20160301.

    Wei W W,Shen F,Zuo P B,et al.2016.Effects of the solar wind background field on the numerical simulation of the Solar Energetic Particle (SEP) transportation.Chinese J.Geophys.(in Chinese),59(3):767-777,doi:10.6038/cjg20160301.

    猜你喜歡
    太陽風磁場
    西安的“磁場”
    當代陜西(2022年6期)2022-04-19 12:11:54
    為什么地球有磁場呢
    應賀共青詩社成立(新韻)
    文脈清江浦 非遺“磁場圈”
    華人時刊(2020年13期)2020-09-25 08:21:42
    多種觀測數據驅動的三維行星際太陽風MHD模擬
    基于ACE飛船觀測的銀河宇宙線與太陽風變化的統計研究
    《磁場》易錯易混知識剖析
    磁場的性質和描述檢測題
    聽,太空的聲音
    2016年春季性感磁場
    Coco薇(2016年1期)2016-01-11 16:53:24
    人人妻人人看人人澡| 在线观看一区二区三区激情| 久久久久人妻精品一区果冻| 少妇被粗大猛烈的视频| 国产精品久久久久久久电影| 国产精品一区二区性色av| 国产精品一区二区在线观看99| 蜜桃在线观看..| 国产又色又爽无遮挡免| 中文资源天堂在线| 久久久久久久久久成人| 夜夜爽夜夜爽视频| 18禁动态无遮挡网站| 韩国高清视频一区二区三区| 啦啦啦啦在线视频资源| 欧美精品亚洲一区二区| 女人十人毛片免费观看3o分钟| 噜噜噜噜噜久久久久久91| 有码 亚洲区| 久久影院123| 搡老乐熟女国产| 久久ye,这里只有精品| 一级毛片电影观看| 伦理电影免费视频| 久久久精品94久久精品| 91久久精品国产一区二区三区| 麻豆国产97在线/欧美| 国产精品欧美亚洲77777| 日产精品乱码卡一卡2卡三| 波野结衣二区三区在线| 日本色播在线视频| 热re99久久精品国产66热6| 欧美日韩一区二区视频在线观看视频在线| 久久这里有精品视频免费| 成人美女网站在线观看视频| 精品人妻一区二区三区麻豆| 涩涩av久久男人的天堂| 国产一区有黄有色的免费视频| 久久精品国产亚洲av天美| 直男gayav资源| 尤物成人国产欧美一区二区三区| 国产久久久一区二区三区| 亚洲国产精品国产精品| 又爽又黄a免费视频| 久久精品国产a三级三级三级| 午夜精品国产一区二区电影| 国产一级毛片在线| 91精品国产九色| 久久精品国产鲁丝片午夜精品| 草草在线视频免费看| 国产午夜精品久久久久久一区二区三区| 亚洲精品日韩在线中文字幕| 在线观看免费日韩欧美大片 | 日韩不卡一区二区三区视频在线| 美女国产视频在线观看| 九九久久精品国产亚洲av麻豆| 日本午夜av视频| 少妇人妻 视频| 自拍欧美九色日韩亚洲蝌蚪91 | 麻豆国产97在线/欧美| 日韩av免费高清视频| 狂野欧美激情性xxxx在线观看| 久久国产亚洲av麻豆专区| 一级毛片我不卡| 校园人妻丝袜中文字幕| 国产亚洲午夜精品一区二区久久| 日本猛色少妇xxxxx猛交久久| 一级二级三级毛片免费看| 国产高清国产精品国产三级 | 欧美97在线视频| 亚洲丝袜综合中文字幕| 伦精品一区二区三区| 人人妻人人澡人人爽人人夜夜| 国产熟女欧美一区二区| av网站免费在线观看视频| 91午夜精品亚洲一区二区三区| 国内少妇人妻偷人精品xxx网站| 三级国产精品欧美在线观看| 天天躁夜夜躁狠狠久久av| 在线观看免费日韩欧美大片 | 美女cb高潮喷水在线观看| 国产亚洲精品久久久com| 韩国av在线不卡| 国产欧美日韩一区二区三区在线 | 18禁裸乳无遮挡动漫免费视频| 欧美精品国产亚洲| 99热网站在线观看| 亚洲国产欧美人成| 激情 狠狠 欧美| 人人妻人人添人人爽欧美一区卜 | 99久久精品国产国产毛片| 成人18禁高潮啪啪吃奶动态图 | 最近2019中文字幕mv第一页| 91精品国产九色| 在线观看人妻少妇| 80岁老熟妇乱子伦牲交| 色5月婷婷丁香| 亚洲,欧美,日韩| 五月开心婷婷网| 日本免费在线观看一区| 天堂中文最新版在线下载| 亚洲av国产av综合av卡| 久久av网站| 一个人看视频在线观看www免费| 黄色配什么色好看| 在线天堂最新版资源| 天天躁夜夜躁狠狠久久av| 亚洲欧美精品专区久久| 久久99蜜桃精品久久| 99久久精品热视频| 五月玫瑰六月丁香| 免费久久久久久久精品成人欧美视频 | 亚洲不卡免费看| 日韩亚洲欧美综合| 人妻少妇偷人精品九色| 久热久热在线精品观看| 欧美日韩亚洲高清精品| 在线观看一区二区三区激情| 久久久欧美国产精品| 久久久久国产网址| 波野结衣二区三区在线| 亚洲激情五月婷婷啪啪| 草草在线视频免费看| 哪个播放器可以免费观看大片| 亚洲精品久久久久久婷婷小说| 一区在线观看完整版| 一本—道久久a久久精品蜜桃钙片| 国产精品欧美亚洲77777| 久久精品久久精品一区二区三区| 男人和女人高潮做爰伦理| 麻豆乱淫一区二区| 99热这里只有精品一区| 亚洲图色成人| 人妻 亚洲 视频| 国产欧美亚洲国产| 国产永久视频网站| 国内揄拍国产精品人妻在线| 色哟哟·www| 麻豆精品久久久久久蜜桃| 成年av动漫网址| 午夜视频国产福利| 中国美白少妇内射xxxbb| 网址你懂的国产日韩在线| 寂寞人妻少妇视频99o| 26uuu在线亚洲综合色| 亚洲内射少妇av| 欧美zozozo另类| 日韩亚洲欧美综合| 亚洲av不卡在线观看| 成人一区二区视频在线观看| 精品少妇黑人巨大在线播放| 亚洲电影在线观看av| 91aial.com中文字幕在线观看| 亚洲一级一片aⅴ在线观看| 99九九线精品视频在线观看视频| 国产免费一区二区三区四区乱码| 97精品久久久久久久久久精品| 亚洲图色成人| 国产毛片在线视频| 肉色欧美久久久久久久蜜桃| 国产精品国产三级专区第一集| 国内揄拍国产精品人妻在线| 搡老乐熟女国产| 男人爽女人下面视频在线观看| 热re99久久精品国产66热6| 精品亚洲成a人片在线观看 | 国产精品免费大片| av在线播放精品| 免费观看无遮挡的男女| a 毛片基地| 欧美高清性xxxxhd video| 久久国产精品男人的天堂亚洲 | 婷婷色av中文字幕| 久久精品国产鲁丝片午夜精品| 国产精品一区www在线观看| 成年人午夜在线观看视频| 妹子高潮喷水视频| 80岁老熟妇乱子伦牲交| 中文资源天堂在线| h日本视频在线播放| 国产女主播在线喷水免费视频网站| 欧美精品人与动牲交sv欧美| 在线观看一区二区三区| 99热网站在线观看| 亚洲精品乱久久久久久| 综合色丁香网| 我要看黄色一级片免费的| 狠狠精品人妻久久久久久综合| 在线播放无遮挡| 在线天堂最新版资源| 亚洲精品日韩av片在线观看| 九色成人免费人妻av| 精品久久久精品久久久| 国产精品久久久久久久久免| 成人综合一区亚洲| 免费大片黄手机在线观看| 午夜日本视频在线| 99久久中文字幕三级久久日本| 国精品久久久久久国模美| 一级毛片我不卡| 狠狠精品人妻久久久久久综合| 亚洲精品一二三| 肉色欧美久久久久久久蜜桃| 午夜福利网站1000一区二区三区| 国产精品成人在线| 日本vs欧美在线观看视频 | 99久久中文字幕三级久久日本| 国产精品一区二区在线不卡| 亚洲欧美清纯卡通| 欧美日韩亚洲高清精品| 99久久精品一区二区三区| 久久女婷五月综合色啪小说| 久久99蜜桃精品久久| 舔av片在线| 一本久久精品| 国产精品国产三级国产专区5o| 精品亚洲成国产av| kizo精华| av在线蜜桃| 菩萨蛮人人尽说江南好唐韦庄| 大码成人一级视频| 亚洲成人手机| 久久久久久人妻| 久久久欧美国产精品| 日本-黄色视频高清免费观看| 日本午夜av视频| av又黄又爽大尺度在线免费看| 99久久综合免费| 丝瓜视频免费看黄片| 久久久久久久久久成人| 日韩av不卡免费在线播放| 日韩三级伦理在线观看| 最近中文字幕高清免费大全6| 免费看av在线观看网站| 18+在线观看网站| 交换朋友夫妻互换小说| av免费在线看不卡| 久久99热这里只频精品6学生| 建设人人有责人人尽责人人享有的 | 91精品国产九色| 国产探花极品一区二区| 午夜福利在线在线| 欧美另类一区| 麻豆成人午夜福利视频| 亚洲国产毛片av蜜桃av| 欧美激情国产日韩精品一区| 黄片wwwwww| 97精品久久久久久久久久精品| 妹子高潮喷水视频| 国产精品99久久久久久久久| 哪个播放器可以免费观看大片| 综合色丁香网| 亚洲av不卡在线观看| 久久鲁丝午夜福利片| 亚洲精品日韩在线中文字幕| 男女国产视频网站| 国产乱来视频区| 亚洲无线观看免费| 少妇人妻一区二区三区视频| 91久久精品国产一区二区三区| 免费黄色在线免费观看| 联通29元200g的流量卡| 国产女主播在线喷水免费视频网站| 国产成人精品婷婷| 日本欧美国产在线视频| 最近最新中文字幕大全电影3| av福利片在线观看| 亚洲综合色惰| 性色avwww在线观看| 少妇精品久久久久久久| 日本免费在线观看一区| 中文精品一卡2卡3卡4更新| 亚洲四区av| 在线观看一区二区三区激情| 男女免费视频国产| 国产成人精品婷婷| 女性被躁到高潮视频| 中文字幕久久专区| 永久免费av网站大全| 有码 亚洲区| 99久国产av精品国产电影| 国产伦精品一区二区三区四那| 天天躁夜夜躁狠狠久久av| 舔av片在线| 特大巨黑吊av在线直播| 女性生殖器流出的白浆| 少妇猛男粗大的猛烈进出视频| 在线看a的网站| 十分钟在线观看高清视频www | 国产 一区精品| 天美传媒精品一区二区| 欧美极品一区二区三区四区| www.av在线官网国产| 男女下面进入的视频免费午夜| 97热精品久久久久久| 国产成人a区在线观看| 日本-黄色视频高清免费观看| 欧美精品国产亚洲| 国产精品久久久久久精品电影小说 | 成人亚洲欧美一区二区av| 爱豆传媒免费全集在线观看| 我要看黄色一级片免费的| 丝瓜视频免费看黄片| 我的老师免费观看完整版| 男女边吃奶边做爰视频| 亚洲欧美一区二区三区国产| av视频免费观看在线观看| 天天躁夜夜躁狠狠久久av| 舔av片在线| 中文字幕制服av| 岛国毛片在线播放| www.色视频.com| 久久久久精品性色| kizo精华| av国产免费在线观看| 亚洲aⅴ乱码一区二区在线播放| 丝袜喷水一区| 久久热精品热| 99九九线精品视频在线观看视频| 欧美成人精品欧美一级黄| 成人亚洲欧美一区二区av| 搡女人真爽免费视频火全软件| 女人十人毛片免费观看3o分钟| 直男gayav资源| 欧美精品一区二区大全| 国产精品欧美亚洲77777| 男人狂女人下面高潮的视频| 国产在线免费精品| 亚洲精品色激情综合| 国产片特级美女逼逼视频| 久久这里有精品视频免费| 国产午夜精品久久久久久一区二区三区| 亚洲aⅴ乱码一区二区在线播放| 精品亚洲成国产av| 国产爱豆传媒在线观看| 欧美性感艳星| 国产91av在线免费观看| 网址你懂的国产日韩在线| 国产亚洲欧美精品永久| 国产精品一二三区在线看| 亚洲经典国产精华液单| 国产精品人妻久久久影院| 亚洲欧美日韩另类电影网站 | 国产真实伦视频高清在线观看| 一区二区三区四区激情视频| 极品少妇高潮喷水抽搐| 精品人妻偷拍中文字幕| 国产高清国产精品国产三级 | 欧美bdsm另类| 美女cb高潮喷水在线观看| 18+在线观看网站| www.色视频.com| 在线免费十八禁| 欧美人与善性xxx| 日韩制服骚丝袜av| 狂野欧美激情性bbbbbb| 日韩国内少妇激情av| 七月丁香在线播放| 少妇人妻精品综合一区二区| 亚洲久久久国产精品| 18禁裸乳无遮挡动漫免费视频| 97精品久久久久久久久久精品| 久久鲁丝午夜福利片| 亚洲国产欧美在线一区| 最近2019中文字幕mv第一页| 国产黄片美女视频| 黑人猛操日本美女一级片| 中文字幕制服av| 高清黄色对白视频在线免费看 | 在线观看国产h片| 新久久久久国产一级毛片| av在线观看视频网站免费| 好男人视频免费观看在线| 秋霞伦理黄片| 草草在线视频免费看| 老司机影院成人| 国产在线视频一区二区| 精品一区二区三卡| 亚洲欧美精品自产自拍| 日韩av免费高清视频| 国产欧美日韩一区二区三区在线 | 校园人妻丝袜中文字幕| 亚洲精品乱码久久久v下载方式| 美女脱内裤让男人舔精品视频| 男女边吃奶边做爰视频| 日本黄色片子视频| 网址你懂的国产日韩在线| 少妇人妻久久综合中文| 久久久色成人| 亚洲色图av天堂| 中国国产av一级| 亚洲人成网站高清观看| 国产一区二区三区av在线| 国内精品宾馆在线| 99re6热这里在线精品视频| 日韩中文字幕视频在线看片 | 日韩欧美 国产精品| 免费观看无遮挡的男女| 男人添女人高潮全过程视频| 最黄视频免费看| 久久国内精品自在自线图片| 亚洲精华国产精华液的使用体验| 日韩精品有码人妻一区| 在线免费观看不下载黄p国产| 伦理电影大哥的女人| 国产精品免费大片| 久久精品国产鲁丝片午夜精品| 国产精品无大码| 国产亚洲精品久久久com| 国产在线视频一区二区| 国产亚洲91精品色在线| 超碰97精品在线观看| 身体一侧抽搐| 天美传媒精品一区二区| 久久国产精品男人的天堂亚洲 | 日韩伦理黄色片| 99热6这里只有精品| 欧美区成人在线视频| 男女国产视频网站| 亚洲av在线观看美女高潮| av在线观看视频网站免费| 日本一二三区视频观看| 国产毛片在线视频| 亚洲美女视频黄频| 久久人人爽人人片av| 美女脱内裤让男人舔精品视频| 欧美日韩视频高清一区二区三区二| 国产免费又黄又爽又色| 少妇被粗大猛烈的视频| 午夜激情福利司机影院| 寂寞人妻少妇视频99o| 国产成人a∨麻豆精品| 欧美变态另类bdsm刘玥| 亚洲综合色惰| 青春草国产在线视频| .国产精品久久| 日韩成人av中文字幕在线观看| 午夜免费观看性视频| 色婷婷av一区二区三区视频| 妹子高潮喷水视频| 男女啪啪激烈高潮av片| 一本久久精品| 99re6热这里在线精品视频| 精品久久国产蜜桃| 久久精品国产a三级三级三级| 欧美日韩在线观看h| 日本欧美国产在线视频| 久久女婷五月综合色啪小说| 激情 狠狠 欧美| 国产亚洲av片在线观看秒播厂| 午夜福利在线在线| av不卡在线播放| 女性生殖器流出的白浆| 亚洲国产精品一区三区| 精品一品国产午夜福利视频| 国产精品国产三级国产av玫瑰| 久热久热在线精品观看| 欧美精品人与动牲交sv欧美| 欧美xxxx黑人xx丫x性爽| 99re6热这里在线精品视频| videos熟女内射| 下体分泌物呈黄色| 亚洲色图av天堂| 18禁在线播放成人免费| 亚洲av成人精品一区久久| 色吧在线观看| 国产 一区 欧美 日韩| 精华霜和精华液先用哪个| 国产精品一二三区在线看| 国产黄色免费在线视频| 国产黄色视频一区二区在线观看| 国产免费一级a男人的天堂| 国产精品嫩草影院av在线观看| 色网站视频免费| 日本av手机在线免费观看| kizo精华| 亚洲国产av新网站| 一级av片app| 国产极品天堂在线| 五月天丁香电影| 国产熟女欧美一区二区| 天天躁日日操中文字幕| 不卡视频在线观看欧美| 激情五月婷婷亚洲| 精品亚洲成国产av| 草草在线视频免费看| 国产精品伦人一区二区| 日韩人妻高清精品专区| 久久国内精品自在自线图片| h日本视频在线播放| 91在线精品国自产拍蜜月| 日韩成人伦理影院| 国产精品欧美亚洲77777| 国产成人a区在线观看| 亚洲国产av新网站| 亚洲精品国产av蜜桃| 蜜臀久久99精品久久宅男| 大香蕉97超碰在线| 波野结衣二区三区在线| 天美传媒精品一区二区| 观看美女的网站| 美女中出高潮动态图| 日韩大片免费观看网站| 久久精品国产鲁丝片午夜精品| 国产成人freesex在线| av一本久久久久| 一级爰片在线观看| 十分钟在线观看高清视频www | 深夜a级毛片| 久久精品夜色国产| 大话2 男鬼变身卡| 2021少妇久久久久久久久久久| 日韩,欧美,国产一区二区三区| 国产国拍精品亚洲av在线观看| 免费看日本二区| 天堂中文最新版在线下载| 18禁裸乳无遮挡动漫免费视频| 亚洲精品国产色婷婷电影| a 毛片基地| 下体分泌物呈黄色| 亚洲成人一二三区av| 女的被弄到高潮叫床怎么办| 精品人妻一区二区三区麻豆| 一区在线观看完整版| 香蕉精品网在线| 欧美日韩精品成人综合77777| 欧美日韩国产mv在线观看视频 | 五月玫瑰六月丁香| av线在线观看网站| 国产精品福利在线免费观看| 国产精品国产三级国产专区5o| 亚洲av国产av综合av卡| 精品人妻熟女av久视频| 国产高潮美女av| 亚洲va在线va天堂va国产| 黑丝袜美女国产一区| 一区二区三区乱码不卡18| 免费av不卡在线播放| 美女内射精品一级片tv| 美女脱内裤让男人舔精品视频| videossex国产| 高清午夜精品一区二区三区| 最后的刺客免费高清国语| 人人妻人人澡人人爽人人夜夜| 99久久综合免费| av福利片在线观看| 久久韩国三级中文字幕| 免费大片黄手机在线观看| 女人十人毛片免费观看3o分钟| 岛国毛片在线播放| 欧美成人a在线观看| 免费av不卡在线播放| 精品久久久精品久久久| 欧美日韩一区二区视频在线观看视频在线| 亚洲欧美清纯卡通| 身体一侧抽搐| 欧美精品亚洲一区二区| 91午夜精品亚洲一区二区三区| 少妇熟女欧美另类| 成人毛片a级毛片在线播放| 一本久久精品| 制服丝袜香蕉在线| 丝袜脚勾引网站| 日日啪夜夜撸| 成人二区视频| 婷婷色综合大香蕉| 久久久久久人妻| 色婷婷av一区二区三区视频| 新久久久久国产一级毛片| 亚洲国产欧美在线一区| 久久国产精品男人的天堂亚洲 | 黄色视频在线播放观看不卡| 大话2 男鬼变身卡| 简卡轻食公司| 日产精品乱码卡一卡2卡三| 韩国高清视频一区二区三区| 国产午夜精品久久久久久一区二区三区| 欧美丝袜亚洲另类| 亚洲精品久久久久久婷婷小说| 超碰97精品在线观看| 少妇人妻一区二区三区视频| 欧美成人午夜免费资源| 99久久精品热视频| videos熟女内射| 国产成人精品婷婷| 午夜视频国产福利| 欧美国产精品一级二级三级 | 女的被弄到高潮叫床怎么办| 精品国产三级普通话版| 老熟女久久久| 国产欧美另类精品又又久久亚洲欧美| 中文字幕久久专区| 日本猛色少妇xxxxx猛交久久| 亚洲人成网站在线观看播放| 久久ye,这里只有精品| 岛国毛片在线播放| 亚洲va在线va天堂va国产| 91久久精品电影网| 韩国av在线不卡| 亚洲成人中文字幕在线播放| 国产在线一区二区三区精| 日本vs欧美在线观看视频 | 国产欧美日韩精品一区二区| 中国国产av一级| 国产精品.久久久| 国语对白做爰xxxⅹ性视频网站| av福利片在线观看| 日韩免费高清中文字幕av| 午夜福利网站1000一区二区三区| 国产精品99久久99久久久不卡 | 美女xxoo啪啪120秒动态图| 国产欧美亚洲国产| 欧美成人a在线观看| 天堂中文最新版在线下载| 久久精品夜色国产|