9512.net
甜梦文库
当前位置:首页 >> >>

Dependence of the liquid-vapor surface tension on the range of interaction a test of the la


Dependence of the liquid-vapor surface tension on the range of interaction: a test of the law of corresponding states
Patrick Gros?ls?

arXiv:0807.0107v1 [cond-mat.stat-mech] 1 Jul 2008

Microgravity Research Center, Chimie Physique E.P. CP 165/62, Universit? e Libre de Bruxelles, Av.F.D.Roosevelt 50, 1050 Brussels, Belgium.? James F. Lutsko Center for Nonlinear Phenomena and Complex Systems CP 231, Universit? e Libre de Bruxelles, Blvd. du Triomphe, 1050 Brussels, Belgium?
(Dated: July 1, 2008)

Abstract
The planar surface tension of coexisting liquid and vapor phases of a ?uid of Lennard-Jones atoms is studied as a function of the range of the potential using both Monte Carlo simulations
? = 2.5 to r ? = 6 and the and Density Functional Theory. The interaction range is varied from rc c

surface tension is determined for temperatures ranging from T ? = 0.7 up to the critical temperature in each case. The results are shown to be consistent with previous studies. The simulation data are well-described by Guggenheim’s law of corresponding states but the agreement of the theoretical results depends on the quality of the bulk equation of state.
PACS numbers: 68.03.Cd,65.20.De

? ?

Electronic address: pgros?@ulb.ac.be Also at Center for Nonlinear Phenomena and Complex Systems CP 231, Universit? e Libre de Bruxelles,

Blvd. du Triomphe, 1050 Brussels, Belgium
?

Electronic address: jlutsko@ulb.ac.be; URL: http://www.lutsko.com

1

I.

INTRODUCTION

One of the most fundamental properties of a ?uid is the surface tension at the liquid-vapor interface. It would seem that such a fundamental property would be an ideal candidate for study via computer simulation. However, the determination of the surface tension from simulation turns out to be frought with di?culties so that even today there is still a substantial amount of e?ort directed towards the development of more reliable algorithms and the re?nement of the reported values even for the paradigmatic case of a simple ?uid modeled with the Lennard-Jones interaction[1]. One of the primary di?culties is that in all simulations, the potential is truncated at some ?nite range and it happens that the surface tension is very sensitive to the value of the cuto?. For that reason, an important part of the development of algorithms has focused on the calculation of the corrections needed to get the in?nite-ranged limit from data obtained using a truncated potential (see, e.g., ref. [1, 2]). This sensitivity is therefore a nuisance when the goal is to get the in?nite range result, but in other ways it can made useful. In particular, one of the important reasons to determine the surface tension from simulation is that it provides a baseline against which theories of inhomogeneous liquids can be tested[3, 4]. For this application, the sensitivity of the surface tension to the range of the potential can be used as a test of the generality of a theory which was probably motivated in the ?rst place by its agreement with some existing simulation data. Furthermore, there has recently been a signi?cant increase in interest in short-ranged potentials in their own right. This is due to the fact that certain complex ?uids, in particular globular proteins, can, in a ?rst approximation, be modeled as a simple ?uid with a very short ranged interaction[5]. It is therefore interesting to study the properties of ?uids with these kinds of interactions and to test that existing theories work in this new domain of interest. For these reasons, we present in this paper a systematic study of the dependence of the surface tension of a Lennard-Jones ?uid as a function of the range of the potential. In this paper, we describe the results of Monte Carlo (MC) simulations of a Lennard-Jones ?uid with the potential truncated at several di?erent points. We have chosen to truncate and shift the Lennard-Jones potential, vLJ (r ), so that the potential used in this work is v (r ; rc ) = vLJ (r ) ? vLJ (rc ) for r < rc and v (r ) = 0 for r ≥ rc . If this shift is not performed, then there is an impulsive contribution to the pressure when atoms move across the r = rc 2

boundary that would have to be taken into account[6]. We do not shift the force, i.e. we
′ ′ do not use v ?(r ; rc ) = vLJ (r ) ? vLJ (rc ) ? (r ? rc )vLJ (rc ) with vLJ (r ) = dvLJ (r )/dr inside the

