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

    Parametric instability of a liquid metal sessile drop under the action of low-frequency alternating magnetic fields*

    2013-06-01 12:29:57LEIZuosheng雷作勝
    水動力學研究與進展 B輯 2013年2期

    LEI Zuo-sheng (雷作勝)

    Shanghai Key Laboratory of Modern Metallurgy and Materials Processing, Shanghai University, Shanghai 200072, China

    SIMAP-EPM Laboratory, Grenoble Institute of Technology/CNRS UMR 5266, France

    E-mail:lei_zsh@staff.shu.edu.cn

    GUO Jia-hong (郭加宏)

    Shanghai Institute of Applied Mathematics and Mechanics, Shanghai University, Shanghai 200072, China

    FAUTRELLE Yves, ERNST Roland, ETAY Jacqueline

    SIMAP-EPM Laboratory, Grenoble Institute of Technology/CNRS UMR 5266, France

    REN Zhong-ming (任忠鳴)

    Shanghai Key Laboratory of Modern Metallurgy and Materials Processing, Shanghai University, Shanghai 200072, China

    Parametric instability of a liquid metal sessile drop under the action of low-frequency alternating magnetic fields*

    LEI Zuo-sheng (雷作勝)

    Shanghai Key Laboratory of Modern Metallurgy and Materials Processing, Shanghai University, Shanghai 200072, China

    SIMAP-EPM Laboratory, Grenoble Institute of Technology/CNRS UMR 5266, France

    E-mail:lei_zsh@staff.shu.edu.cn

    GUO Jia-hong (郭加宏)

    Shanghai Institute of Applied Mathematics and Mechanics, Shanghai University, Shanghai 200072, China

    FAUTRELLE Yves, ERNST Roland, ETAY Jacqueline

    SIMAP-EPM Laboratory, Grenoble Institute of Technology/CNRS UMR 5266, France

    REN Zhong-ming (任忠鳴)

    Shanghai Key Laboratory of Modern Metallurgy and Materials Processing, Shanghai University, Shanghai 200072, China

    (Received November 23, 2012, Revised February 27, 2013)

    A 2-D mathematical model is developed in order to simulate a parametric electromagnetic instability oscillation process of a liquid metal droplet under the action of low frequency magnetic field. The Arbitrary Lagrangian-Eulerian (ALE) method and weak form constraint boundary condition are introduced in this model for implementation of the surface tension and electromagnetic force on liquid droplet free surface. The results of the numerical calculations indicate the appearance of various regimes of oscillation. It is found that according to the magnetic field frequency various types of oscillation modes may be found. The oscillation is originated from an instability phenomenon. The stability diagram of liquid metal droplet in the parameter space of magnetic frequency and magnetic flux density is determined numerically. The diagram is very similar to that found in the so-called parametric instability.

    sessile drop, parametric instability, low-frequency alternating magnetic fields, numerical simulation

    Introduction

    Alternating magnetic fields, including low-, medium- and high-frequency magnetic fields, are usually imposed on liquid metal pool or drops in order to control both internal flow field and free surface behavior. This is a problem of great scientific and technological interest in many processes[1], such as cold crucible[2], electromagnetic levitation[3], welding[4]and electron beam evaporation[5]. In metallurgy and liquid metal refining processes low frequency magnetic field may be used to create free surface agitation which promotes mass transfers in through the interface between the liquid metal and the covering slag or molten salt[1].

    In these applications, electromagnetic force due to interaction between imposing magnetic field and induced electrical currents may cause considerable free surface deformation. It has been observed by Fautrelle et al.[6]that the periphery of the free surface behavior of a thin pool of mercury is determined by the magnetic field amplitude and frequency when it is subjected to a vertical low-frequency magnetic field. They found various regular and irregular free surface patterns among which a structured wave is excited with relatively large magnetic field amplitude. This kind of wave with appearances of polygonal plane shape is considered to be a parametric-type instability caused by the alternating electromagnetic force.

    As an extension of the previous investigation, the present work develops a 2-D mathematical model based on the finite element method able to simulate the parametric instability behavior of liquid metal droplet under a vertical low-frequency magnetic field. This model uses the Arbitrary Lagrangian-Eulerian (ALE) method. The model solves the Navier-Stokes equations, continuity equation and quasi-static magnetic field equations described in the ALE frame. In order to treat surface tension and time-dependent electromagnetic force on the dynamic boundary, a weak form constraint boundary condition is deduced. This model allows direct simulation of transient flow field, magnetic field and free surface deformation of the liquid metal droplet. The free surface oscillation behavior under different magnetic amplitudes and frequencies are examined. Based on the calculation results, the stability diagram of liquid metal droplet in the magnetic field frequency and magnetic flux density parameter space is proposed.

    Fig.1 Sketch of the two-dimensional liquid droplet with free surface submitted to a vertical ac magnetic field

    1. Formulation of the problem

    Consider a thin layer, or a two-dimensional flat liquid metal droplet in the x-yplane with an arbitrarily initial free surface described by a continuous functionΓ= Γ(x, y, t), as shown in Fig.1. The droplet surface area is equal to the one of a circle with prescribed radius of r. It is submitted to a low-frequency AC magnetic field in the vertical direction, or z-direction described asB0= B0cos(2π ft) iz, where B0is the maximum value of the magnetic flux density and is uniform in the whole plane,f is the frequency whilst izis a unit vector along thez-axis. In this case, an electromagnetic force, named as the Lorentz force, will be generated from the interaction between induced electrical currents vectorJ(J being in the xyplane) and applied magnetic field. The electromagnetic force may cause two-dimensional liquid metal motion as well as free surface deformation in thexyplane. In fact, the droplet considered here is an infinite liquid column whose cross section is circular at rest.

    1.1 Magnetic field

    Because here only low-frequency magnetic field (f≤3.5 Hz) is taken into consideration, quasi-static approximation can be introduced and Maxwell-Ampere’s law and Faraday’s law are consequently written as

    where Bz,Hzare respectively the magnetic flux density and magnetic field in the liquid metal, which have only az-direction component. In this case,Bz= Bzizand Hz=Hziz.E is the electric field,uis fluid velocity, and σ,μ0and μrdenote respectively the electric conductivity, permeability in vacuum and relative permeability. Naturally, an electric insulation condition must be satisfied on the free surfaceΓ.

    1.2 Fluid flow equations

    The Navier-Stokes equations govern the velocity and pressure distribution, namely:

    where ρandηare density and dynamic viscosity of the fluid, andpis the pressure field. The body force termFemis the Lorentz force, defined as

    This can be determined from the resolution of magnetic field equations.

    The condition on the free surface separating the liquid pool from the atmosphere is expressed as

    where nis the unit normal directed into the atmosphere. The normal stressf0is caused by curvature-dependent surface tension.

    1.3 Moving mesh and boundary conditions for the fluid flow field

    A commercial finite element method package, the COMSOL Multiphysics?, was used to solve the Eqs.(3)-(8) on a freely moving deformed mesh, which constitutes the fluid domain. By applying the moving mesh (ALE) module in the package, the deformation of the mesh relative to the initial shape of the domain is computed using the Winslow smoothing method[7], in which two elliptic equations are applied for realizing moving mesh and is pre-defined in the COMSOL graphic user interface. On the boundary, the meshes move at the same velocity of the normal component of the fluid flow field. In the interior of the domain the meshes move in a free displacement way.

    The main difficulty of this problem is the implementation of the boundary condition on the dynamic free surface for the fluid flow equations. From the Serret Frenet equation ?t/?s = -Cn , wheresis the arc length,Cthe curvature,tandnthe tangent and normal vectors of free boundaryΓ, and the surface tension can be described as Fs=σs? t/?s, which has an opposite direction to the normal vector with a positive surface tension coefficientσs.

    where u_testand uT_testare respectively the test functions of the velocity components and their associated tangential derivative (See Appendix).

    1.4 Numerical solution

    Several kinds of liquid metal initial shape are taken into account (which will be described later), among which a typical one is an ellipse with long half axis of 30.1 mm and short half axis of 29.9 mm (the area is nearly equal to that corresponding to a circle with radius of 30 mm). In this domain, a mesh with 8 762 triangular elements and 97 973 degrees of freedom was created from unstructured triangular grids with fine elements on the boundary, where the maximum element size is 0.3 mm, and grows to coarse ones in the centre, where the maximum element size is 2.0 mm. Some even finer grids were tried in several cases to make sure the calculation results are grid independent. All the calculations began with the state that the liquid droplet is totally static.

    The full coupled system of unsteady, nonlinear PDEs from Eqs.(3) to (8) was solved by the finite element method using the COMSOL Multiphysics 3.5a. A time-dependent solver was used together with direct (UMFPACK) linear system solver. The absolute and relative tolerance parameters for the ODE solver were respectively 0.001 and 0.01. Computer with dual-core Intel-i7 M620 processors and 4 GB RAM running in windows 7 64-bit system was employed to obtain solutions presented here. Normally, it took one hour to get a one second whole coupled solution of the liquid metal droplet evolution behavior.

    1.5 Model validation

    Firstly, the behavior of a free oscillation liquid mercury droplet is calculated using the mathematical model described above when magnetic field B0is set to zero. The physical properties of mercury are listed in Table 1.

    Table 1 Physical properties of the mercury

    Fig.2 Free oscillation behavior of liquid mercury droplet without magnetic field

    Table 2 Eigen-frequencies fnof a two-dimensional liquid droplet with radius r =30 mm

    2. Calculation results and discussion

    2.1 Typical dynamic oscillation process of a liquid metal droplet

    Beginning from a given initial shape, the transient motion of the liquid mercury droplet can be calculated along with time. Figure 3 shows a typical path in which the droplet evolves to a structured, polygonal shape like wave in the end, where the mode number is 4. In this case, the initial shape is an ellipse with long half axis a =30.1 mm and short half axis b= 29.9 mm, with the maximum value of the magnetic flux density B0=0.014T and frequency f= 1.417 Hz.

    Fig.3 A typical displacements in xdirection of rightmost point under low-frequency magnetic field

    According to Fig.3, the curve of rightmost point displacements in the x direction can be divided into three stages. In the first stage, from the beginning to about 25 s, the liquid mercury droplet oscillates in a mode-2 pattern with a frequency of 0.448 Hz, which can be determined by the FFT analysis. The latter frequency is almost equal to the eigenfrequency of the mode-2. The amplitude decreases slowly. In the second stage, named as a transition stage, there are some disturbances emerging from the symmetric movement and the rightmost point oscillates in a way not as regular as in the first stage. Then, in the third stage, this point oscillates with amplitudes increasing at an exponential rate, and the frequency is 1.417 Hz, which is the same as the magnetic field. Finally a mode-4 drop pattern whose eigenfrequency is close to f is excited.

    Fig.4 Liquid mercury droplet oscillation behavior in one magnetic field cycle

    Figure 4(a) shows the rightmost point displacements in the x direction during one magnetic field cycle for f=1.417 Hz. Figure 4(b) shows the electromagnetic pressure, surface tension and total pressure on the rightmost point along with the time during one magnetic field cycle. One may find that the oscillation frequency of electromagnetic pressure is twice of the magnetic field. The latter result is in agreement with the fact that the electromagnetic forces oscillate at a frequency equal to2f .

    Figure 5 shows the topological shape of the liquid droplet at five typical moments as well as the flow field and total electric current density in the whole domain. In which the surface and streamline mean total current density with a maximum value of 1 989.3 A/m2, the arrow means the liquid metal flow velocity with a maximum value of 23.3 mm/s. It canbe found that the liquid droplet has an obvious mode-4 shape at the moment of 52.26 s. At that moment the electromagnetic pressure on the rightmost point has a positive value, which means that its direction is oriented toward the positive x-axis. But its value is lower than the surface tension toward inside there, so the negative total pressure causes the apexes move inward. At the moment of 52.44 s, the droplet reaches its equilibrium shape, and the electromagnetic pressure changes its direction to toward inside, whilst the velocity of this point reaches its maximum value. The point keeps moving inside to form a concave free surface deformation until the moment of around 52.59 s, when the quadrangle shape is inverted, and at that point both the electromagnetic pressure and surface tension have a direction toward the positive x-axis which causes the point to begin to move toward outside. After reaching its equilibrium position (t=52.77 s) finally it forms an apex at the moment 52.95 s and finishes a movement cycle with a little increased amplitude. After this period the oscillation begins a new cycle.

    Fig.5 Liquid droplet topological shapes at five typical moments, which are shown in Fig.4, as well as the flow field and total current density

    We have applied the same processes under different magnetic flux densities and frequencies. It has been found that at a given magnetic field frequency, there exists a critical magnetic flux density B0_cr. When the imposed magnetic flux density is lower than B0_cr, the liquid droplet admits oscillations which are gradually damped, and finally reaches a static circular state. Some details of the critical magnetic flux density will be given later.

    When the imposed magnetic flux density is higher thanB0_cr, a structured polygonal type pattern is triggered, and it evolves with the same ways as shown in Fig.4. The mode numberndepends on the imposed magnetic field frequency. The oscillation frequencies of these motions are exactly the same as the imposed magnetic field frequencies (half of the Lorentz force). This kind of oscillation is analogous to what happened in the so-called parametric instabilities as pointed out in previous work[11-14]. It must be pointed out that in the polygonal oscillation cases, the oscillation amplitudes increase exponentially after every oscillation cycle, which means the velocity in the liquid droplet keep increasing along with time. So a finer grid is needed in order to get converged solution of liquid metal flow field. However, due to our limited calculation capacity it was not possible to continue the computations for longer time period. Practically we stop the calculation once a well identified mode is triggered, so the maximum deformation of the liquid droplet predicted by numerical simulation is not as high as that observed in experiments.

    2.2 Effects of initial shape on the evolution process

    Several kinds of liquid droplet initial shapes were examined in order to make clear whether the dynamic evolution process depends on the disturbance at the beginning. Figure 6 shows dynamic topological evolutions in one magnetic field cycle of two typical calculation examples when the magnetic flux density is B0=0.012 T and frequency f=2.004 Hz, which is close to the eigen frequency of mode-5.

    One of the liquid droplet initial shapes is an ellipse with long half axisa =32 mm and short half axis b=28 mm. The calculation results show that pentagonal oscillation only begins after 65 s. The other initial shape is a circle with a weak irregular disturbance whose boundary is defined as

    where the disturbance sizes are d =4 mm and radiusa =30 mm.

    Fig.6 Initial shape effects on the evolution process (B0= 0.012 T,f=2.004 Hz)

    The calculation results show that it only need about 10 s before the pentagonal oscillation begins. It is noticeable that the shape of the polygon is nearly the same with different initial shape, even the apex positions are weakly different. Some other initial shapes were examined to make confirmation that the final dynamic oscillation process is independent of initial shape. It is found that it only depends on the magnetic field frequency and the magnetic strength.

    2.3 Effects of magnetic flux density on the oscillation

    amplitude growth rate

    As was described above, once the structured polygonal oscillation modes are excited, the oscillation amplitudes increase exponentially along with time. But under different magnetic flux density the oscillation amplitudes increase at different velocities, which is showed in Fig.7. In that figure the rightmost point displacements in the x direction along with time are still applied when the magnetic frequency f= 1.417 Hz and the magnetic flux density varies from 0.012 T and 0.020 T. It shows that the structured mode-4 oscillations are excited earlier with increasing magnetic field strength. The oscillation amplitudes S are fitted according to an exponentially function as S =exp(k t + q)where the coefficient kis defined as the growth rate of oscillation amplitudes along with time. It is confirmed from Fig.7(b) that the growth rate increases with the magnetic flux density.

    Fig.7 Structured oscillation amplitude growth rate along with time under different magnetic flux densities. Here, magnetic frequency is 1.417 Hz and mode

    2.4 Stability diagram in the f-Bparameter space

    The appearance of non-symmetric oscillation patterns occurs as soon as the magnetic field exceeds a threshold value. That value depends on the magnetic field frequency. We may construct a stability boundary in the parameter space (B0,f). The computed stability boundary is shown in Fig.8. The boundary consists of several V-shape structures, the so-called“tongues” of instability for various non-symmetric modes. Each tongue has its minimum which corresponds to the eigenfrequency of the modes. Such a behavior is encountered in the so-called parametric instability as previous experimental and analytical works have pointed out[12-15], and also observed in related experiments[6]. This kind of behavior is in agreement with the stability diagram of the Mathieu-type equation which for a given mode n of eigenfrequency fnput forth various tongues for fn,fn/2,fn/3, etc.. The strongest instability appears for fn. It is interesting to notice that when the magnetic frequency is between 0.6 Hz and 0.7 Hz, which is in the region near the half of mode-4, i.e.,f4/2, the mode-4 pattern is triggered instead of mode-2 or mode-3. It may be n otice d that the pres ent mode l is able to c atch notonlythemaintonguecenteredaroundfn,butalso weaker tongues of instability localized at fn/2. The latter results indicate that surface tension is well taken into account in the model, since the eigenfrequency are mainly related to surface tension (cf. Eq.(10)).

    Fig.8 Stability boundary of liquid mercury droplet in f-Bparameter space and same calculation examples of various polygonal modes excited in the instability regions

    For a given magnetic frequency, for example, f=1.417 Hz, from Fig.7 it can be seen that in the case ofB0=0.012 T, the liquid droplet oscillates with damping amplitude until 90 s, and when B0>0.014 T the polygonal oscillation modes are triggered in a relatively short time. To find a precise value of B0_crin every given frequency is a time-consuming work, so generally only the range covered the critical point is determined. Here, for example,B0_cris between 0.012 T and 0.014 T whenf=1.417 Hz. All the low limitations of these ranges are determined until the calculation time continues to 120 s and there is no polygonal oscillation modes which are triggered.

    Table 3 Calculation parameters for every case shown in Fig.8

    In the stability diagram of Fig.8, the low limitation points of every frequency are connected by one fine line, which means that when parameter is below this line, the liquid droplet will be stable finally. The upper limit points of every mode number were connected by one bold line, which means that if magnetic flux density is higher than the value on this line, some kind of polygonal modes will emerge. But we have not calculated the case where the magnetic flux densities are far higher than the up limitation points, this stability diagram is only valid when B0<0.1 T.

    Some typical calculation examples of various polygonal modes are also shown in the Fig.8, and the calculation parameters for every case are listed in Table 3.

    Compared with the experimental work[6], the numerical results presented here at least show qualitative agreement, in view of the polygonal oscillation from mode-2 to mode-7 and the structure of stability diagram putting forth tongues of instabilities is localized near the eigenfrequency of each mode. But it must be pointed out that the critical magnetic flux determined by two-dimensional numerical simulation is much lower than in the experiments. That means it is much easier to trigger the non-symmetric mode in the numerical simulation. The possible reasons caused the difference are the following. Firstly, our mathematical model is a two-dimensional one in which there is no gravity, in the experiments gravity is a kind of resistance who tend to restore the droplet to its equilibrium shape during the droplet oscillation. Secondly, in the simulation there is no friction between liquid droplet and the substrate plate and no wetting effects between them either, those are also additional resistance factors. Thirdly, there may be some numerical viscosity during the numerical solution. To improve the results,refiner grid is needed and the corresponding computation is time-consuming, so this work can only give qualitative similarities between the predicted and observed stability diagrams so far.

    In the experiments, there is an unstructured regimes characterized by fingers, ejection and void of liquid metal when the magnetic field is too high. In this numerical simulation we could not take this phenomenon into consideration because it is not suitable to simulate a physical process with topological structure variety by applying the moving mesh technique. This is a problem open for the future work.

    3. Conclusions

    A two-dimensional numerical model based on the finite element method in order to simulate the behavior of liquid metal droplet under a vertical low-frequency magnetic field has been developed. By app lying the ALE method and a weak form constraint boundary condition, the dynamic surface tension and electromagnetic force on dynamic boundary are successfully treated. By using a multi-physics approach, the whole coupled model gives some details on how a liquid two-dimensional droplet is evolving under the action of low-frequency magnetic field.

    The simulation results show that at a given magnetic frequency, when the magnetic flux density is lower than a threshold value, the liquid metal droplet is stable, while when the magnetic flux density is higher than a threshold value, a structured polygonal type oscillation will be triggered, and the growth rate of oscillation amplitude increases with increasing magnetic flux density. The oscillation mode number depends on the magnetic frequency and once a polygonal type oscillation is triggered, the oscillation frequency is the same as the magnetic frequency. This process is independent of initial shape. A stability diagram of liquid metal droplet in the magnetic frequency and magnetic flux density parameter space is proposed. The diagram has a structure similar to the stability diagram of the Mathieu-type equation. All these results prove that the appearance of non-symmetric pattern is a type of parametric electromagnetic instability.

    Acknowledgments

    This work was supported by the China Scholarship Council and Région Rh?ne-Alpes (France) for supporting Lei’s visiting in Grenoble. The hospitality of SIMAP-EPM is gratefully acknowledged as well.

    [1] FAUTRELLE Y., PERRIER D. and ETAY J. Free surface controlled by magnetic fields[J]. ISIJ International, 2003, 43(6): 801-806.

    [2] BOJAREVICS V., ROY A. and PERICLEOUS K. A. Magnetic levitation of large liquid volume[J]. Magnetohydrodynamics, 2010, 46(4): 339-351.

    [3] ETAY J., SCHETELAT P. and BARDET B. et al. Modelling of electromagnetic levitation-consequences on non-contact physical properties measurements[J]. High Temperature Materiales Processes, 2008, 27(6): 439-447.

    [4] ANDREEV O., POTHERAT A. and THESS A. Generation of liquid metal structures of high aspect ratio by application of an ac magnetic field[J]. Journal of Applied Physics, 2010, 107(12): 124093.

    [5] KOCOUREK V., KARCHER C. and CONRATH M. et al. Stability of liquid metal drops affected by a high-frequency magnetic field[J]. Physical Review E, 2006, 74(2): 026303.

    [6] FAUTRELLE Y., ETAY J. and DAUGAN S. Free-surface horizontal waves generated by low-frequency alternating magnetic fields[J]. Journal of Fluid Mechanics, 2005, 527: 285-301.

    [7] HERMANSSON J., HANSBC P. A variable diffusion method for mesh smoothing[J]. Communications in Numerical Methods in Engineering, 2003, 19(11): 897-908.

    [8] WALKLEY M. A., GASKELL P. H. and JIMACK P. K. et al. Finite element simulation of three-dimensional free-surface flow problems[J]. Journal of Scientific Computing, 2005, 24(2): 147-162.

    [9] AZUMA H., YOSHIHARA S. Three-dimensional largeamplitude drop oscillations: Experiments and theoretical analysis[J]. Journal of Fluid Mechanics, 1999, 393: 309-332.

    [10] LAMB H. Hydrodynamics[M]. Cambridge, UK: Cambridge University Press, 1975.

    [11] FAUTRELLE Y., PERRIER D. and ETAY J. Free surface controlled by magnetic fields[J]. ISIJ International, 2003, 43(6): 801-806.

    [12] SHEN C., XIE W. and WEI B. Parametrically excited sectorial oscillation of liquid drops floating in ultrasound[J]. Physical Review E, 2010, 81(4): 046305.

    [13] FAUTRELLE Y., SNEYD A. D. Surface waves created by low-frequency magnetic fields[J]. European Journal of Mechanics B-Fluids, 2005, 24(1): 91-112.

    [14] NOBLIN, X., BUGUIN A. and BROCHARD-WYART F. Triplon modes of puddles[J]. Physical Review Letters, 2005, 94(16): 166102.

    [15] NOBLIN X., BUGUIN A. and BROCHARD-WYART F. Vibrations of sessile drops[J]. The European Physical Journal Special Topics, 2009,166(1): 7-10.

    Appendixes: The free surface boundary condition in weak form of N-S equations

    Consider the N-S equation (A1) with boundary condition (A2) in a two-dimensional domainΩwith the boundary Γ,

    By introducing u _test,v_test, the test functions of velocity components in thexandydirections, the weak form of Eq.(A3) can be deduced by integrating over the domain,

    Consider the first term on the right side, then the x component, for instance, becomes

    By applying the boundary condition (A2) and the Serret Frenet equation ?t/?s = -Cn , we have

    By applying the integrating by parts along the boundary, the first item of right side, the surface tension item gives the contribution,

    In our case, the boundaryΓis a closed curve, so there is no point contribution, the second term of right side equals turns to zero.

    Here uTx _test= tx?u _test/?sis the notation in COMSOL frame, which means thexcomponent of the tangential derivative ofu _testwith t=(tx,ty) the tangential vector.u _testand uTx_testare all predefined in that finite element method package.

    Hence, finally, the wake form of N-S equation (A3) with boundary condition (A2) is,

    The first term on the right side are the pressure imposed on the boundary caused by surface tension. This expression demonstrates the equivalence of Eq.(A9) as claimed. They are implemented by the weak form constraint condition on the boundary setting in the equations in the COMSOL frame.

    10.1016/S1001-6058(13)60367-4

    * Project supported by the National Natural Science Foundation of China (Grant Nos. 51274137, 10872123).

    Biography: LEI Zuo-sheng (1975-), Male, Ph. D., Associate Professor

    GUO Jia-hong,

    E-mail:jhguo@staff.shu.edu.cn

    久久久久久久久中文| 国产私拍福利视频在线观看| 侵犯人妻中文字幕一二三四区| 99精品在免费线老司机午夜| 午夜免费激情av| 黄色视频,在线免费观看| 午夜免费观看网址| 免费一级毛片在线播放高清视频 | 亚洲第一av免费看| 午夜两性在线视频| 在线观看www视频免费| 老熟妇乱子伦视频在线观看| 国产熟女xx| 我的亚洲天堂| 一本综合久久免费| 9191精品国产免费久久| 美女扒开内裤让男人捅视频| 母亲3免费完整高清在线观看| 亚洲精华国产精华精| 久久人妻福利社区极品人妻图片| 999久久久国产精品视频| 精品少妇一区二区三区视频日本电影| 欧美 亚洲 国产 日韩一| 午夜久久久久精精品| 自拍欧美九色日韩亚洲蝌蚪91| 69av精品久久久久久| 久久精品人人爽人人爽视色| 久久精品国产亚洲av高清一级| 99精品欧美一区二区三区四区| 色综合亚洲欧美另类图片| 欧美+亚洲+日韩+国产| 一区二区三区激情视频| 亚洲一区中文字幕在线| 国产视频一区二区在线看| 国产高清视频在线播放一区| 侵犯人妻中文字幕一二三四区| 成人永久免费在线观看视频| 天堂影院成人在线观看| 精品国产国语对白av| 9色porny在线观看| 国产人伦9x9x在线观看| 午夜福利成人在线免费观看| 国产在线精品亚洲第一网站| 成人国语在线视频| 亚洲欧美精品综合久久99| 天堂影院成人在线观看| 国产精品 欧美亚洲| 国产成人av教育| 色av中文字幕| 免费在线观看亚洲国产| 久久久久久久精品吃奶| 国产一卡二卡三卡精品| 久久久国产欧美日韩av| 久久亚洲精品不卡| 日韩中文字幕欧美一区二区| 色尼玛亚洲综合影院| 国产高清有码在线观看视频 | 啦啦啦韩国在线观看视频| 欧美绝顶高潮抽搐喷水| 在线十欧美十亚洲十日本专区| 亚洲自偷自拍图片 自拍| 免费在线观看影片大全网站| 国产又色又爽无遮挡免费看| 99在线人妻在线中文字幕| 亚洲精华国产精华精| 男人操女人黄网站| 51午夜福利影视在线观看| 香蕉丝袜av| 日韩大码丰满熟妇| 欧美丝袜亚洲另类 | 成在线人永久免费视频| 搡老妇女老女人老熟妇| 国产极品粉嫩免费观看在线| 亚洲天堂国产精品一区在线| 99久久99久久久精品蜜桃| 成人国产综合亚洲| 人人妻人人澡欧美一区二区 | svipshipincom国产片| 色播在线永久视频| 日韩免费av在线播放| 国产精品香港三级国产av潘金莲| 淫秽高清视频在线观看| 午夜精品在线福利| 欧美色欧美亚洲另类二区 | 国产三级黄色录像| 一级片免费观看大全| 老司机午夜十八禁免费视频| www.精华液| 久久国产亚洲av麻豆专区| 国产熟女xx| 亚洲专区中文字幕在线| 国产欧美日韩一区二区三| 亚洲性夜色夜夜综合| 精品欧美国产一区二区三| 国产区一区二久久| 真人做人爱边吃奶动态| 久热爱精品视频在线9| 99国产极品粉嫩在线观看| 久热这里只有精品99| 黄色视频,在线免费观看| 午夜福利视频1000在线观看 | 乱人伦中国视频| 欧美黑人精品巨大| videosex国产| 午夜免费激情av| 欧美成人免费av一区二区三区| 国产99久久九九免费精品| 亚洲专区国产一区二区| 精品久久久久久,| 啦啦啦 在线观看视频| 1024视频免费在线观看| 日韩一卡2卡3卡4卡2021年| 国产日韩一区二区三区精品不卡| 两性午夜刺激爽爽歪歪视频在线观看 | 波多野结衣高清无吗| 国产私拍福利视频在线观看| 搡老妇女老女人老熟妇| 国产精品秋霞免费鲁丝片| 日韩大尺度精品在线看网址 | 一级作爱视频免费观看| 欧美 亚洲 国产 日韩一| 国产主播在线观看一区二区| 欧美日韩福利视频一区二区| 国产精品国产高清国产av| 国产一区在线观看成人免费| 在线观看免费日韩欧美大片| 精品久久久精品久久久| 国产一卡二卡三卡精品| 国产亚洲欧美98| 757午夜福利合集在线观看| 少妇被粗大的猛进出69影院| 免费在线观看黄色视频的| 欧美日本视频| 亚洲片人在线观看| 午夜福利视频1000在线观看 | 大型av网站在线播放| 正在播放国产对白刺激| 99久久久亚洲精品蜜臀av| 国产一区二区三区在线臀色熟女| 日韩欧美在线二视频| 国产99白浆流出| 国产精品爽爽va在线观看网站 | 久久午夜综合久久蜜桃| 999精品在线视频| 亚洲av电影在线进入| av片东京热男人的天堂| 夜夜看夜夜爽夜夜摸| 亚洲成av片中文字幕在线观看| 国产国语露脸激情在线看| 欧美绝顶高潮抽搐喷水| 91精品三级在线观看| 性少妇av在线| 曰老女人黄片| 国产三级黄色录像| 久久久久精品国产欧美久久久| 999精品在线视频| 久久中文字幕一级| 九色国产91popny在线| 18禁观看日本| 中亚洲国语对白在线视频| 精品国产亚洲在线| 成在线人永久免费视频| xxx96com| 国产成人系列免费观看| 亚洲色图 男人天堂 中文字幕| 热re99久久国产66热| 亚洲国产看品久久| 人人妻人人澡人人看| 无限看片的www在线观看| 男女做爰动态图高潮gif福利片 | 久久婷婷成人综合色麻豆| 久久久国产精品麻豆| 亚洲欧美精品综合久久99| 久久精品91蜜桃| 国产av精品麻豆| 欧美国产日韩亚洲一区| 日韩免费av在线播放| 两个人视频免费观看高清| 精品第一国产精品| 国产一区在线观看成人免费| 久久久精品国产亚洲av高清涩受| 咕卡用的链子| 日本精品一区二区三区蜜桃| 亚洲av电影不卡..在线观看| 成人国产一区最新在线观看| 1024视频免费在线观看| 中文字幕精品免费在线观看视频| 日韩三级视频一区二区三区| 国产片内射在线| 欧美+亚洲+日韩+国产| 国产精品久久久久久人妻精品电影| 9色porny在线观看| 男女下面进入的视频免费午夜 | 成人国语在线视频| 久久久久久人人人人人| 又大又爽又粗| 超碰成人久久| 国产三级黄色录像| 午夜老司机福利片| 在线观看舔阴道视频| 波多野结衣高清无吗| 亚洲熟妇中文字幕五十中出| cao死你这个sao货| 欧美国产日韩亚洲一区| 亚洲中文字幕日韩| 91精品三级在线观看| 一区在线观看完整版| 精品国产乱子伦一区二区三区| 老司机靠b影院| 欧美激情 高清一区二区三区| 人成视频在线观看免费观看| 免费在线观看黄色视频的| 黄色a级毛片大全视频| 亚洲三区欧美一区| 亚洲第一青青草原| 久久精品成人免费网站| 国内久久婷婷六月综合欲色啪| 亚洲中文字幕日韩| 90打野战视频偷拍视频| 一卡2卡三卡四卡精品乱码亚洲| 亚洲午夜精品一区,二区,三区| 中文字幕人成人乱码亚洲影| 亚洲黑人精品在线| 国产亚洲欧美精品永久| 国产免费男女视频| 精品国产一区二区三区四区第35| 亚洲一码二码三码区别大吗| 久久伊人香网站| 精品不卡国产一区二区三区| 久久精品国产清高在天天线| 午夜免费激情av| 日韩大尺度精品在线看网址 | 好男人在线观看高清免费视频 | 欧美乱色亚洲激情| 97碰自拍视频| 国产私拍福利视频在线观看| 在线观看午夜福利视频| 老司机午夜福利在线观看视频| 久久青草综合色| 在线观看www视频免费| 色综合亚洲欧美另类图片| 亚洲黑人精品在线| svipshipincom国产片| 老司机在亚洲福利影院| 国产高清激情床上av| 夜夜爽天天搞| 国产人伦9x9x在线观看| 亚洲电影在线观看av| 人人妻人人爽人人添夜夜欢视频| tocl精华| 亚洲 欧美一区二区三区| 国产精品久久久久久亚洲av鲁大| 精品国产乱码久久久久久男人| 欧美日韩亚洲国产一区二区在线观看| 美女国产高潮福利片在线看| 婷婷六月久久综合丁香| 亚洲成av人片免费观看| 一级毛片精品| av中文乱码字幕在线| 狠狠狠狠99中文字幕| svipshipincom国产片| 日本黄色视频三级网站网址| 十八禁网站免费在线| 久久久精品国产亚洲av高清涩受| 亚洲国产精品合色在线| 亚洲,欧美精品.| 非洲黑人性xxxx精品又粗又长| 99国产精品一区二区三区| 午夜福利一区二区在线看| 97超级碰碰碰精品色视频在线观看| 欧美日韩黄片免| 免费搜索国产男女视频| 亚洲最大成人中文| 亚洲国产看品久久| 如日韩欧美国产精品一区二区三区| 九色国产91popny在线| 又黄又粗又硬又大视频| 亚洲 国产 在线| 亚洲精品久久成人aⅴ小说| 久久九九热精品免费| 波多野结衣巨乳人妻| 精品一品国产午夜福利视频| 久久久久久久久久久久大奶| 一级毛片高清免费大全| 国产高清激情床上av| 操出白浆在线播放| 麻豆一二三区av精品| 欧美av亚洲av综合av国产av| 久久久久九九精品影院| 在线国产一区二区在线| 欧美黑人精品巨大| 一进一出抽搐动态| 亚洲精品一卡2卡三卡4卡5卡| 国产精品二区激情视频| 极品人妻少妇av视频| 9热在线视频观看99| www.999成人在线观看| 国产单亲对白刺激| 一级黄色大片毛片| 欧美日韩中文字幕国产精品一区二区三区 | 看免费av毛片| av天堂在线播放| 亚洲五月天丁香| 久久青草综合色| 无遮挡黄片免费观看| 久久久久久久久久久久大奶| 看免费av毛片| 久久人人精品亚洲av| 久久精品国产亚洲av高清一级| 制服人妻中文乱码| 人人妻,人人澡人人爽秒播| 亚洲精品久久国产高清桃花| 亚洲av片天天在线观看| 人人妻人人爽人人添夜夜欢视频| 老司机午夜福利在线观看视频| 夜夜看夜夜爽夜夜摸| 狠狠狠狠99中文字幕| 热re99久久国产66热| 88av欧美| 欧美日韩亚洲国产一区二区在线观看| 99久久久亚洲精品蜜臀av| 一级,二级,三级黄色视频| 黑人操中国人逼视频| 精品乱码久久久久久99久播| 成人亚洲精品av一区二区| 久久欧美精品欧美久久欧美| 国产亚洲av嫩草精品影院| 欧美不卡视频在线免费观看 | 久久人人97超碰香蕉20202| 18禁裸乳无遮挡免费网站照片 | 999久久久国产精品视频| 久久亚洲真实| 亚洲午夜理论影院| 99国产精品免费福利视频| svipshipincom国产片| 久久国产乱子伦精品免费另类| 中文字幕人妻熟女乱码| 精品少妇一区二区三区视频日本电影| 久久精品国产亚洲av高清一级| 丝袜美足系列| 国产精品二区激情视频| 丰满人妻熟妇乱又伦精品不卡| 在线视频色国产色| 免费在线观看黄色视频的| 成人特级黄色片久久久久久久| 高清黄色对白视频在线免费看| 精品国产国语对白av| 黄色视频不卡| 91大片在线观看| 精品一品国产午夜福利视频| 女人精品久久久久毛片| 亚洲 国产 在线| 麻豆国产av国片精品| 777久久人妻少妇嫩草av网站| 免费女性裸体啪啪无遮挡网站| 日日爽夜夜爽网站| 他把我摸到了高潮在线观看| 中文字幕色久视频| 乱人伦中国视频| 亚洲人成电影免费在线| 亚洲中文字幕日韩| 亚洲av第一区精品v没综合| 亚洲av片天天在线观看| 亚洲av成人av| 午夜福利成人在线免费观看| 色在线成人网| 亚洲欧美激情综合另类| 欧美亚洲日本最大视频资源| a级毛片在线看网站| 欧美国产精品va在线观看不卡| 淫秽高清视频在线观看| 精品国产乱子伦一区二区三区| 久久精品国产清高在天天线| 一卡2卡三卡四卡精品乱码亚洲| 久久国产乱子伦精品免费另类| 黄色片一级片一级黄色片| 亚洲专区国产一区二区| 91av网站免费观看| 成人欧美大片| 国产精品久久久人人做人人爽| 免费不卡黄色视频| 99久久久亚洲精品蜜臀av| 久久天躁狠狠躁夜夜2o2o| 国产熟女xx| 变态另类丝袜制服| 欧美大码av| 啦啦啦免费观看视频1| 国产极品粉嫩免费观看在线| 女同久久另类99精品国产91| 香蕉丝袜av| 国产激情欧美一区二区| 两个人免费观看高清视频| 国产麻豆成人av免费视频| 精品国产亚洲在线| av网站免费在线观看视频| 成人永久免费在线观看视频| 啦啦啦 在线观看视频| 中文字幕最新亚洲高清| 成人三级黄色视频| 热99re8久久精品国产| 亚洲五月色婷婷综合| 99国产综合亚洲精品| 亚洲 欧美一区二区三区| 18禁国产床啪视频网站| 欧美另类亚洲清纯唯美| 女人高潮潮喷娇喘18禁视频| 国产熟女午夜一区二区三区| 精品不卡国产一区二区三区| 满18在线观看网站| 韩国av一区二区三区四区| 亚洲精品国产精品久久久不卡| 国产成人啪精品午夜网站| 在线视频色国产色| 午夜老司机福利片| 丝袜在线中文字幕| 一边摸一边抽搐一进一小说| 欧美日韩瑟瑟在线播放| 天天添夜夜摸| 国产成人精品久久二区二区免费| 18禁国产床啪视频网站| 亚洲人成伊人成综合网2020| 黄色a级毛片大全视频| avwww免费| 日本一区二区免费在线视频| 精品高清国产在线一区| 黄色女人牲交| www.精华液| 最好的美女福利视频网| 免费在线观看黄色视频的| 黄网站色视频无遮挡免费观看| 亚洲七黄色美女视频| 久久久久久久久中文| 精品国产美女av久久久久小说| 大码成人一级视频| 波多野结衣高清无吗| 国产熟女午夜一区二区三区| 国产99白浆流出| 日韩高清综合在线| 色老头精品视频在线观看| 丝袜美腿诱惑在线| 丁香六月欧美| 天堂动漫精品| 亚洲天堂国产精品一区在线| avwww免费| 老司机午夜福利在线观看视频| 欧美+亚洲+日韩+国产| 国产av又大| 激情视频va一区二区三区| av网站免费在线观看视频| 欧美精品啪啪一区二区三区| 成在线人永久免费视频| 一区二区三区国产精品乱码| av中文乱码字幕在线| 国产欧美日韩综合在线一区二区| 色在线成人网| 首页视频小说图片口味搜索| 欧美+亚洲+日韩+国产| 99在线视频只有这里精品首页| 久久精品亚洲精品国产色婷小说| 19禁男女啪啪无遮挡网站| 久久国产精品人妻蜜桃| 亚洲美女黄片视频| 久99久视频精品免费| 久久久久久人人人人人| 亚洲精华国产精华精| 成年版毛片免费区| 最近最新中文字幕大全免费视频| 国产成人av激情在线播放| 欧美老熟妇乱子伦牲交| 正在播放国产对白刺激| 中文字幕人成人乱码亚洲影| 女人被狂操c到高潮| 波多野结衣av一区二区av| 老司机午夜十八禁免费视频| 亚洲av成人不卡在线观看播放网| 久久久精品国产亚洲av高清涩受| 日本免费一区二区三区高清不卡 | 亚洲精品一卡2卡三卡4卡5卡| 一级片免费观看大全| 激情在线观看视频在线高清| 国产精品av久久久久免费| 男男h啪啪无遮挡| 变态另类丝袜制服| 欧美一级a爱片免费观看看 | 国产亚洲精品一区二区www| 天堂√8在线中文| 一卡2卡三卡四卡精品乱码亚洲| 97人妻精品一区二区三区麻豆 | 满18在线观看网站| 超碰成人久久| 88av欧美| 国产在线精品亚洲第一网站| 国产高清激情床上av| 亚洲avbb在线观看| 欧美在线黄色| 国产亚洲欧美98| 亚洲欧美日韩无卡精品| 最新在线观看一区二区三区| av中文乱码字幕在线| 亚洲成a人片在线一区二区| 国产男靠女视频免费网站| 操美女的视频在线观看| 精品人妻1区二区| 桃色一区二区三区在线观看| 国产野战对白在线观看| 亚洲成国产人片在线观看| 亚洲国产精品sss在线观看| 狂野欧美激情性xxxx| 无遮挡黄片免费观看| 中文字幕色久视频| 日本欧美视频一区| 老司机午夜十八禁免费视频| 国产av在哪里看| 久久草成人影院| 国产在线观看jvid| 欧美大码av| 久久久精品欧美日韩精品| 两性午夜刺激爽爽歪歪视频在线观看 | 97碰自拍视频| 色尼玛亚洲综合影院| 久久久久久久午夜电影| 国产精品99久久99久久久不卡| 99香蕉大伊视频| 大型黄色视频在线免费观看| 老司机午夜十八禁免费视频| 久久久久亚洲av毛片大全| 两性夫妻黄色片| 欧美日韩黄片免| 亚洲色图 男人天堂 中文字幕| 精品国产亚洲在线| 香蕉国产在线看| 中文亚洲av片在线观看爽| 国产精华一区二区三区| 999精品在线视频| 久久性视频一级片| 国产一区二区三区视频了| 国产欧美日韩综合在线一区二区| 天天躁夜夜躁狠狠躁躁| 老司机午夜十八禁免费视频| 国产精品久久久av美女十八| 法律面前人人平等表现在哪些方面| 亚洲片人在线观看| 久久国产精品人妻蜜桃| 午夜免费观看网址| 淫妇啪啪啪对白视频| 久久精品影院6| 精品高清国产在线一区| 久久久精品欧美日韩精品| 女人爽到高潮嗷嗷叫在线视频| 亚洲九九香蕉| 好看av亚洲va欧美ⅴa在| 欧美激情高清一区二区三区| 欧美 亚洲 国产 日韩一| 69av精品久久久久久| 少妇 在线观看| 中国美女看黄片| 欧美成人午夜精品| 久久精品91无色码中文字幕| 侵犯人妻中文字幕一二三四区| 欧美在线一区亚洲| 九色亚洲精品在线播放| 亚洲人成电影观看| 国产午夜福利久久久久久| 久久香蕉激情| 精品人妻1区二区| 人妻久久中文字幕网| 精品人妻在线不人妻| 黄色丝袜av网址大全| 欧美黄色片欧美黄色片| 久久 成人 亚洲| 国内毛片毛片毛片毛片毛片| 曰老女人黄片| 亚洲成av人片免费观看| 亚洲全国av大片| 成人国语在线视频| 亚洲片人在线观看| 国产在线精品亚洲第一网站| 欧美日韩中文字幕国产精品一区二区三区 | 9色porny在线观看| 久久久久久国产a免费观看| 国产伦一二天堂av在线观看| 啪啪无遮挡十八禁网站| 亚洲中文字幕日韩| 一进一出抽搐动态| 亚洲国产毛片av蜜桃av| 亚洲国产中文字幕在线视频| 亚洲激情在线av| 国产激情欧美一区二区| 欧美日韩亚洲综合一区二区三区_| 成年版毛片免费区| 在线观看免费日韩欧美大片| 日本黄色视频三级网站网址| 国产色视频综合| www.自偷自拍.com| 日本撒尿小便嘘嘘汇集6| 国产xxxxx性猛交| 91大片在线观看| 男女床上黄色一级片免费看| 91在线观看av| 日韩av在线大香蕉| 亚洲三区欧美一区| bbb黄色大片| 国产一区在线观看成人免费| 亚洲在线自拍视频| 淫秽高清视频在线观看| 九色亚洲精品在线播放| 国产麻豆成人av免费视频| 久久狼人影院| 日韩大尺度精品在线看网址 | 精品久久久久久成人av| 午夜福利,免费看| 精品久久久精品久久久| 亚洲av熟女| 精品国产亚洲在线| 亚洲色图 男人天堂 中文字幕|