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

    Numerical simulation on dendritic growth of Al–Cu alloy under convection based on the cellular automaton lattice Boltzmann method

    2022-09-24 08:03:54KangWeiWang王康偉MengWuWu吳孟武BingHuiTian田冰輝andShouMeiXiong熊守美
    Chinese Physics B 2022年9期

    Kang-Wei Wang(王康偉) Meng-Wu Wu(吳孟武) Bing-Hui Tian(田冰輝) and Shou-Mei Xiong(熊守美)

    1School of Automotive Engineering,Wuhan University of Technology,Wuhan 430070,China

    2Hubei Key Laboratory of Advanced Technology for Automotive Components,Wuhan University of Technology,Wuhan 430070,China

    3School of Materials Science and Engineering,Tsinghua University,Beijing 100084,China

    Keywords: simulation,cellular automaton,dendritic growth,melt convection

    1. Introduction

    Metal solidification is a complex liquid-solid phase transition process accompanied by solute and temperature changes in both the macro and micro scales.[1,2]Those changes cause natural convection in the solidification process of metal alloys.In addition, there is forced convection induced by external force fields such as electromagnetic stirring or ultrasonic vibration during the metal solidification process.[3-5]The melt flow behavior has a great impact on the dendritic growth,which in turn affects the metal alloy solidification structure and the formation of defects such as solute segregation,shrinkage,and porosity.[6]Therefore,an in-depth understanding and mastery of the dendritic growth law under melt convection could help researchers to give suggestions for controlling the microstructure formation and improving the performance of metal alloys.

    Owing to the rapid development of x-ray radioscopy technique,in situobservation of alloy solidification process becomes possible.[7-12]Shevchenkoet al.[8,9]directly observed the solute segregation during solidification of Ga-In alloys by x-ray radioscopy and studied the effect of natural and forced convection on dendritic growth. Clarkeet al.[12]observed the directional solidification of Al-Cu alloys under different temperature gradients by the synchrotron x-ray technique. They found that melt convection increases the stability of the planner interface at a higher growth velocity. Though the dendritic growth process has already been observed realistically and directly through the experimental method, there are still many limitations. Firstly,the applicable alloy system is limited to a small range. Secondly,it is of a great challenge to construct a precise flow field during metal solidification and observe the dendritic evolution in real time. Finally, there are difficulties in distinguishing the effects of various influencing factors on dendritic growth quantitatively.

    Over the past few decades,numerical simulation has become an effective method to predict microstructure evolution,such as columnar-to-equiaxed transition,[13-15]dendritic motion,[16,17]and porosity formation.[18]The approaches of microstructure simulation include Monte Carlo,[19]phase field,[20,21]and cellular automaton(CA)methods.[22,23]In particular, the CA method can accurately describe the dendrite morphology under the interaction of heat and solute and has the advantage of high computational efficiency with a large calculation domain. Therefore,it has become the mainstream and hotspot in the field of micro-mesoscopic structure simulation. Combining a CA model with momentum and solute transfer,Zhuet al.[24]studied the influence of forced convection on the equiaxed and columnar dendritic morphology of Al-Cu alloy in two dimensions. Subsequently, Sunet al.[25]and Yinet al.[26]developed a cellular automaton-lattice Boltzmann (CA-LB) model, in which the CA method was used to simulate the dendritic growth and the LBM was applied to calculate the flow field, heat and solute transport. They pointed out that the LBM has great advantages, such as high calculation accuracy, superior convergence, and low computational cost. Likewise, Liuet al.[27]used the CA-LBM model to simulate the equiaxed dendritic growth of Al-4.7%Cu alloy and found that the rotation and falling motion significantly affect the dendritic morphology under forced convection. Additionally, Maet al.[28]and Rolchigoet al.[29]applied the CA method to the fields of welding and additive manufacturing,and investigated the influence of process parameters on microstructure evolution in molten pool under forced convection by simulation.

    Up to date, scholars have carried out plenty of research work on microstructure simulation and many valuable conclusions have been drawn in this period. However,existing studies have not perfectly solved the problem of grid anisotropy in the CA method, thus it is difficult to simulate dendritic growth with random orientations. Meanwhile,for the effect of melt convection on dendritic growth,early simulation studies mainly focused on forced convection. In the actual solidification process,there must be more or less temperature or solute concentration difference in the melt. In this case,natural convection is inevitably formed under the action of gravity,and its influence on the solute distribution and dendritic growth cannot be ignored unless its intensity is much less than that of forced convection.

    In this work,a numerical model is developed by coupling the CA method and the LBM to simulate the dendritic growth of metal alloys under convection. An improved decentered square algorithm is embedded into the CA model to eliminate the artificial anisotropy induced by the CA cells. Twoand three-dimensional simulation cases are performed to investigate the effects of both the forced convection and gravitydriven natural convection on the dendritic growth of Al-Cu alloy, based on which the competition growth mechanism of dendrites,the formation of solute plume and interdendritic solute segregation are also discussed.

    2. Model description

    2.1. Microscopic temperature field

    The size unit of the macroscopic grid in the model is mm,but that of the cell in the dendritic growth model isμm. The two orders of magnitude differ greatly. Therefore,the macroscopic data needs to be accurately transmitted to each microscopic cell. In the present work, bilinear interpolation[30]is used to transform the macroscopic temperature data into the microscopic level. Figure 1(a) shows the schematic diagram of the bilinear interpolation method. If the temperature and coordinates of each node in the macroscopic calculation domain are known,the temperature data of the microscopic cell node can be calculated by bilinear interpolation. The calculation method is as follows:

    whereT(P)is the node temperature value of a microscopic element,T(Q) is the node temperature value of a macroscopic element,xandyrepresent the horizontal and vertical coordinates of the node on the macroscopic computing domain, respectively. In the case of three dimensions,select a set of opposite faces,obtain the data of their respective internal nodes through bilinear interpolation, and then interpolate these data again to calculate the temperature value inside the space.

    2.2. The lattice Boltzmann method for the fluid flow and solute transport

    Compared with the traditional computational fluid dynamics (CFD) calculation methods, LBM has the characteristics of a mesoscopic model between the microscopic molecular dynamics model and macroscopic continuum model. It is widely used to simulate fluid flow and solute transport because of its easy handling of complex boundaries,parallel calculation, and simple programming. Based on the LBM, several models have been developed. Among them,the accuracy of the multi-relaxation-time (MRT) model is higher than that of the single-relaxation-time (SRT) model.[31]However, the MRT-LBM model has high model complexity and low computational efficiency, so the SRT-LBM model is preferred in this work. The SRT-LBM equation can be written as[32]

    wherefi(x,t) andfeqi(x,t) are the particle distribution and equilibrium particle distribution functions, respectively;eiis the lattice discrete velocity;τis the relaxation time of the flow field, which is fixed at 1.5 in this work; and Δtis the time step,which can be calculated from the relaxation timeτ. The forcing termFi(x,t)on the right-hand side reads[33]

    wheregis gravitational acceleration,βCis solute expansion coefficient,C0is the reference solute concentration.

    The LBM is similarly applied to the calculation of solute diffusion in the liquid phase region,and its processing method is the same as that of fluid flow calculation. The governing equation of the solute distribution functiongi(x,t)can be written as[34]

    where Δfsis solid fraction increment,kis solute partition coefficient,andCis solute concentration.

    The D2Q9 and D3Q15 models are used to simulate twodimensional and three-dimensional cases, respectively.[35-37]The lattice discrete velocities of the D2Q9 model read

    In the D3Q15 model, the lattice discrete velocities read[38]

    wherecis the lattice velocity, whose value is the ratio of the space step Δxto the time step Δt. The macro parameters are obtained by the corresponding lattice distribution function:

    According to the Chapman-Enskog analysis,[34]the relaxation time of flow fieldτand solute concentration fieldτDcan be calculated by the kinematic viscosityvand solute diffusion coefficientD,respectively,

    The unknown particle distribution function on the boundary node needs to be obtained under boundary conditions, so whether the boundary conditions are reasonable has a significant impact on the simulation.[39]In this study, the nonequilibrium extrapolation method is applied to the inlet and outlet boundaries in fluid flow calculation, while the other boundaries are set as periodic boundaries.[40]For the solute field, all boundaries of the computational domain are set as non-diffusion boundaries, and non-equilibrium extrapolation is used in the calculation. Figures 1(b) and 1(c) show the schematic diagrams related to the boundary conditions under forced and natural convection,respectively. In addition,since the melt cannot pass through the dendrites,the solid-liquid interface is considered to be the no-slip and no-diffusion boundary using the bounce-back rule for flow field and solute field simulation.[41]

    2.3. The cellular automaton model for dendritic growth

    Since simulations on dendritic growth without melt flow based on the CA method were reported in our previous work,the governing equations required for the dendritic growth kinetics are not introduced in this paper. More details can be seen in Ref.[42]. This section focuses on introducing a modified capture rule to simulate dendritic growth with arbitrary orientations.

    The growth orientation of dendrites is generally random during metal alloy solidification. However,due to the artificial anisotropy induced by the CA cells, the simulated dendrites usually grow in the orientations of 0°and 45°regardless of the initial orientation value defined in simulations. An improved decentered square algorithm is adopted in this study to eliminate this artificial anisotropy.The decentered square algorithm was first proposed by Rappaz and Gandin,[43]and then modified by Wanget al.[44]Recently, Zhuet al.[45,46]introduced a geometric factor into the algorithm to ensure the sharp interface of the simulated dendrites. The decentered square algorithm is currently one of the ideal methods for describing the dendritic growth with arbitrary orientation. Figure 1(d)shows the schematic diagram of the decentered square algorithm. When a cell nucleates, its state is transformed from liquid to solid/liquid interface. A father square with the diagonal direction consistent with the crystal orientation is placed in the center. With the dendritic growth driven by the interaction between undercooling and solute diffusion, the solid phase fraction increases,resulting in continuous expansion of the square. When the vertices of the inclined square extend to the surrounding liquid cell, the state of the captured cell is transformed into solid/liquid interface. At the same time,four sub-squares inheriting the growth orientation of the father square are generated with the vertices of the father square as the center.When the solid fraction of the interface cell reaches 1, its state changes to solid, and the corresponding inclined square stops expanding. Additionally, if there are still liquid cells around the solid cell,the liquid cells change to interface cells automatically. The half length of the square diagonalLcan be updated as follows:

    Fig. 1. Schematic diagrams of the issues considered in the model: (a) bilinear interpolation method, (b) boundary conditions under forced convection,(c)boundary conditions under natural convection,(d)decentered square algorithm.

    where Δxis the cell size,the preferential growth orientationθis chosen within-45°to 45°.

    3. Model validation

    The Al-Cu alloy is selected in this work,while the physical parameters required for simulation are listed in Table 1.The melt is supposed to be the incompressible Newtonian fluid under convection. In this section,the side length of the square computational domain is 1 mm. A solid seed is planted in the center of the computational domain and its growth orientation is predefined as 0°, 15°, 30°and 45°, respectively. A cooling rate of 60 K/s is set to the melt of the whole calculation domain.

    For numerical modeling and simulation,grid density has a great impact on the final results. Refining the mesh improves the accuracy of the results, but drastically reduces the computational efficiency. Therefore, it is crucial to select an appropriate number of grids under the premise of balancing computational efficiency and result accuracy. Figure 2 shows the simulated equiaxed dendritic morphology with different growth orientations under different mesh sizes (1 μm and 2μm). Comparing the simulation results, we can realize that the scales of the simulated equiaxed dendrites are basically the same,indicating that the algorithm has already effectively eliminated the impact of grid anisotropy and the present numerical model has the capability and stability to simulate dendritic growth with arbitrary orientations. In addition, the mesh size used in this work has a negligible effect on dendrite growth kinetics, but more details of dendrites, such as secondary and higher-order dendritic arms,can be depicted with a finer mesh size.

    Table 1. Physical parameters of Al-Cu alloy used in the present study.[47,48]

    Fig.2. Simulated equiaxed dendritic morphology with different growth orientations: [(a), (e)]θ =0°, [(b), (f)]θ =15°, [(c), (g)]θ =30°.[(d),(h)]θ =45°. The upper and lower rows show the equiaxed dendritic morphology with mesh sizes of 1μm and 2μm,respectively.

    Fig.3. Simulated equiaxed dendritic morphology with different orientations under forced convection: (a)θ =0°,(b)θ =15°,(c)θ =30°,(d)θ =45°.The forced convection is constructed in the form of a horizontal melt flow from left to right of the calculation domain with a velocity of 0.4 mm/s.

    To verify the feasibility of the established model in simulating dendritic growth under convection,a forced convection is constructed in the form of a horizontal melt flow from left to right of the calculation domain with a velocity of 0.4 mm/s on the basis of the examples in Figs. 2(e)-2(h), and the simulation results are shown in Figs. 3(a)-3(d). It can be seen that the melt flow changes the solute distribution around the equiaxed dendrite,which makes the local undercooling at the dendrite tips different from each other,thereby destroying the symmetry of dendritic morphology. Figure 4 shows the statistical result related to the length of the upstream and downstream primary dendrite trunks under forced convection as indicated in Figs. 3(a)-3(d). The result demonstrates that due to the strong convective transport of solute from the upstream side to the downstream side, the forced convection promotes the growth of dendrite arms at the upstream side and suppresses the growth of downstream dendrite arms. Meanwhile,with the increase of the dendritic orientation angle from 0°to 45°, the four primary dendrite trunks are no longer parallel or perpendicular to the melt flow direction,in the case the length difference between the upstream and downstream primary dendrite trunks gradually decreases.The laws concluded from the simulation results are in accordance with the reality and previously published works. It is confirmed that the CALBM model has the capability of coupling the calculation of temperature field,solute field and flow field.

    Fig. 4. The length of the upstream and downstream primary dendrite trunks under forced convection with a lateral flow velocity of 0.4 mm s as indicated in Figs.3(a)-3(d).

    4. Simulation results and discussion

    4.1. Two-dimensional dendritic growth under forced convection

    Simulation cases are performed by the CA-LBM model to study the multiple equiaxed dendritic growths under forced convection. The computational domain is divided into 600×300 square cells with a cell size of 2μm. To deeply explore the interaction among melt flow, solute distribution and dendritic growth, seeds with random orientation and location are planted in the middle part of the calculation domain,and a lateral melt flow is defined from left to right with a flow velocity of 0.4 mm/s. The cooling rate of the melt in the calculation domain is 60 K/s. Figures 5 and 6 show the simulated multiple equiaxed dendritic growth under forced convection with ten and fifty seeds,respectively.

    It can be noted from Figs. 5(a)-5(d) that the melt flows smoothly through the dendrites at the early solidification stage,bringing the solute from the upstream side of each dendrite to the downstream side, so that the development of the primary and secondary dendrite arms is promoted at the upstream side,while their growth is inhibited on the downstream side. As the solidification proceeds, the dendrites on the left grow rapidly under forced convection and hinder the melt flowing to the right. As indicated in Figs.5(d)and 6(f),the flow velocity in the interdendritic liquid and the right part of the computational domain decreases with the increase of the solid fraction until it approaches zero. In the subsequent stage, the dendrites on the right side grow freely without the action of forced convection. However,the dendrites on the left side are generally still more developed than those on the right side. It can be seen from Fig. 6 that the increase of the seed number makes the interdendritic liquid channels narrower and more complicate.In this case, the blocking effect of dendrites on melt flow is advanced. As shown in Figs. 6(f) and 6(h), the melt flow is already unable to traverse equiaxed dendritic clusters.

    Fig. 5. Simulated multiple equiaxed dendritic growth under forced convection with ten seeds: [(a), (b)]t =0.20 s, [(c), (d)]t =0.45 s, [(e), (f)]t=0.70 s,[(g),(h)]t=0.95 s.The left column shows the solute distribution while the right column shows the flow velocity distribution of the melt.

    Fig. 6. Simulated multiple equiaxed dendritic growth under forced convection with fifty seeds: [(a), (b)] t = 0.20 s, [(c), (d)] t = 0.45 s, [(e), (f)]t=0.70 s,[(g),(h)]t=0.95 s. The left column shows the solute distribution while the right column shows the flow velocity distribution of the melt.

    According to the simulated results in Figs. 5 and 6, it can be found that the dendritic morphology is quite different at different parts of the calculation domain. The dendrites on the left exhibit a highly asymmetric morphology due to a large space offered on the left side of the calculation domain and the action of forced convection. Competitive growth among these dendrites is related to their growth orientations and forced convection, showing a trend growing towards the upstream side. However, the effect of forced convection on the dendritic growth in the middle part and on the right side of the calculation domain is almost negligible. With the obstruction of dendrites on both sides,the dendrites in the central part show a refined and equiaxed dendritic morphology. Similarly,due to a large space offered on the right side of the calculation domain, the dendrites on the right show an asymmetric morphology with the competitive growth between them.

    Figure 7 shows the simulated columnar dendritic growth with and without forced convection. The computational domain is composed of 400×400 square cells with a cell size of 2μm. Nine seeds with random orientation and location are set at the bottom of the calculation domain. With a temperature gradient of 10 K/mm from bottom to top of the calculation domain, the cooling rate of the melt is 60 K/s. It can be noted from Figs.7(a)-7(c)that the competitive growth occurs when columnar dendrites with different growth orientations meet with each other,and the columnar dendrites whose growth orientation is more aligned to the heat flow direction always overcome those dendrites with their growth orientations deviating from the heat flow direction. The columnar dendrites with favorable orientations show a strong growth advantage. Side branches appear on the side of the columnar dendrite trunks,and their growth also reflects a competitive growth trend as the same as that of the columnar dendrite trunks. Additionally,the solute is highly accumulated in the interdendritic liquid.Compared to the case under pure diffusion, the presence of forced convection significantly affects the competitive growth between the columnar dendrites. This can be elucidated as shown in Figs. 7(c) and 7(f) taking the dendrites A and B as examples. It can be seen from Fig.7(f)that the columnar dendrite A with a more favorable growth orientation is blocked by the columnar dendrite B.This is due to the fact that the growth of the columnar dendrite B is promoted under forced convection since its growth direction is more aligned to the upstream side of the melt flow. Consequently,the growth of the columnar dendrite B is dominant and gradually hinders the growth of the columnar dendrite A.

    Fig.7. Simulated columnar dendritic growth under pure diffusion(a)-(c)and forced convection(d)-(f)in the form of a lateral flow with a flow velocity of 0.2 mm/s: [(a),(d)]t=0.25 s,[(b),(e)]t=0.6 s,[(c),(f)]t=0.95 s.

    Simulation cases are conducted to further study the overgrowth mechanism of columnar dendrites under pure diffusion and forced convection. The computational domain consists of 400×300 square cells with a cell size of 2μm. The temperature gradient and cooling rate of the melt are set the same as those in the simulation cases in Fig. 7. Differently, ten seeds are planted at the bottom of the calculation domain with an equal interval, while the first and second halves of the dendrites are named favorably oriented (FO) dendrites and unfavorably oriented(UO)dendrites,respectively according to the heat flow direction. Figure 8 shows the simulated competitive growth of converging columnar dendrites. It can be seen from Figs. 8(a)-8(c) that the FO columnar dendrites always overgrow the UO columnar dendrites during pure diffusioncontrolled solidification, which is consistent with the classical dendritic competitive growth theory proposed by Walton and Chalmers.[49]This overgrowth phenomenon usually occurs because of the overlap of the solute fields between the contacting columnar dendrites. Moreover,the elimination rate of the UO columnar dendrites increases as their growth orientation gradually deviates from the heat flow direction. An unusual competition mechanism appears between converging columnar dendrites with the existence of forced convection.Due to the promotion effect of melt flow on the upstream columnar dendrite arms,the growth of the FO columnar dendrites is hindered by the UO columnar dendrites as shown in Fig. 8(d). As the deviation angle between the growth orientation of the UO columnar dendrites and the heat flow direction increases from 12°to 20°as indicated in Fig. 8(e), the UO columnar dendrites cannot completely block the growth of the FO columnar dendrites. When this deviation angle continuously increases to a specific value such as 30°, the promotion effect of forced convection on the growth of the UO columnar dendrites cannot change their growth weakness,and the growth of those dendrites is inhibited by the FO columnar dendrites as shown in Fig.8(f). That is to say,the competitive growth of the converging columnar dendrites is determined by the interaction between heat flow and forced convection. The simulated results are in accordance with the conclusions drawn by Jaehoonet al.[50]and Pavanet al.[51]

    Fig.8. Simulated competitive growth of converging columnar dendrites under pure diffusion(a)-(c)and forced convection(d)-(f)in the form of a lateral flow with a flow velocity of 0.2 mm/s: [(a), (d)] θFO =10°, θUO =12°, [(b), (e)] θFO =10°, θUO =20°, [(c), (f)] θFO =10°,θUO=30°. Here θFO and θUO are the angles between the columnar dendrite trunks and vertical upward direction while they are specified to be positive clockwise.

    4.2. Two-dimensional dendritic growth under natural convection

    Figure 9 shows the equiaxed dendritic growth in the presence of natural convection, while the cooling rate of the melt is 60 K/s in the whole calculation domain. It can be seen that the natural convection induced by the sinking downwards the rejected heavy Cu solute at the solid-liquid interface changes the solute distribution around the dendrite,thereby destroying the four-fold symmetry of the equiaxed dendritic morphology.This is because with the downward sinking of the solute,it accumulates around the dendrite tips at the downstream side.For the Al-Cu alloy,the solute enrichment will lower the liquidus temperature and consequently reduce the local undercooling of the melt. In this case,the growth of primary dendrite trunk at the downstream side is inhibited, while the growth of the upstream primary dendrite trunk is accelerated to a certain extent.As the solidification proceeds,the heavy solute continues to sink downwards, leading to the formation of solute plume in front of the primary dendrite trunk on the downstream side.Under natural convection as marked by the dashed lines with arrows in Figs.9(b)-9(d),it can be seen that two vortices with opposite rotation directions are formed,which directly have an important influence on the solute distribution.

    Except the highly asymmetric growth of the four primary dendrite trunks, the natural convection also leads to dissimilar growth condition of the secondary dendrite arms as indicated in Fig. 9(d). There are well-developed secondary dendrite arms branching on the upward primary dendrite trunk,while the growth of the secondary dendrite arms branching on the upstream side of the two horizontal primary dendrite trunks is severely inhibited due to the enrichment of the heavy solute.Quite unexpectedly,the secondary dendrite arms branching on the downstream side of the two horizontal primary dendrite trunks grow well,and some of them even grow faster than the downward primary dendrite trunk. This can be ascribed to a severe solute enrichment in front of the downward primary dendrite trunk,while the solute enrichment in front of the secondary dendrite arms is weakened by the two vortices with a strong solute transport capability.

    Fig.9. Simulated equiaxed dendritic growth under natural convection:(a)t=0.08 s,(b)t=0.28 s,(c)t=0.48 s,(d)t=0.68 s.

    Compared with the equiaxed dendritic growth, the gravity-driven natural convection may have a greater effect on the growth of columnar dendrites when they grow parallel to the gravity. Figures 10(a)and 10(b)show the simulated columnar dendritic growth under pure diffusion and gravitydriven natural convection, respectively, with a temperature gradient of 10 K/mm from bottom to top of the calculation domain and a cooling rate of 60 K/s in the melt. With the same calculation time, it can be seen from Fig. 10 that the columnar dendrites under natural convection grow faster than those under pure diffusion. This is because the solute accumulation at the dendrite tips is weakened due to the downward sinking of the heavy solute, which is beneficial to the growth of the columnar dendrite tips. With the sinking downwards of the solute, it can be speculated that the degree of solute segregation increases from top to bottom in the liquid regions between the columnar dendrites.

    Fig.10. Simulated columnar dendritic morphology with the same calculation time t=0.75 s: (a)pure diffusion,(b)natural convection. The columnar dendrites grow in the opposite direction of gravity.

    Keeping the gravitational direction vertical downward,the same number of seeds as in the simulation cases in Fig.10(ten seeds) are planted at the top of the calculation domain.With a cooling rate of 60 K/s in the melt, a temperature gradient of 10 K/mm is set from top to bottom of the calculation domain to ensure that the columnar dendrites grow along the gravitational direction. Figure 11 shows the simulated columnar dendritic growth under natural convection. It can be noted that the rejected heavy solute sinks freely into the liquid phase region,leading to solute accumulation in front of the dendrite tips. With the solute convergence continuously, some solute gradually sinks away from the solute pool and several bulges are formed in front of the solute streams. The expansion of the bulges eventually leads to the formation of chimney-like or mushroom-like solute plumes as illustrated in Fig. 11(c).Compared with the columnar dendritic growth in Fig. 10, it can be seen from Fig.11 that due to a less solute accumulation in the interdendritic liquid,secondary dendrite arms in the direction of gravity are well-developed, resulting in a decrease of the columnar dendritic arm spacing.

    Fig.11. Simulated columnar dendritic growth under natural convection: (a)t =0.25 s,(b)t =0.50 s,(c)t =0.75 s. The columnar dendrites grow along the gravitational direction.

    4.3. Three-dimensional dendritic growth under forced convection

    During metal alloy solidification process, the dendrites actually grow in three dimensions, and the melt flow is also in three dimensions. It can be easily expected that the 3D dendritic growth is more complicate than that in the 2D case,and the 3D simulated result of dendritic growth is much closer to the reality. In this section, simulation cases are performed to investigate the 3D dendritic growth of Al-Cu alloy under forced convection.

    Fig. 12. Three-dimensional equiaxed dendritic growth under forced convection with different inflow velocities: (a)0 mm/s,(b)0.15 mm/s,(c) 0.30 mm/s, (d) 0.45 mm/s. The melt flow direction is from left to right of the calculation domain.

    For the 3D equiaxed dendritic growth, a solid seed is planted in the middle of the calculation domain consisting of 150×150×150 cubic cells with a cell size of 5μm. The temperature of the melt in the calculation domain is assumed to be uniform and it cools down from the liquidus temperature at a constant cooling rate of 10 K/s.Figure 12 shows the simulated 3D equiaxed dendritic growth under forced convection with different inflow velocities. It can be seen that the equiaxed dendritic morphology is symmetrical in three dimensions under pure diffusion. However, this symmetrical feature is broken when forced convection is introduced into the calculation procedure. The asymmetric degree increases with a larger inflow velocity. The 3D simulated result is consistent with that in the 2D case as indicated in Fig. 3. Due to the strong convective transport of solute from the upstream side to the downstream side,the forced convection promotes the growth of the primary dendrite trunk at the upstream side and suppresses the growth of downstream primary dendrite trunk. Here, the solute transport effect of the forced convection can be elucidated quantitatively as shown in Fig.13. Under pure diffusion condition, the solute concentration gradients in front of the dendrite tips on both sides are equal with the same thickness of solute accumulation layer. As the inflow velocity increases,the thickness of the solute accumulation layer in front of the dendrite tip at the upstream side decreases, while it becomes thicker in front of the dendrite tip at the downstream side. As a result,the upstream primary dendrite trunk grows faster than that at the downstream side. It can also be found intuitively from Fig.12 that with the increase of the inflow velocity, the secondary dendrite arms branching on the two vertical primary dendrite trunks grow well on the upstream side,and they gradually occupy a dominant position in the competition growth with the secondary dendrite arms branching on the horizontal primary dendrite trunk at the upstream side.

    Fig.13. Solute distribution along the centerline from left to right of the calculation domain as indicated in Fig.12.

    For the 3D columnar dendritic growth,the calculation domain is divided into 150×40×200 cubic cells with a cell size of 5 μm. Three seeds are planted at the bottom of the calculation domain with a uniform interval. With a temperature gradient of 10 K/mm from bottom to top of the calculation domain, the cooling rate of the melt is set to be 10 K/s. Figure 14 shows the simulated 3D columnar dendritic growth with and without forced convection. Under pure diffusion condition, the columnar dendrites grow straightly from bottom to top of the calculation domain,accompanied by solute rejection at the solid-liquid interface. With a similar solute distribution in front of the columnar dendrite tips,the three columnar dendrites grow shoulder to shoulder as shown in Figs. 14(a) and 14(d). When forced convection is introduced, a conspicuous phenomenon can be observed in Figs.14(b),14(c)and 14(e),14(f)that the upstream secondary dendrite arms branching on the three columnar dendrite trunks are more well-developed than those on the downstream side. Meanwhile,the growth of the three columnar dendrites from left to right is gradually accelerated. Here,the liquid region between the left and middle columnar dendrites is defined as regionL,and the liquid region between the middle and right columnar dendrites is defined as regionR.It can be noted from Fig.14(f)that the solute concentration in regionLis significantly higher than that in regionR.This can be attributed to the blocking effect of columnar dendrites on melt flow. Since the left columnar dendrite is close to the inlet of melt flow,the intensity of melt flow is relatively strong, bringing the solute at the upstream side to regionL.However,due to the obstruction of the middle columnar dendrite, the melt flow intensity is weakened to a certain extent.In such a case,the solute in regionLno longer easily diffuses into regionR.

    Fig. 14. Three-dimensional columnar dendritic growth under forced convection with different inflow velocities: [(a), (d)] 0 mm/s, [(b), (e)]0.5 mm/s,[(c),(f)]1.0 mm/s. The melt flow direction is from left to right of the calculation domain. (a)-(c)The simulated columnar dendritic morphology. (d)-(f)The corresponding solute distribution at the middle cross section parallel to the front view.

    5. Conclusions

    In this study,a CA-LBM model is developed to simulate the equiaxed and columnar dendritic growth of Al-Cu alloy in both two and three dimensions. The effects of forced and natural convection on the dendritic growth are studied. The following conclusions can be drawn:

    (i) By introducing an improved decentered square algorithm into the CA method to eliminate the artificial anisotropy induced by CA cells,the present numerical model has the capability of simulating dendritic growth with arbitrary orientations.

    (ii) During multiple equiaxed dendritic growth under forced convection, the increase of solid fraction makes the interdendritic liquid channels narrower and more complicate.The blocking effect of dendrites on melt flow is advanced with a larger number of seeds. The growth of columnar dendrites is promoted when their growth orientation is aligned to the heat flow direction and upstream side of forced convection. The competitive growth of the converging columnar dendrites is determined by the interaction between heat flow and forced convection.

    (iii)Gravity-driven natural convection leads to the highly asymmetric growth of equiaxed dendrites. With the sinking downwards of the heavy Cu solute, severe solute segregation forms in the interdendritic liquid when the columnar dendrites grow in the opposite direction of gravity. Conversely,chimney-like or mushroom-like solute plumes are formed in the melt in front of the columnar dendrites when they grow along the gravitational direction.

    (iv) The present numerical model has the capability of simulating dendritic growth of Al-Cu alloy under convection in three dimensions.The 3D simulated result is consistent with that in the 2D case. More details on dendritic growth are revealed by 3D simulations.

    Acknowledgements

    The authors thank Dr. Zhipeng Guo at Tsinghua University (China) and Dr. Ang Zhang at Chongqing University(China)for their helpful discussion.

    Project supported by the National Natural Science Foundation of China (Grant No. 51805389), the Key R&D Program of Hubei Province, China (Grant No. 2021BAA048),the 111 Project (Grant No. B17034) and the fund of Hubei Key Laboratory of Advanced Technology for Automotive Components,Wuhan University of Technology(Grant No.XDQCKF2021011).

    少妇高潮的动态图| 波多野结衣巨乳人妻| 成年版毛片免费区| 亚洲乱码一区二区免费版| 九色成人免费人妻av| 色播亚洲综合网| 亚洲精品456在线播放app | 中亚洲国语对白在线视频| 又爽又黄无遮挡网站| 亚洲av中文字字幕乱码综合| 97碰自拍视频| 岛国视频午夜一区免费看| 亚洲aⅴ乱码一区二区在线播放| 午夜激情福利司机影院| 免费av观看视频| 国产精品久久久久久久电影 | 看免费av毛片| 色哟哟哟哟哟哟| 久久久久久久久大av| 免费人成视频x8x8入口观看| 丰满人妻熟妇乱又伦精品不卡| 内射极品少妇av片p| 最好的美女福利视频网| 亚洲av一区综合| 免费在线观看亚洲国产| 亚洲一区二区三区不卡视频| 美女大奶头视频| 99国产综合亚洲精品| 成人性生交大片免费视频hd| 色综合婷婷激情| 在线观看av片永久免费下载| 欧美日韩中文字幕国产精品一区二区三区| 精品熟女少妇八av免费久了| 18+在线观看网站| 久久精品人妻少妇| 听说在线观看完整版免费高清| 国产美女午夜福利| 国产av麻豆久久久久久久| 国产极品精品免费视频能看的| 91麻豆精品激情在线观看国产| 中亚洲国语对白在线视频| 蜜桃久久精品国产亚洲av| 岛国在线观看网站| 丁香欧美五月| 一个人观看的视频www高清免费观看| 成年女人毛片免费观看观看9| 国产伦在线观看视频一区| 亚洲精品在线观看二区| 一个人免费在线观看电影| 夜夜夜夜夜久久久久| 欧美高清成人免费视频www| 日本免费一区二区三区高清不卡| 亚洲国产欧美人成| 久久久国产成人精品二区| 国产熟女xx| 午夜福利在线在线| 一夜夜www| 成年女人毛片免费观看观看9| 久久精品国产亚洲av涩爱 | 日本a在线网址| 国产单亲对白刺激| 一个人看的www免费观看视频| 少妇丰满av| 欧美zozozo另类| 欧美另类亚洲清纯唯美| 亚洲中文字幕日韩| 国产一区二区激情短视频| 99精品久久久久人妻精品| 亚洲在线观看片| 国产一级毛片七仙女欲春2| 99热只有精品国产| 禁无遮挡网站| 在线国产一区二区在线| 午夜福利欧美成人| 18禁裸乳无遮挡免费网站照片| 国产av一区在线观看免费| 亚洲专区国产一区二区| 日本精品一区二区三区蜜桃| 中文字幕精品亚洲无线码一区| 午夜日韩欧美国产| 国产一区在线观看成人免费| 69av精品久久久久久| 成熟少妇高潮喷水视频| 欧美区成人在线视频| 深夜精品福利| 欧美色欧美亚洲另类二区| 在线观看舔阴道视频| 亚洲美女黄片视频| 亚洲精品久久国产高清桃花| 国产欧美日韩精品一区二区| 久久国产乱子伦精品免费另类| 香蕉丝袜av| 亚洲精品456在线播放app | 国产伦一二天堂av在线观看| 一级a爱片免费观看的视频| 成人午夜高清在线视频| 国产综合懂色| 熟女电影av网| 国产一区在线观看成人免费| www.999成人在线观看| 日本免费a在线| 18禁美女被吸乳视频| 亚洲黑人精品在线| 两个人看的免费小视频| 不卡一级毛片| 欧美一区二区亚洲| 一卡2卡三卡四卡精品乱码亚洲| 国语自产精品视频在线第100页| 高清日韩中文字幕在线| 99精品久久久久人妻精品| 中文字幕av在线有码专区| 国产男靠女视频免费网站| 男女之事视频高清在线观看| 欧洲精品卡2卡3卡4卡5卡区| 99久久久亚洲精品蜜臀av| 久久午夜亚洲精品久久| 午夜免费男女啪啪视频观看 | 色尼玛亚洲综合影院| 全区人妻精品视频| 亚洲精品国产精品久久久不卡| 51国产日韩欧美| 国产伦一二天堂av在线观看| 久久久成人免费电影| 午夜福利视频1000在线观看| 亚洲国产日韩欧美精品在线观看 | 真实男女啪啪啪动态图| 每晚都被弄得嗷嗷叫到高潮| 亚洲狠狠婷婷综合久久图片| 亚洲第一欧美日韩一区二区三区| 看免费av毛片| 免费在线观看影片大全网站| 精品乱码久久久久久99久播| 日韩亚洲欧美综合| 高潮久久久久久久久久久不卡| 国产一级毛片七仙女欲春2| 亚洲美女视频黄频| 我要搜黄色片| 国产精品一区二区免费欧美| 亚洲乱码一区二区免费版| 丝袜美腿在线中文| 国产激情偷乱视频一区二区| 日本 av在线| 日本一二三区视频观看| 看片在线看免费视频| 久久人人精品亚洲av| 母亲3免费完整高清在线观看| 黑人欧美特级aaaaaa片| 波多野结衣高清无吗| 搡老妇女老女人老熟妇| 国产成人av激情在线播放| 大型黄色视频在线免费观看| 男女视频在线观看网站免费| 可以在线观看毛片的网站| 亚洲真实伦在线观看| 欧美激情在线99| 观看免费一级毛片| 久久欧美精品欧美久久欧美| 欧美xxxx黑人xx丫x性爽| 国产探花在线观看一区二区| 国产一区在线观看成人免费| av女优亚洲男人天堂| 91在线观看av| 在线免费观看的www视频| 操出白浆在线播放| 亚洲欧美日韩东京热| 高清日韩中文字幕在线| 99热6这里只有精品| 少妇熟女aⅴ在线视频| 亚洲,欧美精品.| 热99re8久久精品国产| 成人高潮视频无遮挡免费网站| 久久久久久久久大av| 国产美女午夜福利| 国内精品一区二区在线观看| 国产精品三级大全| 久久伊人香网站| 12—13女人毛片做爰片一| 在线观看美女被高潮喷水网站 | 成人欧美大片| 热99在线观看视频| 日本撒尿小便嘘嘘汇集6| 国产精品98久久久久久宅男小说| 亚洲av电影在线进入| 校园春色视频在线观看| 国产三级中文精品| 欧美黑人巨大hd| 可以在线观看毛片的网站| 美女免费视频网站| av女优亚洲男人天堂| 国产精品乱码一区二三区的特点| 国产精品一区二区三区四区久久| 夜夜夜夜夜久久久久| 女同久久另类99精品国产91| 国产精品精品国产色婷婷| 成人亚洲精品av一区二区| 亚洲精品在线美女| 欧美在线一区亚洲| 99久久99久久久精品蜜桃| 久久婷婷人人爽人人干人人爱| 国产精品嫩草影院av在线观看 | 99热这里只有精品一区| 好男人电影高清在线观看| 亚洲欧美日韩卡通动漫| 国产亚洲精品av在线| 男女做爰动态图高潮gif福利片| av女优亚洲男人天堂| 久久久久国产精品人妻aⅴ院| 天堂√8在线中文| 欧美中文日本在线观看视频| 大型黄色视频在线免费观看| 俺也久久电影网| 欧美乱码精品一区二区三区| 天堂√8在线中文| 亚洲av不卡在线观看| 色哟哟哟哟哟哟| 亚洲七黄色美女视频| 99精品在免费线老司机午夜| 国产精品久久久久久精品电影| 国产成人aa在线观看| 国语自产精品视频在线第100页| 国产精品日韩av在线免费观看| 久久伊人香网站| 神马国产精品三级电影在线观看| 一本一本综合久久| av福利片在线观看| 久久久精品欧美日韩精品| 久久久久精品国产欧美久久久| 久久久久亚洲av毛片大全| 18禁裸乳无遮挡免费网站照片| 亚洲无线观看免费| 哪里可以看免费的av片| 尤物成人国产欧美一区二区三区| 亚洲国产精品sss在线观看| tocl精华| 麻豆成人av在线观看| 叶爱在线成人免费视频播放| 精品久久久久久久久久免费视频| 久久婷婷人人爽人人干人人爱| 亚洲欧美激情综合另类| 亚洲成人免费电影在线观看| 亚洲黑人精品在线| 国产探花极品一区二区| or卡值多少钱| 窝窝影院91人妻| 精品国产三级普通话版| 国产精品久久久久久久久免 | 亚洲成人中文字幕在线播放| 亚洲人与动物交配视频| 色哟哟哟哟哟哟| 丝袜美腿在线中文| 搞女人的毛片| 乱人视频在线观看| 波多野结衣巨乳人妻| 国产精品一及| 亚洲国产中文字幕在线视频| 国产真实伦视频高清在线观看 | 国产成人影院久久av| 国产亚洲精品一区二区www| 岛国视频午夜一区免费看| 欧美成人免费av一区二区三区| 午夜福利高清视频| 久久精品91无色码中文字幕| 婷婷亚洲欧美| 99精品在免费线老司机午夜| 18禁国产床啪视频网站| 亚洲人成网站在线播| 日韩欧美国产一区二区入口| 在线观看免费午夜福利视频| 午夜福利在线观看免费完整高清在 | 午夜a级毛片| 午夜激情福利司机影院| 午夜福利高清视频| 麻豆国产av国片精品| h日本视频在线播放| 亚洲狠狠婷婷综合久久图片| 国产激情偷乱视频一区二区| 久久精品影院6| 久久久久性生活片| 老汉色∧v一级毛片| 亚洲五月天丁香| 熟妇人妻久久中文字幕3abv| 听说在线观看完整版免费高清| 日本在线视频免费播放| 国产爱豆传媒在线观看| 一区二区三区免费毛片| 中文资源天堂在线| 精品国产美女av久久久久小说| 热99re8久久精品国产| 精品一区二区三区视频在线观看免费| 日韩高清综合在线| 国产黄片美女视频| 99久久无色码亚洲精品果冻| 天天一区二区日本电影三级| 欧美性猛交╳xxx乱大交人| 淫妇啪啪啪对白视频| 国产精品香港三级国产av潘金莲| 成年女人毛片免费观看观看9| 亚洲 国产 在线| 久久久色成人| 色哟哟哟哟哟哟| 亚洲精品久久国产高清桃花| 三级男女做爰猛烈吃奶摸视频| 国产精品精品国产色婷婷| 国产精品电影一区二区三区| 亚洲国产欧洲综合997久久,| 夜夜躁狠狠躁天天躁| 一进一出好大好爽视频| 国产高清videossex| 少妇丰满av| 国产成人欧美在线观看| 亚洲国产精品sss在线观看| 99热6这里只有精品| 亚洲成人久久性| 18禁裸乳无遮挡免费网站照片| 欧美不卡视频在线免费观看| 欧美激情在线99| 熟妇人妻久久中文字幕3abv| 一个人免费在线观看电影| 日韩欧美在线乱码| 久久久精品欧美日韩精品| 日日摸夜夜添夜夜添小说| 中文字幕人妻丝袜一区二区| 国产成人av激情在线播放| www.色视频.com| 内地一区二区视频在线| 久久婷婷人人爽人人干人人爱| 色av中文字幕| 日本与韩国留学比较| 三级毛片av免费| 久久久久久人人人人人| 国产美女午夜福利| 99久久久亚洲精品蜜臀av| 亚洲最大成人手机在线| 天天一区二区日本电影三级| 国产成人福利小说| 欧美中文综合在线视频| 一区二区三区激情视频| 国产一区二区在线av高清观看| 日韩亚洲欧美综合| 可以在线观看的亚洲视频| 伊人久久精品亚洲午夜| 亚洲av免费在线观看| 成人亚洲精品av一区二区| 日本五十路高清| 青草久久国产| 三级男女做爰猛烈吃奶摸视频| 一个人免费在线观看的高清视频| 人妻久久中文字幕网| ponron亚洲| 久久久久久久精品吃奶| 久久久久久久久中文| 3wmmmm亚洲av在线观看| 午夜福利在线观看吧| 午夜视频国产福利| 成人三级黄色视频| 男女之事视频高清在线观看| 欧美中文综合在线视频| 又紧又爽又黄一区二区| 黄片小视频在线播放| av女优亚洲男人天堂| 国产亚洲精品综合一区在线观看| 免费看a级黄色片| 听说在线观看完整版免费高清| 精品日产1卡2卡| 精品福利观看| 久久久久久久午夜电影| e午夜精品久久久久久久| 又紧又爽又黄一区二区| 天天添夜夜摸| 怎么达到女性高潮| 国产真人三级小视频在线观看| 久久九九热精品免费| 特级一级黄色大片| 美女高潮喷水抽搐中文字幕| 免费看日本二区| 三级国产精品欧美在线观看| 一夜夜www| 美女被艹到高潮喷水动态| 精品国产三级普通话版| 97超级碰碰碰精品色视频在线观看| 国产欧美日韩一区二区三| 久久久久久久久中文| 国产成人av教育| 深夜精品福利| 亚洲五月婷婷丁香| 婷婷丁香在线五月| 亚洲最大成人中文| 久久亚洲精品不卡| a级毛片a级免费在线| 成人三级黄色视频| 淫秽高清视频在线观看| 欧美日韩福利视频一区二区| 小蜜桃在线观看免费完整版高清| 黄片小视频在线播放| 内地一区二区视频在线| 少妇人妻一区二区三区视频| 亚洲av电影不卡..在线观看| av中文乱码字幕在线| 精品无人区乱码1区二区| 不卡一级毛片| 麻豆成人午夜福利视频| 国产精品亚洲一级av第二区| 色播亚洲综合网| 又黄又粗又硬又大视频| 欧美在线一区亚洲| 欧美日韩综合久久久久久 | 欧美不卡视频在线免费观看| 欧美成人性av电影在线观看| 国产乱人伦免费视频| 51午夜福利影视在线观看| 一个人看视频在线观看www免费 | 成年女人看的毛片在线观看| 国产精品嫩草影院av在线观看 | 99精品在免费线老司机午夜| 国产一区二区亚洲精品在线观看| 波多野结衣巨乳人妻| 久久久久久久午夜电影| 香蕉久久夜色| 岛国在线免费视频观看| 国产精品99久久99久久久不卡| 免费搜索国产男女视频| 国产精品98久久久久久宅男小说| 青草久久国产| 亚洲专区中文字幕在线| 国产97色在线日韩免费| 日韩 欧美 亚洲 中文字幕| 欧美成人免费av一区二区三区| 五月伊人婷婷丁香| 国内久久婷婷六月综合欲色啪| 免费av毛片视频| 九色国产91popny在线| 欧美黑人欧美精品刺激| 成人精品一区二区免费| 国产精品女同一区二区软件 | 小蜜桃在线观看免费完整版高清| 亚洲av一区综合| 日本熟妇午夜| 欧美黑人巨大hd| 久久香蕉精品热| 午夜精品在线福利| 91字幕亚洲| 国产亚洲欧美98| 美女免费视频网站| 成人18禁在线播放| 十八禁网站免费在线| 成人性生交大片免费视频hd| 99久久99久久久精品蜜桃| 国产精品久久久人人做人人爽| 非洲黑人性xxxx精品又粗又长| 国产老妇女一区| or卡值多少钱| 观看免费一级毛片| 亚洲自拍偷在线| 毛片女人毛片| 国产精品乱码一区二三区的特点| 岛国视频午夜一区免费看| 最新在线观看一区二区三区| 亚洲 欧美 日韩 在线 免费| 一个人免费在线观看电影| 亚洲不卡免费看| 中文字幕av成人在线电影| 18美女黄网站色大片免费观看| 看黄色毛片网站| 亚洲国产欧洲综合997久久,| 法律面前人人平等表现在哪些方面| 国产精品综合久久久久久久免费| 欧美xxxx黑人xx丫x性爽| 国产精品野战在线观看| 日本 av在线| 成年免费大片在线观看| 色哟哟哟哟哟哟| 1000部很黄的大片| 国内毛片毛片毛片毛片毛片| 人妻丰满熟妇av一区二区三区| 亚洲在线观看片| av天堂在线播放| 国产亚洲av嫩草精品影院| 最近最新中文字幕大全电影3| 老汉色∧v一级毛片| 男插女下体视频免费在线播放| 亚洲av二区三区四区| 亚洲精品色激情综合| 99热这里只有精品一区| 久久人人精品亚洲av| 久久九九热精品免费| 变态另类成人亚洲欧美熟女| 日日干狠狠操夜夜爽| 少妇的逼水好多| 变态另类丝袜制服| 一本精品99久久精品77| 精品无人区乱码1区二区| 757午夜福利合集在线观看| 露出奶头的视频| 国产在线精品亚洲第一网站| 一区福利在线观看| 国产中年淑女户外野战色| 麻豆成人av在线观看| 亚洲熟妇熟女久久| 色视频www国产| 中文亚洲av片在线观看爽| 久久香蕉国产精品| 免费在线观看成人毛片| 欧美日韩综合久久久久久 | 俺也久久电影网| 欧美高清成人免费视频www| 亚洲国产精品999在线| 亚洲欧美日韩高清专用| 波野结衣二区三区在线 | 老司机午夜十八禁免费视频| 国产一级毛片七仙女欲春2| 亚洲国产精品成人综合色| 91九色精品人成在线观看| 99在线人妻在线中文字幕| 又爽又黄无遮挡网站| 人人妻人人澡欧美一区二区| 国产淫片久久久久久久久 | 看黄色毛片网站| 欧美一区二区国产精品久久精品| 国产色爽女视频免费观看| 国产成年人精品一区二区| 中文字幕久久专区| 午夜免费观看网址| 中文在线观看免费www的网站| 国产精品免费一区二区三区在线| 午夜久久久久精精品| 国产一区二区在线观看日韩 | 精品福利观看| 99久久无色码亚洲精品果冻| 国产探花极品一区二区| 男女视频在线观看网站免费| 免费av不卡在线播放| 久久6这里有精品| 国产午夜精品论理片| 天堂√8在线中文| www.熟女人妻精品国产| 亚洲成人中文字幕在线播放| 亚洲精品456在线播放app | 在线a可以看的网站| 18禁黄网站禁片午夜丰满| 久久精品夜夜夜夜夜久久蜜豆| 国产精品电影一区二区三区| 亚洲国产精品成人综合色| 精品人妻偷拍中文字幕| 国产精品 欧美亚洲| 女生性感内裤真人,穿戴方法视频| 少妇熟女aⅴ在线视频| 亚洲国产色片| 日本 av在线| 中文字幕av在线有码专区| 最近最新免费中文字幕在线| 此物有八面人人有两片| 亚洲第一电影网av| 淫秽高清视频在线观看| 中文亚洲av片在线观看爽| 成人鲁丝片一二三区免费| 美女高潮喷水抽搐中文字幕| 久久精品亚洲精品国产色婷小说| 亚洲成av人片免费观看| 日韩成人在线观看一区二区三区| 日韩大尺度精品在线看网址| 亚洲欧美日韩卡通动漫| 国产爱豆传媒在线观看| 亚洲国产中文字幕在线视频| 久久久久久国产a免费观看| 在线观看66精品国产| 久99久视频精品免费| 久久精品夜夜夜夜夜久久蜜豆| 看黄色毛片网站| 国产成人影院久久av| 免费看a级黄色片| 亚洲人成网站高清观看| 天美传媒精品一区二区| 老汉色∧v一级毛片| 亚洲男人的天堂狠狠| 成人特级av手机在线观看| 欧美区成人在线视频| 97超级碰碰碰精品色视频在线观看| 国产精品免费一区二区三区在线| 日韩国内少妇激情av| 精品一区二区三区视频在线 | 看免费av毛片| 欧美zozozo另类| 成人国产一区最新在线观看| 亚洲精品乱码久久久v下载方式 | 网址你懂的国产日韩在线| 国产美女午夜福利| 热99在线观看视频| 久久精品国产综合久久久| 久久久久国内视频| 亚洲精品一区av在线观看| 精品久久久久久久久久免费视频| 欧美成人性av电影在线观看| 97碰自拍视频| 日本一本二区三区精品| 99国产精品一区二区蜜桃av| 亚洲av成人不卡在线观看播放网| 高清在线国产一区| 国产精品女同一区二区软件 | 成人特级黄色片久久久久久久| 性欧美人与动物交配| 国产av不卡久久| av天堂中文字幕网| 欧美3d第一页| 成年女人永久免费观看视频| 一级黄色大片毛片| 欧美中文日本在线观看视频| 国产国拍精品亚洲av在线观看 | 岛国在线免费视频观看| 男人舔奶头视频| 久久国产精品人妻蜜桃| 少妇裸体淫交视频免费看高清| 麻豆久久精品国产亚洲av| 久久伊人香网站| 亚洲狠狠婷婷综合久久图片| 麻豆久久精品国产亚洲av|