cuto?, as is usually done in molecular dynamics simulations to avoid impulsive forces: our potential is truncated and shifted but the force is not shifted. This choice was made in order to allow for comparison with previous MC studies. In the simulations a slab of liquid is bounded on both sides by vapor. The surface tension is determined using the method of Bennett[6, 7] as there seems to be some evidence that this method is more robust than other commonly used techniques[8]. It is often the case that the quantity of interest is the surface tension for the in?nite-ranged potential. Since simulations almost always make use of truncated potentials, various techniques have been developed to approximate the so-called long range corrections, i.e. the di?erence between quantities calculated with the truncated potential and the in?nite-ranged quantities[2]. We do not include any such corrections here since our goal is actually to study the truncated potentials. Thus, each value of the cuto? de?nes a di?erent potential with its own coexistence curve and thermodynamics. In Section II, we present the simulation techniques used in our work. Section III contains a discussion of our results including a comparison to previous work. Since one of the motivations for this work is to provide a baseline for testing theories of the liquid state, we illustrate this by comparing our results to Density Functional Theory (DFT) calculations and by testing the law of corresponding states. We give our conclusions in the last Section.

II.

SIMULATION METHODS

Simulations are performed with a standard Metropolis Monte-Carlo algorithm (MC-NVT) for a system of N = 2000 particles of mass m at temperature T in a volume V = Lx Ly Lz where Lx , Ly , and Lz are the dimensions of the rectangular simulation cell. Periodic boundary conditions are used in all directions. Particles interact via the Lennard-Jones potential, vLJ (r ) = 4? σ r
12

?

σ r

6

(1)

which is truncated and shifted so that the potential simulated is ? ? v (r ) ? v (r ) : r < r LJ LJ c c v (r ) = ? 0 : r ≥ rc 3

(2)

where rc is the cuto? radius. Each simulation starts from a rectangular box (Lx = Ly = L, Lz = 4 L) ?lled with a dense disordered arrangement of particles (ρ? ≡ ρσ 3 = 0.8) surrounded along the z-direction by two similar rectangular boxes containing particles in a low density state (ρ? ? 0.01) . The total simulation box has sides of length Lx = Ly = 9.15σ and Lz = 109.63σ . The liquid ?lm located in the middle of the box has a thickness ?z ? 27σ so that the two interfaces do not in?uence each others. The system is ?rst equilibrated during 5 × 105 Monte-Carlo cycles (one cycle = N updates) after which the positions of the particles are saved every 20 cycles during 5 × 105 cycles. This ensemble of 2.5 × 104 con?gurations is used to compute the density pro?le and the surface tension by the Bennett’s method. Although several methods are available for the computation of the surface tension, the Bennett’s approach has been chosen because of its accuracy[8]. In the Bennett’s method the calculation of the surface tension follows from the de?nition γ = ?F ?A (3)
N,V,T

where F is the free energy and A is the area of the liquid-vapor interface. In its implementation the method requires that one performs two simulations: one for system 0 of interface area A0 , and another for system 1 of interface area A1 = A0 + ?A. In this work ?A/A = 5 × 10?4 . The free energy di?erence ?F between the two systems is evaluated by the method of acceptation ratio which starts with the computation of ?E01 = E01 ? E00 which is the di?erence between E00 , the energy of a con?guration of system 0, and E01 , the energy of a new con?guration obtained from the previous one by rescaling the positions of the particles[6, 9] : x′ = x (A1 /A0 ) 2 , y ′ = y (A1 /A0 ) 2 , and z ′ = z (A0 /A1 ). Similarly one computes ?E10 = E10 ? E11 obtained from a con?guration of system 1 following an inverse rescaling of the positions. ?F is obtained by requiring that f (?E01 ? ?F ) =
n0 n1
1 1

f (?E10 + ?F )

(4)

where

n0

(

n1 )

is a sum over the con?gurations of systems 0 (1), and f (x) = (1 +

exp(βx))?1 . Then, taking into account the fact that the system contains two ?at interfaces, the value of the surface tension is given by γ = ?F/(2?A) .

4

1

rc =6

*

γ? 0.5 rc =2.5
*

0 0.6

0.8

T*

1

1.2

FIG. 1: (Color online)The surface tension as a function of temperature for two di?erent cuto?s. The open circles are our data, the ?lled circles are from Duque et al[11], the squares are from Mecke et al[1], the diamonds are from Poto? et al[12] and the triangles are from [10]). Note that the Mecke and Poto? data both include long-ranged corrections. III. A. RESULTS Comparison to previous results

