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

    On numerical modeling of low-head direct chill ingot caster for magnesium alloy AZ31

    2014-04-21 02:45:09*
    Journal of Magnesium and Alloys 2014年4期

    *

    Department of Mining&Materials Engineering,McGill University,M.H.Wong Building,3610 University Street,Montreal,QC H3A 0C5,Canada

    On numerical modeling of low-head direct chill ingot caster for magnesium alloy AZ31

    Mainul Hasan*,1,Latifa Begum

    Department of Mining&Materials Engineering,McGill University,M.H.Wong Building,3610 University Street,Montreal,QC H3A 0C5,Canada

    A comprehensive 3D turbulent CFD study has been carried out to simulate a Low-Head(LH)vertical Direct Chill(DC)rolling ingot caster for the common magnesium alloy AZ31.The model used in this study takes into account the coupled laminar/turbulent melt f l ow and solidif i cation aspects of the process and is based on the control-volume f i nite-difference approach.Following the aluminum/magnesium DC casting industrial practices,the LH mold is taken as 30 mm with a hot top of 60 mm.The previously verif i ed in-house code has been modif i ed to model the present casting process.Important quantitative results are obtained for four casting speeds,for three inlet melt pouring temperatures(superheats)and for three metal-mold contact heat transfer coeff i cients for the steady state operational phase of the caster.The variable cooling water temperatures reported by the industry are considered for the primary and secondary cooling zones during the simulations.Specif i cally,the temperature and velocity f i elds,sump depth and sump prof i les,mushy region thickness,solid shell thickness at the exit of the mold and axial temperature prof i les at the center and at three strategic locations at the surface of the slab are presented and discussed.

    Low-head DC caster;Magnesium alloy AZ31;3D CFD modeling;Turbulent melt f l ow;Solidif i cation;Mushy region thickness;Sump prof i le

    1.Introduction

    Magnesium and many of its alloys are the lightest among all engineering metallic materials in use today.When compared with aluminum,magnesium is about 36 pct lighter, while it is about 78 pct lighter compared to steel[1,2].Magnesium and its alloys have better desirable properties,such as specif i c strength,stiffness,castability,damping capacity, machinability,shielding capacity,recyclability,etc compared to the most widely used materials like steel and aluminum alloys.Demand for magnesium alloys is increasing each year particularly in automotive,electronics and aerospace industries where weight-saving engineering solutionsare continuously being searched for reducing the cost for energy consumption.Among the many casting routes,today about 20 pct of magnesium alloys is cast via the direct chill casting (DCC)process[3].This process allows casting a large volume of magnesium alloys in an energy eff i cient manner.One of the major drawbacks in using this technology for magnesium alloys is that it costs more compared to steel and aluminum alloys[4].It is to be stressed here that DCC process for magnesium varies signif i cantly from plant to plant and unlike DCC of aluminum,the magnesium casting plants in this regard keep their own successful casting practices quite guarded.

    The direct chill cast magnesium alloys suffer from a number of unwanted defects,such as hot tears,cold cracks, center cracks,surface folds,etc.Another technical problem is that modeled magnesium alloys have a low volumetric heatcapacity and as a result they freeze quickly making the start-up quite diff i cult.Since magnesium is highly reactive with oxygen hence during the casting operations care is to be taken so that such contact does not take place and in this regard during the DCC operation the top part of the melt delivery system is kept under pressure using a mixture of CO2(90%by volume) and SF6(10%by volume)gases[5].Also,one needs to be very careful while casting these alloys using the DCC process, since they may cause a hydrogen explosion if the alloys come accidentally in contact with the primary or secondary cooling water of the process.

    The vertical DCC process,which is used for aluminum, magnesium and a few other non-ferrous metals like zinc, copper,etc,is quite complex due to the interactions of various parameters,likealloy composition,thethermo-physical properties of the alloy,casting speed,pouring temperature of the melt,mold length and cross-section,cooling water temperature and its f l ow rate in the mold and post mold region,etc [6].Though this process has seen many improvements over the past several decades but due to the complex interactions of the above mentioned factors,the successful operation of this process still depends upon the prior experience of the operators gained in earlier years in the plant.

    Using a suitable melt delivery system,in a conventional vertical DCC process superheated molten magnesium alloy is supplied from the top to a water-cooled open-ended mold. Before the start of casting,the open bottom end of the mold is kept plugged by a bottom block(also called a starter block). When the supplied molten metal level in the bottom-blocked mold reaches a certain level and the solidif i ed layers form against the top of the bottom block and the mold faces,the partially solidif i ed ingot sitting on the top of the starter block is slowly lowered down towards the casting pit below by the f i tted hydraulic ram on the bottom block.During the lowering of the block-ingot assembly by the ram its speed is increased gradually until a constant casting speed is reached.As the partially solidif i ed ingot containing liquid metal within it exits the bottom of the mold it is further cooled by directly spraying mold water in the form of jets or streams from the bottom of the mold onto the outside surface of the solid shell.The water stream from the impingement region then forms a continuously falling f i lm down the slab.The indirect cooling of the metal that takes place in the mold is called primary cooling where only 10-15%of the heat content of the ingot is removed by the cooling water circulating through the various passages within the mold body[3].Through the direct cooling of the partially solidif i ed ingot by the water from the mold in the post mold region(also called secondary cooling region), the rest(85%-90%)of the total heat content of the ingot is removed[3].The process is initially unsteady.However, depending on the casting speed,the thermo-physical properties of the magnesium alloy and the cooling conditions in the mold and post mold regions after about 0.5 m-1.0 m length of casting the process reaches steady state[3].In a normal slab casting operation after casting of about 8-10 m in length,the melt delivery to the mold is stopped and the casting operation is suspended.The cast is then taken out from the water containing casting pit for further downstream processing.The cast magnesium slabs are usually rolled to produce plates, sheets,foils and strips for various end applications while the billets are extruded and forged to form bars,rods,pipes,tubes, etc.for various practical uses.

    Due to the customers'demands for the higher quality products from the castings,the major manufacturers of the aluminum and magnesium castings have adopted a number of technological changes in their DCC practices with a view to improving the cast quality.Because of the interactions of several interlinked parameters in this process it is diff i cult to prevent the formation of defects by controlling one parameter at a time.The use of a low-head mold with an active mold length between 25 and 30 mm has been seen to suppress the formation of an air gap in the mold.And as a result,the reheating of the shell is minimized which signif i cantly reduces exudation,surface macro-segregation and surface cracks,etc. In this regard,the casting company Wagstaff,Inc.(Spokane Valley,WA 99206,USA),a major supplier of DC casting machines,developed a low-head composite(LHC)ingot casting technology in 80s which the company since has been marketing to the various aluminum and magnesium casting industries under their trade name ‘LHC?’[7].In comparison with the conventional DCC mold(which is usually between 70 and 80 mm in length)the LHC mold has been shown to have the following advantages;it produce casts with a better surface quality,it offers safer operations,it minimizes scalp rates,it allows higher casting speeds and results in a longer mold life. The latter improvements have led to the better cast quality for a lower cost.In this regard,the interested readers can consult the recently published book on DCC process[2]for the historical developments,the modern versions of this process as well as for the various experimental and modeling studies related to casting of aluminum and magnesium alloys.It is interesting to note that unlike aluminum alloys,there is not enough modeling studies available in the literature concerning the present problem.This can be verif i ed by looking at the cited references of Hu et al.[8]who very recently,modeled the 3-D temperature f i eld for a DC castslab of 800 mm × 300 mm cross-section of AZ31 alloy for casting speeds ranging from 20 to 70 mm/min all for a f i xed pouring temperature of 600°C.They have reported the sold shell thickness,and the temperature contours and temperature gradients at the wide symmetry plane for various mentioned casting speeds.Although the paper stated that a thermo-f l uid analysis was carried out but no supporting results have been reported,nor have they provided any explanation on how the melt f l ow analysis was carried out.

    To the best knowledge of the authors,no one has ever modeled the 3-D coupled laminar/turbulent melt and solidif ication heat transfer which take place in a low-head vertical DCC slab caster for any magnesium alloy and certainly not for the alloy AZ31 modeled in this study.Since casts of AZ31 slabs of various sizes are produced by the DCC industries,the objective of the present investigation is to carry out a comprehensive 3-D numerical simulation study of the L-H DCC process for the alloy AZ31 in order to gain somefundamental understanding of the solidif i cation behavior of this alloy and also to provide some useful guidelines and offer some concrete suggestions to the DCC operators to improve their casting campaigns with regard to the productivity and the quality of the casts from the said alloy.

    2.Mathematical model

    The system modeled is shown schematically in Fig.1 which represents a low-head rectangular-shaped industrialsized caster with an open-top melt feeding system.Since the modeled geometry has a two-fold symmetry,to save the computational time and cost only the top left-hand quarter of the slab is investigated.A f i xed domain having dimensions of 2500 mm × 865 mm × 330 mm inx,y,andzdirections, respectively with its origin located at the central top part of the ingot with the axial coordinate in the direction of gravity is modeled.The melt of the alloy AZ31 is considered to have been supplied at an uniform temperature and velocity along the entire top cross-section of the mold.The solid-shell which forms in the mold and in the sub-mold regions is assumed to move in the axial direction at the selected but constant casting speed.

    2.1.Test material

    Fig.1.Schematic of the low-head DC caster and the computational domain for an open-top melt feeding scheme.

    As stated earlier,the common magnesium alloy AZ31 was selected to predict the solidif i cation behavior of this alloy in a vertical direct chill slab casting process.The chemical composition limits of this alloy are listed in Table 1[9]and its thermo-physical properties are summarized in Table 2[9].In this magnesium alloy,aluminum and zinc are the two major alloying elements along with other minor elements like manganese,silicon,copper,etc.This alloy possesses good roomtemperature strength and ductility,has good corrosion resistance property and can be welded.The alloy AZ31 f i nds applications in many useful products such as aircraft fuselages, cell phones and laptop cases,intricate components of automobiles,vibration tests equipment,inspection gauges,speaker cones,etc.[9].

    2.2.Assumptions and simplif i cations

    In order to simplify the modeling effort of the present complex industrial slab DCC problem,the following realistic assumptions were invoked:

    [1]A f i xed-coordinatesystem (Eulerian approach)was employed in the simulation.

    [2]Local thermodynamic equilibrium during solidif i cation was assumed to prevail.Heat f l ow was taken to be very fast and every point reached equilibrium with its neighboring points instantaneously.

    [3]A mushy-f l uid solidif i cation model was assumed and pore formation was ignored.

    [4]Flow in the mushy region was modeled similar to a f l ow through a porous medium.

    [5]Molten aluminum was assumed to behave as an incompressible Newtonian f l uid and turbulence effects were approximated using the popular two equationk- ε model of turbulence.

    [6]Any heat released due to solid-solid transformation was not taken into consideration.The formation of intermetallic compounds during solidif i cations was ignored.

    [7]Evolution of latent heat in the solidif i cation domain was not inf l uenced by the microscale species transport,that is, the liquidus and solidus temperatures were f i xed.

    [8]Because of the lack of experimental results,the variation of liquid fraction in the mushy zone was assumed to be a linear function of temperature.

    [9]The thermo-physical properties of magnesium alloy were invariant with respect to the temperature and there was no viscous dissipation effect.The thermal buoyancy effect wasincorporated in the momentum equationsby employing the Boussinesq approximation.

    [10]Only the steady-state phase of the caster's operation was modeled without taking into account the transient start-up condition.

    2.3.Numerical formulation

    2.3.1.Governing equations for turbulent melt velocity and enthalpy

    The 3-D general mean transport equations for an incompressible laminar/turbulent melt velocity and enthalpy in arectangular domain can be expressed in a general form of a dimensional partial differential equation.The conservation equation in the general form at steady state in Cartesian tensor notation can be written as follows:

    Table 1Chemical composition limits of AZ31.

    Table 2Thermo-physical properties of AZ31.

    where,ρ is the melt density anduiis the velocity component in thexi-direction.For continuity equation Φ is equal to 1.For momentum equations it represents the velocity componentsui, while for the energy equation it represents temperatureT,and for the turbulent kinetic energy and its rate of dissipation it representskand ε,respectively. ΓΦis the effective diffusion coeff i cient of Φ andSΦis the source term of Φ.The gravity direction is taken in the positivex-direction.As the gravity term in thex-momentum equation may have an effect on the melt f l ow,the density in the buoyancy force terms is assumed to be dependent on the temperature f i eld based on the Boussinesq approximation.To calculate the turbulence quantities in the liquid pool and mushy region the low-Reynolds number version of the populark-ε eddy viscosity concept proposed by Launder and Sharma[10]were used.The details derivation of the above general equation and turbulence model can be found in Begum[11].

    The enthalpy-porosity technique was adopted to account for liquid-mushy and mushy-solid phase changes in the domain, in which a single set of transport equations was appropriately applied for the liquid,mushy and solid regions.In order to model the f l ow in the mushy zone where 0 <fl< 1,the Darcy's law of porous media based on the Carman-Kozeny equation was adopted[12,13].Equation(1)and the associated boundary conditions were suitably non-dimensionalized to suppress the numerical instabilities which often occur during the solution process of the dimensional discretized equations due to the large change in the values of the common variable Φ.

    The following dimensionless variables were employed to non-dimensionalize the governing coupled partial differential equations and the associated boundary conditions in order to ascertain the non-dimensional parameters which govern this process:

    whereDis the hydraulic diameter of the caster,uinis the inlet velocity and ΔHfis the latent heat of solidif i cation.All the appropriate variables of the transport equations in the nondimensional form are available in Table 3 in Ref.[14].The non-dimensional form of the boundary conditions is also given in Eqs.(10)-(16)in Ref.[14].In order to keep the article short,the non-dimensional transport equations and boundary conditions are not reported here.The interested readers can consult the recently published PhD dissertation of Begum [11],which can be easily obtained free of charge by Googling the title,for further details.

    The values of the effective heat transfer coeff i cient along the length of the DC slab caster were varied and were taken from the very well cited article by Vreeman and Incropera [15],which are given below:

    Adiabatic section:

    Mold-metal contact region:

    Secondary cooling zone:

    Table 3Values of the parameter used in open-top melt delivery arrangement.

    wherex1=90 mm,x2=100,x3=140 mm and γmax=20.0 W/m2K, γflim=10.0 W/m2K.In the present study,three different cooling water temperatures,namely 10, 20 and 30°C,respectively for the mold,impingement and streaming regions have been considered to ref l ect the reality of the thermal conditions of the cooling water used by the L-H DCC industry[16].A complete description of the solution methods employed in this work is provided in details in Ref. [11].

    Before carrying out extensive parametric studies the grid independency tests were f i rst performed using three sets of successively ref i ned meshes(coarse,moderate,and f i ne).The relative difference in the local surface heat f l ux between the coarse and f i ne grids was within 4.0%,indicating that the coarse grid was suff i cient for predicting the results of engineering accuracy.It is to be stressed here that the local surface heat f l ux is the most sensitive quantity with regard to grid in any heat transfer problem and this is particularly true for turbulent convective heat transfer problems.A comparison (see Ref.[11])of the surface temperature prof i les showed less than 1%variation between the f i ne and course grid systems. Therefore,the coarse grid (60 × 42 × 24)was used for all production runs reported in this work.The details of the predicted results along with the comparisons can be found in Ref.[11].

    A code validation test was also carried out by simulating the solidif i cation of a rolling ingot for AA-3104 alloy,using the identical experimental casting conditions and set-up for the vertical DCC process performed by Jones et al.[17].The present in-house developed code showed an acceptable agreement with the experimental data.These results were presented in Ref.[11]and are not given here to avoid duplication.

    3.Results and discussion

    A detailed study of all casting parameters namely,casting speed,inlet melt superheat,HTC,cooling water f l ow rate, mold dimensions,etc.will certainly generate too many results for an arbitrarily selected 3-D geometry.Since the main purpose of this work is to show the inf l uence on the solidif i cation behavior of the most important casting parameters hence an industrial-sized caster with a low-head DC caster was selected for modeling studies by varying three process parameters, namely,casting speed,inlet melt superheat,and HTC at the metal-mold contact region.A set of twelve separate cases for the AZ31 alloy which is listed in Table 3 were selected to convey the most important casting information.

    3.1.Temperature and velocity f i elds for cases(1-4)of Table 3

    Fig.2.3-D surface plots of temperature f i eld for four cases(1-4).

    The surface temperature distributions in f l ooded format are presented in Fig.2(i-iv)for the range of casting speed from 60 to 180 mm/min representing cases(1-4)of Table 3.These numerically predicted temperature data clearly show the progression of the solidif i cation process for the above mentioned cases.The solidif i cation range of magnesium alloy AZ31 is 55°C which is an intermediate freezing range alloy.The present results indicate that the mushy region,which isbounded by the liquidus(630°C)and solidus(575°C)isotherms,is expanded as the casting speed is increased due to the stronger thermal convective effects.Furthermore,a classical parabolic-shaped solidif i cation front is seen to have been developed particularly if one looks from the central axis at each casting speed.Since near the wall the melt loses superheat and as a result the effect of thermal convection reduces compared to thermal conduction effect as the melt moves from the wall region to the central region.The parabolic shape of the solidif i cation prof i le is the manifestation of the above fact. As the casting speed increases the solidif i cation front as well as the other isotherms are becoming steeper and are moving downward,which has already been found by several other authors[18,19].

    Fig.3.Longitudinal 2-D view of temperature contours and velocity vectors for case 1.

    Fig.2(i-iv)shows only the temperature f i elds for the two symmetric planes(x-yandx-zplanes)and for the top free surface.In order to assess the melt f l ow and heat transfer along the longitudinal planes of the domain,the 2-D projections of the temperature and velocity f i elds at the wide symmetry plane (z=0)and parallel to the wide symmetry plane at two inward locations,viz.,z=62.5 mm,andz=312.5 mm,for two cases(1 and 4)are selected.The right hand panels ofFigs.3(a)and 4(a)are plots for the velocity vectors at wide symmetry plane(z=0)for two casting speeds of 60 and 180 mm/min.As mentioned earlier,in this study it has been assumed that the melt is delivered uniformly at a specif i ed casting speed across the entire top cross-section.Due to this simple melt feeding system,from these f i gures it can be clearly seen that all of the melt f l ows downward.Along the narrow face a part of the vertical melt that are adjacent to the wall is diverted by the developed solidif i ed shell,and then they move along the solidif i cation front.A comparison between these two f i gures shows that the magnitude of the resultant velocity vector increases with the increase of the casting speed.The above observation is due to the conservation of mass into the domain.

    The right hand panels of Figs.3(b)and 4(b)are plots for the velocity f i elds on the longitudinal planes at an inward distance ofz=62.5 mm for the same casting speeds.As depicted in those f i gures,the developing solidif i ed shell is still obstructing the corner stream.The right hand panels of Figs.3(c)and 4(c) are showing the velocity f i elds atz=312.5 mm.Notice that the plane atz=312.5 mm is nearer to the rolling face.The difference among the three planes is that forz=312.5 mm the f l ow is moving downward uniformly without any obstruction.This is because on the third plane the melt is almost solidif i ed and the magnitude of the velocity vectors there is slightly higher than the corresponding casting speed,which signif i es that the melt is accelerated in this plane with respect to the imposed casting speed on the solidif i ed shell.

    Fig.4.Longitudinal 2-D view of temperature contours and velocity vectors for case 4.

    The heat transfer rate increases as one proceeds from the heated ingot center toward the cooled wall.The numerically predicted temperature contours at the three longitudinal planes give a better visualization of this trend.The location of a f i xed isotherm is seen to have been lifted upward progressively at two inward planes atz=62.5 mm andz=312.5 mm with respect to the wide symmetry plane atz=0 mm.The left hand panels of Figs.3(a)and 4(a),Figs.3(b)and 4(b),and Figs.3(c) and 4(c)also show the impact of the change in casting speed on the temperature f i elds.The results in these f i gures depict that all the isotherms have moved downward and the width of the mushy region is increasing with increasing casting speeds.

    Among the various possible options to represent the solidif i cation heat transfer phenomena,the 2-D transverse crosssectional planes(y-zplane)are chosen to get a vivid picture of the solidif i cation process happening within the ingot.The liquidus and solidus isotherms are plotted in Fig.5(a)-(d)for four cases(1-4).A total of seven transverse cross-sections starting from the top up to an axial distance of 700 mm have been selected.The numerical results show that a nearly uniform thick solid shell on the narrow and wide faces is formed and as one moves downward the extent of the solid shell increases progressively.Since more heat is being extracted from the ingot through the mold and by the chilled water jets both of the latter factor are contributing in increasing the thickness of the solid-shell along the cast direction.It appears from the progression of the solid-shell thickness that an almost uniform rate of heat extraction is taking place from both the sides.At the corner region heat isbeing extracted from both the rolling and narrow faces,this has resulted in a higher rate of heat extraction at that location compared to the other places.As a result,an almost roundshaped solid layer and mushy zone have formed there.The location of the solidus isotherm with the decrease in the casting speeds has shifted more toward the ingot center over each cross-section compared to the higher casting speed, which can be clearly seen in the above f i gures.The primary reason in the reduction of heat transfer with increasing casting speed is that for a higher casting speed the amount of incoming melt as well as the energy content of the melt are both higher.

    Fig.5.Contours of solidus and liquidus temperatures at various transverse cross-sectional planes for four cases(1-4).

    3.2.Quantitative analysis forΔT=32°C

    3.2.1.Sump depth and mushy layer thickness

    The sump depth and mushy thickness are indirect indicators of the quality of the cast.The latter two quantities were obtained by post processing the predicted temperatures from the CFD model.In this regard,the numerically predicted sump depth and the mushy layer thickness at the ingot center for four cases(1-4)are provided in Table 4.The present authors are surprised by the fact that in none of the vast literature on DCC process has reported the quantitative values of the sump depth and mushy thickness for any alloy.A closer look at this table shows that a deeper sump and a thicker mushy layer have developed at a higher casting speed.To get a comparative picture of the effect of the casting speed on sump depth and mushy thickness the lowest casting speed of 60 mm/min was taken as the base case for relative comparisons.The sump depth is found to be approximately 487.803 mm from the top of the mold for a withdrawal speed of 60 mm/min.The relative difference in sump depth for higher casting speeds of 100,140, and 180 mm/min,is respectively about 79.72%,197.47%, 308.75%higher.For the casting speed of 60 mm/min,the mushy thickness is found to be 318.091 mm,whereas for casting speeds of 100 mm/min,140 mm/min,and 180 mm/ min,the mushy thickness increases by 80.52%,215.91%,and 333.614%,respectively compared to the casting speed of 60 mm/min.The reason for the proportionate increase of the mushy layer thickness is diff i cult to explain,since a complex interaction between the f l ow and thermal f i eld prevails in that region.

    3.2.2.Predicted shell thickness

    As the casting speed inf l uences the growth rate of the solid shell,the predicted shell thickness from the narrow side at the wide symmetry plane and from the wide side at the narrowsymmetry plane at the exit of the mold are plotted in the form of a bar chart in Fig.6 in order to get a quantitative picture of the differences among the results for four cases(1-4).The numerical value of the predicted shell thickness(in mm)is provided inside each bar of the chart of Fig.6.For both locations,with the increase of the casting speed the shell thickness decreases which is consistent with the physics of the problem.The rate of decrease of the shell thickness is smaller for the higher casting speeds(>100 mm/min)compared to the lower casting speeds(<100 mm/min).The shell thickness results further show that the rate of growth of the shell at the center of the narrow side is greater compared to the center of the rolling face for all four simulated casting speeds.The reason behind this is the proximity of the middle of the narrow side from the slab corner(region of high heat extraction) compared to the middle of the wide slab side.

    Table 4Sump depth and mushy thickness in mm at the ingot center for cases(1-4).

    3.2.3.Predicted strand surface temperature

    Fig.6.Shell thickness at the middle of the slab faces at mold exit vs.casting speed for cases(1-4).

    Fig.7.Variations of surface temperature along the axial direction of the strand at four locations of the caster for:(a)us=60 mm/min(case-1)(b) us=180 mm/min(case-4).

    Fig.7(a)and(b)illustrates the temperature distributions along the caster for two cases(1 and 4)which correspond to the casting speed of 60 and 180 mm/min,respectively.It is to be realized that the surface temperatures were not known initially.In order to obtain the surface temperature at each grid location,a heat balance on the control volume corresponding to the grid point was done after the converged solution was obtained.A total of four strategic locations are selected in this study,namely:(a)center;(b)mid-distance of the wide face; (c)mid-distance of the narrow face;and(d)corner point of the caster to portray the surface temperatures.These locations are thought to be the most probable sensitive zone for the developments of the compressive and tensile stresses.The four locations are given on the top of Fig.7(a)and(b).Except for the temperature at the ingot central region,the temperatures are seen to follow the same pattern at three other locations. The temperature drops sharply from the top of the mold up to the initial part of the cooling water streaming zone.i.e.in the region within 60 mm-122 from the top surface.Beyond that region the surface temperatures are seen to rebound due to the imposition of a lower heat transfer coeff i cient and higher cooling water temperature to ref l ect the industrial operating conditions(refer to Section 2.3.1).It is to be noted that in the DCC process the cooling water for the impingement region and for the free streaming region comes from the bottom of the mold and as a results the water temperatures in the two regions increase due to the gain of heat from the ingot.The surface temperatures increased by about 7°C within the length of about 20 mm(142 mm from the top)in the streaming zone and then they decreased gradually up to the end of the simulated casting length.The rebound of the surface temperatures could promote cold crack along with other surface cracks.The corner temperatures in both the f i gures show relatively lower values with respect to the other locations due to the simultaneous extractions of heat through the wide and narrow sides of the slab.The temperatures at the ingot center in the cast direction are the highest compared to the temperatures at three other locations as expected.A comparison of the two f i gures clearly show that for a higher casting speed the surface remains overheated for the same axial location compared to the lower casting speed.This trend is visible in all the four ingot locations.The reason behind this is the fact that the predictions were carried out for the lower and higher casting speeds by invoking the same convective heat transfer boundary conditions without considering the enhanced cooling requirements for a higher casting speed.

    3.3.Effect of primary coolant heat transfer coeff i cient

    There are serious disagreements among the researchers working in this f i eld regarding the values of the effective heat transfer coeff i cients in the mold region.From the literature one gets conf l icting values of the HTCs in the metal-mold contact region.Even the well refereed literature on this issue has reported a wide range of values for the HTCs.In the metal-mold contact region,a value of HTC ranging from 1000 to 5000 W/ (m2-K)has been reported in the well-cited literature[2,3].In order to ascertain the impact of the heat transfer coeff i cient (HTC)at the mold-metal contact region,for an inlet superheat of 32°C,three HTCs values of 750,1500 and 3000 W/m2K were selected for the two casting speeds namely,100 mm/min and 180 mm/min.Fig.8 reports the values of the shell thickness against the HTCs in the form of a bar chart for the above mentioned two casting speeds.

    The shell thickness is estimated from the temperature prof i les at the mold exit at an axial distance of 90 mm from the top surface on the wide symmetry plane from the narrow face. It can be seen that for an increase of HTC the magnitude of the shell thickness increases.Similar trends are found for both casting speeds.For a f i xed HTC,the shell thickness decreases with casting speed as noticed earlier.For a lower casting speed of 100 mm/min,an increase of the HTC from 750 to 1500 W/ (m2-K)has caused an enhancement of 6.67%in solid-shell thickness and for the increase of HTC from 1500 to 3000 W/(m2-K)the increment is almost the same and is approximately 6.79%.For a higher casting speed of 180 mm/ min,the corresponding increase in the shell thickness is about 37.45%,and 28.64%,respectively.Hence,it is clear that the effect of lower HTCs(from 750 to 1500 W/m2-K)on shell thickness is more pronounced at a higher casting speed.

    Fig.8.Shell thickness at the middle of the narrow face at mold exit vs. effective HTC(W/m2-K)for six cases(2,4,and 5-8).

    Fig.9.Shell thickness at the middle of the narrow face at mold exit vs.melt superheat(°C)for six cases(2,4,and 9-12).

    3.4.The effect of inlet melt superheat(ΔT)on shell thickness

    The inlet melt superheat is usually f i xed by the casting operator in advance of the casting campaign.If a too low superheat is used there could a possibility of the melt-freeze in the launder-trough assembly during delivery of the melt into the mold which in turn could lead to the suspension of the casting process.On the contrary,if a very high inlet melt superheat is imposed,there could be the possibility of the formation of inter-metallic compounds which would adversely affect the quality of the cast.To study the effect of inlet melt superheat,a total of three superheats,namely 16,32 and 64°C were selected for each of the two casting speeds of 100 and 180 mm/min.The results of the shell thickness at the widesymmetry plane from the narrow face at the exit of the mold are presented in Fig.9 in the form of a bar chart for easy understanding.It can be seen from the above f i gure that for an increase of the superheat,the magnitude of the shell thickness decreases as technically expected.The latter observation is true for both casting speeds.For a lower casting speed of 100 mm/min an increase of the superheat from 16 to 32°C has caused in the reduction of shell thickness of only about 1.45%, but for the increase of superheat from 32 to 64°C the decrease of shell thickness is relatively greater and is about 5.41%.On the contrary,for a higher casting speed of 180 mm/min,the corresponding decrease in shell thickness is about 3.68%and 26.35%,respectively.The predicted results thus indicate that at a higher casting speed the change in superheat from 32 to 64°C causes a greater impact on the solidif i cation process.

    4.Conclusions

    After investigating the effects of the casting speed,inlet melt superheat and effective heat transfer coeff i cient at the metal-mold contact region through the 3-D turbulent CFD modeling study of a low-head,vertical DC slab casting process for the magnesium alloy AZ31,the following conclusions were drawn:

    [1]For the range of casting speeds from 60 to 180 mm/min, for the superheat range of 16°C-64°C and for the metalmold effective heat transfer coeff i cient ranging between 750 W/m2K to 3000 W/m2K the sump depth and the mushy thickness at the center of the ingot do remain conf i ned within the simulated axial length of 2500-mm of the caster.

    [2]A shallower sump and a lower depth of the melt pool at the center of the caster are obtained for the lower casting speeds compared to the higher casting speeds.For the standard cases of 1-4 of Table 3,the sump depth(SD in mm)and the melt pool depth(MP in mm)can be represented respectively by the following linear equations with regard to the casting speed (usin mm/min): SD = -325.41 + 12.74us(R2= 0.993)and MP=-264.07+9.02us(R2=0.99)

    [3]The vertical extent of the mushy region(which is bounded by the liquidus and solidus temperatures of the alloy)at the ingot center increases with the increase of the casting speed.For the standard cases of 1-4 of Table 3,the mushy thickness(MT in mm)can be represented by the following linear equation with regard to the casting speed(usin mm/min):MT=-264.07+9.02us(R2=0.99)

    [4]The solid shell thickness at the exit of the mold decreases with the increase in the casting speed which also resulted in the increase of the surface temperatures of the strand.

    [5]The shell thickness at the mold exit decreases with the increase of the pouring temperature(melt superheat)and the reduction in shell thickness is more signif i cant for the higher range of the superheat,particularly at the higher casting speed.

    [6]For a constant cooling water temperature in the mold, with a higher imposed HTC at the metal-mold contact region the shell thickness at the mold exit is greater.The enhancement of the shell thickness is more pronounced for the lower HTCs at higher casting speeds compared to the lower casting speeds.

    [7]The higher casting speeds have resulted in the stronger thermal convective f l ows in the melt pool and in the mushy region compared to the lower casting speeds.The implication of this is that a greater number of broken dendrites can be expected for the higher casting speeds for this predominantly forced convective melt f l ow system and thus the macro-segregation patterns within the cast will be different for higher casting speeds compared to the lower casting speeds.

    Acknowledgments

    This work is partially supported from the National Sciences and Engineering Research Council(NSERC)of Canada Discovery Grant RGPIN48158 awarded to M.Hasan of McGill University,Montreal,for which authors are grateful.

    Nomenclature

    A Darcy coeff i cientap,anb,bcoeff i cients in the discretized governing equationsc1,c2,cμempirical constants for low Reynolds number modelcp specif i c heatC morphology constantD hydraulic diameter of the caster(the entire crosssection of the ingot)Dk extra dissipation term ink-equationEε extra generation term in ε-equationf1,f2,fμempirical constants used in low-Re version ofk-? modelsfl liquid fractionfs solid fractionG production term in turbulent kinetic energy equation Gr Grashof number Grm modif i ed Grashof number(Gr/Ste)h sensible heath∞ ambient enthalpyH total heat(sensible and latent)k turbulent kinetic energyP hydrodynamic pressure Pe Peclet number Pr laminar Prandtl number Re Reynolds number Ret turbulent Reynolds number based on the turbulent quantitiesS source termSΦ source term associated with ΦSk* non-dimensional source term associated with buoyancy term ink-equationSε* non-dimensional source term associated with buoyancy term in ε-equation Ste Stefan numberT temperatureT′ f l uctuation of temperatureTin inlet temperatureTl liquidus temperatureTs solidus temperatureTsurf slab surface temperatureui velocity component in thei-thdirection;corresponding touui time-average velocity component in thei-th directionu′i l uctuation of velocity in thei-th directionuin inlet velocityus casting speedU,V,Wnon-dimensional form of theu,vandwvelocitiesUs non-dimensional form ofusx axial directiony horizontal direction parallel to the wide facez horizontal direction parallel to the narrow faceX,Y,Znon-dimensional form ofx,y,zi,j,k coordinate direction f

    Greek symbolsΔH nodal latent heat ΔHf latent heat of solidif i cation Γφ diffusion coeff i cient associated with Φ ρ alloy density μl laminar viscosity μt turbulent viscosity or eddy diffusivity μe effective viscosity equal to μl+ μtΦ generalized dependent variable Γeff effective diffusivity ? rate of energy dissipation ?in inlet rate of energy dissipation γ effective convective heat transfer coeff i cient σk, σ? turbulence model constants σt turbulent Prandtl numberSuperscripts* non-dimensional variables - time-averaged variables′fl uctuation of variables.

    [1]H.Watari,T.Haga,N.Koga,K.Davey,J.Mater,Process.Technol. 192-193(2007)300.

    [2]J.F.Grandf i eld,D.G.Eskin,I.Bainbridge,Direct-chill Casting of Light Alloys:Science and Technology,John Wiley and Sons,Inc.,Hoboken, New Jersey,2013,p.242,334.

    [3]E.J.F.R.Caron,M.A.Wells,Metall.Mater.Trans.B 40B(2009)585.

    [4]P.W.Baker,P.T.McGlade,Magnesium Direct Chill Casting:a Comparison with Aluminium,Light Metals,TMS,Warrendale,PA,2001,pp. 855-862.

    [5]S.S.Park,W.J.Park,C.H.Kim,B.S.You,N.J.Kim,JOM 61(8)(2009) 14.

    [6]D.G.Eskin,Physical Metallurgy of Direct Chill Casting of Aluminum Alloys,Advances in Metallic Alloys,CRC Press,Taylor&Francis Group,Boca Raton,FL,2008,pp.120-122.

    [7]B.Rinderer,P.Austen,A.Tuff,Casthouse Modif i cations for Improved Slab Quality,Light Metals,The Minerals,Metals&Materials Society, Warrendale,PA,2003,pp.1-6.

    [8]W.Hu,Q.Le,Z.Zhang,L.Bao,J.Cui,J.Magnesium Alloys 1(2013) 88-93.

    [9]http://asm.matweb.com and www.aircraftmaterials.com/data/aluminum/ 5052/html.

    [10]B.E.Launder,B.I.Sharma,Lett.Heat Mass Transfer 1(1974)131.

    [11]L.Begum,3-D Transport Phenomena in Vertical Direct Chill Casting Processes,Ph.D.thesis,Dept.of Mining and Materials Engineering, McGill University,Montreal,Quebec,Canada,2013.

    [12]V.R.Voller,C.A.Prakash,Int.J.Heat Mass Transfer 30(1987)1709.

    [13]S.H.Seyedein,M.Hasan,Int.J.Heat Mass Transfer 40(1997)4405.

    [14]L.Begum,M.Hasan,Int.J.Heat Mass Transfer 73(2014)42.

    [15]C.J.Vreeman,F.P.Incropera,Int.J.Heat Mass Transfer 43(2000)687.

    [16]Personal communications,http://www.Wagstaff.com.

    [17]W.K.Jones Jr.,D.Xu,J.W.Evans,W.E.Williams,D.P.Cook,Light Met. (1999)841-845.

    [18]M.R.Amin,D.Greif,Int.J.Heat Mass Transfer 42(1999)2883.

    [19]M.R.Amin,J.Thermophys.Heat Transfer 14(2)(2000)170.

    Received 18 September 2014;accepted 20 November 2014 Available online 29 December 2014

    *Corresponding author.Tel.:+1 514 398 2524,+1 514 924 8542(mobile); fax:+1 514 398 4492.

    E-mail addresses:Mainul.hasan@mcgill.ca(M.Hasan),Latifa.begum@ mail.mcgill.ca(L.Begum).

    Peer review under responsibility of National Engineering Research Center for Magnesium Alloys of China,Chongqing University.

    1ASME Member(#1964840).

    http://dx.doi.org/10.1016/j.jma.2014.11.008.

    2213-9567/Copyright 2014,National Engineering Research Center for Magnesium Alloys of China,Chongqing University.Production and hosting by Elsevier B.V.All rights reserved.

    Copyright 2014,National Engineering Research Center for Magnesium Alloys of China,Chongqing University.Production and hosting by Elsevier B.V.All rights reserved.

    在线免费观看不下载黄p国产 | 国产精品影院久久| 国产精品久久电影中文字幕| 中国美女看黄片| 午夜亚洲福利在线播放| 成在线人永久免费视频| 国产私拍福利视频在线观看| 亚洲精品在线观看二区| 国产午夜福利久久久久久| 国内精品一区二区在线观看| 小蜜桃在线观看免费完整版高清| 成人亚洲精品av一区二区| 国产欧美日韩一区二区精品| 观看美女的网站| 热99在线观看视频| 久久久久久九九精品二区国产| 不卡av一区二区三区| 一本综合久久免费| 观看免费一级毛片| 欧美一级毛片孕妇| 夜夜夜夜夜久久久久| 99热6这里只有精品| 久久九九热精品免费| 亚洲av电影在线进入| 伊人久久大香线蕉亚洲五| 三级国产精品欧美在线观看 | 成人特级黄色片久久久久久久| 国产精品永久免费网站| 很黄的视频免费| 成人高潮视频无遮挡免费网站| 制服人妻中文乱码| 熟女少妇亚洲综合色aaa.| 三级毛片av免费| 嫁个100分男人电影在线观看| 国产精品自产拍在线观看55亚洲| 男女做爰动态图高潮gif福利片| 午夜视频精品福利| 亚洲 欧美 日韩 在线 免费| 欧美日本视频| 国产人伦9x9x在线观看| 狂野欧美白嫩少妇大欣赏| a在线观看视频网站| 极品教师在线免费播放| 亚洲精品美女久久久久99蜜臀| 久久人人精品亚洲av| 国产极品精品免费视频能看的| 日韩精品中文字幕看吧| 啦啦啦韩国在线观看视频| www.熟女人妻精品国产| 久久久国产成人免费| 欧美xxxx黑人xx丫x性爽| a级毛片a级免费在线| 黄色女人牲交| 亚洲自拍偷在线| 日韩欧美 国产精品| 夜夜看夜夜爽夜夜摸| 变态另类成人亚洲欧美熟女| 一个人看视频在线观看www免费 | 成人亚洲精品av一区二区| 国产麻豆成人av免费视频| 成年人黄色毛片网站| 69av精品久久久久久| 女生性感内裤真人,穿戴方法视频| av国产免费在线观看| 91av网一区二区| 一级a爱片免费观看的视频| 午夜影院日韩av| 亚洲成人中文字幕在线播放| 午夜a级毛片| 99re在线观看精品视频| 中文字幕人成人乱码亚洲影| 日韩 欧美 亚洲 中文字幕| 男人舔奶头视频| 嫩草影视91久久| 日韩人妻高清精品专区| 国产精品久久久久久亚洲av鲁大| av天堂在线播放| 天天躁狠狠躁夜夜躁狠狠躁| 男人舔奶头视频| av国产免费在线观看| 欧美在线黄色| 国产高清激情床上av| 亚洲国产看品久久| 在线永久观看黄色视频| 国产精品免费一区二区三区在线| 美女高潮喷水抽搐中文字幕| 18禁国产床啪视频网站| 在线免费观看的www视频| 国产精品免费一区二区三区在线| 在线免费观看的www视频| 五月玫瑰六月丁香| 日韩欧美三级三区| 国产精品爽爽va在线观看网站| 制服丝袜大香蕉在线| 欧美激情在线99| 69av精品久久久久久| 久久亚洲真实| 69av精品久久久久久| 久久久久九九精品影院| 欧美乱妇无乱码| 中国美女看黄片| 99精品久久久久人妻精品| 可以在线观看的亚洲视频| 国产美女午夜福利| 在线看三级毛片| 色视频www国产| 色综合欧美亚洲国产小说| 法律面前人人平等表现在哪些方面| 亚洲人成网站高清观看| 精品一区二区三区av网在线观看| 国产欧美日韩精品一区二区| 啦啦啦韩国在线观看视频| 啦啦啦韩国在线观看视频| 搞女人的毛片| 亚洲国产高清在线一区二区三| 国产精品国产高清国产av| 国产精品一区二区三区四区免费观看 | 欧美大码av| 在线播放国产精品三级| 每晚都被弄得嗷嗷叫到高潮| 999精品在线视频| 日韩免费av在线播放| 日韩成人在线观看一区二区三区| av国产免费在线观看| 国产高清激情床上av| 国产精品一区二区三区四区久久| 黄色视频,在线免费观看| 精品国产乱码久久久久久男人| 制服丝袜大香蕉在线| 少妇裸体淫交视频免费看高清| 久久久国产精品麻豆| 搡老妇女老女人老熟妇| 欧美xxxx黑人xx丫x性爽| 中亚洲国语对白在线视频| 久久亚洲真实| 国产精华一区二区三区| 三级毛片av免费| 欧美黄色淫秽网站| 精品久久久久久久久久久久久| 日韩欧美在线乱码| 国产精品一区二区精品视频观看| 一区二区三区激情视频| ponron亚洲| 亚洲国产精品sss在线观看| 国产免费av片在线观看野外av| 淫秽高清视频在线观看| 麻豆av在线久日| 亚洲欧美日韩高清专用| 国产精品亚洲av一区麻豆| 亚洲人成网站高清观看| 91九色精品人成在线观看| 一个人看的www免费观看视频| 999久久久国产精品视频| 热99在线观看视频| 日韩三级视频一区二区三区| 欧美性猛交╳xxx乱大交人| 一个人免费在线观看电影 | 中文资源天堂在线| 日韩国内少妇激情av| 黄色成人免费大全| 国产三级中文精品| 色吧在线观看| 欧美日韩乱码在线| 午夜福利视频1000在线观看| 欧美日韩亚洲国产一区二区在线观看| 精品一区二区三区av网在线观看| 91在线观看av| 国产午夜福利久久久久久| 神马国产精品三级电影在线观看| 国产野战对白在线观看| 久久中文看片网| 日本a在线网址| 1024手机看黄色片| 99久国产av精品| 国产伦精品一区二区三区四那| 成人国产综合亚洲| 久久久久亚洲av毛片大全| 国产欧美日韩一区二区三| 免费观看人在逋| 好男人在线观看高清免费视频| 亚洲精品456在线播放app | 一级毛片高清免费大全| av福利片在线观看| 在线视频色国产色| 亚洲18禁久久av| 69av精品久久久久久| 在线免费观看不下载黄p国产 | 99久久99久久久精品蜜桃| 中文亚洲av片在线观看爽| 国产精品影院久久| 亚洲av免费在线观看| 免费大片18禁| 欧美乱妇无乱码| 99久国产av精品| 国产精品电影一区二区三区| 亚洲专区字幕在线| 免费在线观看日本一区| 亚洲黑人精品在线| 午夜精品久久久久久毛片777| 99热这里只有精品一区 | 欧美三级亚洲精品| 人人妻人人澡欧美一区二区| 久久久久久国产a免费观看| 久久亚洲真实| 国产成人aa在线观看| 国产精品一区二区免费欧美| 岛国在线观看网站| 亚洲真实伦在线观看| 少妇丰满av| 国产欧美日韩精品亚洲av| 又大又爽又粗| 巨乳人妻的诱惑在线观看| 村上凉子中文字幕在线| 哪里可以看免费的av片| 天天躁日日操中文字幕| 免费看日本二区| 久久伊人香网站| 国产亚洲精品综合一区在线观看| 亚洲最大成人中文| av天堂在线播放| 亚洲精品美女久久久久99蜜臀| a级毛片a级免费在线| 99国产精品一区二区蜜桃av| 午夜精品久久久久久毛片777| 人妻久久中文字幕网| 麻豆一二三区av精品| 在线a可以看的网站| 亚洲一区二区三区色噜噜| 好男人电影高清在线观看| 亚洲精品在线观看二区| 97超视频在线观看视频| 国内毛片毛片毛片毛片毛片| 99在线视频只有这里精品首页| 国产综合懂色| 日本成人三级电影网站| av欧美777| 亚洲性夜色夜夜综合| 90打野战视频偷拍视频| 两性夫妻黄色片| 法律面前人人平等表现在哪些方面| 不卡av一区二区三区| 亚洲av成人不卡在线观看播放网| 国产一区二区在线观看日韩 | 999久久久国产精品视频| 在线永久观看黄色视频| 欧美zozozo另类| 99热6这里只有精品| 真人做人爱边吃奶动态| 三级国产精品欧美在线观看 | 欧美大码av| 国产精品一区二区精品视频观看| 一区福利在线观看| 黄频高清免费视频| 午夜激情欧美在线| 观看免费一级毛片| 亚洲国产高清在线一区二区三| 一区福利在线观看| 91在线精品国自产拍蜜月 | 国产淫片久久久久久久久 | 激情在线观看视频在线高清| 久久精品人妻少妇| 长腿黑丝高跟| 久久午夜亚洲精品久久| av黄色大香蕉| 一个人看的www免费观看视频| 亚洲欧美日韩无卡精品| 久久久久久大精品| 别揉我奶头~嗯~啊~动态视频| 国产高清三级在线| 欧美日韩综合久久久久久 | 神马国产精品三级电影在线观看| 亚洲熟妇中文字幕五十中出| 99在线视频只有这里精品首页| 国产精品电影一区二区三区| 三级男女做爰猛烈吃奶摸视频| 国产又色又爽无遮挡免费看| 日本免费a在线| 最新在线观看一区二区三区| 成人精品一区二区免费| 亚洲一区二区三区色噜噜| 国产三级黄色录像| 久久这里只有精品19| 不卡av一区二区三区| 亚洲片人在线观看| 日本a在线网址| 欧美3d第一页| 亚洲中文日韩欧美视频| 嫩草影视91久久| 久久久久久久久久黄片| 99riav亚洲国产免费| 老司机午夜福利在线观看视频| 在线观看舔阴道视频| 久久国产乱子伦精品免费另类| 欧美性猛交╳xxx乱大交人| 欧美黄色淫秽网站| 欧美乱码精品一区二区三区| 日韩欧美 国产精品| 精品乱码久久久久久99久播| 黄片大片在线免费观看| 男女那种视频在线观看| 动漫黄色视频在线观看| 岛国视频午夜一区免费看| 久久精品91无色码中文字幕| 99热这里只有精品一区 | 国产精品影院久久| 亚洲熟妇熟女久久| 最近最新免费中文字幕在线| 色综合欧美亚洲国产小说| 色老头精品视频在线观看| 美女午夜性视频免费| 两个人视频免费观看高清| 精品无人区乱码1区二区| 日日夜夜操网爽| 午夜福利在线观看免费完整高清在 | 国产精品亚洲美女久久久| 欧美+亚洲+日韩+国产| 免费高清视频大片| 亚洲精品久久国产高清桃花| 国产精品久久电影中文字幕| 欧美激情久久久久久爽电影| 国产伦精品一区二区三区视频9 | 黄频高清免费视频| 热99re8久久精品国产| 99久久精品一区二区三区| 日本免费一区二区三区高清不卡| www.999成人在线观看| 男插女下体视频免费在线播放| 国产aⅴ精品一区二区三区波| 熟女电影av网| 亚洲精品色激情综合| 青草久久国产| 日日摸夜夜添夜夜添小说| 午夜福利在线在线| 偷拍熟女少妇极品色| 成人亚洲精品av一区二区| 亚洲精品一卡2卡三卡4卡5卡| 国产蜜桃级精品一区二区三区| 亚洲成av人片免费观看| 免费观看精品视频网站| 伊人久久大香线蕉亚洲五| 亚洲av电影在线进入| 成年女人毛片免费观看观看9| 国产黄片美女视频| 一区福利在线观看| 嫩草影院入口| 中文字幕最新亚洲高清| 亚洲电影在线观看av| 亚洲精品粉嫩美女一区| 夜夜夜夜夜久久久久| 757午夜福利合集在线观看| 久久亚洲精品不卡| 这个男人来自地球电影免费观看| 精品国产乱子伦一区二区三区| av在线蜜桃| 国产一区二区三区在线臀色熟女| 成人鲁丝片一二三区免费| 亚洲欧美精品综合一区二区三区| 色吧在线观看| 在线观看免费午夜福利视频| 国产亚洲精品av在线| 美女午夜性视频免费| 欧美xxxx黑人xx丫x性爽| 国产 一区 欧美 日韩| 日韩人妻高清精品专区| 日韩av在线大香蕉| 亚洲在线自拍视频| 欧美+亚洲+日韩+国产| 国产黄色小视频在线观看| 日本撒尿小便嘘嘘汇集6| 真实男女啪啪啪动态图| 1000部很黄的大片| 亚洲激情在线av| 欧美在线一区亚洲| 一本精品99久久精品77| 好男人电影高清在线观看| 我要搜黄色片| 丰满人妻熟妇乱又伦精品不卡| 国产精品综合久久久久久久免费| 精品久久久久久久毛片微露脸| 婷婷丁香在线五月| 欧美三级亚洲精品| 国产精品日韩av在线免费观看| 精品99又大又爽又粗少妇毛片 | 嫩草影视91久久| 制服丝袜大香蕉在线| 久久热在线av| 精品电影一区二区在线| 国产成人精品无人区| 国内精品久久久久精免费| 美女大奶头视频| 一本久久中文字幕| 真实男女啪啪啪动态图| 久久久久国产精品人妻aⅴ院| 色尼玛亚洲综合影院| 男人的好看免费观看在线视频| 一夜夜www| 欧美日韩中文字幕国产精品一区二区三区| 丰满人妻一区二区三区视频av | 在线观看免费视频日本深夜| 国产成年人精品一区二区| 非洲黑人性xxxx精品又粗又长| 国产精品98久久久久久宅男小说| 国产精品爽爽va在线观看网站| 一本精品99久久精品77| 午夜免费成人在线视频| 亚洲aⅴ乱码一区二区在线播放| 怎么达到女性高潮| 国产男靠女视频免费网站| 中国美女看黄片| 国产熟女xx| 天堂√8在线中文| 一卡2卡三卡四卡精品乱码亚洲| 午夜福利在线观看免费完整高清在 | 一进一出好大好爽视频| 中文字幕精品亚洲无线码一区| 老司机在亚洲福利影院| 亚洲性夜色夜夜综合| 天堂影院成人在线观看| 亚洲精品在线观看二区| 中亚洲国语对白在线视频| 一级毛片精品| 1024香蕉在线观看| 首页视频小说图片口味搜索| 在线观看66精品国产| 久久中文字幕人妻熟女| 级片在线观看| 村上凉子中文字幕在线| 国产成人影院久久av| 亚洲午夜理论影院| 色吧在线观看| 性欧美人与动物交配| 夜夜看夜夜爽夜夜摸| 国产精品精品国产色婷婷| 黑人欧美特级aaaaaa片| 亚洲国产精品sss在线观看| 99在线视频只有这里精品首页| 网址你懂的国产日韩在线| 亚洲精品美女久久久久99蜜臀| 婷婷丁香在线五月| 在线观看66精品国产| 俄罗斯特黄特色一大片| 高潮久久久久久久久久久不卡| 亚洲自拍偷在线| 美女cb高潮喷水在线观看 | 亚洲精品美女久久av网站| 欧美成人免费av一区二区三区| 一进一出抽搐gif免费好疼| 性色av乱码一区二区三区2| 亚洲在线自拍视频| 国产精品久久久人人做人人爽| 老熟妇仑乱视频hdxx| 每晚都被弄得嗷嗷叫到高潮| 69av精品久久久久久| 亚洲国产日韩欧美精品在线观看 | 99精品欧美一区二区三区四区| 狂野欧美激情性xxxx| 大型黄色视频在线免费观看| 18美女黄网站色大片免费观看| 欧美乱码精品一区二区三区| 男女床上黄色一级片免费看| 成人亚洲精品av一区二区| 久久精品人妻少妇| 美女cb高潮喷水在线观看 | 香蕉丝袜av| 啦啦啦观看免费观看视频高清| 亚洲片人在线观看| 久久精品国产99精品国产亚洲性色| 精品电影一区二区在线| 给我免费播放毛片高清在线观看| 手机成人av网站| 亚洲中文av在线| 亚洲av成人精品一区久久| 国产精品av视频在线免费观看| 精品一区二区三区四区五区乱码| 国产一区二区三区视频了| 日日干狠狠操夜夜爽| 国产亚洲精品av在线| 床上黄色一级片| 亚洲在线自拍视频| 久久99热这里只有精品18| 中出人妻视频一区二区| 日韩人妻高清精品专区| 久久性视频一级片| 免费av不卡在线播放| 国产成人精品久久二区二区91| 亚洲午夜理论影院| 亚洲精品国产精品久久久不卡| 夜夜躁狠狠躁天天躁| 色吧在线观看| www.熟女人妻精品国产| 欧美在线黄色| 成年免费大片在线观看| 亚洲专区中文字幕在线| 成年人黄色毛片网站| 小说图片视频综合网站| a级毛片在线看网站| 日韩国内少妇激情av| 偷拍熟女少妇极品色| 真人做人爱边吃奶动态| 一二三四社区在线视频社区8| 中文在线观看免费www的网站| 国产主播在线观看一区二区| 99久久精品国产亚洲精品| 久久精品人妻少妇| av女优亚洲男人天堂 | 国产精品 国内视频| 国产精品乱码一区二三区的特点| 精品久久久久久久毛片微露脸| 日韩大尺度精品在线看网址| 香蕉av资源在线| 蜜桃久久精品国产亚洲av| 久久久国产欧美日韩av| 国产成人aa在线观看| 999久久久精品免费观看国产| 亚洲五月天丁香| 国产精品国产高清国产av| 叶爱在线成人免费视频播放| 波多野结衣巨乳人妻| 长腿黑丝高跟| 亚洲欧美日韩高清在线视频| 丰满的人妻完整版| 日韩欧美在线二视频| 国产不卡一卡二| 亚洲欧美精品综合久久99| 久久久久亚洲av毛片大全| 国产成人av教育| 一进一出抽搐gif免费好疼| 亚洲成av人片免费观看| 国产精品久久久久久久电影 | 欧美xxxx黑人xx丫x性爽| 日韩精品中文字幕看吧| 一夜夜www| 国产精品99久久久久久久久| 极品教师在线免费播放| 日韩欧美在线乱码| 中文字幕熟女人妻在线| 亚洲 欧美一区二区三区| 天堂动漫精品| 久久久精品大字幕| 色吧在线观看| 最新美女视频免费是黄的| 18禁国产床啪视频网站| 亚洲精品乱码久久久v下载方式 | 国产又色又爽无遮挡免费看| 两人在一起打扑克的视频| 好男人电影高清在线观看| 国产欧美日韩一区二区精品| 国产午夜福利久久久久久| 禁无遮挡网站| 国产精品久久久久久久电影 | 青草久久国产| 日本黄色视频三级网站网址| 免费看十八禁软件| 国产视频一区二区在线看| 午夜福利在线观看吧| 久久国产精品影院| 亚洲精品久久国产高清桃花| 亚洲欧美精品综合一区二区三区| 久久亚洲真实| 日本成人三级电影网站| 国产成人影院久久av| 国产极品精品免费视频能看的| 禁无遮挡网站| 国产亚洲av高清不卡| 18禁裸乳无遮挡免费网站照片| 99视频精品全部免费 在线 | 中文资源天堂在线| 午夜亚洲福利在线播放| 国产亚洲精品一区二区www| 我要搜黄色片| 国产伦在线观看视频一区| 黄色视频,在线免费观看| 国产精品99久久久久久久久| 成人午夜高清在线视频| 激情在线观看视频在线高清| 天天躁日日操中文字幕| 欧美日韩瑟瑟在线播放| 久久亚洲真实| 亚洲精品在线美女| 亚洲色图 男人天堂 中文字幕| 国产极品精品免费视频能看的| 久9热在线精品视频| 一区二区三区高清视频在线| 午夜激情福利司机影院| 午夜福利18| 99久久综合精品五月天人人| 在线观看一区二区三区| 制服人妻中文乱码| 国产三级黄色录像| 欧美高清成人免费视频www| 国产伦人伦偷精品视频| 欧美精品啪啪一区二区三区| 亚洲国产欧美一区二区综合| 国产亚洲精品久久久久久毛片| 欧美色视频一区免费| 久久久久亚洲av毛片大全| 又爽又黄无遮挡网站| 国产真人三级小视频在线观看| www.熟女人妻精品国产| 免费高清视频大片| 黑人操中国人逼视频| 国产成人精品无人区| 99视频精品全部免费 在线 | 成人特级黄色片久久久久久久| 别揉我奶头~嗯~啊~动态视频| 女人高潮潮喷娇喘18禁视频| 亚洲avbb在线观看| 精品国产三级普通话版| 久久热在线av| 久久香蕉国产精品| 超碰成人久久| 成人无遮挡网站| 99久久久亚洲精品蜜臀av| 精品一区二区三区四区五区乱码| 亚洲欧美日韩无卡精品|