Numerical Simulation Of Three-Dimensional Unsteady Flow Past Ice Crystals

15 SEPTEMBER 19972263WANG AND JIFIG. 2. The general configuration of the theoretical problem. Icecrystals are assumed to fall with their long axes in the horizontaldirection.FIG. 1. The three types of ice crystals considered in this study.by finite circular cylinders), hexagonal ice plates, andbroad-branch crystals. Figure 1 shows a schematicsketch of these three types of crystals. The quantity arepresents the ‘‘radius’’ of the ice crystals as defined inthe figure. We shall also assume that these ice crystalsfall with their broad dimensions oriented in the horizontal direction, which is known to be the common fallorientation of many medium-sized ice crystals (Pruppacher and Klett 1978). It is known that large ice crystals also exhibit zigzag fall attitude, but this is not simulated here due to the limitation of computer resource.The schematic configuration of the theoretical problemconsidered here is shown in Fig. 2.To facilitate the numerical analysis, we first introducethe dimensionless variablesxx9 5 ,aP9 5V9 5P,rzV z 2V,zV zRe 5t9 5V50at the surface of the ice crystal,V 5 1·e z at infinity,and(4)(5)where ez is the unit vector in the general flow direction.In real numerical computations, of course, the domainis always finite and the condition (5) can only be takento mean that the velocity is constant at an outer boundary that is sufficiently far away from the crystal. It isdifficult at present to determine from purely theoreticalgrounds how far the distance should be in order to be‘‘sufficiently far.’’ We did this by trial and error, andthe outer boundary is considered far enough when thecomputed results do not change within a few percentsas we move the boundary out. Similar treatment wastzV z,a2zV za,n(1)where x (or y, z) is one of three Cartesian coordinates,V is the fluid velocity, V the free-stream velocity equalto the terminal fall velocity of the ice crystal, P thedynamic pressure, and n the kinematic viscosity of thefluid. Here, Re is the Reynolds number relevant to theflow. All primed quantities are nondimensional. Usingthese dimensionless variables, we can write down thenondimensional Navier–Stokes equation and the continuity equation as (after dropping the primes)]V2 21 V· V 5 2¹P 1¹ V,]tRe ·V 5 0.(2)(3)The ideal boundary conditions appropriate for thepresent problems areFIG. 3. The initial perturbation imposed on the steady flow fieldin order to generate time-dependent flow behavior. The magnitude ofperturbation in region A and B is 30% of the free-stream velocity,but the directions of the perturbations are opposite.

2264JOURNAL OF THE ATMOSPHERIC SCIENCESVOLUME 54TABLE 1. Outer boundaries of the computational domains for thethree crystal cases. The radius is 1.BoundariesColumnarcrystalUpstreamLateral gonal 0915.27condition at the downstream boundary is replaced by aweaker condition, ]V/]z 5 0.The pressure field can be determined from the Navier–Stokes equation at all boundaries except at thedownstream boundary where the condition ]P/]z 5 0is used. Since we are dealing with unsteady flow here,we also need initial conditions to close the equations.The initial conditions are P 5 0 and V 5 1·ez everywhere at t 5 0 except at the surface of the crystal. Thecondition on the surface is V 5 0 (nonslip condition)at all t.c. Generation of unsteady flow featuresFIG. 4. The nonuniform grid used for numerically solving the Navier–Stokes equations for flow past a columnar ice crystal. (a) Broadside view. (b) End view.done for all outer boundaries. Table 1 shows the locations of the upstream, downstream, and lateral boundaries for determining the numerical flow fields for thethree cases.While condition (5) is approximately valid in the upstream and lateral boundaries, it is usually not valid inthe downstream boundaries. This is due to the fact thatat higher Reynolds number ranges as investigated here,the shedding of eddies may occur. The disturbances often propagate downstream for a long distance. Thus theAlthough the Navier–Stokes equation (2) is writtenas a time-dependent equation, this does not mean thatthe computational results will always result in timedependent flow features such as the shedding of eddies.Indeed, Dennis and Chang (1970) has shown that forflow past two-dimensional cylinders starting with symmetric initial conditions, the eddy shedding does notoccur even at high Reynolds numbers. In order to generate these time-dependent, or unsteady, features, it isnecessary to implement an asymmetric initial perturbation field. There are many ways of implementing thisperturbation. For example, Braza et al. (1986) achievedthis on a two-dimensional cylinder by performing a rotation of the cylinder along its axis. In the present study,we achieve this by implementing a velocity perturbationof magnitude 0.3V in the downstream region immediately behind the crystal to the steady-state solutionsas shown in Fig. 3. The directions of the perturbationare opposite to each other in regions A and B, so as toform a shear along the central plane of the flow. As weshall see later, at high enough Reynolds numbers, thisperturbation will generate a periodic eddy-shedding pattern in the simulated flow. On the other hand, the perturbation will be damped out in a short time if the Reynolds number is low.4. The numerical schemeTo solve Eqs. (2) and (3) with the appropriate initialand boundary conditions, we adopt a numerical approach utilizing the finite difference method. It is alsonecessary to set up a mesh grid. Due to the more complicated shapes of these ice crystals, it is decided thatthe simplest way to set up grids is to use the Cartesian