Our results for the surface tension as a function of the cuto? are given in Table 1. Note that all quantities are reported in reduced units so that the reduced temperature is T ? = T /?,
? the reduced cuto?s are rc = rc /σ and the reduced surface tension is γ ? = γσ 2 /?. In Fig. 1 ? we show our results for cuto?s of rc = 2.5 and 6.0 compared to the MC data of Haye and

Bruin[10] for the shorter cuto? and to the MD data of Duque et al[11] (who appear to shift the forces) and Poto? et al[12] and Mecke et al[1]. The latter two are shown even though they include long-ranged corrections. Our data are seen to be very consistent with the MC data obtained without long-ranged corrections and to lie slightly below the corrected data, as expected.

5

B.

Comparison to DFT

In Fig. 2, we compare our results to the predictions of a recently proposed Density Functional Theory model[4]. The DFT requires knowledge of the bulk equation of state and the ?gure shows results using two di?erent inputs: the 33-parameter equation of state of Johnson, Zollweg and Gubbins (JZG)[13] and ?rst order Barker-Henderson thermodynamic perturbation theory[14, 15]. Both versions of the theory are in good qualitative agreement with the data, showing the decrease in surface tension as the range of the potential decreases. It might be thought that use of an empirical equation of state should automatically give superior results to an approximation, like thermodynamic perturbation theory, but this is not necessarily the case since the equation of state is ?tted to data for the in?nite-ranged potential. The ?nite cuto? is accounted for using simple mean-?eld corrections[6, 13] and these become increasing inaccurate as the cuto? becomes shorter and, for ?xed cuto?, as the ?uid density becomes higher. The latter condition means, in the present context, increasing inaccuracy as the temperature decreases. Both of these trends are con?rmed by the ?gure. The decrease in accuracy with decreasing cuto? can be seen in the fact that the critical point (corresponding to the temperature at which the surface tension extrapolates to zero) is less accurately estimated for the smaller cuto?s than for the larger cuto?s. The perturbation theory, on the other hand, takes the cuto? into account more accurately and consistently so that no strong change in accuracy is expected as the cuto? decreases. However, the theory itself is expected to be less accurate for higher densities so again a drop in accuracy with decreasing temperature would be expected and that is indeed seen in the ?gure. Furthermore, perturbation theory is in general going to be inaccurate near the critical point as it does not take into account renormalization e?ects which tend to lower the critical point. These e?ects are less pronounced for shorter-ranged potentials and indeed the perturbation theory seems more consistent with shorter-ranged potential.

C.

Corresponding states

The principle of corresponding states is a generalization of the results of the van der Waals equation of state[16]. The idea is that the properties of simple liquids should be universal functions of the state variables, density and temperature, scaled to the critical

6

1

1

0.5

0.5

0 0.6

0.8

1

0 0.6

0.8

1

FIG. 2: (Color online)Comparison of our simulation data to a DFT model[4]. The panel on the left shows the results of the theory using an empirical equation of state while the results on the right were obtained using thermodynamic perturbation theory.The lines are ordered from the smallest
? = 2.5, 3, 4, 6, ∞. cuto? (lowest lines) to the largest cuto? (highest lines) and were calculated for rc

The data is represented by circles(2.5), squares(3), diamonds(4), ?lled diamonds (4 - larger system) and triangles(6).

point. In this section, we test this hypothesis by applying it to the surface tension. The ?rst step is therefore to determine the critical temperatures and densities of the various truncated potentials. Since the theoretical calculations require as input an equation of state, the critical points are easily determined. To determine them from the simulations, we took ?ve independent averages over 5000 con?gurations and ?tted the density pro?les in each case to a hyperbolic tangent and then from these extract the coexisting vapor and liquid densities at each temperature. The ?ve values obtained at each temperature were averaged and the variance used as an estimate of the errors in the values. The critical temperature was then estimated by using the lowest order renormalization group (RG) result,
1 (ρ 2 l

? ρv ) = A(Tc ? T )0.325 [17, 18]. In some applications[19], higher order terms

are included but we did not feel that the accuracy of our data warranted use of anything

7

TABLE I: Surface tension determined from simulation as a function of temperature for di?erent cuto?s. Temperature 0.70 0.72 0.75 0.80 0.85 0.90 0.95 1.0 1.05 1.1 1.15 1.2
a b

