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

    Exploration of the coupled lattice Boltzmann model based on a multiphase field model: A study of the solid–liquid–gas interaction mechanism in the solidification process

    2024-03-25 09:33:08ChangShengZhu朱昶勝LiJunWang王利軍ZiHaoGao高梓豪ShuoLiu劉碩andGuangZhaoLi李廣召
    Chinese Physics B 2024年3期
    關(guān)鍵詞:李廣

    Chang-Sheng Zhu(朱昶勝), Li-Jun Wang(王利軍), Zi-Hao Gao(高梓豪),3,Shuo Liu(劉碩), and Guang-Zhao Li(李廣召)

    1College of Computer and Communication,Lanzhou University of Technology,Lanzhou 730050,China

    2State Key Laboratory of Gansu Advanced Processing and Recycling of Non-Ferrous Metal,Lanzhou University of Technology,Lanzhou 730050,China

    3Institute of Modern Physics,Chinese Academy of Sciences,Lanzhou 730000,China

    Keywords: multiphase field model,lattice Boltzmann model(LBM),Shan-Chen multiphase flow,solidification organization

    1.Introduction

    With the development of the field of materials science and engineering,the study of the behavior and interaction mechanisms of complex multiphase fluids during solidification has become increasingly important.Bubbles and dendrites are two common microstructures in the solidification process,and their morphological evolution and interactions during solidification have an important impact on the properties and structure of materials.[1]Therefore,an in-depth understanding of the interaction mechanism between bubbles and dendrites as well as the influence of various factors on the solidification tissue morphology is of great theoretical significance and application value.[2,3]

    In order to investigate the growth evolution mechanism of solidified tissues as well as the bubble formation mechanism,scholars have conducted a large number of experimental studies.However, the opacity of metallic materials largely limits researchers’ real-time observation of tissue evolution during alloy solidification.With the rapid development of computer technology,numerical simulation has become an effective way to solve such complex problems.[4]

    In recent years, scholars have widely applied the phasefield method(PFM),level-set method(LSM),and cellular automata(CA)methods in microporosity simulations.The goal of these simulations is to study the generation and evolution of microporosity in order to obtain simulation results that couple the real-time mesoscopic and microscopic scale microstructure morphology.Among them, the phase field method[5,7](PFM) is a simulation method that describes the kinetic evolution process of the system based on thermodynamics, and it has many advantages over other methods, such as the ability to accurately simulate a variety of complex microstructures without having to track the complex solid-liquid interface.[8,9]In addition, multiphase flow systems play an important role in natural and industrial processes, and the lattice Boltzmann method[10](LBM), as a mesoscopic simulation method between macroscopic continuum simulation and microscopic molecular dynamics simulation, has been developing rapidly in the field of computational fluid dynamics(CFD)and shows great potential.[11]Currently,many scholars have successfully simulated complex multiphase flow systems using LBM.[12]

    Meidaniet al.[13]developed a two-dimensional/threedimensional phase field model to simulate the growth process of irregular pores confined by the dendrite lattice using a multiphase field approach,and investigated the effect of factors such as the solid-phase fraction on the pores morphology.However, the model did not consider the effect of dendrite growth.Zhanget al.[14]developed a solid-liquid-gas multiphase field lattice Boltzmann model for the study of complex multiphase flow interactions during solidification of pure metals,and they developed a mixed-phase-field lattice Boltzmann method as well as a parallel adaptive lattice refinement algorithm for simulating the interactions between the gas bubbles and the solid-growth fronts during the solidification process and successfully solved the solid growth and bubble evolution problems using the phase field method.Sunet al.[15]successfully simulated the formation of pores and microstructure evolution of binary sub-eutectic Al-Si-Mg alloys during solidification using a two-dimensional multicomponent multiphase metatronic cellular automata model,considered the solute and hydrogen partitioning at the solid-liquid interface, and simulated the growth of dendrites and solute migration using the CA-FDM method.Guet al.,[16]on the other hand,used meta cellular automata to simulate the coupled process of hydrogen porosity evolution and dendrite growth during the solidification of the alloy,calculated the diffusion of solute and hydrogen using the finite difference method,observed the competing growth and porosity between columnar dendrites and isometric dendrites, and explored the effect of the cooling rate on the radius of the porosity; it is worth noting that the model can only describe the static state of the gas bubbles during the solidification process.In addition, Zhanget al.[17]proposed a coupled LB-CA-FD model for simulating the formation of hydrogen pores during the solidification of binary subeutectic aluminum alloy dendrites.They used the LBM simulation to calculate the hydrogen transport as well as the nucleation,growth and motion of bubbles.

    It should be noted that the meta cellular automata (CA)method is mostly used in the above studies.As for the cellular automata (CA) method compared with the phase field(PF) method, the PF method can more realistically simulate the complete process of dendrite growth.[18]Moreover,due to the relevant simulation techniques,there are almost no studies on the interaction between dendrites with different orientations and gas bubbles during the directional solidification process.Therefore,in the face of the above research gap in this field,a multiphase field simulation coupled lattice Boltzmann method(PF-LBM)is proposed in this study for simulating the distribution mechanism of bubbles and solutes at the solid-liquid interface during solidification of SCN-0.24wt.%butanedinitrile alloys,the interactions between dendrites and bubbles,as well as the effects of different temperature gradients, anisotropic strengths,and tilting angles on the solidified organization.The model combines the advantages of multiphase field simulation of dendrite growth and Shan-Chen multiphase flow-based LBM simulation of bubble motion,and is able to describe the multiphase flow phenomena during solidification more comprehensively.Through this study, we reveal the important physical mechanisms of bubble nucleation,growth and motion evolution and interaction with dendrite growth during solidification.These research results are of great significance for the in-depth understanding of the behavior of complex multiphase flow during solidification, and provide new ideas and theoretical basis for the preparation and property optimization of materials.

    2.Multiphase field model and Shan–Chen model for LBM multiphase flow

    2.1.Multi-phase field mathematical model

    A multi-phase field mathematical model is a physical model used to describe a system of matter with different phases.In this model, matter is viewed as consisting of multiple phases, and each phase has different characteristics and properties.These phases can be transformed through interactions and form specific spatial and temporal distributions throughout the system.The multi-phase field mathematical model provides an effective tool for describing the phase transition process.Its emergence effectively solves the problem of describing the interaction and evolution laws of multiple phases(different states)in matter, and has important applications and significance in multiple fields.In this study,a quantitative multi-order parametric phase field model for dilute binary alloys developed by Ofori-Opokuet al.[19]is used,where the following multi-phase field equations are used in this paper:

    Assuming negligible solute diffusion in the solid phase,the temperature is replaced by an approximation of the condensation temperatureT=Tref+G(x-Vpt),whereTrefis the reference temperature and the position above the solidification axis is denoted byx.τis the interfacial relaxation time,Wis the thickness of the diffusive interfacial layer,φiis a whole representing the phase field variables,Uis the liquid flow velocity,kis the solute partition coefficient at the solid-liquid interface.The strength of the interfacial energy anisotropy is given by

    whereas=1-3ε4,θare the direction of the thermal gradient and the angle perpendicular to the grain interface,respectively,andΘi0is the top angle of the first grain interface.

    Thec0l=c0/kis the concentration of the liquid portion of the solid-liquid interface, the coupling constantλrepresents the width of the diffusion interface,d0=a1W0/λis the solute capillary length,τ0is the interfacial relaxation time,x,yare the horizontal and vertical axes respectively,Gis the temperature gradient at the solid-liquid interface.ε4is the anisotropic strength of the interfacial free energy,and the solute saturation in Eq.(3) dimensionlessUis obtained by solving the solute conservation equation with the solute equation

    Due to the influence of the interface width selection, the scales simulated by the phase field method model are relatively small, which is difficult to realize in the real process, so the thin interface theory is used to analyze the interface so that the selection of the interface width is greater than the width of the actual interface.And when the width of the interface increases, solute shedding occurs, and this occurrence introduces a degree of computational error.In order to minimize the error generated in this part will be introduced in the solute equation the inverse solute trap flow proposed by Karmaet al.[19]in the anti-solute trapping flow.τ0=a1a2W30/(d0Dl),

    wherea1,a2take the results of the derivation of the asymptotic analysis of thin interfaces by Karma and Rappel,[19,21]a1=0.8839 anda2=0.6762.When the sample is fixed and cooled at a given temperature gradientGat cooling rateR,the rate of movement of the isotherm is called the pulling speedVp=R/Gand the inverse solute trap flowJatin Eq.(3)is

    2.2.Shan–Chen model for LBM multiphase flow

    The Shan-Chen model is a multiphase flow model in the lattice Boltzmann method(LBM)used to simulate phase separation phenomena in complex fluid systems.LBM is a numerical simulation method based on the discrete grids,which approximately discretizes a continuous medium into discrete spatial and time step grids, and describes fluid dynamics behavior using the microscopic velocity distribution function of particles.LBM has the advantages of strong parallelism and ease of handling complex geometric shapes and boundary conditions.The main use of different mutual characteristics can be simulated by adjusting model parameters to handle various multiphase flow situations such as gas-liquid, solid-liquid,and so on.The following are the equations of the LBM multiphase flow Shan-Chen model used in this article.The Boltzmann evolution equation for LBM multiphase flow lattice can be expressed as

    Regarding the particle equilibrium state distribution function calculation expression:ρσ(r,t)anduσ(r,t)represent the density and velocity of the fluid in theσcomponent, respectively,ωiis the weight coefficient, in the D2Q9 model,the value ofωiis taken as follows:ω0=4/9,ω1-4=1/9,ω5-8=1/36,eirepresents the velocity of the particle in theidirection of motion,andeiis the configured velocity

    The value of pressureρ(r,t) at the lattice point of the fluid phase can be obtained from the following equation of state:

    In the formula,ρ(r,t) for the density of the fluid, when different phases are transformed into each other, will cause the density of the particles change, such as gas to liquid conversion with the increase of the gas phase liquid phase will be gradually reduced,the following for the gas-liquid phase particles density changes in the relationship is

    whereρgandρldenote the particle densities of the gas and liquid phases of the fluid, respectively, Δφis the amount of conversion from the liquid phase to the gas phase, which is very much related to the amount of conversion of the liquid phase to the solid phase.Δfsdenotes the amount of conversion of the particles from the liquid phase to the solid phase,andAis the coefficient,which is taken to be the value of 10-6in this paper.Coupling the multiphase field model with the LBM multiphase flow Shan-Chen model, whereCandDare the hydrogen concentration as well as diffusion coefficients,respectively, andKis the hydrogen partition coefficient, the evolution of bubbles and dendrites growth can be solved coupled by the following equations:

    2.3.Physical parameters

    SCN-0.24wt.% butanedinitrile alloy was used as the object of study,the grid spacing of the computational range was set tox=1W0,and the time step was selected ast=0.002τ0.Periodic boundary conditions were applied to the LBM model for simulation of the periodic fluid behavior during the interval of solidified tissue calculation in order to facilitate a better study of the phenomena between different states.The dendrite tilt angle was measured according to the angle of rotation in the counterclockwise direction.The physical parameters used in this paper are shown in the table.

    Table 1.Calculation parameters used in the simulation process.

    2.4.Validation of the PF-LBM model

    whereGσadsis the coefficient of interaction between fluidsolid, which determines the strength of the force between the fluid particles and the solid phase,s(r+eiΔt,t)eiis the labeling function, which is used to label the solid-liquid phase composition of the current lattice node,with a value of 1 when it is the solid phase,and a value of 0 when it is the liquid phase.

    In order to verify the correctness of the PF-LBM model,the simulation of the wetting phenomenon between the solidliquid-gas three-phase is used to verify our model.Figure 1 shows the schematic diagram when the bubble is in a stable state after contact with a smooth solid-phase surface.At the solid-liquid-gas three-phase contact point,the angleθf(wàn)ormed between the horizontal line of the solid-gas interface and the tangent line of the gas-liquid interface is the contact angle.The contact angle can be calculated by the following expression:

    whereγlg,γls, andγsgare the gas-liquid, solid-liquid and solid-gas interfacial tensions,respectively.

    Fig.1.Relationship between dendrites and bubbles at different temperature gradients.

    There is a relationship between the contact angleθwhen it has three specific values and the flow-fluid interaction coefficientGcin Eq.(13)and the flow-solid interaction coefficientGadsin Eq.(23)

    whereψis a function that calculates the effective density of the fluid with the expression of Eq.(14).

    Equation(25)can be understood as follows:whenθ=0°,i.e., there is no gas-phase portion above the unit circle at the contact point O in the schematic diagram of Fig.1,so only the solid-phase portion and the liquid-phase portion remain in the unit circle.At this point the density of the phase above the unit circle at the contact point O is only the density of the liquid phase.Whenθ=180°,it follows from the same reasoning that the density of the phase above the unit circle at the contact point O at this point is only the density of the gas phase.Whenθ=90°, at this point, the upper left part of the unit circle at point O has only the liquid phase portion,and the upper right part has only the gas phase portion, so the density above the unit circle at point O is half of the sum of the density of the gas phase and the density of the liquid phase.

    Fig.2.Relative errors of the three special contact angles with respect to the contact angles obtained by simulation.

    Fig.3.Morphology of the bubbles on the surface of the smooth solid phase: (a)Gads=-4.77,(b)Gads=-3.90,(c)Gads=-3.03,(d)Gads=-2.16,(e)Gads=-1.29,(f)Gads=-0.42.

    In the cases ofGc=-6.0 andGc=-6.5, respectively,theGadstheoretical values of the three special angle values are calculated,and then theGadstheoretical values are substituted into the model for calculation, and the obtained contact angle simulation values are compared with the three special contact angle values,respectively,to calculate the relative errors.Equation(26)is used to calculate the size of the obtained contact angle simulation values and format(27)is used to calculate the contact angle simulation values and special contact angle relative errors,respectively,

    In Eq.(26),His the value of the length of the vertical line from the highest point of the bubble to the surface of the solid phases,andLis the value of the bottom of the contact between the bubble and the surface of the solid phases, i.e., the value from point O to point A in the schematic diagram of Fig.1.In Eq.(27),ERis the relative error,θsis the simulated angle of the obtained contact angle,andθtis the theoretical value of the three special contact angles.The simulation region is set to be 120×120,the initial density of the liquid phaseρl=2.8,the initial density of the gas phaseρg=0.07,the relaxation time of both phases is 1.0, and the initial radius of the bubbleRis set to be 15.The bottom boundary condition of the model is used to simulate the solid phase,and the rest of the boundaries are used for the periodic condition.From Fig.2,it can be seen that the relative error does not exceed 2%,which indicates the correctness of our model.Figure 3 shows the different contact angle topographies resulting from the contact of a bubble with a smooth solid phase surface under theGc=-6.0 condition.The red region in the figure is the bubble and the surrounding red-blue region is the liquid phase.From the figure, it can be seen that the contact angle obtained from the simulation is gradually increasing as the fluid-solid interaction coefficientGadsincreases,i.e.,the wettability between the fluid phase and the solid phase is gradually decreasing.

    3.Results and discussion

    3.1.Interaction between dendrites and bubbles at different temperature gradients

    In the numerical simulation of directional solidification of alloys,numerous physical parameters have a significant impact on the morphology changes of dendrites, and dendrite morphology is of great significance for studying the interaction between dendrites and bubbles.Therefore, statistical characterization of some interesting quantities and exploring the internal mechanisms of research from different appearances are issues that we are more concerned about.Among them, temperature gradient and anisotropic strength are the two main factors that affect the morphology of dendrites.

    Figure 4 demonstrates the morphological changes of the interaction between dendrites and bubbles under different temperature gradient conditions.At the beginning of the simulation, we set up a region as a square of 800×800 grids,and three directionally growing dendrite nuclei were set up in the region at the beginning of solidification.With time, the dendrites began to grow slowly and the solid-liquid interface gradually advanced upward.The density of the liquid at the left bottom of the region is lower than that of the liquid at the right bottom, and the high-density region is gradually transformed into a solid state with the prolongation of solidification time.According to the law of conservation of energy,the density of low density areas will decrease.When the liquid density drops to the critical value of the solid-liquid interface,the hydrogen concentration in the dendrites gap will exceed the hydrogen concentration at the initial solidification,resulting in the formation of bubbles.Due to the tendency of high-density areas to promote bubble formation,bubbles are more likely to form and develop within the area.This process indicates that in the early stage of solidification, dendrites and bubbles will grow together, and high-density liquid regions are conducive to the formation of bubbles.

    At low temperature gradients, as shown in Figs.4(a1)-4(a3), the solute content around the solid-liquid interface is higher,resulting in a relatively low concentration of the liquid phase, which generates some small bubbles in the dendrites interstices.Due to the hydrogen concentration difference,the growth of dendrites simultaneously inhibits the gas bubbles from entering the liquid phase, and the solute concentration around the solid-liquid interface is higher than that elsewhere,which provides the conditions for the gas bubbles to be transformed.Bubbles between dendrites have an important effect on dendrite growth.Gas bubbles start in a dispersed state,but they fuse together due to the pressure difference at the solidliquid interface.Thus, bubble formation promotes dendrite growth,and in the middle stage of solidification,the high concentration of solute inhibits the growth of secondary dendrite arms.

    As the temperature gradient is further increased,as shown in Figs.4(b1)-4(b3), the gas bubbles absorb hydrogen atoms around the surrounding liquid-phase dendrites during growth,leading to overflow,enlargement and fusion.As the secondary dendrite arms grow and the bubbles enlarge and fuse,the liquid and gas generate different pressures, creating a pressure difference in the interstitial space.This causes the bubbles to move irregularly between the dendrites, being squeezed from the secondary dendrite arms into the interstices of the dendrite arms,causing the bubbles to deform,or fuse with another bubble.For example,the bubble M1 in Figs.4(b1)-4(b3)initially appears within the secondary dendrite arm, and as solidification proceeds,the bubble increases in size and passes through the secondary dendrite arm to the primary dendrite arm gap M2.The increase in size and deformation of the bubbles is accompanied by the growth of the dendrites, resulting in the formation of bubble M3,which takes on different shapes.

    Fig.5.Evolution of bubble mean radius with time step for temperature gradients of G = 0.0002 K/μm, G = 0.0004 K/μm, and G =0.0006 K/μm.

    Fig.6.Evolution of the total fraction ratio between bubbles and solids for temperature gradients of G=0.0002 K/μm,G=0.0004 K/μm,and G=0.0006 K/μm.

    When the temperature gradient increases to 0.0006, as shown in Figs.4(c1)-4(c3),the initial bubble N1 grows in the interstices of the dendrites,and the co-growth of the dendrites and bubbles leads to the extrusion and deformation of the bubble at N2, forming two parts.As the dendrite grows over a larger area, the bubble at N3 is more significantly extruded and deformed.Over time,the growth of dendrites and bubbles compete together for growth space, and the growing bubbles absorb solute excretion,contributing to the growth of dendrites and the development of secondary dendrite arms.At the same time, there is also competition between bubbles, while some bubbles grow,other bubbles may disappear or fuse,and large bubbles are more likely to aggregate with surrounding small bubbles to form irregular large bubbles.Figure 5 shows the variation curves of the average radius of bubbles with time step for temperature gradientGat 0.0002, 0.0004 and 0.0006, respectively.It can be seen that the curve fluctuates more whenG=0.0002.As the temperature gradientGreaches 0.0006,the bubble radius becomes smaller and has better continuity.Figure 6 shows the variation curves of the total fraction ratio between bubbles and solids for temperature gradientGat 0.0002, 0.0004 and 0.0006, respectively.The curves in the figure show that the bubbles are formed gradually during the solidification organization process, the radius of the bubbles is getting bigger and bigger, and the dendrites do not stop growing because of the appearance of the bubbles,but are still growing.The curve shows a monotonous upward trend, and the bubble continues to get bigger and bigger accompanied by a decrease in the ratio of bubble to solid, and finally until stabilized, which shows that the bubble has a great role in promoting the solidification of dendrites.

    3.2.Interaction between dendrites and bubbles at different anisotropic strengths

    Under different anisotropic strength conditions, we perform simulations of dendrite growth and stomatal evolution modeling to quantitatively explore the effect of anisotropic strength on the interaction between oriented columnar dendrites and bubbles during solidification.We set the computational range as a square area of 800×800 grids,each of which represents a computational grid.The temperature gradient in the whole area is assumed to be uniform and the pulling speed is kept constant.The anisotropy intensity is increased from 0.02 to 0.06 to observe the evolution process of dendrites and bubbles.

    At the beginning of the solidification simulation, we set up three initial nuclei in the calculation area.As the solidification process advances, the number of dendrites starts to increase,and it can be clearly observed from Fig.7 that the oriented dendrites and pores change with time, and the primary dendrites grow gradually and interact with the bubbles.From Figs.7(a1)-7(a3), it can be seen that when the anisotropic strength is small, the primary dendrites are thicker and generate a large number of secondary dendrite arms.At the same time,solutes accumulate in the gap between the primary dendrite and the secondary dendrite arms, and the growth of the secondary dendrite arms is slowed down by the accumulation of solutes.When the dendrite grows to a critical value, hydrogen bubbles emerge from the low end of the dendrite(e.g.,at L1) and then continue to increase until the local hydrogen concentration falls below the hydrogen saturation or the dendrite grows around it hindering the growth of the bubbles.And the lush secondary dendrite arms keep colliding with the bubbles, forming a squeeze on the bubbles, and then the bubbles deform,and part of the bubbles will be divided into two parts by the secondary dendrite arms.Some other bubbles are compressed by the secondary dendrite arms, which leads to the fusion phenomenon, as in Figs.7(a1)-7(a3), the bubbles at L1 gradually migrate and fuse together (at L2).With the increase of the time step,the volume fraction of microhydrogen bubbles increases, and the pores also increase, and when the pressure difference in dendrites is larger than the critical value of the bubbles it will divide the bubbles into two parts (e.g.,at L3),and the density of the bubbles is not much affected by the increase of the time step.It is worth noting that the solutes released from the solid phase in the dendrites gap cannot be completely discharged and accumulate in the dendrites gap resulting in the inability of the secondary dendrites arms to continue to develop,and the effect of the secondary dendrites arms on the bubble will gradually slow down.

    With the further increase of the anisotropic strength,the primary dendrite spacing becomes correspondingly larger and the secondary dendrite arms gradually decrease whenε4=0.04.At the same moment, compared with the smaller anisotropy strength,the interaction between the dendrites and the bubbles slows down,and the deformation rate of the bubbles decreases significantly.As shown in Figs.7(b1)-7(b3),due to the content of hydrogen concentration of bubbles in the dendrite organization growth process is less than the concentration in the dendrite, bubbles are crowded by the dendrite piled together, such as M1, the dendrite’s continuous growth process, the density of the dendrite is greater than that of the hydrogen bubbles that is divided into two sections, such as M2, with the increase of the time step, the dendrite channel bubbles accompanied by dendrites competing for the growth of the dendrite infinite close to the connected together(e.g.,at M3).

    Fig.7.Relationship between dendrites and bubbles under different anisotropy.

    Whenε4=0.06,the overall morphology of the dendrite becomes elongated, only part of the dendrite has weak secondary dendrite arms,and the surface of the primary dendrite arms of the rest of the dendrite is nearly smooth.As shown in Figs.7(c1)-7(c3),only a small number of secondary dendrite arms are in contact with the bubbles, and the primary dendrite arms of the rest of the dendrites have less contact with the bubbles,and the interaction between the two is very weak.Among them, the dendrite at N1 will be close to the bubble during the growth process,and when the bubble grows to the tip of the dendrite,it will be interspersed by the tip of the dendrite through the bubble to divide the original bubble into two bubbles of different sizes (as shown at N2).Until the liquid phase that grows to the solid-liquid interface contains more hydrogen than in the hydrogen bubbles, it engulfs the small bubbles,causing them to disappear(shown at position N3).

    As shown in Fig.8,we count the average radii of the bubbles under different anisotropic strength conditions.It is worth noting that, when the anisotropy strength is increasing, along with the primary dendrite arms become finer and their number increases,and the surface of the dendrite is also gradually smooth, the secondary dendrite arms are reduced, and the interstitial channel between the dendrite arms is wider, and the dendrite gradually reduces the pressure on the bubbles,resulting in the bubbles from more continuity to the existence of a separate dendrites channel, which is less prone to the splitting and fusion phenomena.It can be found that the larger the anisotropy strength is,the more stable the curve fluctuation expressed by the mean radius of bubbles is, and also due to the reduction of fusion phenomenon of different bubbles,thus the mean radius of bubbles is also relatively reduced.

    Fig.8.Variation curves of bubble mean radius with time step for anisotropy strengths of ε4=0.02,ε4=0.04,and ε4=0.06.

    3.3.Interaction of tilted dendrites with gas bubbles

    The interrelationship between dendrites and gas bubbles during directional solidification at an initial angle of dendrite growthΘi0=70°is illustrated in Figs.9(a1)-9(a3).As the solidification process progresses,due to the relatively low solubility of hydrogen atoms in the solid phase, the hydrogen atoms are repelled by the solid-liquid interface and thus undergo migration into the liquid phase.When the hydrogen concentration in the liquid phase exceeds a critical value,bubbles begin to form and start generating from the interstices of the dendrites and move irregularly,in which case the morphology of the dendrites will not be changed by the deformation of the bubbles.Bubbles show irregular motion patterns as they are captured by the solid phase.The coupled growth of bubbles and dendrites leads to the appearance of micropores in the solidified tissue.It is well known that the defects in the microstructure due to bubble generation are mainly due to the trapping of bubbles at the solid-liquid interface.The interaction between bubbles and dendrites directly determines the formation of defects in the microstructure.In Fig.9(a1) at L1 the bubble does irregular movement between the dendrites,due to the solid-liquid interface and the different densities of the hydrogen bubbles,the bubble is divided into two sections by the dendrites while passing through the dendrites, and the bubble absorbs the hydrogen atoms from the dendritic interstices into itself during the growth of solidified tissues,so that the two sections of the differently shaped bubbles undergo deformation and are fused together (e.g., Fig.9(a2)).With the continuous growth of the dendrite,hydrogen atoms will be excluded by the solid-liquid interface,which will lead to the reduction of hydrogen atoms in the liquid phase,along with the growth of bubbles, the microstructure of the dendrite grows relatively fast, when the bubble encounters the dendrite during the growth process,its growth will be inhibited,and in the dendrite gap in the new position with other bubbles connected together to form a worm-like bubble as shown in Fig.9(a3).It is worth noting that due to the tilted growth of the dendrite,the growth direction of the secondary dendrite arm on the left side of the dendrite is closer to the direction of the temperature gradient, resulting in the secondary dendrite arm on the left side of the dendrite being relatively more developed,which will be more helpful for the interaction between the dendrite and the bubbles,resulting in bubbles being more susceptible to deformation and the phenomena of multi-bubble fusion and bubble splitting.

    Figures 9(b1)-9(b3) reflect the interaction between dendrite and bubble when the initial angle of dendrite growth is 70°and 90°,respectively.Dendrite A is grown along the direction of temperature gradient, and the preferential orientation angle of dendrite B is 70°.Dendrite A and dendrite B show convergence in the left layer region and divergence in the right side, and the interstitial channels of converging bicrystals are narrower than those of parallel dendrites, and the interstitial channels of diverging bicrystals are wider than those of parallel dendrites.Bubble growth exists in the interstices of the dendrites,and the continuous growth of the dendrites with tilt angles of 70°and 90°creates the convergence phenomenon,and the bubbles are located above the dendrites at this moment(as in Fig.9 M1).Subsequently, with the continuous growth of the dendrites,the bubble grows to be in infinite contact with the converging dendrites on both sides(as in Fig.9 M2).

    When the dendrite contacts the bubble, the gap between the converging dendrites decreases, the bubble is compressed and deformed, and finally the bubble is squeezed by the converging dendrites in three directions in the shape of a flower(as in Fig.9 M3).In the space of divergent twin crystals on the right,the bubble(as in Fig.9 N1)is sandwiched in a narrow channel and squeezed into a strip shape by the secondary dendrite arms.As the dendrite continues to grow, the angle between the two dendrites grows wider and wider, giving the bubbles enough space to migrate in motion.The bubbles pile up into the wide channel (as in Fig.9 N2).When the angle is large enough,the bubbles are extruded from the narrow channel to the uppermost position to fuse with other bubbles to form larger bubbles (as in Fig.9 N3).It can be seen that when the bubbles located in the gap between the dendrites interact with the dendrites at different inclination angles, when the bubbles are located between the converging dendrites,the bubbles will be squeezed more obviously, and the interaction between the dendrites and the bubbles will be more intense.When the bubbles exist in the space formed by the divergent bicrystals,due to the gap between the bicrystals is expanding,and the bubbles will have a narrow channel through the divergent tendency of the dendrites to be squeezed into the larger dendritic gaps to form large bubbles.

    Figure 10 shows the curves of the average radius of bubbles with time when the growth angle of both kinds of dendrites is 70°.It can be seen that when the growth direction of the dendrites is the same, the gap between the dendrites is more uniform,the radius of the bubbles is in a constant state,and their areas are roughly the same, which is smoother and less volatile in the curve.Figure 11 reflects the curve of the average radius of the bubbles as a function of the time step when the bicrystal forms a pinch point.When the dendrites perpendicular to the 70°dendrites to form the angle, such as Fig.11,AΘi0=70°dendrites gap from wide to narrow the formation of convergent dendrites, the bubble will be squeezed from large to small.The radius of the curve tends to decrease from large to small downwards, and vice versa the formation of evanescent dendrites bubbles within the narrow channels of the dendrites.As the dendrites diverge, the gap suddenly becomes larger, and under pressure the bubbles rush out and become larger in the narrow channel.

    Fig.10.Curve of average radius of bubbles versus time for growth angle of 70° for both dendrite AΘio and dendrite BΘio.

    Fig.11.Curve of average radius of bubbles versus time when dendrite A and dendrite B are formed at angles of AΘi0=70° and BΘi0=90°.

    Fig.12.Relationship between bubble growth rate and time step for different dendrite growth angles.

    Figure 12 examines the bubble growth rate curves of dendrite A and dendrite B under different growth angles.When the growth angle of both dendrite A and dendrite B is 70°,the gap channel between dendrites is narrower,the size of bubbles is relatively uniform, the bubble movement is in a relatively uniform speed, and the bubble growth rate curve is relatively smoother.When the growth angle of dendrite A is 70°, the growth angle of dendrite B is 90°, the bubble growth rate at the beginning of solidification is in a smooth state, with the formation of dispersive bicrystals, the solid-liquid interface and the pressure difference between the hydrogen pore lead to the bubbles be rushed out to the wider channels to migrate a little faster.

    4.Conclusion

    A multiphase field model coupled with an LBM model is successfully proposed.For in-depth investigation of the complex phenomenon of solid-liquid-gas interactions between the three phases of an alloy during solidification.With this model,we are able to simulate the distribution mechanism of bubbles and solutes at the solid-liquid interface to explore the interactions between dendrites and bubbles, as well as the effects of different temperatures,anisotropies,and tilt angles on the solidification organization.In this study, we have drawn some interesting conclusions and findings through an in-depth analysis of the relevant physical parameters and elucidated the physical mechanisms of bubble nucleation, growth and kinematic evolution during solidification as well as the interactions with the growth of dendrites.

    First, our simulation results show that during the directional growth of columnar dendrites, bubbles will precipitate out at the bottom of the dendrite gap, and then gradually rise due to the differences in solid-liquid densities and differential pressures, and this process involves interactions between the bubbles and the dendrites,such as extrusion,spillage,fusion,and disappearance.All these phenomena have a significant effect on the solidification organization.

    Secondly,in cases where the gap between dendrite channels is wide,bubbles have more opportunities to fuse and form irregularly larger bubbles.At the dense channels,the bubbles are deformed by the extrusion of the dendrites.This suggests that the channel structure has an important influence on the evolution and morphology of bubbles,which is of great engineering significance in the actual solidification process.

    Finally,our study also reveal that with the divergence and convergence of dendrites in the simulation area, bubbles will precipitate out under the compression and diffusion of dendrites, and also lead to physical phenomena such as bubble fusion and overflow.This phenomenon manifests the dynamic complexity of the solidification process and provides a strong theoretical basis for better understanding and controlling the material solidification process.

    This study has successfully investigated in depth the interaction between bubbles and dendrites during the solidification process through the application of the PF-LBM model,which provides important insights into our understanding of the physical mechanisms involved in the solidification process of materials.These research results not only help to optimize the solidification of industrial arts of materials, but also provide strong support for related engineering applications and useful information for further development in the field of materials science and engineering.In the future,we can continue to extend this model to deepen our understanding of the microscopic phenomena in the solidification process and to promote further research in related fields.

    Acknowledgments

    Project supported by the National Natural Science Foundation of China (Grant Nos.52161002, 51661020, and 11364024), the Postdoctoral Science Foundation of China(Grant No.2014M560371), and the Funds for Distinguished Young Scientists of Lanzhou University of Technology of China(Grant No.J201304).

    猜你喜歡
    李廣
    念(外一首)
    How to Arouse Middle School Students’ Interest in English Study
    速讀·中旬(2018年7期)2018-08-17 07:22:00
    李廣的故事
    隨機(jī)微分方程的樣本Lyapunov二次型估計(jì)
    An improved potential field method for mobile robot navigation①
    飛將軍李廣
    小伙伴們一起跑
    水中的精靈
    喙嘴龍
    花生
    大型av网站在线播放| 日韩精品中文字幕看吧| 国内精品久久久久精免费| 国产精品永久免费网站| 久久伊人香网站| 人人澡人人妻人| 一区在线观看完整版| 午夜福利在线观看吧| 免费在线观看黄色视频的| 欧美黄色淫秽网站| 黑人欧美特级aaaaaa片| bbb黄色大片| 国产xxxxx性猛交| 欧美亚洲日本最大视频资源| 一级毛片高清免费大全| 亚洲性夜色夜夜综合| 精品久久久久久,| 纯流量卡能插随身wifi吗| 成年版毛片免费区| 日韩国内少妇激情av| 怎么达到女性高潮| 看黄色毛片网站| 丝袜美足系列| 日本一区二区免费在线视频| av视频在线观看入口| 久久中文字幕人妻熟女| 一个人观看的视频www高清免费观看 | 成熟少妇高潮喷水视频| 亚洲精品国产一区二区精华液| 精品久久久久久成人av| 亚洲精品国产一区二区精华液| 国产区一区二久久| 美女 人体艺术 gogo| 午夜福利成人在线免费观看| 国产成人av激情在线播放| 国产激情欧美一区二区| 法律面前人人平等表现在哪些方面| 十分钟在线观看高清视频www| 这个男人来自地球电影免费观看| 亚洲熟女毛片儿| 成人国产综合亚洲| 国产午夜福利久久久久久| 国产男靠女视频免费网站| 国产不卡一卡二| 一区福利在线观看| 亚洲国产看品久久| 久久人人97超碰香蕉20202| 亚洲人成网站在线播放欧美日韩| 午夜日韩欧美国产| 久久 成人 亚洲| 成年女人毛片免费观看观看9| 国产区一区二久久| 狂野欧美激情性xxxx| 国产亚洲欧美98| 亚洲欧美精品综合久久99| 麻豆一二三区av精品| 少妇 在线观看| 无遮挡黄片免费观看| 日韩三级视频一区二区三区| 亚洲色图 男人天堂 中文字幕| 国产成人欧美在线观看| 咕卡用的链子| 国产欧美日韩一区二区三区在线| 老司机靠b影院| 国产一级毛片七仙女欲春2 | 免费女性裸体啪啪无遮挡网站| 亚洲久久久国产精品| 精品一品国产午夜福利视频| 婷婷精品国产亚洲av在线| 欧美成人午夜精品| 一级毛片女人18水好多| 搡老熟女国产l中国老女人| 又黄又粗又硬又大视频| 免费无遮挡裸体视频| 不卡一级毛片| 色老头精品视频在线观看| 侵犯人妻中文字幕一二三四区| 免费在线观看日本一区| 黄色丝袜av网址大全| 日本在线视频免费播放| 麻豆一二三区av精品| 91在线观看av| 国产人伦9x9x在线观看| 亚洲精品国产色婷婷电影| 国产区一区二久久| 性少妇av在线| x7x7x7水蜜桃| 国产国语露脸激情在线看| 成年人黄色毛片网站| 欧美激情高清一区二区三区| 国产又爽黄色视频| 国产精品香港三级国产av潘金莲| 亚洲一码二码三码区别大吗| 精品日产1卡2卡| 一本大道久久a久久精品| 久久久久亚洲av毛片大全| 日韩视频一区二区在线观看| 老司机午夜福利在线观看视频| 久久精品国产综合久久久| 亚洲精品一区av在线观看| 欧美激情高清一区二区三区| 级片在线观看| 成人国产一区最新在线观看| 99国产精品一区二区蜜桃av| 日韩欧美国产在线观看| 妹子高潮喷水视频| 免费高清视频大片| 欧美日韩福利视频一区二区| 国产成+人综合+亚洲专区| 丁香欧美五月| 这个男人来自地球电影免费观看| 国产成人av激情在线播放| 欧美日韩黄片免| 美女高潮到喷水免费观看| 美女 人体艺术 gogo| 久久国产亚洲av麻豆专区| 亚洲精品在线美女| 久久久精品欧美日韩精品| 成人国语在线视频| 欧美精品亚洲一区二区| 亚洲,欧美精品.| a在线观看视频网站| 国产精品一区二区免费欧美| 黑人巨大精品欧美一区二区蜜桃| 久久精品人人爽人人爽视色| 国产97色在线日韩免费| 在线观看免费视频日本深夜| 久久久久久国产a免费观看| 精品无人区乱码1区二区| 亚洲午夜理论影院| 亚洲片人在线观看| 男人的好看免费观看在线视频 | 欧美乱妇无乱码| 欧美最黄视频在线播放免费| 精品人妻1区二区| 久久中文字幕一级| 九色国产91popny在线| 亚洲国产欧美日韩在线播放| 久久精品人人爽人人爽视色| 久久精品国产99精品国产亚洲性色 | 日本黄色视频三级网站网址| 国产精品精品国产色婷婷| 成熟少妇高潮喷水视频| 国产精品综合久久久久久久免费 | 久久草成人影院| 国产精品98久久久久久宅男小说| 久久精品国产亚洲av香蕉五月| 高清黄色对白视频在线免费看| 日韩欧美在线二视频| 麻豆av在线久日| 亚洲第一电影网av| 欧美日韩黄片免| 久久久精品欧美日韩精品| 人人妻人人澡人人看| 久久久久精品国产欧美久久久| 搞女人的毛片| 男女下面进入的视频免费午夜 | 色在线成人网| 韩国av一区二区三区四区| 每晚都被弄得嗷嗷叫到高潮| 亚洲av美国av| 亚洲精品久久国产高清桃花| 丰满人妻熟妇乱又伦精品不卡| 中文字幕久久专区| 9191精品国产免费久久| 精品卡一卡二卡四卡免费| 女性被躁到高潮视频| 天天躁狠狠躁夜夜躁狠狠躁| 国产色视频综合| 色播在线永久视频| 一区福利在线观看| 极品人妻少妇av视频| 欧美 亚洲 国产 日韩一| 自线自在国产av| 久久人妻福利社区极品人妻图片| 在线播放国产精品三级| 欧美日本亚洲视频在线播放| 午夜福利视频1000在线观看 | 久久亚洲精品不卡| 69精品国产乱码久久久| 久久精品亚洲熟妇少妇任你| 亚洲男人天堂网一区| 人人妻人人澡人人看| 欧美绝顶高潮抽搐喷水| 国产一区二区三区综合在线观看| 可以免费在线观看a视频的电影网站| 久久久久国产精品人妻aⅴ院| 国产蜜桃级精品一区二区三区| 亚洲成av人片免费观看| 两性午夜刺激爽爽歪歪视频在线观看 | 青草久久国产| 黄色a级毛片大全视频| 亚洲伊人色综图| 亚洲欧美一区二区三区黑人| 淫妇啪啪啪对白视频| 欧美日本亚洲视频在线播放| 久久香蕉激情| 女同久久另类99精品国产91| 欧美日韩乱码在线| 国产1区2区3区精品| 日韩欧美国产一区二区入口| 亚洲五月天丁香| 欧美最黄视频在线播放免费| 亚洲视频免费观看视频| 成年女人毛片免费观看观看9| 精品不卡国产一区二区三区| 久久中文字幕人妻熟女| 窝窝影院91人妻| 在线观看一区二区三区| 夜夜夜夜夜久久久久| 亚洲三区欧美一区| 一个人免费在线观看的高清视频| 法律面前人人平等表现在哪些方面| 老汉色∧v一级毛片| 俄罗斯特黄特色一大片| 中文字幕最新亚洲高清| 国产亚洲精品第一综合不卡| 久久香蕉激情| 久久国产精品影院| 国产高清有码在线观看视频 | 变态另类丝袜制服| 老司机午夜福利在线观看视频| 欧美成人午夜精品| 黄频高清免费视频| 日韩大尺度精品在线看网址 | 熟女少妇亚洲综合色aaa.| 亚洲人成网站在线播放欧美日韩| 国产精品自产拍在线观看55亚洲| 亚洲欧洲精品一区二区精品久久久| 午夜精品国产一区二区电影| 中亚洲国语对白在线视频| 麻豆一二三区av精品| 国产精品久久久人人做人人爽| bbb黄色大片| 侵犯人妻中文字幕一二三四区| 色综合站精品国产| 人妻丰满熟妇av一区二区三区| 村上凉子中文字幕在线| 黄色片一级片一级黄色片| 国产精品98久久久久久宅男小说| 丝袜美腿诱惑在线| 91精品国产国语对白视频| 久久狼人影院| 无人区码免费观看不卡| 国产亚洲av高清不卡| 久久中文字幕人妻熟女| 亚洲中文日韩欧美视频| 亚洲精品美女久久久久99蜜臀| 久久性视频一级片| 9热在线视频观看99| 久久精品影院6| www.精华液| 神马国产精品三级电影在线观看 | 18禁美女被吸乳视频| 日日夜夜操网爽| 亚洲片人在线观看| 久久久久国产一级毛片高清牌| 在线观看www视频免费| 嫩草影院精品99| 久久人人爽av亚洲精品天堂| 亚洲片人在线观看| 欧美色欧美亚洲另类二区 | av福利片在线| 中文亚洲av片在线观看爽| 久久人妻福利社区极品人妻图片| 国产激情欧美一区二区| 国产精品二区激情视频| 天堂影院成人在线观看| 国产高清激情床上av| 国内精品久久久久久久电影| 国产欧美日韩一区二区三区在线| 一本久久中文字幕| 欧美精品亚洲一区二区| 午夜亚洲福利在线播放| 国产亚洲精品综合一区在线观看 | 巨乳人妻的诱惑在线观看| 国产熟女午夜一区二区三区| 一二三四在线观看免费中文在| 久久久精品欧美日韩精品| 高潮久久久久久久久久久不卡| 国产一区在线观看成人免费| 成人国语在线视频| 欧美+亚洲+日韩+国产| 青草久久国产| 日韩精品中文字幕看吧| 亚洲熟妇中文字幕五十中出| 大陆偷拍与自拍| 亚洲中文av在线| 久久久久国内视频| 黄色毛片三级朝国网站| 久久精品91蜜桃| 天天一区二区日本电影三级 | 久久久久久久午夜电影| 亚洲人成伊人成综合网2020| 老司机深夜福利视频在线观看| 最好的美女福利视频网| 国产精品久久久久久人妻精品电影| 午夜福利成人在线免费观看| 一区二区三区激情视频| 在线天堂中文资源库| 国产精品秋霞免费鲁丝片| av电影中文网址| 欧美av亚洲av综合av国产av| 久久精品91无色码中文字幕| 中文字幕久久专区| 91成人精品电影| 最近最新中文字幕大全免费视频| 国产av一区在线观看免费| 老司机靠b影院| 久久人妻熟女aⅴ| 精品无人区乱码1区二区| 国产人伦9x9x在线观看| 一区二区三区高清视频在线| 久久香蕉国产精品| 很黄的视频免费| 欧美日韩瑟瑟在线播放| 婷婷六月久久综合丁香| 一级作爱视频免费观看| 国产熟女xx| 欧美在线黄色| 这个男人来自地球电影免费观看| 中文字幕最新亚洲高清| 国产午夜福利久久久久久| 国产欧美日韩精品亚洲av| 免费不卡黄色视频| 一边摸一边抽搐一进一出视频| 欧美午夜高清在线| 一进一出抽搐动态| 97人妻天天添夜夜摸| 精品国产美女av久久久久小说| 国产麻豆69| 999久久久精品免费观看国产| 高清黄色对白视频在线免费看| 亚洲av成人一区二区三| 露出奶头的视频| 免费一级毛片在线播放高清视频 | 中文字幕高清在线视频| 啦啦啦 在线观看视频| 怎么达到女性高潮| 黄色视频不卡| 日日夜夜操网爽| 深夜精品福利| svipshipincom国产片| 午夜福利欧美成人| 亚洲av电影不卡..在线观看| 国产主播在线观看一区二区| 久久久久精品国产欧美久久久| 国产真人三级小视频在线观看| 午夜老司机福利片| 亚洲精品一区av在线观看| 亚洲人成伊人成综合网2020| 国产男靠女视频免费网站| 国内精品久久久久久久电影| 国产亚洲精品第一综合不卡| 老司机午夜福利在线观看视频| 波多野结衣高清无吗| 国产精品美女特级片免费视频播放器 | 久久久久久大精品| 一级毛片高清免费大全| 久久久久久国产a免费观看| 99精品在免费线老司机午夜| 久久久国产精品麻豆| 淫妇啪啪啪对白视频| 深夜精品福利| 欧美一级a爱片免费观看看 | 欧美在线黄色| √禁漫天堂资源中文www| 日韩成人在线观看一区二区三区| 日韩一卡2卡3卡4卡2021年| www.自偷自拍.com| 精品一区二区三区四区五区乱码| 伊人久久大香线蕉亚洲五| 国产欧美日韩综合在线一区二区| 国产精品综合久久久久久久免费 | av免费在线观看网站| 国产精品 国内视频| 欧美国产精品va在线观看不卡| 国产一区二区三区在线臀色熟女| 色综合婷婷激情| 久久久久久大精品| 在线观看舔阴道视频| 亚洲精品国产精品久久久不卡| 精品一区二区三区av网在线观看| 美女高潮喷水抽搐中文字幕| 国产激情欧美一区二区| 国产一区二区三区在线臀色熟女| 欧美乱码精品一区二区三区| 亚洲七黄色美女视频| 两性午夜刺激爽爽歪歪视频在线观看 | 此物有八面人人有两片| 日本撒尿小便嘘嘘汇集6| 国产精品久久久久久亚洲av鲁大| 在线观看免费午夜福利视频| 韩国精品一区二区三区| 操美女的视频在线观看| 中文字幕人成人乱码亚洲影| 在线观看免费午夜福利视频| 欧美精品亚洲一区二区| 又黄又粗又硬又大视频| 99国产精品99久久久久| 天天躁狠狠躁夜夜躁狠狠躁| 久久精品91无色码中文字幕| 国产熟女xx| 精品乱码久久久久久99久播| 精品日产1卡2卡| 黄色 视频免费看| 黄频高清免费视频| 久久久久久久久免费视频了| 国产精品香港三级国产av潘金莲| 一级作爱视频免费观看| 国产乱人伦免费视频| 淫秽高清视频在线观看| 午夜久久久久精精品| 亚洲国产精品999在线| 久久香蕉激情| 97人妻天天添夜夜摸| 91成年电影在线观看| 国产精品一区二区免费欧美| 免费少妇av软件| 久久久精品国产亚洲av高清涩受| 老司机在亚洲福利影院| 国产精品精品国产色婷婷| 啦啦啦韩国在线观看视频| 国产精品二区激情视频| 99国产精品99久久久久| 午夜亚洲福利在线播放| 老熟妇乱子伦视频在线观看| 黄片大片在线免费观看| 亚洲成人久久性| 国产成人av教育| 久久人人精品亚洲av| 亚洲av电影不卡..在线观看| 99热只有精品国产| 丁香六月欧美| 一夜夜www| av免费在线观看网站| 久久性视频一级片| 亚洲av成人不卡在线观看播放网| 国产成人av教育| 在线视频色国产色| 国产片内射在线| 在线免费观看的www视频| 久久热在线av| 亚洲熟妇熟女久久| 香蕉久久夜色| 啪啪无遮挡十八禁网站| 亚洲一区二区三区色噜噜| 黄色片一级片一级黄色片| 真人一进一出gif抽搐免费| 亚洲欧美精品综合一区二区三区| 欧美日韩乱码在线| 国语自产精品视频在线第100页| 国产精品久久视频播放| 午夜视频精品福利| 亚洲一区高清亚洲精品| 精品国产美女av久久久久小说| 日韩欧美在线二视频| 男女床上黄色一级片免费看| 黑人操中国人逼视频| 国产精品自产拍在线观看55亚洲| 亚洲国产日韩欧美精品在线观看 | 长腿黑丝高跟| 国产熟女xx| 人成视频在线观看免费观看| 黄色成人免费大全| 欧美日韩福利视频一区二区| 97人妻精品一区二区三区麻豆 | 亚洲欧洲精品一区二区精品久久久| 久热这里只有精品99| 日韩欧美国产在线观看| 十八禁网站免费在线| 韩国av一区二区三区四区| 免费搜索国产男女视频| 老司机靠b影院| 夜夜夜夜夜久久久久| 久热这里只有精品99| 黄色a级毛片大全视频| 免费高清在线观看日韩| 免费在线观看视频国产中文字幕亚洲| 最近最新中文字幕大全免费视频| 黑人巨大精品欧美一区二区mp4| 一本久久中文字幕| 午夜免费成人在线视频| 欧美中文综合在线视频| 啦啦啦 在线观看视频| 久久精品亚洲精品国产色婷小说| 久久久久久久精品吃奶| 国产成人欧美| 97超级碰碰碰精品色视频在线观看| 精品人妻1区二区| 熟女少妇亚洲综合色aaa.| 99精品欧美一区二区三区四区| 中文字幕av电影在线播放| 99国产综合亚洲精品| 欧美 亚洲 国产 日韩一| 一个人观看的视频www高清免费观看 | 国产真人三级小视频在线观看| 亚洲第一电影网av| 国产三级黄色录像| 欧美不卡视频在线免费观看 | 天天添夜夜摸| 在线观看66精品国产| 国产精品美女特级片免费视频播放器 | 久久久国产精品麻豆| 国产精品野战在线观看| 欧美在线一区亚洲| 亚洲国产欧美日韩在线播放| 女人精品久久久久毛片| 久热这里只有精品99| 黄色毛片三级朝国网站| 免费女性裸体啪啪无遮挡网站| 一边摸一边抽搐一进一出视频| 高清毛片免费观看视频网站| 亚洲少妇的诱惑av| 多毛熟女@视频| 人人妻人人爽人人添夜夜欢视频| 国产精品电影一区二区三区| 99国产综合亚洲精品| 成人国产综合亚洲| 91国产中文字幕| 久久人人97超碰香蕉20202| 欧美一级毛片孕妇| 长腿黑丝高跟| 欧美激情久久久久久爽电影 | 妹子高潮喷水视频| 国产熟女午夜一区二区三区| 日本免费一区二区三区高清不卡 | 亚洲国产看品久久| 免费久久久久久久精品成人欧美视频| 在线视频色国产色| 国产aⅴ精品一区二区三区波| 19禁男女啪啪无遮挡网站| 日韩国内少妇激情av| 午夜精品在线福利| 欧美在线黄色| 成人欧美大片| 一级片免费观看大全| 欧美日韩黄片免| 十分钟在线观看高清视频www| 精品人妻在线不人妻| 看免费av毛片| 黄频高清免费视频| 每晚都被弄得嗷嗷叫到高潮| 熟女少妇亚洲综合色aaa.| 久久香蕉国产精品| 成年版毛片免费区| 青草久久国产| 国产又爽黄色视频| 天堂影院成人在线观看| 免费看美女性在线毛片视频| 国产午夜福利久久久久久| 最新在线观看一区二区三区| 18禁国产床啪视频网站| 丝袜美足系列| 国产精品久久视频播放| 国产男靠女视频免费网站| 久久人妻av系列| 亚洲狠狠婷婷综合久久图片| 欧美精品亚洲一区二区| 亚洲人成电影免费在线| 国产精品香港三级国产av潘金莲| 免费在线观看日本一区| 黄色视频,在线免费观看| 黄色成人免费大全| 国内精品久久久久久久电影| 国产一卡二卡三卡精品| 久久精品国产综合久久久| 一区二区三区精品91| 亚洲自拍偷在线| av天堂久久9| 国产成人精品在线电影| 搞女人的毛片| 久9热在线精品视频| 日本五十路高清| 国产精品永久免费网站| 俄罗斯特黄特色一大片| 亚洲精品国产区一区二| 自拍欧美九色日韩亚洲蝌蚪91| 亚洲精品美女久久久久99蜜臀| 欧美乱色亚洲激情| 一本久久中文字幕| 天天躁狠狠躁夜夜躁狠狠躁| 亚洲色图综合在线观看| 大型av网站在线播放| 日本免费一区二区三区高清不卡 | 村上凉子中文字幕在线| 啦啦啦观看免费观看视频高清 | 国产日韩一区二区三区精品不卡| 国产成人精品在线电影| 亚洲欧美日韩另类电影网站| 国产麻豆69| 色综合亚洲欧美另类图片| 精品久久久精品久久久| 国产欧美日韩一区二区精品| 婷婷六月久久综合丁香| 国产亚洲精品一区二区www| 99在线视频只有这里精品首页| 精品久久久久久,| 国产激情欧美一区二区| 这个男人来自地球电影免费观看| 国产精品综合久久久久久久免费 | 999精品在线视频| bbb黄色大片| 男人舔女人下体高潮全视频| 9191精品国产免费久久| 亚洲精品久久国产高清桃花| 国产在线观看jvid| 成人精品一区二区免费| 在线观看免费日韩欧美大片| 欧美av亚洲av综合av国产av| 久久 成人 亚洲| 亚洲黑人精品在线| 涩涩av久久男人的天堂| 99精品在免费线老司机午夜| 久久久久国内视频|