15 SEPTEMBER 19972265WANG AND JITABLE 2. Dimensions of columnar ice crystals treated in thepresent study. Units: dimensionless.ReDiameter (d)Length 2.66TABLE 3. Dimensions of hexagonal ice plates treated in the presentstudy. Units: dimensionless.coordinate system. In order to prescribe the inner boundary conditions with adequate precision, the grid spacingnear the crystal surface has to be small. On the otherhand, the grid spacing far from the crystal can be largerto save computing time. This results in nonuniform gridsused in the present study, as shown in Fig. 4.As indicated before, the primitive velocity formulation of the Navier–Stokes equation is adopted for thisstudy. The velocity at each time step is obtained by apredictor–corrector method. First, the velocity predictorV* is determined by solving the equationV* 2 V n2 2 n1 (V n· )V n 5¹ V,DtRe(6)where V is the velocity solved at time step n and Dtis the time increment. The pressure at time step n 1 1is then given byn¹ 2 P n11 5 ·V*Dt(7)(e.g., Peyret and Taylor 1983). Finally, the velocity attime step n 1 1 is determined byV n11 2 V*5 2 P n11 .Dt(8)The scheme of velocity interpolation at each time stepis the modified QUICK (Quadratic Upstream Interpolation for Convective Kinematics) scheme with secondorder accuracy developed by Leonard (1979) and extended by Davis (1984) and Freitas et al. (1985). Sincewe use the nonuniform mesh grid here, the uniform gridformulation given by Leonard (1979) cannot be used.Instead, we adopted the formulation of Freitas et al.(1985) for nonuniform grid but modified to suit ourthree-dimensional flow cases. The details of the formulation and the stability criterion are given in Ji andWang (1990, 1991).The Poisson equation for pressure, Eq. (7), is solvedby the standard successive over-relaxation (SOR) method as described in Peyret and Taylor (1983) and Anderson et al. (1984).The time step Dt used in the integration varies fromReDiameter (d)Thickness .02880.015 to 0.03, depending on the local grid spacing, suchthat the stability criterion is satisfied. The smallest gridspacing was Dx equals 0.0775. The largest grid size usedwas 59 3 75 3 89. Naturally, a larger grid size willresult in better accuracy but will increase the computing time considerably. The typical computing time for10 000 time steps is on the order of a few hours on aCRAY X/MP computer. The computation on a CRAY-2computer is somewhat faster. It appears that the SORscheme in solving the pressure equation is the mainbottleneck of the computation. The grid size used in thisstudy represents a compromise between accuracy andavailable computing resources.5. Results and discussionThe size, aspect ratios, and Reynolds numbers of thecolumns, hexagonal plates, and broad-branch crystalsare listed in Tables 2, 3, and 4, respectively. Their dimensions are chosen to overlap those adopted by someprevious work (Schlamp et al. 1975; Pitter et al. 1973,1974; Pitter 1977; Miller and Wang 1989) so that theresults can be compared. In the following we shall discuss the results for each crystal type separately.a. General features of the flow fields around fallingcolumnar ice crystalsAs mentioned earlier, this type of crystal is approximated by a circular cylinder of finite length. Becauseof the finite length, the cylinder, as well as the flow fieldTABLE 4. Dimensions of broad-branch crystals treated in thepresent study. Units: dimensionless.ReDiameter (d)Thickness 0.0750.070.04570.0400.0330.0300.0260.0235

2266JOURNAL OF THE ATMOSPHERIC SCIENCESVOLUME 54FIG. 5. Streak pattern, or ‘‘snapshot’’ field, of massless tracer particles for flow past an ice column. (a) Broad-side view, Re 5 40. (b)End view, Re 5 40. (c) Experimental photograph of a falling short cylinder at Re 5 40. (d) Broad-side view, Re 5 70. (e) End view, Re5 70. (f) Experimental photograph of a falling short cylinder at Re 5 70 (photo courtesy of Dr. K. O. L. F. Jayaweera).