? = 2.5a rc

? = 3a rc

? = 4a rc

? = 4b rc

? = 6b rc

0.584 (27) 0.561 (26) 0.511 (20) 0.421 (19) 0.315 (13) 0.228 (13) 0.181 (11) 0.106 (7)

0.770 (21) 0.726 (25) 0.698 (28) 0.608 (16) 0.480 (12) 0.384 (15) 0.314 (12) 0.234 (8) 0.156 (11) 0.111 (16) 0.074 (10) 0.054 (6)

0.964(46) 0.899(22) 0.825(31) 0.748(25) 0.633(24) 0.542(16) 0.443(22) 0.348(11) 0.258(23) 0.200 (12) 0.108 (14) 0.067 (15)

0.914(30)

1.070(13) 1.034(7) 0.959(8)

0.736 (32)

0.847(13)

0.484(33)

0.638(8)

0.313 (59)

0.438(15)

0.261(13)

0.063 (19)

0.105(5)

using approximately 2000 atoms. using approximately 8000 atoms.

but the lowest order function. The critical density was then estimated using the law of the
1 rectilinear diameter, 2 (ρl + ρv ) = ρc + B (Tc ? T )[16]. There are again higher order corrections

to this formula which can be calculated using RG methods, but for the reasons just given, we have not attempted to include them. The results of these ?ts are illustrated in Fig. 3 and summarized in Table 2. The largest errors in this procedure are in the determination of the critical density. Figure 4 shows the surface tensions, as determined from simulation and theory using the JZG equation of state, scaled to the critical density and temperature as a function of distance from the critical temperature. Despite the wide range of cuto?s and the mixture of data from simulations and theory, it is nevertheless seen that the data do in fact obey the law of corresponding states to a good approximation. However, the same scaling of the theoretical calculations using the equation of state from thermodynamic perturbation theory, shown in
? Figure 5, does not give a single curve. While the data for the shorter cuto?s, rc = 2.5 and

8

0.5

0.4 (ρL-ρv)/2

0.3

0.2

0.1

0.8

1 * T

1.2

FIG. 3: (Color online) The ?t of the di?erence in liquid and vapor densities, as determined from simulation (symbols), to the RG functional form (lines). The data are shown as circles (Rc = 6.0), open squares (Rc = 4.0, larger cell), ?lled squares (Rc = 4.0, smaller cell), diamonds (Rc = 3.0) and triangles (Rc = 2.5).
? rc = 3 appear to coincide, the data for the larger cuto?s does not. This appears to be due,

at least in part, to the fact that the estimate of the critical density as a function of the cuto? calculated using the perturbation theory is not monotonic (see Table II) which is at odds with the quantities as determined from simulation which clearly are monotonic in the cuto?.

IV.

CONCLUSIONS

We have presented our determination of the liquid-vapor surface tension in a LennardJones ?uid as a function of the range of the potential. The data give a systematic picture of the variation of surface tension with the cuto? and are in agreement with previous studies. It is hoped that this can serve as a useful benchmark for the development of theories of inhomogeneous liquids. Indeed, the results were compared here to calculations made using a

9

TABLE II: The critical points for the LJ potential truncated at di?erent values. The theoretical values were determined using the empirical JZG equation of state[13] (JZG) and the ?rst order Barker-Henderson perturbation theory (BH).
? Rc

MC Tc? ρ? c 0.31 (9)a 0.31 (7)a 0.31 (5)a 0.30 (6)b 0.32 (9)b 0.317
c

Theory - JZG Tc? 1.04 1.15 1.25 — 1.29 1.311 ρ? c 0.26 0.28 0.32 — 0.35 0.351 Tc?

Theory - BH ρ? c 0.325 0.342 0.341 — 0.341 0.312

2.5 3.0 4.0 4.0 6.0 ∞
a b

1.10 (1)a 1.18 (1)a 1.26 (1)a 1.25 (2)b 1.30 (2)b 1.31
c

1.18 1.27 1.34 — 1.38 1.40

using approximately 2000 atoms. using approximately 8000 atoms. c From ref. [20].

