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

    Dynamic development model for long gap discharge streamer-leader system based on fractal theory

    2022-08-29 00:42:50JianghaiGENG耿江海QuanshengWANG王權(quán)圣XiaominLIU劉曉敏PingWANG王平FangchengLYU律方成andYujianDING丁玉劍
    Plasma Science and Technology 2022年9期
    關(guān)鍵詞:方成王權(quán)王平

    Jianghai GENG (耿江海), Quansheng WANG (王權(quán)圣) ,Xiaomin LIU (劉曉敏), Ping WANG (王平), Fangcheng LYU (律方成) and Yujian DING (丁玉劍)

    1 Hebei Provincial Key Laboratory of Power Transmission Equipment Security Defense, North China Electric Power University, Baoding 071003, People’s Republic of China

    2 China Electric Power Research Institute, Beijing 100192, People’s Republic of China

    Abstract Establishing a long air gap discharge model considering the streamer-leader transition and randomness of the discharge path is of great significance to improve the accuracy of discharge characteristic prediction and optimize external insulation design. Based on fractal theory and thermal ionization theory of streamer-leader transition, this work establishes a dynamic development model for the long air gap discharge streamer-leader system, which includes streamer inception, streamer development, leader inception, development of streamer-leader system and final jump. The positive discharge process of a 3 m rod plate is simulated to obtain the fractal distribution of the discharge path and the law of leader development for comparison with the discharge test results.The results show that the simulation model is similar to test results in the development characteristics of leader path distribution,each stage time and leader velocity.Finally,a simulation calculation of a 50%breakdown voltage of the rod-plate gap and ball-plate gap is carried out, with results fairly consistent with test data, proving the effectiveness and practicality of the model.

    Keywords: long air gap, fractal theory, thermal ionization, streamer-leader system, simulation model

    1. Introduction

    As the core part of the converter station of a DC transmission project,the valve hall is a closed space with multiple voltage levels and complex air gap types. The selected gap distance determines the operational reliability of related equipment,but many long gap discharge characteristic tests with huge costs are required. Therefore, predicting discharge characteristics through a long air gap discharge simulation model is helpful in order to save test costs, optimize insulation design and improve the economic efficiency and safety of a DC transmission project [1-3].

    Aiming at the physical mechanism of each long gap discharge stage, domestic and foreign scholars put forward simplified models of each stage based on a large number of test results. Bondiout and Gallimberti proposed a positive rod-plate gap discharge model to describe each stage of long gap discharge [4] and calculated the streamer current using the equivalent electron avalanche model, Goeliany et al calculated the space charge quantity using the distortion of space potential and simulated the gradual development of a streamer-leader system [5]. Subsequently, Becerra and Cooray optimized the charge distribution in the streamer zone and simulated the space charge using several ring charges with different radiuses [6]. In [7], the critical volume method is adopted to consider the dispersion of streamer inception and a space charge model considering the variable form of the streamer zone is established based on the charge simulation method, and a separately developed discharge path is simulated. The above models have realized simulation of the whole discharge process by calculating key physical parameters in each discharge stage but ignored the influence of discharge channel branch distribution and space charge distribution.However,the discharge channel branch distribution,especially the branch structure of the streamer zone, will affect distribution of the electric field in space [8], affecting the development of the leader channel and resulting in errors between simulation results and test values.

    In order to simulate the branching phenomenon in the streamer discharge process, domestic and foreign scholars have established various streamer discharge models based on fractal theory. Among them, the most classic fractal models are the DBM model [9] and WZ model [10], considering the influence of electric field distribution and random factors on streamer development. Reference [11] applied fractal theory to the simulation of a long gap discharge path and assumed that the streamer and leader are conductive channels with different conductivities. Reference [12] proposed a streamer growth probability model combined with the streamer discharge theory and fractal theory to simulate the space charge accumulation process with time in the streamer zone. Reference [13] presented an advanced 3D numerical model of lightning development using fractal theory, which has features such as bidirectional propagation and simultaneous growth of multiple branches. In a subsequent work, [14]studied the progression of negative lightning stepped leader including stepping and branching.This kind of model focuses on branch development characteristics of the simulation path and distribution of the predicted discharge path. It does not include the physical process of development of the streamerleader system.Therefore,it is difficult to predict the discharge characteristics.

    References [15, 16] proposed a transport model of discharge in a thunderstorm cloud, in which the probability of the emergence of streamers is considered in a new method,and the leader channel is formed by combining the currents of tens of thousands of streamers. This model combines the advantages of the two models mentioned above, and its morphology and electrical parameters of the calculated discharge tree agree with the existing lightning-discharge data.However,the simulated lightning discharge is different to the typical air gap discharge in transmission engineering in terms of gap size and voltage type.

    Therefore, in view of the shortcomings of the models above, this work establishes a dynamic fractal development model of long gap discharge, which includes streamer inception, streamer development, leader inception, development of a streamer-leader system and final jump based on the basic principles of fractal theory and thermal ionization theory of streamer-leader transition. This aims at the physical mechanism and random development characteristics of discharge paths in the long gap discharge process.First,a 3 m rod-plate gap discharge process was simulated and a long gap discharge test was carried out. Then, the simulated discharge path distribution and characteristic parameters of each stage were verified using the test results. Finally, the U50%breakdown voltages of the rod-plate gap and ball-plate gap were calculated,respectively,based on the model in this work,and the effectiveness and practicality of the model were verified.

    2. Dynamic fractal development model for long gap discharge

    Long gap discharge includes ionization and energy conversion processes with randomness. Based on the above characteristics, this work establishes a model for the long gap discharge development process, which considers the dispersion of streamer inception time delay according to critical volume theory,simulates the fractal development of streamers and leader development path based on fractal theory, and determines streamer-leader transition and development based on thermal ionization theory, to realize a whole process simulation of dynamic fractal development for long gap discharge.

    For the rod-plate gap structure shown in figure 1, the air gap is divided into M×N square grids with side length h,the rod electrode length is set as L and the gap distance is set as H in the fractal development model. The spatial potential φ of each grid point is obtained using the finite difference method:the potential of a grid point outside the discharge channel conforms to the results obtained from the Laplace equation;the potential inside the channel conforms to the results obtained from the Poisson equation:

    where ρ is the charge density inside the channel and ε0is the vacuum dielectric constant. The boundary conditions are set as follows: the potentials of the rod electrode and plate electrode are the first boundary conditions, and the electrode potential is equal to the applied voltage U0(t);the boundary in contact with the air in the solution area is the second boundary condition.

    2.1. Streamer inception

    Under the action of positive switching impulse voltage, the inception time delay of streamer t includes the boosting time delay t0required for boosting the critical electric field Ecof streamer inception and the statistical time delay tsrequired for the emergence of effective free electrons, among which the statistical time delay is one of the main factors leading to the dispersion of discharge inception.

    Figure 1. Grid division diagram of the rod-plate gap.

    According to the classic streamer theory, free electrons collide with air molecules in the process of moving towards the electrode under the action of electric field force in a strong electric field area where the ionization coefficient of electron collisions near the electrode head is greater than the attachment coefficient, resulting in an initial electron avalanche. When the positive charge generated by the initial electron avalanche is enough to support the formation of a secondary electron avalanche and realize the self-sustaining development of the discharge,the streamer incepts;that is,the boosting time delay t0.

    where Ncritis the critical positive charge quantity, that is,Ncrit=0.55×108; α is the ionization coefficient of electron collision and β is the electron attachment coefficient, both of which are functions of electric field intensity. R0is the boundary of the collision ionization zone, where α(Ec)=β(Ec), and Ecis the critical field intensity of the collision ionization zone.

    According to the critical volume theory of streamer inception,when the electric field intensity of the electrode head reaches E0,effective free electrons will be generated in a critical volume before streamer inception. At ts, the occurrence probability of effective electrons is proportional to the critical volume and voltage change rate, the statistical streamer inception time delay tscomplies with Rayleigh distribution,and its probability density function is as follows [17]:

    where k is the function of voltage change ratio and critical volume.Since the radius of the electrode will affect the critical volume,k finally represents the influence of the voltage change rate and electrode shape on the occurrence frequency of effective electrons, and its empirical formula can be obtained by fitting the experimental data.

    After the boosting time delay t0is obtained using formula(2), the statistical time delay tsis calculated according to the probability distribution of streamer inception, and the streamer inception time delay t is obtained.

    2.2. Streamer development

    After the inception of initial streamers, a point in the critical volume of the electrode rod head is randomly selected as the development point of the initial streamers.

    Figure 2.Schematic diagram of the streamer growth probability model.

    When calculating the development in the step m+1 of the streamers, the grid points surrounding them in which the electric field intensity is greater than Ecare selected as the possible development points in step m+1 for all spatial points that have become streamers, as shown in figure 2.

    Assuming that the time tirequired from the emergence of secondary electrons from a point to be developed (i, j-1) to the formation of a streamer channel with the developed point(i, j) satisfies the Weibull distribution, the probability density of the time required for each development step of the streamer p(t) is,

    where r(E) is a streamer growth probability function. Then,the probability of the next formation time tiof the streamer p(E) is,

    Assuming p(E)=ξ, ξ is a random number from 0 to 1 that complies with uniform distribution, tican be expressed as,

    According to the Biller model [18],

    where τ is the formation time constant and τ=2 μs is taken in this simulation, E is the electric field intensity at the point to be developed, η is the development probability index representing the relationship between discharge development probability and electric field intensity (its value is related to voltage type and electrode shape) and η=2 is taken in this simulation.

    When there are n possible points to be developed in the streamer zone,the formation time ti(0<i≤n)of the points to be developed is calculated in turn, and the point to be developed with the shortest formation time is set as the streamer development point in step m+1 and the shortest formation time is the time step tm+1of this step.

    The long gap discharge test shows that the average electric field intensity Esin the streamer zone is about 400-500 kV m-1,so the electric potential of each point in the streamer zone can be obtained using formula (8), and the calculated electric potential is set as the first boundary condition:

    where 1≤i≤m+1, d is the streamer length developed in step i+1 and φ1is the applied voltage U0(t) to the electrode rod.

    The Poisson equation is solved by the finite difference method to obtain the spatial potential distribution, and the charge quantity and current of the streamer channel are obtained by combining the current continuity equation and Ohm’s law. After discretization, the relationship between potential and charge is obtained as follows [12]:

    where tm+1is the development time of the streamer in step m+1, σ is the conductivity of the streamer channel and σ=0.025 S m-1is taken in this simulation.

    The total development time tstreamerin the streamer zone is,

    In the streamer development process, when there are no points to be developed around all developed points, it is considered that the streamer stops developing, and if the electric field around the developed points still does not meet the requirements of critical field intensity Ecwith the rise in voltage, it is considered that the discharge is over.

    2.3. Leader inception

    According to the thermal ionization theory for streamer-leader transition, when the streamer develops rapidly, a quantity of free electrons generated by the streamer branches enters the electrode along the common origin of the streamer branches.In the streamer root transformation zone, part of the energy generated by elastic and inelastic collisions between electrons and gas molecules is directly used for temperature rise,and another part of the energy stored in the form of electronic excitation energy, rotational kinetic energy, translational kinetic energy and vibrational kinetic energy is used to transfer the gas molecules, both of which jointly contribute to the temperature rise of gas molecules [19]. When the temperature of the streamer root reaches the critical temperature Tc=1500 K,negative ions decompose under the action of the thermal ionization mechanism,forming a large number of electrons and neutral molecules,making the conductivity in the streamer root transformation zone increase sharply [20]. Therefore, the streamer root is transformed into the initial leader.

    In the simulation, the vibrational kinetic energy Wv,translational kinetic energy Whand translational kinetic energy temperature Th[4] of the streamer root are calculated according to Gallimberti’s energy balance equation after each step of streamer development. When the temperature Th>1500 K,it is considered that the leader is generated and the streamer branches disappear.

    where k is the Boltzmann constant,nhis the molecular density of gas, fe, ft, frand fvare the distribution coefficients of electronic excitation energy,translational kinetic energy,rotational kinetic energy, and vibrational kinetic energy, respectively. In the simulation, fe+ft+fr=0.03 and fv=0.97, E is the field intensity at the streamer root and I is the current intensity flowing through the streamer root, τvtis the conversion time from vibrational kinetic energy to translational kinetic energy(the simulation uses τvt=5 μs), r is the initial radius of the streamer branch root(300 μm is taken in the simulation),εvis the energy required for vibration excitation of nitrogen molecules (0.28 eV is taken in the simulation) [19].

    Due to the inhibitive effect of distorted electric fields caused by the residual positive charge in space, the initial pilot cannot be developed immediately after being generated.At this time,the field intensity of the leader head is the sum of the background field intensity and the distorted field intensity caused by the space charge. In order to consider recovery of the electric field of the leader head, the conditions of initial leader inception are determined using formula (2). When the electric field of the leader head is recovered enough to meet the requirements of the field strength of streamer inception,the leader can be developed further.

    2.4. Streamer-leader development system

    After leader inception,the leader head generates a new streamer zone, guided by the head streamer branch; the streamer-leader system develops continuously. In the development process,electrons generated in the streamer zone transfer energy to neutral molecules in the transformation zone. When the temperature of the gas molecules reaches critical temperature Tc,it is considered that the main streamer channel is transformed into a new leader.The field intensity in front of the channel increases,which further promotes development of the streamer zone.Therefore, the streamer-leader system can develop on its own.

    In the simulation, the development process of the leader head streamer zone is consistent with section 2.2. The streamer root develops from the head of the new leader section.It is the first condition for temperature to reach Tcin order to achieve leader development.

    When selecting the leader channel, this model uses the leader head as the starting point,then selects the path with high charge density in the streamer zone as the leader development channel, and notes that the length of the ith leader microelement segment converted from the streamer in the path is li,and the total length of each newly developed leader segment is L=∑li. When the amount of charge qLrequired for the development of the leader channel per unit length is shown in a formula(14),it meets the second condition ΔQ≥LqL,that is,the amount of charge generated in the streamer zone is sufficient to support generation of this section of the leader into leader development. If the condition is not met, the streamer region continues to develop until a sufficient amount of charge is generated to support leader generation.where Kgis the leader shape coefficient,which is taken as 0.38 in simulation,tstreameris obtained using formula(8),andΔThis the difference between thermal ionization temperature and room temperature.

    The electron charge entering the leader head is the total amount of ionized electrons in the streamer zone, and the current Ilinjected into the leader head is expressed as,

    In the streamer-leader system development process, the leader channel is regarded as a conductor with a voltage drop,and the streamer zone is a weak plasma channel with a multibranch structure. In the simulation, the leader channel radius is 1.0 mm and the average field intensity Elin the leader channel is 0.5 kV cm-1.

    The leader head potential can be expressed as,

    The calculation steps of potential in the streamer zone are consistent with section 2.2.

    Considering the influence of current,the velocity v of the leader is expressed as [21],

    Then, the growth time tstepof each step of the streamerleader system can be calculated:

    where tleaderis the leader development time.

    2.5. Final jump

    During the development of the streamer-leader system, a large number of free charges accumulate in the air gap, and the streamer area gradually increases until the streamer extends to the plate electrode at a certain moment,the final jump of the leader occurs,then the gap breaks down, and the simulation process stops.

    Figure 3.Flow chart of leader initiation.

    After the streamer develops to the plate electrode, the leader channel and streamer are transformed into an electric arc channel with higher conductivity and higher temperature within a very short time. Therefore, only the development time of the streamer region tstreameris calculated.

    2.6. Simulation procedure

    The leader inception simulation procedure is shown in figure 3.First, the simulation applies operational impulse voltage U0(t) to the electrode and solves the potential and electric field distribution.Then, the streamer inception conditions are determined through formula(2)and the streamer inception time delay is calculated in combination with formula (3). According to the fractal development model of a streamer, the voltage boundary conditions are updated, the simulation streamer develops step by step, and the accumulation amount and development time of each step are recorded. After streamer development stops, the translational kinetic energy temperature of the streamer root is solved according to the energy balance equation of streamer root. If the temperature reaches the critical temperature of 1500 K, the streamer will be transformed into the initial leader.As the voltage continues to rise,the initial leader head electric field is recalculated, and the leader starts to develop when the electric field restores to the critical field intensity E0at the start of the streamer.

    Figure 4. Flow chart of streamer-leader development.

    After leader inception, the development of the streamerleader system is simulated in each step, including the fractal development of the streamer zone and leader transformation.Then, the physical parameters are calculated, such as the potential of the leader head and development time, as shown in figure 4. When the leader head streamer zone develops to the plate electrode, the leader has a final jump, followed by gap breakdown. During the simulation, if the streamer zone stops developing, it is considered to be gap tolerance.

    3. Model simulation results and verification

    According to the above simulation flow,the discharge process of a 3 m rod-plate gap was simulated under operational impulse voltage.The rod electrode length L was set to be 1 m,and the curvature radius of the head was set to be 2 cm while the gap distance H was 3 m and a standard operating shock wave with the peak value of 1000 kV was applied to the rod electrode. Typical simulation results are shown in figure 5.

    In this simulation, the inception time of the initial streamer is t1=20.32 μs, in which the statistical time delay ts=5.43 μs. At t2=26.57 μs, the initial streamer stops developing and turns into a leader. However, due to the shielding of space charge, the leader develops again at t3=33.64 μs and moves to the plate electrode. The leader head streamer reaches the plate electrode at t=114.43 μs,and then gap breakdown occurs.

    3.1. Analysis of discharge characteristics

    To verify the effectiveness of the simulation model,we built a rod-plate long air gap discharge observation system,as shown in figure 6. In the system, a high-speed camera is used to photograph the optical shape of the positive leader discharge channel, and a photomultiplier tube (PMT) is used to synchronously measure the instantaneous optical power in the discharge process, and each discharge stage is determined through the optical image and amplitude of optical power. In the test, a 6000 kV/450 kJ impulse voltage generator is used to apply the operational impulse voltage with the waveform of 250/2500 μs to the rod electrode, and the measurement accuracy of the capacitor voltage divider is 0.5%. The rod electrode is a 6 m long metal rod with a curvature radius of 2 cm. The ground electrode is a galvanized iron plate 11 m×11 m. In the test, three model H10721-210 PMTs are used,with a sampling frequency of 159 MHz(response speed of 6.29 ns)The PMTs are respectively placed at points 1.3 m from the gap axis, 1, 1.6 and 2.6 m from the ground. A Photron Fastcam SA5 high-speed camera is used, and the photographing speed is 7.5×104fps (photographing time interval is 13.33 μs). In order to achieve synchronous triggering of observation equipment, a Teledyne Lecroy 3024 oscilloscope is used to collect the voltage signal from the capacitor voltage divider as the trigger source. When the amplitude exceeds the trigger level, the oscilloscope synchronously collects the voltage signal, and the instantaneous optical power signal triggers the high-speed camera to photograph, and the delay time of each item of equipment meets the synchronous triggering condition [22].

    U50%is applied to the 3 m rod-plate gap for 40 repeated tests, and a typical test result is shown in figure 7.

    According to figure 7(a), the simulation results of leader discharge have high similarity with the optical images obtained from the test, which shows that the model in this paper is effective for predicting the discharge path.

    By comparing the high-speed image with the instantaneous optical power, each stage of the test discharge process is divided. PMT is more sensitive to the optical signal of the discharge process and higher in time resolution. Therefore,each stage of discharge is divided by the change in optical power in different stages.

    It can be seen from figure 7 that the initial streamer incepts at t1=11.38 μs, the optical power collected by PMT1 jumps with an amplitude of 200.8 pW,then the optical power gradually decreases and enters the dark period at t2=14.01 μs. At t3=19.47 μs, the optical power jumps again, indicating generation of the leader, and the high-speed camera captures the leader image at 28.04 μs. With development of the leader head, the optical power collected by PMT1 increases rapidly,exceeding the measured amplitude at t4=32.64 μs.It is considered that the leader has reached the position of PMT1.The development position and stage of the leader head can be determined by three saturation moments of PMT optical power. In the leader development stage, the leader develops 1 m downwards through T1=65.61 μs, and the average axial velocity is 1.52 cm μs-1. At the final jump stage, the leader reaches PMT3 through T2=65.61 μs, and the average velocity is 25.1 cm μs-1. Based on the optical observation results of the high-speed camera, the comparison between the position of the leader head in the 3 m gap and the simulation results is shown in figure 8.According to the above analysis, the simulation and test have similar results in leader path distribution, each stage time and leader development.

    Figure 5. Typical process simulation chart of leader discharge development in 3 m rod-plate gap.

    Figure 6. Layout of the long gap discharge observation system.

    Table 1.Comparison of the simulated results and experimental values for 3 m gap discharge.

    Figure 7.(a) Optical image of the discharge process and (b) change in optical power during discharge.

    Under the same simulation parameters, the leader development simulation is carried out many times and compared with the 40 test results. The comparison of physical macroscopic characteristics in the discharge process of a 3 m rod-plate gap is shown in table 1.

    From the table, it can be seen that the rod-plate gap has strong randomness and dispersion in discharge,but on the whole,the simulation results have good agreement with test results in leader development characteristics, such as breakdown time,leader velocity, charge accumulation, discharge current, etc,and the model has high rationality, in which the statistical time delay is 1.8-24.6 μs,duration of the dark period is 4.2-17.8 μs,leader velocity is 1.4-2.0 cm μs-1, final jump velocity is 14.5-32 cm μs-1, charge density of the leader channel is 42-57 μC m-1and leader discharge current is 0.5-1.4 A.

    Figure 8. Streamer-leader development position.

    Figure 9.Comparison of breakdown voltage.

    3.2. Breakdown voltage prediction

    In order to further verify the effectiveness of the simulation model and achieve prediction of breakdown voltage of the classic external insulation gap in power transmission and transformation projects,the U50%breakdown voltage of rod-plate gap and large-size ball-plate gap is simulated and calculated in this work, and compared with test results and the empirical formula[24] based on the gap coefficient method proposed by Paris. In the simulation, the gap distances with different lengths are set,250/2500 μs standard operation impulse voltage with positive polarity is applied to the rod electrode with a radius of 2 cm and the ball electrode with a radius of 1 m. The different gaps are simulated 30 times using a boosting method, respectively. The results are shown in figure 9.

    It can be seen from figure 9 that for the simulation of the rod-plate gap,the simulation value is similar to the test value and empirical formula. With the gap length increasing to 10 m,the maximum error reaches 5.5%;for the ball-plate gap,and the error between the simulation value and the test value is 4.6%-7.1%. For simulation of long gaps in large electrodes, the simulation value is generally higher than the test value.This is mainly because in practical simulation,in order to save calculation time and running memory,the subdivision size of the mesh is enlarged, which leads to a decrease in accuracy. However, the error is still within an acceptable range.This further illustrates the rationality and practicality of this model for the simulation of long gap discharge since the model not only simulates the physical process of discharge,but also predicts the breakdown voltage of the insulation gap.

    4. Conclusion

    Aiming at two characteristics, namely randomness and leader transformation of the long gap discharge path, this work establishes a dynamic development model for a long gap discharge streamer-leader system based on fractal theory and thermal ionization theory to provide a new path to study the randomness of long gap discharge paths and the law of discharge development.

    (1) The streamer-leader system in the model illustrates the development of a streamer zone and generation of space charge in the form of fractal development, which facilitates analysis of the guiding mechanism of streamer branch to leader path in the discharge process.(2) Through the model established in this work, the leader discharge process of a 3 m rod-plate gap is simulated to obtain the dark period time of 4.2-17.8 μs, breakdown time of 96.3-115.6 μs,leader velocity of 1.4-2.0 cm μs-1and final jump velocity of 14.5-32 cm μs-1. In addition,the rationality of the simulation results is verified through test results.(3) According to the model in the paper, the U50%breakdown voltage of the rod-plate gap and ball-plate gap is simulated and calculated, and the calculation results are highly consistent with test results and empirical formulas.This shows the practicality of the model in this paper in predicting the breakdown voltage of typical gaps in power transmission and transformation projects.

    Acknowledgments

    This work is supported by the Natural Science Foundation of Hebei Province (No. E2021502025).

    ORCID iDs

    Quansheng WANG (王權(quán)圣) https://orcid.org/0000-0001-9321-3758

    猜你喜歡
    方成王權(quán)王平
    篆刻作品欣賞
    法國舊制度時期的王家騎士團(tuán)與王權(quán)
    我眼中的太陽
    Transmission-type reconfigurable metasurface for linear-to-circular and linear-to-linear polarization conversions
    篆刻作品欣賞
    吐蕃王權(quán)研究海外學(xué)術(shù)史鉤沉①
    淺析赫梯國王維護(hù)王權(quán)統(tǒng)治的監(jiān)管措施
    厚德載物 方成大道
    如果歷史是一群喵
    王萬奇:吃苦耐勞方成大事
    金色年華(2016年7期)2016-02-28 01:39:29
    а√天堂www在线а√下载| 亚洲男人的天堂狠狠| 国产激情欧美一区二区| 日韩中文字幕欧美一区二区| 欧美黄色片欧美黄色片| 一本一本综合久久| 国产精品永久免费网站| 又紧又爽又黄一区二区| 97人妻精品一区二区三区麻豆 | 成人三级做爰电影| 日本熟妇午夜| 亚洲av成人不卡在线观看播放网| 国产一区二区激情短视频| 狂野欧美激情性xxxx| 精品久久久久久久久久久久久 | videosex国产| 久久久精品国产亚洲av高清涩受| 国产熟女xx| 成人国产一区最新在线观看| 女性生殖器流出的白浆| 女警被强在线播放| 757午夜福利合集在线观看| 黑丝袜美女国产一区| 亚洲人成电影免费在线| 一区二区三区国产精品乱码| 亚洲av成人不卡在线观看播放网| 老汉色av国产亚洲站长工具| 国产国语露脸激情在线看| 欧美一级a爱片免费观看看 | 国产av不卡久久| 男人的好看免费观看在线视频 | 久久香蕉精品热| 日韩欧美一区视频在线观看| 久久亚洲真实| 中文字幕最新亚洲高清| 在线观看www视频免费| 一进一出好大好爽视频| 给我免费播放毛片高清在线观看| 精品国产乱子伦一区二区三区| 又紧又爽又黄一区二区| 久久久久久久久久黄片| 欧美乱妇无乱码| 一本精品99久久精品77| 欧美av亚洲av综合av国产av| 亚洲av第一区精品v没综合| 欧美中文日本在线观看视频| 曰老女人黄片| 国产精品免费一区二区三区在线| 好男人在线观看高清免费视频 | 久久人妻福利社区极品人妻图片| 国产视频内射| 久久国产精品人妻蜜桃| 老司机深夜福利视频在线观看| 婷婷精品国产亚洲av| 免费在线观看影片大全网站| 国产爱豆传媒在线观看 | 国产在线观看jvid| 19禁男女啪啪无遮挡网站| 这个男人来自地球电影免费观看| 免费在线观看黄色视频的| 日韩欧美在线二视频| 午夜日韩欧美国产| 又黄又爽又免费观看的视频| 久久国产精品影院| 999精品在线视频| 亚洲一区高清亚洲精品| 午夜久久久在线观看| 国产精品野战在线观看| 欧美日韩瑟瑟在线播放| 老汉色∧v一级毛片| 国产成人系列免费观看| 十八禁网站免费在线| 国产精品久久久久久人妻精品电影| 女警被强在线播放| 亚洲va日本ⅴa欧美va伊人久久| 丝袜人妻中文字幕| 欧美国产日韩亚洲一区| 久久精品91蜜桃| 欧美黑人欧美精品刺激| 国产国语露脸激情在线看| 色综合婷婷激情| 国产精品综合久久久久久久免费| 欧美日韩中文字幕国产精品一区二区三区| 久久草成人影院| 日本熟妇午夜| 欧美黄色片欧美黄色片| 精品久久久久久久人妻蜜臀av| 亚洲 欧美 日韩 在线 免费| 美女高潮到喷水免费观看| 国产国语露脸激情在线看| 国产又爽黄色视频| 俺也久久电影网| 人人澡人人妻人| 色哟哟哟哟哟哟| 啦啦啦 在线观看视频| 一进一出抽搐动态| 在线国产一区二区在线| 亚洲国产精品sss在线观看| 日日干狠狠操夜夜爽| 美女扒开内裤让男人捅视频| 国产亚洲精品久久久久久毛片| 欧美激情 高清一区二区三区| 午夜精品在线福利| 国产亚洲精品综合一区在线观看 | 亚洲狠狠婷婷综合久久图片| 女人被狂操c到高潮| 少妇熟女aⅴ在线视频| avwww免费| 又大又爽又粗| 亚洲在线自拍视频| 日本五十路高清| 欧美激情高清一区二区三区| 欧美成人免费av一区二区三区| 婷婷精品国产亚洲av| 岛国视频午夜一区免费看| 久久久国产成人免费| 成人精品一区二区免费| 美女大奶头视频| 国产精品久久久久久人妻精品电影| 国产人伦9x9x在线观看| 999精品在线视频| 亚洲无线在线观看| 国产成人啪精品午夜网站| 国产视频一区二区在线看| 成年人黄色毛片网站| 久久 成人 亚洲| 男女午夜视频在线观看| 中国美女看黄片| 淫妇啪啪啪对白视频| 亚洲专区中文字幕在线| 免费搜索国产男女视频| 国产精品电影一区二区三区| 97人妻精品一区二区三区麻豆 | 亚洲自偷自拍图片 自拍| 日韩成人在线观看一区二区三区| www.www免费av| 怎么达到女性高潮| 午夜免费成人在线视频| 免费看十八禁软件| 精品久久久久久久人妻蜜臀av| 亚洲成人精品中文字幕电影| 十八禁人妻一区二区| 久久青草综合色| 日本五十路高清| cao死你这个sao货| 日日摸夜夜添夜夜添小说| 波多野结衣高清无吗| 狠狠狠狠99中文字幕| 久久久久精品国产欧美久久久| 999精品在线视频| 亚洲狠狠婷婷综合久久图片| 色尼玛亚洲综合影院| 999精品在线视频| 女性生殖器流出的白浆| 日本一区二区免费在线视频| 久久青草综合色| 亚洲国产高清在线一区二区三 | 99久久精品国产亚洲精品| 国产精品99久久99久久久不卡| 男人舔奶头视频| 少妇裸体淫交视频免费看高清 | 两个人看的免费小视频| 国产av在哪里看| 免费搜索国产男女视频| 国产精品电影一区二区三区| 一a级毛片在线观看| 久久久久九九精品影院| 国产亚洲精品一区二区www| 久久99热这里只有精品18| 1024手机看黄色片| 欧美不卡视频在线免费观看 | 1024手机看黄色片| 欧美激情久久久久久爽电影| 国产黄a三级三级三级人| 国产97色在线日韩免费| 亚洲精品国产区一区二| 免费看美女性在线毛片视频| 亚洲欧洲精品一区二区精品久久久| 午夜视频精品福利| 成人午夜高清在线视频 | 亚洲人成77777在线视频| 午夜福利成人在线免费观看| 精品电影一区二区在线| 日本 欧美在线| 国产激情欧美一区二区| 哪里可以看免费的av片| 99国产精品一区二区三区| 好男人电影高清在线观看| 久久精品国产清高在天天线| 精品久久久久久久人妻蜜臀av| 欧美绝顶高潮抽搐喷水| 大型黄色视频在线免费观看| 亚洲精品美女久久久久99蜜臀| 欧美色视频一区免费| 美国免费a级毛片| 婷婷精品国产亚洲av在线| 欧美日韩亚洲综合一区二区三区_| 一进一出抽搐动态| av电影中文网址| 久久精品夜夜夜夜夜久久蜜豆 | 国产主播在线观看一区二区| 久久人妻福利社区极品人妻图片| 两人在一起打扑克的视频| 国产99久久九九免费精品| 亚洲美女黄片视频| 最近在线观看免费完整版| 欧美国产精品va在线观看不卡| 99热只有精品国产| 亚洲第一欧美日韩一区二区三区| 日本免费一区二区三区高清不卡| 亚洲一码二码三码区别大吗| 色综合婷婷激情| 精品国产美女av久久久久小说| 非洲黑人性xxxx精品又粗又长| 激情在线观看视频在线高清| av有码第一页| www.熟女人妻精品国产| 自线自在国产av| 日日干狠狠操夜夜爽| 看免费av毛片| 久久精品国产清高在天天线| 成年女人毛片免费观看观看9| 亚洲精品在线美女| 欧美黑人巨大hd| 日日爽夜夜爽网站| 国内揄拍国产精品人妻在线 | 妹子高潮喷水视频| 亚洲av熟女| 亚洲 欧美一区二区三区| 午夜福利欧美成人| 最新在线观看一区二区三区| 亚洲成人久久性| 人成视频在线观看免费观看| 中文字幕高清在线视频| 欧美黄色片欧美黄色片| 欧洲精品卡2卡3卡4卡5卡区| 天天躁狠狠躁夜夜躁狠狠躁| 久99久视频精品免费| 亚洲精品一卡2卡三卡4卡5卡| 亚洲无线在线观看| 精品国内亚洲2022精品成人| 国产私拍福利视频在线观看| 亚洲精品国产一区二区精华液| www.自偷自拍.com| 好男人在线观看高清免费视频 | 一级作爱视频免费观看| 国产视频一区二区在线看| 免费观看人在逋| 成人一区二区视频在线观看| 草草在线视频免费看| 三级毛片av免费| 女同久久另类99精品国产91| 99久久无色码亚洲精品果冻| 欧美 亚洲 国产 日韩一| 亚洲熟妇中文字幕五十中出| 日韩欧美国产在线观看| 久久国产精品男人的天堂亚洲| av福利片在线| 999久久久精品免费观看国产| 亚洲精品中文字幕一二三四区| 日日摸夜夜添夜夜添小说| 一进一出好大好爽视频| 国产激情久久老熟女| 亚洲av成人不卡在线观看播放网| 欧美日韩一级在线毛片| 欧美日韩黄片免| 亚洲av熟女| 91在线观看av| 老汉色av国产亚洲站长工具| www.熟女人妻精品国产| 在线观看一区二区三区| 成年女人毛片免费观看观看9| ponron亚洲| 十八禁人妻一区二区| 757午夜福利合集在线观看| 91成年电影在线观看| 91老司机精品| 丰满的人妻完整版| 亚洲avbb在线观看| 丁香欧美五月| 夜夜夜夜夜久久久久| 亚洲成人久久性| a级毛片在线看网站| 国产成+人综合+亚洲专区| 国内揄拍国产精品人妻在线 | 搞女人的毛片| 免费高清视频大片| 老司机午夜福利在线观看视频| 成人av一区二区三区在线看| 午夜精品在线福利| 国产一区二区三区视频了| 色精品久久人妻99蜜桃| 日本熟妇午夜| 搡老妇女老女人老熟妇| 成人永久免费在线观看视频| 久久香蕉国产精品| 日日爽夜夜爽网站| 亚洲精品粉嫩美女一区| 啦啦啦观看免费观看视频高清| 少妇被粗大的猛进出69影院| 美女扒开内裤让男人捅视频| av在线播放免费不卡| 国产区一区二久久| 午夜影院日韩av| 亚洲成人久久爱视频| 欧美日韩黄片免| 久久久久国内视频| 精品少妇一区二区三区视频日本电影| 欧美成人一区二区免费高清观看 | 9191精品国产免费久久| 欧美性猛交黑人性爽| 午夜免费鲁丝| 长腿黑丝高跟| 亚洲国产精品合色在线| 男女之事视频高清在线观看| 日韩国内少妇激情av| 国产精品精品国产色婷婷| 日日爽夜夜爽网站| 亚洲人成伊人成综合网2020| www.自偷自拍.com| 免费在线观看亚洲国产| 精品人妻1区二区| 国产国语露脸激情在线看| 亚洲一区二区三区色噜噜| 又大又爽又粗| 午夜两性在线视频| 国产精品久久视频播放| 国产精品亚洲一级av第二区| 日韩大码丰满熟妇| 亚洲熟女毛片儿| 久久精品成人免费网站| 成年版毛片免费区| 久久精品91无色码中文字幕| 国产精品免费视频内射| 久久性视频一级片| 欧美中文综合在线视频| 丁香六月欧美| 欧美黄色片欧美黄色片| 国产亚洲欧美98| 国产亚洲精品综合一区在线观看 | 免费观看人在逋| 高潮久久久久久久久久久不卡| 欧美在线黄色| 亚洲成国产人片在线观看| 悠悠久久av| 日韩大码丰满熟妇| 麻豆成人av在线观看| 欧美性猛交黑人性爽| 久久热在线av| 一区二区日韩欧美中文字幕| 婷婷丁香在线五月| 亚洲美女黄片视频| 亚洲一区高清亚洲精品| 欧美+亚洲+日韩+国产| 99在线视频只有这里精品首页| 一级a爱视频在线免费观看| 欧美激情高清一区二区三区| 看黄色毛片网站| 美女高潮喷水抽搐中文字幕| 黄色丝袜av网址大全| 日韩精品免费视频一区二区三区| 一区二区三区精品91| 欧美性长视频在线观看| 精品免费久久久久久久清纯| 又黄又粗又硬又大视频| 日本免费a在线| 在线国产一区二区在线| 啦啦啦 在线观看视频| 成人18禁高潮啪啪吃奶动态图| 久久婷婷人人爽人人干人人爱| 亚洲第一欧美日韩一区二区三区| 久久国产精品男人的天堂亚洲| 久久久久久亚洲精品国产蜜桃av| 亚洲美女黄片视频| 老司机福利观看| 国产精品久久久av美女十八| 91麻豆av在线| 在线观看一区二区三区| 国产成人啪精品午夜网站| 久久国产精品男人的天堂亚洲| 久久草成人影院| 色在线成人网| 成人三级做爰电影| av在线播放免费不卡| 亚洲一区高清亚洲精品| 午夜福利欧美成人| 午夜福利成人在线免费观看| 亚洲专区中文字幕在线| 久久精品成人免费网站| 久久香蕉精品热| 18禁美女被吸乳视频| 黄片播放在线免费| 久久中文看片网| 黄色丝袜av网址大全| 午夜免费观看网址| 国产在线精品亚洲第一网站| 十八禁网站免费在线| 人人澡人人妻人| 亚洲熟妇中文字幕五十中出| 成人一区二区视频在线观看| 久久久久久久午夜电影| netflix在线观看网站| 淫秽高清视频在线观看| 搡老熟女国产l中国老女人| 嫩草影院精品99| 精品国产国语对白av| 嫩草影院精品99| 成人永久免费在线观看视频| 免费人成视频x8x8入口观看| 黄色成人免费大全| 桃红色精品国产亚洲av| 亚洲精品中文字幕在线视频| 嫁个100分男人电影在线观看| 在线观看免费视频日本深夜| 两性夫妻黄色片| 999久久久精品免费观看国产| 最近最新免费中文字幕在线| 欧美激情久久久久久爽电影| av超薄肉色丝袜交足视频| 国产精品98久久久久久宅男小说| 久久性视频一级片| 91字幕亚洲| 亚洲人成伊人成综合网2020| 美女免费视频网站| 午夜成年电影在线免费观看| 少妇被粗大的猛进出69影院| 男女那种视频在线观看| 久久精品影院6| 99国产精品99久久久久| 99国产综合亚洲精品| 久久精品亚洲精品国产色婷小说| 国产成人精品无人区| 免费无遮挡裸体视频| 最近最新免费中文字幕在线| 999精品在线视频| 欧美日韩一级在线毛片| 国产av在哪里看| 亚洲一区高清亚洲精品| 欧美黑人精品巨大| 久久热在线av| 亚洲av成人av| 人人妻人人看人人澡| 极品教师在线免费播放| 国产亚洲精品综合一区在线观看 | 久久中文字幕一级| 中文字幕精品免费在线观看视频| 一个人观看的视频www高清免费观看 | 精品高清国产在线一区| 成人亚洲精品一区在线观看| 一级毛片高清免费大全| 精品国产乱码久久久久久男人| 国产精品国产高清国产av| 最近最新免费中文字幕在线| 悠悠久久av| 国产日本99.免费观看| 亚洲国产精品999在线| 亚洲精品国产一区二区精华液| 国产一卡二卡三卡精品| 麻豆一二三区av精品| av在线天堂中文字幕| 男女午夜视频在线观看| 国产亚洲欧美98| 亚洲精品久久成人aⅴ小说| 久久久久精品国产欧美久久久| 免费在线观看完整版高清| 免费在线观看日本一区| 亚洲精品美女久久久久99蜜臀| 丁香六月欧美| 亚洲精品中文字幕一二三四区| 欧美日韩精品网址| 欧美乱妇无乱码| 91字幕亚洲| 亚洲第一欧美日韩一区二区三区| 宅男免费午夜| 真人做人爱边吃奶动态| 免费女性裸体啪啪无遮挡网站| 免费无遮挡裸体视频| 欧美乱妇无乱码| 久久欧美精品欧美久久欧美| 亚洲欧美精品综合久久99| 婷婷六月久久综合丁香| 夜夜躁狠狠躁天天躁| 久久久国产精品麻豆| 国产免费男女视频| 久久精品成人免费网站| 久久精品人妻少妇| 精品午夜福利视频在线观看一区| 91麻豆精品激情在线观看国产| 精品国内亚洲2022精品成人| 性色av乱码一区二区三区2| 免费无遮挡裸体视频| 真人一进一出gif抽搐免费| 波多野结衣高清无吗| 91麻豆av在线| 国产精品,欧美在线| 99国产精品99久久久久| 日韩av在线大香蕉| 少妇粗大呻吟视频| 婷婷丁香在线五月| 天堂动漫精品| 久久久久免费精品人妻一区二区 | 成人精品一区二区免费| 免费看美女性在线毛片视频| 一级a爱视频在线免费观看| 在线观看日韩欧美| 搡老岳熟女国产| 欧美久久黑人一区二区| 国产一区二区三区在线臀色熟女| 日韩精品中文字幕看吧| 91国产中文字幕| 日韩精品青青久久久久久| 久久精品国产99精品国产亚洲性色| 国产蜜桃级精品一区二区三区| 中文字幕av电影在线播放| 免费女性裸体啪啪无遮挡网站| 熟女电影av网| 一进一出抽搐动态| 国产又色又爽无遮挡免费看| 脱女人内裤的视频| 搡老妇女老女人老熟妇| 成人特级黄色片久久久久久久| 日本撒尿小便嘘嘘汇集6| 久久婷婷人人爽人人干人人爱| 久久 成人 亚洲| 日韩av在线大香蕉| 精品电影一区二区在线| xxxwww97欧美| 宅男免费午夜| 女人爽到高潮嗷嗷叫在线视频| 这个男人来自地球电影免费观看| 成人国产一区最新在线观看| 一卡2卡三卡四卡精品乱码亚洲| 在线免费观看的www视频| 亚洲国产欧美网| 母亲3免费完整高清在线观看| 亚洲国产欧美日韩在线播放| 精品日产1卡2卡| 久久99热这里只有精品18| 丁香六月欧美| av超薄肉色丝袜交足视频| 最近最新中文字幕大全免费视频| 999久久久精品免费观看国产| 欧美黄色淫秽网站| 欧美黑人欧美精品刺激| av电影中文网址| 又大又爽又粗| 大型黄色视频在线免费观看| 精品久久蜜臀av无| 免费av毛片视频| 热re99久久国产66热| 啪啪无遮挡十八禁网站| 久久久久久九九精品二区国产 | 在线十欧美十亚洲十日本专区| 又大又爽又粗| 午夜影院日韩av| 欧美久久黑人一区二区| 搡老妇女老女人老熟妇| 悠悠久久av| 亚洲avbb在线观看| 亚洲精品国产精品久久久不卡| www日本在线高清视频| 国产av又大| 搡老岳熟女国产| 国产激情欧美一区二区| 最近最新中文字幕大全免费视频| 国产精品综合久久久久久久免费| 一卡2卡三卡四卡精品乱码亚洲| 久久久久久国产a免费观看| 亚洲色图av天堂| 欧美乱妇无乱码| 巨乳人妻的诱惑在线观看| 50天的宝宝边吃奶边哭怎么回事| 精品免费久久久久久久清纯| 日韩欧美 国产精品| av片东京热男人的天堂| 18禁国产床啪视频网站| 亚洲午夜理论影院| 丝袜人妻中文字幕| 2021天堂中文幕一二区在线观 | 俄罗斯特黄特色一大片| 精品国产超薄肉色丝袜足j| 亚洲欧美精品综合久久99| 女生性感内裤真人,穿戴方法视频| 国产精品98久久久久久宅男小说| aaaaa片日本免费| 亚洲男人天堂网一区| 两性午夜刺激爽爽歪歪视频在线观看 | 欧美成人性av电影在线观看| 悠悠久久av| 免费看日本二区| 老汉色av国产亚洲站长工具| 欧美性长视频在线观看| 免费在线观看黄色视频的| 两个人免费观看高清视频| 女人爽到高潮嗷嗷叫在线视频| 在线观看免费日韩欧美大片| 啦啦啦免费观看视频1| 窝窝影院91人妻| av超薄肉色丝袜交足视频| 国产三级在线视频| 我的亚洲天堂| 69av精品久久久久久| 成人国产一区最新在线观看| 最近最新中文字幕大全电影3 | 亚洲av中文字字幕乱码综合 | 欧美又色又爽又黄视频| 精品久久久久久久久久免费视频| 一夜夜www| 一进一出抽搐gif免费好疼| 亚洲国产精品久久男人天堂| 丝袜美腿诱惑在线| 欧美亚洲日本最大视频资源| 日本三级黄在线观看| 欧美激情高清一区二区三区| 老熟妇仑乱视频hdxx|