15 SEPTEMBER 1997WANG AND JI2267FIG. 6. Drag coefficients for flow past cylinders of various d/l ratios. The solid line and thetriangles are for infinitely long cylinders. It is seen that as the Re increases, the cylinder becomeslonger and the drag coefficient becomes closer to that of the infinitely long cylinder.FIG. 7. Velocity fields of flow past a hexagonal ice plate at Re 5 2. The cross-section location is indicated by the dashed line in theupper right corner of the figures. Vectors represent projections of the 3D vectors onto that cross-sectional plane.

2268JOURNAL OF THE ATMOSPHERIC SCIENCESFIG. 8. Velocity fields of flow past a hexagonal ice plate at Re 5 20.The cross-section location is indicated by the dashed line in the upperright corner of the figures. Vectors represent projections of the 3D vectorsonto that cross-sectional plane.VOLUME 54

15 SEPTEMBER 19972269WANG AND JIwhen the Reynolds number increases to 70, the flowpast the same cylinder developed a full-fledged periodiceddy shedding. This indicates that the current schemeis capable of simulating the eddy shedding phenomenoncorrectly. Thus no perturbation was necessary for computing the flow fields with Re # 50.Since many of the features of the flow past finitecylinders have been discussed in one of our previouspapers (Ji and Wang 1991), only a brief summary andsome highlights are given here. It is sufficient to saythat the numerical results reproduced the features of flowfields observed in the laboratory experiments of Jayaweera and Mason (1965), such as the pyramidal wakeregion for Re # 40 and the shedding of eddies for flowin greater Reynolds number range. Figure 5 gives apictorial comparison between the computed masslesstracer streaks and experimental photographs of Jayaweera and Mason (1965). It is seen here that there aregood similarities between the two sets of pictures.Figure 6 shows the comparison between the computed drag coefficients with those obtained by othertheoretical and experimental results. The drag coefficient is defined asCD 5FIG. 9. The unsteady flow field past a hexagonal ice plate at Re 5140.around it, is no longer cylindrically symmetric. The dimensions and aspect ratios of the cylinders chosen forthe computation are shown in Table 2 and are the sameas those given by Schlamp et al. (1975) for Reynoldsnumbers between 0.2 and 20, and by Jayaweera andMason (1965) for Reynolds numbers 40 and 70. Severalhigher Reynolds numbers cases were also computed forthe purpose of checking, but the details of these willnot be discussed here. The aspect ratios of the cylindersspecified by Schlamp et al. (1975) are taken from theactual samples whose diameter-length relations were reported by Auer and Veal (1970). For the higher Reynolds number cases, the Jayaweera and Mason’s (1965)cylinder dimensions are also used here because theirresults are the only experimental data available for verification purpose for free-falling finite cylinders. In allcases, the ice columns become longer as compared tothe diameters as the Reynolds number increases.Previous experimental studies of both two- and threedimensional flow past circular cylinders indicated thatthe flow remains steady up to Re ø 50 (e.g., Kovasznay1949; Jayaweera and Mason 1965). In Ji and Wang(1990, 1991), it was shown that even with the 0.3V perturbation, the periodic shedding of eddies did notdevelop in the simulated flow field when Re 5 40. Instead, it just produced a transient disturbance that quickly dissipated in about 45 time steps. On the other hand,D,rV 2 a(9)where D is the drag force and a is one-half of the crosssectional area of the cylinder normal to the flow direction. Obviously, the drag coefficients of the present results differ from the results for infinite long cylinders.The difference is greater for smaller Reynolds number.This is due to the fact that the dimensions of columnsfor lower Reynolds numbers are such that their shapesdiffer more from an infinitely long cylinder. On the otherhand, columns for higher Reynolds numbers are closerto the shape of an infinitely long cylinders, hence theirdrag coefficients are closer to each other. This impliesthat the theoretical results obtained by previous investigators regarding the behavior of columnar ice crystalsare probably reasonable for the case of larger ice columns but not for short columns. The collection efficiencies of small droplets by columns and ventilationfactors of falling ice columns computed based on ournew scheme did show such trends and will be reportedelsewhere. The discrepancies exist mainly for Re , 10.It is seen here that the coefficients are practically thesame as that of infinitely long cylinders for Reynoldsnumber . 10. The drag coefficients calculated here canbe fitted by the expressionlog 10 CD 5 2.44389 2 4.21639A 2 0.20098A 21 2.32216A 3 ,(10)log 10Re 1 1.0.3.60206(11)whereA5