recently developed DFT and the strengths and weaknesses of the theory are evident: while it gives a good semi-quantitative estimate of the surface tension for all cuto?s, errors on the order of 10% are present indicating that further improvement is possible. We have also tested the law of corresponding states by showing our results from both simulation and theory scaled to the critical density and temperature. For the simulation data and the theoretical calculations based on an empirical equation of state, the law of corresponding states appears to be obeyed. However, the calculations based on the equation of state from ?rst order perturbation theory do not appear to scale well at all. This failure appears to be due to poor behavior of the critical density as a function of the cuto? and indicates that care must be exercised before using the law of corresponding states to extrapolate calculations.

Acknowledgments

This work was supported by the European Space Agency under contract number ESA AO-2004-070 and by the projet ARCHIMEDES of the Communaut? e Fran? caise de Belgique

10

2

2

1.5 γ
??

1.5

1

1

0.5

0.5

0 0

0.2 0.4 1-T/Tc

0 0

0.2 1-T/Tc

0.4

FIG. 4: (Color online) The scaled surface tension, γ ?? ≡

γ 2/3 , Tc ρc

as a function of distance from the

? = ∞), critical temperature. The left panel includes the theoretical curves, shown as full line (Rc ? = 8), dashed line (R? = 6), dash-dot line (R? = 4), dash-dot-dot line (R? = 3), dotted line (Rc c c c ? = 2.5), and the simulation data, shown as circles (R? = 6), ?lled squares and line+circles (Rc c ? = 4, 2000 atoms), open squares (R? = 4, 8000 atoms), diamonds (R? = 3) and triangles (Rc c c ? = 2.5). The right hand panel shows only the data from simulation as well as the estimated (Rc

error. In both cases, the thin line is a best ?t to all of the data (theory and simulation) of the form
? (1 ? T /T ) with γ ? = 3.41. γ ? = γ0 c 0

(ARC 2004-09).

[1] M. Mecke, J. Winkelmann, and J. Fischer, J. Chem. Phys. 107, 9264 (1997). [2] D. Duque and L. F. Vega, J. Chem. Phys. 121, 8611 (2004). [3] K. Katsov and J. D. Weeks, J. Phys. Chem. B 106, 8429 (2002). [4] J. F. Lutsko, J. Chem. Phys. 128, 184711 (2008). [5] P. R. ten Wolde and D. Frenkel, Science 77, 1975 (1997). [6] D. Frenkel and B. Smit, Understanding Molecular Simulation (Academic Press, Inc., Orlando,

11

2

1.5 γ
??

1

0.5

0 0

0.1

0.2

0.3 1-T/Tc

0.4

0.5

FIG. 5: (Color online) The scaled surface tension, γ ?? ≡

γ 2/3 , Tc ρc

as a function of distance from

the critical temperature as calculated using the perturbative equation of state displayed as a full
? = ∞), dashed line (R? = 6), dash-dot line (R? = 4), dash-dot-dot line (R? = 3), and line (Rc c c c ? = 2.5). line+circles (Rc

FL, USA, 2001). [7] C. H. Bennett, J. Comput. Phys. 22 (1976). [8] J. R. Errington and D. A. Kofke, J. Chem. Phys. 127 (2007). [9] E. Salomons and M. Mareschal, J.Phys.: Condens.Matter 3 (1991). [10] M. J. Haye and C. Bruin, J. Chem. Phys. 100, 556 (1994). [11] D. Duque, J. C. P` amies, and L. F. Vega, J. Chem. Phys. 121, 11395 (2004). [12] J. J. Poto? and A. Z. Panagiotopoulos, J. Chem. Phys. 112, 6411 (2000). [13] J. K. Johnson, J. A. Zollweg, and K. E. Gubbins, Molecular Physics 78, 591 (1993). [14] J. A. Barker and D. Henderson, J. Chem. Phys. 47, 4714 (1967). [15] J.-P. Hansen and I. McDonald, Theory of Simple Liquids (Academic Press, San Diego, Ca, 1986). [16] E. A. Guggenheim, J. Chem. Phys. 13, 253 (1945). [17] J. C. Le Guillou and J. Zinn-Justin, Phys. Rev. Lett. 39, 95 (1977).

12

[18] J. C. Le Guillou and J. Zinn-Justin, Phys. Rev. B 21, 3976 (1980). [19] H. Okumura and F. Yonezawa, J. Chem. Phys. 113, 9162 (2000). [20] J. P? erez-Pellitero, P. Ungerer, G. Orkoulas, and A. D. Mackie, J. Chem. Phys. 125, 054515 (2006).

13



更多相关文章:
Evaporation Preempts Complete Wetting_图文.pdf
of the liquid-vapor surface tension and the ...range of any intermolecular forces present, then ...(3) keeping terms with dependence on L , a ,...
CaF2-CaO-Al2O3-MgO-SiO2渣系表面张力计算模型.pdf
CaF2-CaO-Al2O3-MgO-SiO2渣系表面张力计算模型_电子/电路_工程科技_专业资料...dependence oftheSUrfaee tensionsinmolten slagsystemsas reproducedbythepresent...
...simulation study of the dimethyl sulfoxide liquidvapor_....pdf
the temperature dependence of surface tension, and the dependence of ...simulations on DMSO liquidvapor system in the temperature range of 298 ...
Surface tension effects on particle collection efficiency_....pdf
the dependence of particle collection efficiency on the surface tension of ...fliquid interaction time is greater than the time required for the vapor ...
Section 3 Electro-capillary curve_图文.ppt
to measure For liquid metal, it can be measured...dependence of surface tension on the applied ...Aliphatic compounds adsorb in a range of about 1...
Surface tension of liquid AlCuAg ternary alloys.pdf
Due to the high vapour pressure of silver, the...dependence of the surface tension on temperature ...decreasing cB in the range of cB \ 10 at.%....
Equation to estimate the surface tensions.pdf
tension (γms ) between the liquid metal and ...temperature dependence of the surface tension (dγ...range of stainless steel compositions in terms of...
...Calculation of the Dependence on Gravity of the Surface T_....pdf
of the Surface Tension and Bending Rigidity of a Fluid Interface’ In a ... calculation of the gravity dependence of properties of a liquid-vapor ...
Wettability and surface tension of fluorite_图文.pdf
of interaction; one is apolar (by Lifshitz-van...surface and their dependence on the drying ...liquid tension, yr is the liquid surface tension...
sk专业英语.doc
Surface tension can be modified to some extent ...of the dependence of electrode melting rate on ...within the range of low-current welding processes...
Molecular Dynamics Simulations and Surface Tension ....pdf
dependence (beyond 0.5 M) of surface tension ...toward the vapor phase with respect to that of ... and chloride in a wide range of concentrations....
物理化学(应化专业)教学大纲_New.doc
dependence of reaction enthalpies Chapter 3 The ...liquid surface Surface Tension Fundamental equations ...The interaction of light and matters; ...
RE2O3-MgO-SiO2(RE=La,Nd,Sm,Gd和Y)熔体表面张力的计算模型.pdf
The temperature and composition dependence of the surface tensions in molten RE2O3?MgO?SiO2 slag systems was reproduced by the present model using surface...
...of a liquid with zero surface-tension....pdf
Dependence of the liquid... 13页 免费 Stability in the Stefan ... ...behavior in a granular jet Emergence of a liquid with zero surface-tension...
Calculating Models on the Surface Tension of CaO-FeO-SiO2 ....pdf
molecule coexistence theory(IMCT)of slag structure and Butler’s equation.The temperature nda composition dependence of the surface tensions in molten CaO--...
...as a demonstration showing the effect of surface tension_....pdf
The Dependence of Surface tension on circumference of green balloons 160 140...The surface tension is most commonly referred to a property of a liquid. ...
MEASUREMENT OF DENSITY, SOUND VELCICITY, SURFACE TENSION, AND....pdf
SURFACE TENSION, AND VISCOSITY OF FREELY ...liquid range is comparable or even larger than ...The time dependence of the oscillations was ...
Behavior of the Dripping Faucet over a Wide Range of the Flow....pdf
Behavior of the Dripping Faucet over a Wide Range of the Flow Rate_专业...Schelly [2] showed that the dependence of the surface tension and of the...
Capillary-gravity waves The effect of viscosity on the wave ....pdf
the liquid density, γ the liquid-air surface tension, and g the ...The above equation may also be written as a dependence of wave speed c ...
MC2014-Chapter4_图文.ppt
and the range of the force is only several ... surface tension, viscosity and other physical and...Two groups are miscible both in liquid and solid...
更多相关标签:

All rights reserved Powered by 甜梦文库 9512.net

copyright ©right 2010-2021。
甜梦文库内容来自网络,如有侵犯请联系客服。zhit325@126.com|网站地图