2270JOURNAL OF THE ATMOSPHERIC SCIENCESVOLUME 54FIG. 10. Velocity fields of flow past a broad-branch ice crystal at Re 5 2. The cross-section location is indicated by the dashed line inthe upper right corner of the figure.This formula is valid within the range 0.2 , Re , 100.It fits the computed data to within a few percent. Notethat in reality the drag coefficient is also a function ofthe aspect ratio of the cylinder, which is not explicitlyrepresented in Eq. (10); hence, strictly speaking, this fitis only applicable to those cases indicated in Table 2.But judging from the smooth behavior of this relation,we feel that it is probably applicable to columnar crystals with dimensions satisfying Auer and Veal’s (1970)relations and with flow Reynolds numbers in the aforementioned range. It would be desirable to find a relationof CD as a function of the aspect ratio. However, morecalculations are needed to establish this relation. Thisis not done for the present study due to the constraintof computing resources.b. General features of the flow fields around fallinghexagonal ice crystalsPioneering numerical work on the flow fields aroundfalling ice plates were performed by Pitter et al. (1973),who used thin oblate spheroids to approximate planarice crystals. In the present study, we use the actual hexagonal shape to model the ice plates whose dimensionsare given in Table 3. The corresponding range of theReynolds numbers is from 1 to 120. However, additionalcases were also computed as needed to demonstrate theflow fields.Examples of steady flow fields around hexagonalcrystals are shown in Figs. 7 and 8, which representcases of Reynolds numbers 2 and 20, respectively. Theseare steady flow cases. The flow fields look similar tothose obtained by Pitter et al. (1973).The flow field of Re 5 1 (figure not shown) does notindicate the existence of standing eddies. But there arealready standing eddies formed in the wake region ofthe crystal at Re 5 2. This is consistent with Pitter etal. (1973), who indicated that the eddies start to appearat Re 5 1.5. As expected, the eddies become larger athigher Reynolds numbers.Experiments of Willmarth et al. (1964) showed thatat Re 100, eddy shedding occurs in the downstreamof a falling disk. Such unsteady behavior can be simulated using the same technique as we did for the finitecylinders. Figure 9 shows the simulated unsteady flowfields for flow past hexagonal plates at Re 5 140. Theinitial perturbation introduced (after the steady-state solution has been obtained) was again 0.3V , which has

15 SEPTEMBER 1997WANG AND JI2271FIG. 11. Velocity fields of flow past a broad-branch ice crystal at Re 5 20. The cross-section location is indicated by the dashed line inthe upper right corner of the figure.been proven to be adequate for kicking up the shedding.It can be seen that the flow field is obviously asymmetricdue to the shedding. Detailed analysis of how sheddingstarts has not been done yet, but it is expected that theshedding would start at a particular corner and the pointof detachment would rotate around on the plate.It must be stressed here that the flow fields describedabove are computed assuming that the plate position isfixed with respect to the incoming air flow; that is, theangle between the c axis (normal to the plate basal surface) of the plate and the general flow is kept at 908.In reality, falling plates are known to perform zigzagmotion, which implies that the angle is not constant butis actually a function of time. In order to simulate suchzigzag motions, one has to use very small time stepsfor adequate accuracy. Due to the constraint of computing resources, these cases are not yet simulated herebut are currently being studied by us. However, it isexpected that the above results should give good approximations, especially when the variations of the angles are not large.There seems to be no experimental measurementsavailable for flow properties past hexagonal plates. Willmarth et al.’s (1964) results for flow past circular platesare the closest cases for the comparison purpose. Buthere the comparison is difficult to make because theaspect ratios of the computed and experimental resultsare different. For the same reason the comparison between our present results and those of Pitter et al. (1973)is also difficult to make. The aspect ratio of the platescalculated here varies with the Reynolds number, whereas the thin oblate spheroids in Pitter et al. (1973) havefixed aspect rati

Corresponding author address: Dr. Pao K. Wang, Department of Atmospheric and Oceanic Sciences, University of Wisconsin-Mad-ison, 1225 W. Dayton Street, Madison, WI 53706. E-mail: pao@windy.meteor.wisc.edu of ambient water vapor toward them. The flow fields created by their motion influence the vapor density gra-