A biviscous modified Bingham model of snow avalanche motion by Jimmie Duane Dent A thesis submitted in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY in Civil Engineering Montana State University © Copyright by Jimmie Duane Dent (1982) Abstract: Contained in this thesis is a description of snow avalanche dynamics. This information is gained from the researches of others and from an experimental program that is described within. Tests were carried out in which 2.1 m^3 volumes of snow were decelerated to a stop from 18 m/s on a flat level snow covered runnout. Data concerning leading edge position versus time and flow velocity versus depth were collected. Based upon this information a biviscous modified Bingham flow model is proposed and developed. This flow model consists of a constitutive relation that is linear between the stress and velocity gradients over two distinct regions. There is a region of high viscosity, for low stresses, where little deformation will take place and a second region with a much lower viscosity that operates for stresses above a specified level. This model, in the form of a finite difference computer code, is used to simulate the results of the snow flow tests. In carrying out this modeling, values of the parameters that appear in the biviscous flow law are evaluated.  A B!VISCOUS MODIFIED BINGHAM MODEL OF SNOW AVALANCHE' MOTION by ' - ' JIMMIE DUANE DENT A. thes&s submitted in partial fulfillment of the requirements.for the.degree of DOCTOR OF PHILOSOPHY in Civil Engineering Approved: Head, Major ^aduateDean MQNTANA STATE UNIVERSITY; Bozeman.-, Montana June, 19.82 iii • ' ACKNOWLEDGMENT This work was made possible through grants from the U. S. Depart­ ment of Agriculture, Forest Service (cooperative agreement Nos. Rm-16-778-CA, Rm-16-884-CA, Rm-81-165-CA) and from the Montana, State ,University Computing Center. The support from these.agencies, and in particular Mr. Mario Martinelli, Jr. of the Forest Service and Range Experiment Station and Mr. Robert Motsch and Mr. Phillip Adolf of the M.S.U. Scientific Subsystems computing facility is gratefully appreciated Also this work would not have been possible without the facilities of the Civil Engineering and Engineering Mechanics Department and Mr. Theodore Williams. The author would especially like to thank Dr. Theodore Lang for his help and guidance in the preparation of this thesis, for the many hours spent in the field shoveling snow and for the many hours spent listening and counseling. To all those people who are not mentioned here, but who have participated in some way in this research, I would also like to say thank you. TABLE OF CONTENTS Page M Q N T A /........... if ACKNOWLEDGMENTS . ...................... i;Li LIST OF FIGURES.................... vi ABSTRACT.......... ix CHAPTER I. INTRODUCTION ........................................ I II. A DESCRIPTION OF AVALANCHE MOTION............. 3 Observations ...................................... 3 Measurements............................... .. . . 4 III. HISTORICAL ANALYSIS OF AVALANCHE MOTION. . . . . . . . n Open Channel Flow................................. 12 Center of Mass Motion..................... . . . . 17 Linear Viscous Model ............................. 19 Bingham Model.................. 20 Granular Theories. .... ......................... 24 IV. FIELD TEST PROGRAM ....................... 26 The Test Facility........................... 26 Observations and Measurements............ 29 Window T e s t s ............................... 42 Discussion of Test Results........ '............. 42 V. B!VISCOUS BINGHAM M O D E L ................. 47 Linear Viscous Model ............................. 47 S M A C ............................. 49 Linear Viscous Model Results ..................... 56 Bingham Model................... 58 Biviscous Bingham Model............ 60 BVSMAC ........................................ .. . 68 BiviscouS. Modeling Results. . . . . . . ........ .84 VTable of Contents (Continued) VI. SUMMARY AND CONCLUSIONS. ........... . LIST OF REFERENCES............................. APPENDICES. . . . . ........................... A. BVSMAC Fortran Computer Code......... B. Time Sequence Particle Plots of BVSMAC Calculation of Flowing Snow Test. . . 2 2 2 I. T0/p= 2.2 m /s , V = .002 m /s, 2. T V 1 = 0.10 m /s ............. 3.0 m^/s^, V = .010 m^/s. 0/p V = 0.20 m /s o/p V 1.0 m2/s2, V '0.10 m /s .010 m /s, C. I. V = .001 2. V = .010 3. V = .100 4. V = .200 Page 96 101 106 106 . 132 Viscous Deformation of a Rectangular Block . . . 2 144 154 164 165 168 172 177 vi LIST OF FIGURES Page 1. Profile of an avalanche derived from pressure sensors (Schearer 1980)............................. ■ 9 2. Shear stress of granular snow on an over-riding surface-roughened steel plate versus speed of the plate. (Lang et al. 1981)......................... 23 3. Schematic of field test facility..................... 27 4. Photograph of field test facility................... 28 5. Glass window used to measure velocity profile . . . . . 30 6. Dump geometry, test No. 3-2-23-80 ................... 31 7. Power Cloud . ................. ...................... 33 8. Dump geometry, test No. 2-2-23-80 .................. 35 9. Leading edge position versus time, test No. 2-2-23-80 ......................... .. . . . 36 10. Debris distribution, test No. 2-2-23-80 ............ 37 11. Debris distribution, test No. 3-2-23-80 ............ 38 12. Dump geometry, test No. 1-2-25-80 (Block Test). . . . 40 13. Debris distribution, test No. 1-2-25-80 ............. 41 14. Cross section of debris, chunk test 1-2-25-80 . . . . 43 15. Cross section of debris, chunk test 1-2-25-80 . . . . .44 16. Data points of particle speed versus depth of flow. .. - 4 5 17. 8MAC computing g r i d ............................. 52 18. . G'MftC method of specifying the tangential boundary condition at a rigid wall........................... 54 19. One dimensional form of the Bingham and Biyiscous constitutive laws............................... 63 Vii List of Figures (Continued! Paige 20. Mohr's circles for the rotational.transformation . of the deyiatoric stresses and the velocity gradients................ 67 21. Flow between parallel plates; geometry stress and velocity gradient.............. 71 22. BVSMAC particle plot sequence for the solution to the flow of a Bingham material between parallel plates; T /p = 0.30, V = 0.20, V 1 = 2 . 0 .......... 74 23. Steady state velocity profile, comparison between analytical and numerical solutions ................. 78 24. BVSMAC particld plot, flow between parallel plates at t = 1.0 s; TqZP = 0.30, V = 0.20, V = 4.0. . . . 79 25. BVSMAC particle plot, flow between parallel plates at t = 1.0 s; Tq/P = 0.60, V = .20, V 1 •= 2.0 . . . . 80 26. Constitutive relations for the Biviscous model simulations of the flow between parallel plates .......................................... .. . 81 27. BVSMAC particle plot, flow between parallel plates t = I. s; TqZP = 0.60,V = 0.10, V 1 = 2 . 0 .......... 83 28. Input flow configuration for BVSMAC simulation of snow test 2-2-23-80........ ......................... 86 29. Constitutive relations that adequately modeled the leading edge motion of snow test 2-2-23-80 ........ 88 30. Leading edge position versus time.; comparison between experiment and BVgMAC „ SEUIV R 2.2 in Zs ; v = .002 m /s# V' - 0.10 Itf1Zs)........................... - . .. • . 90 31. Velocity profile comparison between snow test and BVSMAC calculations for various.combinations of BVSMAC parameters .................................... 91 viii List of Figures (.Continued) Page 32. Final debris depth profile; comparison between experiment and' BVSMAC model..................... . . 93 ix ABSTRACT Contained in this thesis is a description of snow avalanche dyna­ mics. This information is gained from the researches of others and from an experimental program that is described within, Tests were carried out in which 2.1 m .volumes of, snow were decelerated to,a stop from'18 m/s on a flat level.snow - covered, runnout. Data concerning, leading edge position versus .time and .flow,velocity versus depth were.collected. Based upon this information a biviscous modified Bingham flow model is proposed and developed. This flow model consists of a con^ stitutive relation that is linear between the stress and velocity gradients over two distinct regions. There is a.region.of high vis­ cosity, for low stresses, where little deformation will take place and a second region with a much lower viscosity that operates for stresses above a specified level. This model, in the form of a finite difference computer code, is used to simulate the results of the snow flow tests. In carrying out this modeling, values of the parameters that appear in the biviscous flow law are evaluated. CHAPTER I INTRODUCTION ' Snow avalanches represent a significant hazard in the mountainous regions of the world. Every year many people are killed and millions of dollars of damage are done by avalanches. The United States alone, averages seven deaths and $300,000 in property destruction annually. CPerla and Martinelli, 1975). The European Alpine countries, because of their much more densely settled mountain regions have a much more severe problem. Single avalanches have taken as many as 96 lives in the United States and an avalanche consisting of ice, snow,and earth completely buried a village in the Peruvian Andes in 1970, killing over 20,000 people. As activity increases in the mountainous regions of the world, the potential for avalanche destruction mounts. It thus becomes increasingly important to be able to identify avalanche hazards. Avalanche paths and their runout zones must be mapped in those regions people frequent, so that construction in hazardous areas can be mini­ mized. In those instances where building in an avalanche path cannot be avoided, for example, roads and bridges, engineering design criteria must be established to minimize potential danger. This can take the form of construction of defensive sttructures for protection, or designing the primary structure itself to withstand.anticipated ava­ lanche forces. In either case, it is necessary to be able to predict avalanche motion and the forces that the avalanche would exert on the 2structure under consideration. It is these concerns that motivate the study of avalanche dynamics. Through an understanding of the mechanics of flowing snow, predictive models may be developed that will allow engineers to build safely in avalanche hazard areas. Contained in this thesis is a description of avalanche mechanics. Much of this information is observational in nature, since great dif­ ficulty is encountered in attempting to instrument and measure full scale avalanche motion. From this understanding of avalanches, a flow model is proposed and developed. This model, based upon the classical equations of continuum mechanics, may be used to predict avalanche motion and impact forces. Analytical data on moving snow are scarce; therefore, to test the model, it was necessary to institute a field program for the collection of data. The field data were used to help verify and then to calibrate the flow model. This provided a source of analytic, as well as qualitative information on flowing snow that has led to a better understanding of the mechanics of snow flow. The model may now be used to predict avalanche motion and impact forces. CHAPTER IJ A DESCRIPTION OF AVALANCHE MOTION ,Observations A snow avalanche can range from a few snow particles-.rolling and bouncing down a slope to a massive assemblage of ice, snow, and earth., moving with enough velocity to destroy virtually anything in its path.. The majority of avalanches are of the former type; harmless sluffs of small volumes of snow that occur as new snow is deposited on mountain slopes. Large avalanches occur much less often and only achieve de- 3 structive potential when the volume of snow reaches about 1000 m . These destructive avalanches typically start high on a mountain side in snow catchment areas called starting zones. As snow is deposited on these 30 to 45 degree slopes, the snow load eventually exceeds the strength of the bonds holding the snow in place. Volumes of snow up- 5 3wards of 10 m can then be released at one single time. The released snow may be made up,of new cohesionless snow that has recently been deposited, in which case the resulting avalanche is a loose snow ava­ lanche. If the snow has resided on the ground for a period of time, it will have sintered into a solid matrix. When it then starts to move, this slab avalanche initially moves as a large rigid block that im­ mediately begins to break into pieces. Finally, if the snow pack has warmed and there is free water in the snow when.it is released, a wet snow or slush avalanche results. Combinations and variations of these ' 4avalanches, of course, also may be observed. Once initiated the ,avalanche gains speed as it flows down the avalanche track. The blocks that may have been present at release begin to break up. They tumble and collide with each other rounding corners, and reducing their sizes. The smaller snow particles can be thrown into the air by collisions with other particles, creating a dust cloud that may completely surround the flow. It is even possible for loose snow avalanches to become entirely airborne, with all the snow being air. entrained so that there is no central flowing core. However, the majority of avalanches consist of a dense flowing core that contains most of the momentum and is responsible for the high damage potential. This core is surrounded by a large diffuse dust cloud that makes obser­ vation of the flowing core difficult. As.a large avalanche continues to move down the mountain its velocity may reach upwards of 50 m/s. Typical flow speeds however, are more in the range of 20-40 m/s. The avalanche keeps flowing until the slope angle decreases to a point where gravity can no longer sustain the motion or until some obstacle is encountered. .Measurements The description.of avalanche motion.involves specifying such parameters as velocity, density, and volume, as fundamental.variables. These quantities must be known as functions of space and time, to totally define the flowing system. But, because of the formidable nature of 5snow .avalanches and their hostile environment, even the simplest measurements have proven very difficult to perform.. As a result, much of the information about avalanches has been gathered., by examining the debris after it has. stopped. Mass estimates, type of snow, and particle sizes can.all be ascertained. The slope profile can be measured and simple energy analysis used to determine average friction values and average flow velocities (Korner 1981). Data.collection on moving avalanches was at first limited to eye witness accounts and their estimates of speed and length of time that the avalanche ran. A few avalanches were timed to yield average velocity values. Only in recent years have attempts been made to collect time dependent data. These efforts have mostly centered on the measurement of avalanche velocities, of which most measurements have been only of the leading edge of the avalanche. For example, Shoda (1966), LaChapelle and Lang (1979), and Martinelli et al. (1980) have all used motion picture photography to gather this information. They filmed avalanches as they descended a slope and then analyzed the film to determine the number of frames necessary for the front of the avalanche to traverse known distances between landmarks- This is a tedious procedure that is frought with many difficulties. Among them are parallax errors, badly spaced or poorly defined landmarks, and obscurement of the flowing mass by the powder dust cloud. Photogrammetric techniques eliminate the parallax and landmark difficulties at the expense of more effort for 6data collection and reduction. Bruikhanoy (.1967). in the Soviet Union and Brugnot (1980) in France have successfully employed this method tp measure the time-dependent leading edge velocity for several full scale ,avalanches. There is only in existence then a very limited, number of data sets providing avalanche leading edge velocities. Measurements of flow velocities in avalanches at points other than the leading edge are extremely rare. Some measurements have been made by Schearer (1973), Schearer and Salway (1980) and.Shimizu et al. (1973, 1975, 1977, 1981). These researchers mounted sets of pressure transducers in avalanche paths. By correlating the pressure signals observed in separate sensors and the time difference.between the sig­ nals, internal velocities could be calculated. Only sporadic measure­ ments could be made, however, due to the difficulty in finding cor­ relations between sensor signals. A potentially powerful technique for measuring velocities lies in the use of microwave radar. Gubler (1980) has used this technique to measure leading edge velocities and is currently experimenting with the use of active and passive sensors imbedded in the avalanche flow. It is hoped that development of this system will lead to the capability of measuring internal flow velocities. Because velocity data is so hard to gather by field measurements, physical modeling has provided another means of obtaining useful ava­ lanche data. Powder avalanches in particular have lent themselves to 7this type of analysis. In scale modeling, dynamic similarity must be maintained and boundary conditions faithfully reproduced. Dynamic similarity is achieved by matching certain characteristic dimensionless . numbers that involve the flow parameters. in the case of powder ava-.\ . Ianche modeling it is the densimetric Froude number, which is the ratio of inertial forces to gravitational forces, that must be matched. In addition, the Reynolds number, which is the ratio of inertial to vis­ cous forces must be maintained high so that viscous effects are negli­ gible. Hopfinger and To chan-Danguy (1975) and Tochan-Danguy and Hopfinger (1977) were able to match the Froude number of powder ava­ lanches by using the gravity flow of salt solutions. They constructed a three meter long inclined tank to allow a salt brine to flow under a less dense ambient fluid. The qualitative results of this model strong­ ly resemble the motion of the natural event. Quantitative data, however, was limited. As sparse as data is for velocities, it is much better than the data available for snow density within the flowing avalanche. Direct density measurements have been very difficult to make. Eyber-Berard et. al. (1980) attempted to use gamma rays to obtain some average density measurements.in.flowing snow, but they obtained.little.useful data. Indirect evidence based on flow depths, has also been .,used to estimate the average snow density in an avalanche. These estimates range from 60-90 kg/rn^ for dry. snow avalanches to 300-400 kg/m^ for wet snow 8avalanches, Schaerer. (.1975) . Schaerer and Selway (1980) concluded from evidence gathered from their pressure gauges located at Rogers Pass in British Columbia "that well-devloped dry snow avalanches... consist of three stratified components: dense flowing snow at the bottom, light flowing snow, and powder snow" at the top (Fig. I). They go on to say that the dense zone was 0.5 to 1.2 m deep and formed the principle mass of the avalanche. The light flowing snow, with an average density from 10 to 30% of that of the dense zone, consisted of a mixture of powder and lumps of snow of up to 60 mm in diameter. Their speculation is that this I to 5 m layer of light snow is the result of lumps being tossed upwards from the surface of the dense flowing snow by the turbu­ lent motion at that surface. In the runout zone the light flow is often observed to separate from the dense flow either at bends in the avalanche track or during rapid deceleration. The third component of the avalanche reported by Schaerer and Salway is the powder cloud which. could not be detected by their instruments, but was observed visually. Systematic direct measurements of snow density in avalanches will probably have to await further development of remote sensing.techniques. There is, in general, a lack of documentation of avalanche flow. Very little time-dependent data is available and virtually no systematic ■measurements have been made. Also lacking is experimental .data under controlled flow conditions, detailed data.that could.be used.for evalua­ tion of theoretical models. Therefore, it was necessary to initiate 9POWDER SNOW LIGHT FLOW DENSE FLOW SCALE m 0 5 10 15 20 FRONTAL SPEED 27.7 m s Figure I. PROFILE OF AN AVALANCHE DERIVED FROM PRESSURE SENSORS (SCHAERER 1981). 10 a field test program fpp.the express purpose of documenting snow flowing under controlled conditions. This effort is described In Chapter XV. The next chapter gives an historical review of the develop­ ment of snow avalanche theory. CHAPTER III HISTORICAL ANALYSIS OF AVALANCHE MOTION If one we^e to try to characterize snow avalanches, in general terms, a description might read: 'A snow avalanche is the transient, 3-dimensional motion of a variable mass system made up of. a non-rigid, non-rotund,non-uniform ,assemblage of cohesive granular snow fragments flowing down a non-uniform slope of varying surface resistance.1 The driving force in an avalanche is the force of gravity, but the resist­ ive forces are dependent upon the interactions between the individual snow fragments. Even if the form of the interaction forces were known, the large number and highly variable nature of the particles would make analysis unmanageable, except in some statistical sense. Given, then, that the exact mechanical solution to the avalanche motion problem may be impossible, it is also unnecessary. The motion of individual snow fragments is only of limited interest. For most purposes the time and i ■ special transients are small. A possible exception would be the impact of an individual snow clod upon an area of similar or smaller size. In general, it is the overall motion - the average motion, in some sense - of the avalanche that is of interest. Approximations to the mechanics of the fragment interactions must be made so that analysis can be carried out. in this way, some appropriate average velocity may be determined from which, runout distances and impact forces can be calcu­ lated. It is generally true that the more realistically the interaction forces are modeled, the more accurate and general will be the analytical 12 solution. ■ It Ihas-Ilbeeri the, objective, of the avalanche dynamicigt to formulate avalanche flow models .simple enough, to solve, yet complicated enough, to include as much of the mechanical behavior of the .avalanche as possible. Open Channel Flow Among the first to present an analytical formulation.for.avalanche flow was A. Voellmy in his classic monograph presented in.1955. Voellmy based his analysis on the equations developed to describe open- channel fluid flow. He extended the open-channel hydraulics approach, with all its implicit assumptions, to the flow of snow on a mountain slope. Voelimy, and many others since him, have based their theories on the observation that snow avalanches seem to "flow". That is, as the material deforms, it conforms to the slope, and is often seen to flow around obstacles. This motion is much like the movement of water along a streambed - hence the analogy with open channel hydraulics. Implicit in the assumptions of this approach is that flowing snow behaves as a fluid. In most continuum mechanical theories, fluids are considered to be continuous, homogeneous, and isotropic. That is, any granular nature of the material is. ignored, and spatial variations in the fluid are neglected. Furthermore, hydraulics deals with materials, that can be considered incompressible. The degree to which, these considerations influence, snow avalanche movement determines, the extent to which an incompressible fluid model can give useful results. The 13 incompressible fluid model.does,however, provide ..for simple analysis. Using a continuum mechanical theory, the avalanche flow problem is reduced to finding point values for the stresses within the material. This can be done most easily by using a momentum-balance approach. The stresses must then be related to.the deformation by a constitutive equation that takes into account the properties of the fluid. The open channel hydraulics analysis that Voellmy applied, contains, in addition to the fluid assumptions, approximations related to fluid motion in an open channel. Chief among these approximations is that the flow is steady: that is, that there is no time-dependence to the motion. It is assumed that the avalanche quickly accelerates on a constant slope to a terminal velocity and thereafter exhibits steady motion. The momentum-balance equation is then commonly integrated over the slope-normal cross sectional area so that only the boundary stress­ es need be considered. The dependent variable becomes the average velocity of the material in cross-section. In addition, a uniform motion is usually assumed, so that the flow variations along the slope can be neglected. Under these assumptions, the momentum balance equa­ tion becomes the equilibrium condition, T = pgh sin© C3.ll in which T is the resistive shear stress at the fluid.boundary, p is the average fluid.density, g is the acceleration of gravity, h is the 14 flow depth., and Q is the slope angle. In order to derive an expression for the velocity of the flow, an equation; relating, the velocity to the stresses on the boundary must, be introduced. For this constitutive. equation Voellmy used the Prantl mixing length approximation CPao 1961)., which, for turbulent fluid flow, in which the velocity and stress­ es at any point fluctuate randomly, yields the following boundary drag: T = Ic1V2 (.3.2) However, Voellmy also noted that a minimum slope angle was required for an avalanche to start, so he proposed to modify this equation by adding a constant. 2 ' T = k-v + TJghp cos0 . ,(3.3) J- >• The form of the additional term in this equation is the same as that of a conventional coulomb friction force in which y is the coefficient of friction. The addition of the friction term was found to be necessary, because of the cohesion or locking property that snow exhibits; a minimum stress-level must be exceeded before deformation occurs. Using equation 3.3 for the shear stress at the avalanche snow surface inter­ face, the momentum balance equation may be solved for the flow velo­ city yielding V h(sin0 - pcos 9). (3.4). where the drag coefficients and the density have been combined into a single parameter £. The flow velocity calculated by this equation is the average steady state or terminal velocity that an avalanche is 15 supposed to attain undet uniform conditions on a constant slope. (Jn the case of channellized .flow, the height of. the avalanche, h, is often replaced by the hydraulic radius, R, which is defined as- the ratio of the cross-sectional, area to the wetted perimeter.} To calculate runout distances, the avalanche path....mus.t be divided into three parts: the starting zone, where the avalanche, is initiated and accelerates; the avalanche track, with a constant slope which allows steady-state flow equations, to be applied; and the runout.zone, where the slope of the track has decreased to allow, the avalanche, to dec­ elerate to a stop. If it is assumed that the avalanche enters a con­ stant-slope runout zone with the velocity calculated from the steady- state conditions of the track flow, the runout distance S, can be found from a simple energy-balance equation. The equation Voellmy derived is ' o P S = v /[2g(ycos0 -tan0) + v g/^h^] (3.5) in which h is the mean deposition depth. By applying these equations • m to avalanche paths, avalanche runout may be calculated. The usual procedure is to subjectively break the avalanche track into its three sections, the important point being the transition point into the runout zone. The steady velocity is calculated, from the average slope of the track and then the runout distance is computed from the transi­ tion point. The result o£ the analysis is highly dependent upon exactly how the path is divided and where the transition point is chosen. 16 Through experience, individuals have become very proficient at choosing these points, and this method has yielded good results. As the ava­ lanche paths become complicated, the application,of Voellmy1s equations become more difficult. Also, the method is found to.be difficult for an inexperienced person to successfully apply. The material parameters appearing in Voellmy's equations ,are the turbulent friction coefficient, and p, the dry friction coefficient. The.value of the coefficient C is dependent upon.snow type and motion; . generally the heavier the snow and the slower, the motion the smaller the turbulent coefficient. As speed increases and the snow becomes more powdery B, becomes larger.. The coefficient C was said to be in the range 2 400-600 m/s by Voellmy. Later investigators, Mears (1976), Schaerer (1975), Leaf and Martinelli (1977), Buser and Fruitiger (1980) , and 2 Martinelli et al. (1980), have allowed C to become as large as 1600 m/s for powder avalanches, and have recommended average values of about 1400 , 2 .m/s . Bucher and Roch (1946) were among the first to publish data for the dry friction coefficient, p. They found that for hard, wet snow sliding slowly over a hard, wet snow surface, the friction coefficient ranged from .I to .4. Inaho (19411 slid blocks of dry snow over, a granular .. surface and found,}4 to range from .4 to .6. Reimgactner [1977,1 measured values for p from .22 to .39 bn a small test slide at WeissfuhjocIi., and Sommerhalder (1972) computed a maximum friction coefficient of 0.5 and ■17 a,n average value of Q.32 from data Qn snow flowing oyer' gnowsbeds. Using data reduced from 16 mm movie film, Martinelli et. al (19.80). . computed the friction coefficient for a wet slab avalanche in the runout zone to be from .30 to .35. Center of Mess Motion Salm (1966). expanded upon Voellmy1 s. theoretical, treatment' and derived a time-dependent equation for the velocity of the mass-center of an avalanche. His resistive, force was a Taylor series in velocity,- truncated after the third term. R ~ Ro + Ri i i + r 2 f r . ( 3 - 6 ) The constant term (Rq) represents a dry friction force, can be attributed to linear viscosity, and the third term is due to turbulent friction and the ploughing of the avalanche into the stationary snow. Salm neglected the second term with respect to the third for high speed flowing avalanches, assuming that the energy dissipated by internal friction (R^ term) is small compared to the turbulent energy dissipation (the Rg term). The resulting expression for maximum velocity on a constant slope then took the same form as Voellmy1s equation. Salm, also used energy considerations to derive an equation for runout dis­ tance. It is a more comprehensive treatment than Voellmy1s and included the possibility of artificial obstacles. The equation is.considerably more complicated but it does not contain the troublesome singularity that appears in Voellmy1s equation. 18 Mellor in 1968 presented a,n analysis similar to that of Salm (19661 except that he considered the possibility of having a variable mass system. The mass was allowed to vary linearly with the avalanche velocity. m = m (I + kv) (3.7)" Q . Here m is the initial avalanche mass and k is an entrainment co- o efficient. This is an empirical relation based upon a correlation between avalanche mass and velocity. The dynamic equations that result from this analysis are cumbersome and "afford little insight into the effects of entrainment as they stand." (Mellor, 1968). Perla et al. (.1980) also developed a center-of-mass system of equations similar to these models. They included a mass-entrainment term proportional to the square of the speed in the dissipative forces. 2 This term combined with other v terms into a single term with the coefficient D, resulted in an equation of motion that was solved successively over small intervals along the slope. In each interval the slope and friction coefficients were considered to be constant. The expression for the maximum velocity over any interval has the same form as Voellmy1s equation for terminal velocity, V = /-^ (sinO - J-lcos 0) (3.8).max V D in which m/P replaces Yoellmy'.5 £h. In addition, at the juncture of the two segments, a correction for momentum loss due to the slope transition is included. ValIues for the parameters ]_t and m/D for this model are found in the same way that y and £ are found for Voellmy1s model. The paper5 by Perla et al. (.1980). and. BakkeWi et al. (1980). discuss and give examples for these parameter's. In another work, Perla (.1980) analyzed the effects of snow en­ trainment that occur at a constant rate per unit length, of avalanche travel. "Initially entrainment and drag cooperate to oppose accelera­ tion, but if entrainment continues the mass increase will eventually work against the drag and velocity will increase. However, ....the velocity boost is only significant given continual entrainment on very Ipng paths." Linear Viscous Model . A major drawback to the mass-center approach to avalanche dynamics is the lumping of the avalanche mass at a single point. An avalanche is an extended mass of variable densty and changing velocity. The use of the center-of-mass equations does not provide any insight into the spatial properties of the flow. The flow depth, for example, may appear as a parameter, but it is never calculated. Tq account for spatial characteristics, a general stress-deformation relation must be applied over the entire flowing mass. One of the simplest relations is the continuum mechanical representation of a linear, viscous fluid, where the flow in two dimensions, is governed by the equation; 20 Here, u and. v are the slope parallel and slope normal velocity compo­ nents , and x and y are the slope parallel and slope normal coordinates. s The constant of proportionality y, between the shear stress, T, and the velocity gradient term is called the dynamic coefficient of viscosity. In addition the kinematic coefficient of viscosity is defined as the dynamic viscosity divided by the mass density, p, of the fluid. When this constitutive relation is substituted into the equation of motion, which is derived from momentum balance considerations, a non­ linear partial differential equation results - the Navier-Stokes Equa­ tion. For the flow geometries encountered in avalanche motion only computer solutions to this equation can be found, and then only in two dimensions and for incompressible flow. The method of solution and the use of this model is detailed in Chapter 5. The principle parameter for linear viscous fluid flow, the Newtonian viscosity, has never been measured directly for moving snow. Estimates -5 of this parameter range from the kinematic viscosity of air, 1.5 x. 10 2 , (Perla, 1980), to the viscosity associated with creep motion 63.4 2 m /S, (Shen and Roper, 1970). Bingham Model . The linear viscous flow model provides a method for calculating the spatial characteristics of flowing snow. However, it is questionable whether flowing snow obeys a linear constitutive equation. The exact dependence of the internal resistive forces on velocity gradient is not 21 known. Many have speculated that at low speeds, the forces may be linearly viscous and that at higher speeds the forces, take on a velocity squared dependence due to turbulent motion (VoelImy 1955, Perla, 1980). Thus at least for flows at low velocity a linear viscous model would seem to be appropriate. However, snow exhibits the granular property of requiring a minimum stress to initiate deformation. Regions of material below that minimum stress level do not deform, they translate as rigid bodies. This is seen clearly in the initiation phases of an avalanche where large chunks of material move as blocks until the internal stress­ es become large enough to break them apart. The linear viscous equations cannot model this behavior since deformation occurs until all internal stresses are relieved. A block of linear viscous material initially at rest on a horizontal surface will deform until it has completely col­ lapsed to zero height and infinite surface extent. The simplest mathematical model that exhibits the locking property observed in snow and other granular materials is the Bingham model. In the Bingham model no deformation occurs until the shear stress in the material exceeds a specified value AB ), then deformation proceeds as a linear viscous material at a ■rate proportional, to the amount that the stress exceeds T , The constitutive equation for a Bingham fluid then consists of two parts. For T < T — + — " Q (.3.10) o 3y $x 22 and For T > Tq,W T - Tg C3.1U Substitution of this constitutive relation into the equations of motion again yields a partial differential equation. The details of the com­ puter method of solution of this equation is also included in Chapter 5. Direct measurements of the parameters in the Bingham model, the shear stress (Tq) and viscosity coefficient, (y), have not been carried out. However, Bucher and Roch (1946), in measuring the frictional resistance of hard wet snow for speeds between 0.2 and 2.4 m/s, found that the linear fit to their data yielded a constant of proportionality of 475 N- ^ ec . If it is assumed that there was a 2 mm layer of granu- 3 lated snow of density 300 kg/m between the sliding surfaces and that the velocity gradient was linear in this region, then the viscosity in 2 this layer would be about .003 m /s. Similar tests by Lang et al. ■ (1981), using hard sintered snow over the velocity range 0.1 to 0.25 m/s 2 yielded a viscosity coefficient of .004 m /sec and a EU value of 540 2 N/m (.Figure 2) . These values are for a very narrow range of slow speeds and will most likely differ at higher speeds,, but.do serve as possible ojrder-of-magnitude estimates. The viscosity oj; fluidized snow has been measured in Japan by Maeno. and Nishimura (1979), and Maeno et al. (1980). By passing air upward through a snow column they were able to suspend the snow and create a S H E A R S T R E S S ■ 800 23 I 700 600 500 400 SNOW DENSITY: 300 kgm"3 TEMPERATURE: -18°C NOMINAL SNOW GRAIN SIZE: 0.4 mm RELATIVE HUMIDITY: 88% OVER-BURDEN PRESSURE OF PLATE: 1170 Nm'2 300 200 - 100 - 0 0 .05 10 __ L .15 20 25 Fiyure 2. SHEAR STRESS OF GRANULAR SNOW ON AN OVER-RIDING SURFACE-ROUGHENED STEEL PLATE VERSUS SPEED OF THE PLATE. (LANG ET. AL. 1981). ■24 fluidized bed. "The general behavior of the fluidized enow resembled that of a liquid." ■(Maeno et, al. 1980). The dynamic viscosity was measured by standard techniques for a fully fluidized bed and was found -3 2 to be about .3 to 10 N-sec/m , varying with the particle . size and the fluidization velocity. When this number is divided by the density of the snow, values of the kinematic viscosity are found to be in the range 5 x 10 ^ m^/s to .001 m^/s. As the fluidization air velocity is de­ creased below the velocity necessary for complete fluidization, the measured viscosity coefficients rise sharply. At the minimum possible air flow velocity, below which the incomplete fluidized bed becomes unstable and it becomes impossible to measure the viscosity, the mea- 2 2 surements yield results in the range .0001 m /sec to .001 m /sec. Granular Theories Since flowing snow is made up of granular materials, some of the work in grain flow may eventually be used in snow avalanche models. Mears (.1980) has. tried to apply some of the results of Bagnold1 s (1954, 1956) work on cohesionless grain flow to snow. Additional references in the granular material area are articles by Goodman and Cowin (1971, 1972) entitled "A Continuum Theory For Granular Materials" and "Two problems in. the Gravity Flow of Granular Materials, a monograph by Savage C1979)., "The Gravity Flow of Cohesionless Granular Materials in Chutes- and Channels", an article by Nedderman and Laohakul (1980), "The 25 Thickness of the Shear Zone of Flowing Granular Materials" an<3 a paper entitled "A Micropolar Continuum Theory for the Flow of Granular Mat- . erials" by Kanatani (1979). Constitutive equations for granular mat­ erials are formulated in these papers and the rheology of. grain flow is probed by experimental as well as theoretical techniques. Only a par­ tial list of the work in grain flow has been provided, other investi­ gations either completed or in progress may also have relevance to snow flow. It remains for further research to make use of these theories and formulate better mathematical models for flowing snow. Before much progress can be made in this area however, better data sets must be gathered to check the models against. The next chapter gives the re­ sults of one such data collection effort. CHAPTER JV FIELD TEST PROGRAM The Test Facility A 30P slope near the Bridger Bpwl Ski Area in Southwestern Montana was. chosen as the location of an experimental test facility. This facility was used to run a number of controlled flow tests of avalan- . ching snow. In these experiments, volumes of prepared snow up to 2.7 m^ were released. This snow descended a 30° slope in a polyethelene lined, one meter radius, semi-circular cross-section track. The poly­ ethelene was utilized to minimize friction to achieve maximum snow velocities. The snow mass accelerated down the length of the track and exited onto a flat-level runnout area of packed snow. Using track lengths of 60 m, generated flow speeds of up to 20 m/s at the entrance to the runnout area. Upon entering the level runout zone the snow mass decelerated and came to rest. The runout was marked off in one meter increments and a 16 mm movie camera was used to record this phase of the motion. It was the deceleration phase of the snow motion that constituted the snow flow test. A schematic and photograph of the test facility is shown in Figures 3 and 4. The snow used in the tests ranged in density from 250 to 300 kg/m3. The ambient temperatures on days of testing were always below O0C, with the skies overcast. The snow was released from a polyethelene lined plywood box by unlocking a hinged door. The mass of show then 27 DUMP BOX HINGED DOOR POLYETHELENE LINED SEMI-CIRCULAR TRACK POWDER CLOUD ABATEMENT RUNOUT ZONE PACKED SNOW CARPET SURFACE Figure 3: SCHEMATIC OF FIELD TEST FACILITY. 28 Figure 4: PHOTOGRAPH OF FIELD TEST FACILITY. accelerated as a semi-rigid body. At the entrance to the runout zone a piece of carpeting was installed just above the track, in an attempt to reduce the dust cloud that always accompanied the flowing snow. Depth gauges at the entrance to the runout zone and along the track were used to -monitor the depth of flow, fn addition, the final debris pile was cross-sectioned and the depth profile was measured. An attempt also was made to measure the velocity profile as a function of depth at one.point in the decelerating flow. This was accomplished by placing a glass window in the Runout zone and filming the flow as it passed in front of the window. Marker particles and dye were injected into the snow just prior to reaching the window to en­ hance visiblity. A photograph of this.facility is shown in Figure 5. Observations and Measurements A number of tests were run on this experimental setup during the winter of 1979-80. For each of these tests the dump box was filled with sifted snow, snow chunks, or combinations of both. In addition, powdered dye was added to the snow at specific locations in the release dump. One such release setup is shown in Figure 6. Upon release, the snow mass accelerated, down the polyethelene track and was observed to spread slowly longitudinally. This spreading appeared as crevices formed in the snow. Deformation was also observed, caused.mostly by irregularities in the track. As the mass spread, its depth diminished correspondingly. 29 Figure 5: GLASS WINDOW USED TO MEASURE VELOCITY PROFILE. 31 2.5 1.15.ANGULAR SNOW BLOCK (10 cm max dimension) MIXED GRANULARIZED SNOW AND SNOW BLOCK DIRECTION OF RELEASE Figure 6: DUMP GEOMETRY, TEST NO. 3-2-23-80. 32 The flow eventually peached speeds of 15-20 m/s and then exited onto the runout area where it decelerated. . The variations, in maximum speed were due mainly to the varying condition of the polyethylene lined slope. The sun shining on the plastic caused.the.snow beneath it to melt, and this material then refroze at night. Over a period of several weeks this process increased local irregularities in the track which then retarded the flow. As the speed of the flow reached about 10 m/s, a ,powder cloud was observed to form. This made observation of the main flowing mass difficult. Because of this cloud, the depth of flow could not be measured, but was estimated to be about 30 cm at the entrance to the runout area. Deceleration of the snow continued in the . runout area until the flow finally stopped. Total runnout dis­ tances ranged from 14-21 m depending upon the inflow velocity. This motion was recorded on 16 mm motion picture film which was then analyzed Again, the powder cloud obscured the motion of all but the very front of the flowing snow (Figure 7). Analysis of the film provided only position data on the leading edge as a function of time. Upon coming to rest, the debris depth distribution was measured as a function of position. The maximum debris depth was.generally 20 - 25 cm and the length of the debris was usually about 10 m. Beyond 10 m there was a minimal amount of snow of less than 5 cm depth deposited all the way back to the runout entrance. The dump dimensions, runout position 33 Figure 7: POWPER CLOUD. 34 versus time, end debris distribution for a typical test ate. shown in Figures 8,9 and 10. This particular test will be used as the prototype for calibration of the flow models discussed in Chapter 5. The debris distribution for.each test was cross-sectioned to determine the locations of the dyes that were inserted in the snow prior to.release. Jt was generally found that the dye maintained the same relative position in the debris that existed in the . released snow. For example, the debris dye locations for the release shown in Figure 7 are plotted in Figure 11, In this figure the distribution of the different colored dyes in the debris pile is shown where concentrations were greatest. The red dye which started at the back and bottom of the dump ended up with some dye on the surface of the debris near the 15.0m position and then penetrated into the snow until maximum concentration was observed at station 14.0 near the bottom of the debris. The blue dye, which started at the bottom front of the released snow, ended up as a ribbon of snow dyed blue I to 5 cm thick at the bottom of the debris.extending from the front to the rear of the entire debris pro­ file. There was also a concentration of blue dye near the bottom front of the debris, near Station 18. The spreading of dye, initially at the front edge of the flow, to the entire length of the runnout zone was observed in other tests as.well. In addition, a thin layer about 1.0 cm in thickness of snow dyed blue was left trailing the debris all the way back to the entrance of the runout zone. 35 0.75 GRANULARIZED AND SMALL CHUNK SNOW SIFTED SNOW DIRECTION OF RELEASE Figure 8: DUMP GEOMETRY, TEST NO, 2-2-23-80, PO SI T I O N (m ) 36 18.0 16.0 14.0 12.0 10.0 TIME (S) Figure 9. LEADING EDGE POSITION VERSUS TIME FOR TEST NUMBER 2-2-23-80 I 37 TOP VIEW CROSS SECTION „ * * . X M * ’ X I X * f X X I 3 K X Xr V : X 6 3 10 I2 14 16 IE3 * 1 20 POSITION Cm) Figure 10: DEBRIS DISTRIBUTION TEST NO. 2-2-23-80 38 TOP VIEW Individual Blocks. 3,01 BROWN DYE BLUE DYE GREEN DYE RED DYE (Surface to YELLOW DYE Bottom) CROSS SECTION I I * X X [ : I X X < X X X : K < X * X 6 8 10 12 14 l b Iij 2D POSITION Cm) Figure 11: DEBRIS DISTRIBUTION TEST NO. 3-2-23-80. 39 For the majority of tests, including those discussed above, the release mass was made up of snow sifted through a 6 mm wire mesh to at least a depth of 30 cm. The upper half of the release.mass was then made.up of granularized.snow with some interspersed blocks, of up to 5 cm maximum dimension. In another series of tests, considerable effort was taken to fill the release box with .individual, show blocks with a minimum of interstitual granular snow. Large blocks of snow were taken from the mid-winter sintered snowpack and cut by a shovel into smaller blocks with maximum dimensions between 10 and 20 cm.j These I blocks were individually placed in the dump box, with the estimated granular snow residue being less than 5% of the total snow content. Overall dimensions of the release dump and location of local concentra­ tions of dye for a typical chunk test are shown in Figure 12. Outflow of the main body of snow for this test extended to 14.0 m, with individ­ ual snow blocks projecting out to 18.0 m (Figure 13). The movement of dye was similar to that observed in other tests. In this case, the brown and blue dyes remained localized on the surface, the green dye stayed in the trailing section and the red dye extended as a :ribbon over the length of the debris at the bottom. The yellow block placed on the surface remained on the surface, of the debris .after the test. Ip examining the. content of the debris pile, it was found that no solid blocks extended down to the runnout surface, and the snow.at the bottom of the debris was all granular. Blocks.that were deposited on the 40 GREEN DYE ON BACK FACE ANGULAR SNOW BLOCK (10 to 20 cm max dimension) DIRECTION OF RELEASE Figure 12: DUMP GEOMETRY, TEST NO. 1-2-25-80 (Block Test). 41 TOP VIEW INDIVIDUAL BLOCKS xxx DEPTH INCLUDING BLOCKS- O O O GRANULAR SNOW DEPTHCROSS SECTION X X X I * I x X : X : * " e ' » . * * O ----- '---- I----- 1----- POSITION Cm) Figure 13: DEBRIS DISTRIBUTION TEST NO. 1-2-25-80 42 surface retained their angular shapes, while blocks.that were, found interspersed with granular snow buried in the. debris were small and rounded showing considerable wear and abrasion. Figures .14 and 15 are photographs.of the debris crossrsection showing this snow distribution. The average depth of the surface of the debris, including blocks, was recorded and then the blocks were removed and the depth of the granularized snow was measured. This information is plotted in Figure. 13. Window Tests In a number of granular snow tests, the motion of the flowing snow was photographed from behind a glass window ihseb.ted in the wall near the entrance of the runout area. This enabled the measurement of the flow velocity as a function of depth in the flow. Figure 16 shows the resulting velocity profile as measured from several tests. These measurements were only taken 0.1 to 0.2 seconds into the flow past the window to minimize the effects of deceleration of the snow. These results show a strong velocity gradient near the lower surface and a nearly uniform velocity for the upper 75% of the flowing snow. The dashed line in the figure is an exponential curve fit to the data. Discussion,of. Test Results In examining the test results, several conclusions can be. made. First, there is little mixing or relative motion in the upper depths of the flow. Dye that started on the surface in tests was found to still 43 Figure 14: CROSS SECTION OF DEBRIS, CHUNK TEST 1-2-25-80.. 44 Figure 15: CROSS SECTION OF DEBRIS, CHUNK TEST 1-2-25-80. DE PT H (m ) 45 0.35 0.30 0.25 0.20 7.0(1 0.15 0.10 VELOCITY (ms 1J Figure 16: DATA POINTS OF PARTICLE SPEED VERSUS DEPTH OF FLOW 46 be on the surface when■the flow came, to rest. There was. Very little penetration of the dye into the body of the flow despite an increase in the surface area. Also, the angular snow chunks that were near the surface of the flow did not significantly abrade. The measured velo­ city profile also shows Very little Velocity gradient in. the upper regions of the flow. This indicates that there were large areas of the flow where little deformation was taking place. A second conclusion has to do with the motion of the flow over the stationary packed snow surface. It is inferred from the dye tests that there is deposition of material from the bottom front of the flow onto the entire runout area. In test No. 2^2-23-80 red dye initially at the bottom front of the flowing snow mass was deposited in a thin layer 1-5 cm in thickness at the bottom of the entire length of the debris pile. The red dye extended all the way back to the entrance to the runnout area. From this evidence/ and from the shape of velocity profile measured near the lower boundary, it is strongly suggested that a no-slip boundary condition exists at the interface between the flowing Show and the runout Surface. The very large Velocity gradient in this lower region then results, as seen in the block tests, in the grShUlari^ zation.of the snow located there. This produces a very uniform material in which the majority of. deformation takes place. It can be speculated that the material parameters describing avalanche flow.and used in avalanche modeling may not Vary greatly between avalanches of similar types. CHAPTER V.' B!VISCOUS BINGHAM MODEL Linear viscous Model The first attempt, to model snow, avalanches, employed a linear viscous model. This formulation was chosen because it is the simplest mathematical description, of a fluid that retains internal dissipation. It was also assumed that the flow was incompressible, which.simplified the continuity equation, and eliminated energy considerations. With these assumptions, only numerical methods could be employed to solve the resulting equations and then only in two dimensions. The governing equations for two-dimensional linear viscous flow begin with the constitutive relationships where T ,0 , and 0. are the shear stress and the nominal stresses, u xy x ' y and v are the velocity components in the x and y directions, p is the velocity gradients. Next, equations 5.1-5.3 are substituted into Cauchy's equations of motion; which, are derived from momentum balance, considerations • (See Malvern. 1969,. Chapter 5, for example!. In two (5.1) (5.2) (5.3) hydrostatic pressure, and p is the dynamic coefficient of viscosity. As can be .seen these equations linearly relate the stresses to the 48 dimensions,. Caucli^Vs ^ eguetions. pf motion:are .r-ti-ffv-. , _xi 3 y 3 y + pb CS. 41. (5.5), where b and b. are the components of the body force per unit mass and x y p is the mass density. Upon substitution the following, equations result: P af"" & - S - P - 32v ^ 3x3y + ^ x (5.6) ^dv 3YM ;A 3p , 9^v , ,,2 3 2u - , (5.7) P .dt = ^ 3y2 ^ »x2 v .. 3x3y y These equations can be simplified somewhat using the incompressibility condition. This condition is derived from conservation of mass con­ siderations and is commonly called the continuity equation. For two- dimensional incompressible flow, the continuity equation takes the form I f -.o (5.8) It is also common to write the final flow equations in terms of the convective time derivative rather than.the Lagrangian time derivative. The relationship between the.two derivatives can be expressed as: 4u dt dv dt ■ 3u ■ 3t ■ W . 9t + U- 3v + v- 3y 9y CS. 91 ■ (5,101 49 where is the Lagrangian derivative and is the.Eulerian deri­ vative . Substituting 5.8 - 5.10 into 5.6 and 5.7 and dividing by p yields: CS.Ill 9v 3t 3v 9y 36 3y CS. 121 where V/ the kinematic viscosity, has been substituted for p/p and (J) = p/p. These equations are called the Navier-Stokes equations, and as can be seen are non-linear partial differential equations for the velocity components u and v. Analytic solutions to these equations exist only for very simple flow geometries. Time dependent, non- uniform, flows with a free surface must be solved numerically. SMAC A computer code to solve the non-steady, two-dimensional Navier- Stokes equations for incompressible flow was developed at Los Alamos Scientific Laboratory by Anthony Amsden and Francis Harlow. (1970)• The code, called SMAC, for "simplified marker, and cell" method, uses finite difference .techniques to solve equations 5.11. and. 5,. 12, The marker and cell method utilizes a set of marker.particles, to delineate the location of the..fluid. These dimensionless and massless particles do not participate in any of the velocity calculations. They are 50 simply advected with the calculated fluid velocities,, and serve to. locate the position of the fluid, and in particular the fluid surface. In' solving the.velocity equations, a two-step process is employed. This is necessary, because,of the non-linear character of the .equations. The first step is a calculation of tentative velocities at the advanced time step. This calculation is carried out by solving equations' 5.11 and 5.12 subject to an arbitrary pressure field, . The correct pressure at the free surface is, however, specified. This imposed boundary condition is the satisfaction of the normal stress require­ ments, which for surfaces of small curvature can be approximated by, Bu - 2V - ^ = ° (5.13) (Hirt and Shannon, 1968)., where n and u . are the surface normal directionn and surface normal velocity. Upon specifying the correct velocity , boundary conditions, it is insured that the tentative advanced time velocity field will contain the correct fluid vorticity, defined by the equation, W = V_ x v ■ (5.14) Here V x is the vector curl operator. Physically the vorticity is the local rate of rotation of. the fluid. At this point in the. calculation the.tentative velocity field.will not satisfy the incompressibility condition. The. second step in' the calculation, is to modify the., tentative, velocity field to.'satisfy the incompressibility condition without I51 disturbing the already correct vorticity- distribution. This is... accomplished through the use of a potential function, f , such..that the spatial gradients of ip become .the velocity changes necessary to insure incompressibility. .Since V. xV}p = 0, where V is the vector, gradient operator, does not change the fluid vorticity already!, calculated. In order to carry•out the first step in this calculation, equations 5.11 and 5.12 are rewritten in a conservative form and'then finite differenced. The differencing is carried out over a rectangular grid superimposed on the flow domain. Figure 17 shows the locations of the various variables and parameters as defined on this superimposed grid. Simple centered differencing is employed .for.most of the special .deriva­ tives and the time derivative is pure forward differencing; Applying the correct velocity boundary conditions involves satisfying the shear stress condition at the free surface , 3u 9u n + ^ = b , (5.15)3m ' 3n where m refers to the free surface tangential direction. Equation 5,15 is also used to satisfy the normal and tangential boundary conditions ' on the rigid walls. The.tangential velocity condition for ordinary viscous fluids is a no-slip condition where.the tangential.velocity must vanish..at the. wall. This condition is applied by specifying the ■ tangential velocity in the finite.difference grid.cell outside the • , \ 52 X B B B B B B B X B E E E E E E E B B h - ''S-., E E E E B B F F F ' W /Z E E E B B F F F F s \ E E B X B B B B B B B X i j+ 1/2 1-1/2, i+1/2,j+1/2 i+1/2,j i+1/2,j-1/2 CELL i,j B = BOUNDARY CELL F = FULL CELL S = SURFACE CELL E = EMPTY CELL Figure 17: SMAC COMPUTING GRID fluid domain adjacent to the wall. This velocity is chosen to be the negative of the velocity in the fluid cell it is next to. Using a linear velocity interpolation between the tvyo cells yields zero velo­ city at the wall (Figure.18}. Provision can also be made for a free- slip boundary condition for cases in which the boundary la'yer at the wall may be too small to resolve, with the computing mesh. This is done by specifying the velocity in the boundary cell to be the same as that in the adjacent fluid cell (Figure 18}. The boundary layer effect may also be applicable to the free surface shear stress condition, in which case equation 5.15 must be modified. The free surface shear stress condition then becomes, 53 S F = 0 and = o (5.17) A friction type boundary condition can also be applied using this pro­ cedure. The boundary cell velocity may be set to produce any viscous shear stress at the wall that may be desired. This shear stress may be determined from any given friction condition. The second phase of the calculation is to modify the tentative velocity field to satisfy the incompressibility condition. By ex­ pressing this modification, as the gradient of a potential function, the vorticity field remains unchanged. The potential function then is found fpom ft solution of Poisson's- equation = D (5.17} 54 FLOW CELL WALL NO SLIP CONDITION VELOCITY GRADIENT Ay FLOW CELL WALL VELOCITY GRADIENT BOUNDARY CELL FREE SLIP CONDITION J - L FLOW CELL wall WALL VELOCITY GRADIENT BOUNDARY C wall FRICTION BOUNDARY CONDITION Figure 18: SMAC METHOD OF SPECIFYING THE TANGENTIAL BOUNDARY CONDITION AT A RIGID WALL. 55 ■ 2 3 '3' where V . is the operator — - + — ; 9 x 9 y that is found in the tentative velocity field 2 # and D is the divergence, V + Iy- If the tentative velo­ city field is designated v, then the final velocities v.. afe given fay V = V + Vlp . CS.181 Equation 5.17 is derived fay considering the incompressibility condition which can be written V • V = 0 (.5,191 When equation 5.18 is substituted into this relation the following equation results: V . V + V 2Ip = 0. (5.20) This equation is rewritten to correspond to Equation 5.17. Equation 5.17 is finite differenced and solved, for ip using an iteration procedure. The final velocities are then calculated according to Equation 5.18. The final feature in the SMAC code is the movement of the marker particles. These particles are moved using.the velocity of the com­ puting cells in which they reside. Simple forward differencing of the position-velocity relationship results in an equation for this move­ ment. The computing cells are then identified if they contain parti­ cles and whether, or.not they are then adjacent to a cell.which, does not contain.particles. Full cells are cells.containing particles,, and . 56 which, have no empty cell neighbors; they are treated, to the full compu­ tational system. Surface cells contain particles, but have at least one empty cell adjacent to them. In these cells the free surface boundary ■ conditions must be applied. Linear Viscous.Mode! Results' In modeling the snow flow tests described in Chapter 4, the linear viscous program SMAC proved to be inadequate as originally written. It was found that the single viscosity parameter in the governing equations was insufficient to simulate the observed snow motion. Using a no-slip boundary condition, a low viscosity was needed to model the correct runout distance, and a high viscosity was necessary to give the correct debris distribution. This contradiction was resolved by assuming the formation of a friction boundary layer at the bottom interface between ■■ ■ the flowing snow and the packed snow runout. This layer was treated as a friction boundary condition as described in the previous section. As a first guess, the friction force exerted on the flowing snow was assumed to have the form of a Coulomb, or dry friction force; F = a N (5.21) where N is the normal force and a is the coefficient of friction. How­ ever, it was quickly found- that for reasonable values of the friction coefficient (see Chapter 2)., this constant force was unable to bring the computed snow flow to a stop soon enough.. It was .decided that rather 57 than go to unreasonable and unphysical values.of cx, a second friction force would be tried. This friction .force was..,assumed, proportional to the slip velocity, Vg, i.e . F = g Bo SAzAYYe where g is the constant of proportionality. Modeling with, this fric­ tion force showed good high speed correlation between calculations and data for the leading edge motion,. but proved to be insufficient at low- speeds . Finally, the preceding two friction forces were combined into a single friction term containing both coefficients a and $. F = CXN + 3vs (5.23) One of these coefficients was eliminated by writing it in terms of the other, utilizing data from a steady-state flow configuration. This again yielded a one parameter friction model. By adjusting this re­ maining parameter, a good fit to the velocity versus time curve for the snow tests was achieved. This was accomplished for a value of 0.5 for. the coefficient a, and a corresponding calculated $ coefficient of .050. One other form of the friction force was investigated. It was made up of a term proportional to the square of the slip velocity in addition to the Coulomb,friction term, i.e. 2 F = OtN + gVg (5.24) Again, g was calculated as .a function of a using equilibrium data. JFor a value of a = 0.6, and g = 0.075, this friction force also adequately modeled the leading edge motion of the snow tests. For the flow speeds ,encountered in the snow tests, it was not possible to.choose which of the two friction forces provided more accurate simulation.. In order for that distinction to be made data for higher speed snow flows must be collected. However, based on historical models and terminal velocity data, arguments can be made that the velocity-squared force is more appropriate. Equilibrium considerations using this type of friction force yield a terminal velocity expression identical to those used by Voellmy (1955),. Salm (1966), Perla (1980), and others as discussed in Chapter 3. Ths'final debris distribution calculated by the linear viscous model proved to be.principally a function of the viscosity coefficient. It was found that for this phase of the modeling, a kinematic viscos- ' 2ity of. 0.04 m /s yielded the best fit to the data. A detailed com­ parison between model and experiment, as well as further discussion,. is^ contained in the article by bent and Lang (1981). Bingham Model . A. major failure, of ..the linear viscous model for flowing snow is .■ its response, to low stresses. In the modeling just described, it was necessary to halt the computations at the point where the leading edge of the flow fell below an arbitrarily small velocity. To continue the 59 calculations beyond that point, would have eventually allowed.the fluid to deform until the.material depth was reduced to zero and the horizon­ tal dimension had become.infinite. This motion exemplifies one obvious difference between flowing snow and a linear viscous fluid. Snow is seen to come to rest with., a finite depth, where the stress-state is non­ zero. This indicates a threshold stress state in snow that must be overcome in order for deformation to take place. This property is due in part to the cohesion of the individual snow particles, but is mainly a result of the granular nature of the material. The simplest continuum mechanical model to exhibit this locking property is called a Bingham fluid. The constitutive equation for a Bingham fluid is made up of two parts. First, if the stress is below a threshold value, Tq , no deformation takes place. Second, if the stress is above this value, deformation takes place, and is proportional to the amount that the stress exceeds T . In equation form the constitutiveo relation is expressed by: 1 if T < T (5.25) oy ox o and ' i£ T - T° t5-261 Well k n o w n S i n g h a m type m a t e r i a l s include, paints, greases, concrete, a n d toothpaste. In addition to modeling the locking property of flowing snow, a Bingham model contains the necessary stress-deformation components to 6Q model the boundary, layer that.was treated as a friction force, in the linear viscous model. In part, this is due to the additional parameter, Tq , involved in .the Bingham, equations. But also the very nature of the Bingham model, being physically more accurate, allows a more realistic representation of the motion, and as will be seen by the results, provides a,very good fit to the data. Biviscous Bingham Model The implementation of the Bingham model proved to be a difficult task. Of primary importance to this model is the location of.the yield surface that separates the two flow regimes. On one side of this sur­ face the material is. locked and behaves as a rigid body. On the other side, the constitutive equation, when substituted into Cauchy1s equa­ tions, provides the Navier-Stokes equations. In both cases, existing methods easily allow the governing equations to be integrated. Un­ fortunately, calculating the location of T = EU is not a simple pro- cedure,- since the constitutive equation does not define the stress when it is less than Tq. The surface for which T - EU cannot be approached from below. The surface can be found in the deforming region; it is the location at which the velocity gradient first goes .to zero- ‘ Using this surface and applying the governing velocity equations over the two defined regions will.result in an updated velocity calculation. How­ ever, the location of the surface T = T will not have changed. Thiso surface can be relocated by considering the equations of motion and 61 calculating the stress- state ,from ■ the acceleration. Howyer/ the two equations of motion in two dimensions contain three,unknown stress components. The constitutive equation in the deforming ,region must be . invoked to provide additional relationships. When the equations.of , • ■ . '• motion and the constitutive equation are solved for the.stresses, a second order partial differential equation results., This equation, along with the two partial differential equations of motion, then comprise the system of equations that must be solved for the ‘ location of the surface where'T = Tq . This procedure represents a formidable problem, far beyond the original task of intergrating the Navier-Stokes equation. For this reason, an alternate formulation to perform essen­ tially the same task was considered. This new approach allowed small deformations to take place accord­ ing to a linear viscous flow law in the locked portion of the flow. The viscosity used in this region is taken high enough so that the resulting deformation can be neglected relative to deformations outside the region. The small deformations and linear visCpus flow.law allow calculation of stress values from the constitutive equation... Location of the yield surface for which T ^ EU is then easy to find. Outside the small-deformation region*, as in the pure Bingham model, . the.flow is still linear viscous, but. with.a different viscosity. This two . ■ . viscosity.system was dubbed the.bivisdous modified. Bingham model. .62. Figure 19 shows a one-dimensional characterization of this flow law. along with the'corresponding pure Bingham relationship. The mathe­ matical representation of the Mviscous Bingham model is: T “ I-v t I i + S ' 1 for T S T0 and T = TQ f | f '+ I i I " ^ ^ f ° r T > To CS,271 CS. 281 where y' and y are the viscosities in the two flow regions. The term T / y' in the second equation is a small.correction to account for the fact that the velocity gradient is not zero at the point where the stress becomes equal to EU (see Figure 19). In order to irorm a general constitutive equation,the yield surface calculation T = EU must not be dependent upon the coordinate system used. Since the shear stress is a function of the orientation of the applied coordinate system, a more general form of the stress must, be found to compare with the constant Tq. The invariants of the stress tensor can be used. The second or third invariant, which happen to be equal in magnitude in two dimensions, are appropriate. In addition, the pressure can be eliminated from the analysis by considering only the deviatoric stresses, defined in two-dimensions as CT = O1x - (7 (5.29) o; - Oy - 9 (5.30) 63 BINGHAM MATERIAL B!VISCOUS MATERIAL Figure 19: ONE DIMENSIONAL FORM OF THE BINGHAM AND B!VISCOUS CONSTITUTIVE LANS. ■64 ' T ' = T' xy xy where o' is the average normal stress I 15,311 2 % + CTy) (.5,32).. The second inyarient of the stress deviator then, takes the form IXT* <7 '<7 ’ + T ,2x y '.xy CS. 331 Substituting into this relationship the definitions, 5,29. - 5,32, yields a rotationally invariant function of the stresses O , a and . x y Txy: mmEm a m nY. E - , t (5.34) This equation is used as the yield condition in the Biviscous Bingham model. ' IIT , = To (5.35) This analysis follows the recommendations of Hohenemser and Prager (1932) as presented in Malvern (1968), and is the equivalent of the I i ■ Treska yield criterion found in solid mechanics. The total two-dimen­ sional system of equations defining the Biviscous Bingham model can be written: when W T, < T0 = TJ' + 15,36) 9X' = P + 2 "E 'CS,37). Cy = P + 2 ^ ^ CS.38I I■ 65 when' IIm' < T T o f «o + 2 « # - ^ > 0 + 2 ^ . CS. 391 CS.40) CS. 411 Where $ and t are the.shear and normal stress values in the original coordinate system, which satisfy equation 5.35. Implementation of these equations is dependent upon finding values of II^11 at each calculational point in the fluid domain. To accomplish this, the stresses at each point are required. Since it is not known a priori which of the two forms of the constitutive equation to apply a condition on the velocity gradients is derived. This is done by noting that, due to the form of the constitutive equation, the com­ ponents of the stress deviator, 0^,0^* , T ^ 1 , and the components of the velocity gradients, 2 g— , 2 (g^ + g^ ) , transform in the same way for coordinate rotations. Therefore, the coordinate system - g ' 3v that maximizes T'_, will also maximize the term SNi • + ^rl . For this particular coordinate system, 0^' = 0, 0 1 = II 1 = T ^ ,so that equation 5.35 becomes T xymax 3v Sb - *' o, and I ^xy-m9lx ^ To The fluid will then be in the large deformation region if (5,421 . In' tfe - &Sy ^ 3x^maxI > To' (5.431 Calculation of the maximum value of is accomplished using the 66 rotational transformation equations, in two-dimensions, a Mohr's circle construction provides a convenient method to carry out these calculations (See Figure 20.).. Thus: ' 3u + 3v 2 ^ y 9 x ^ max. , 9u 9v. 2 5? + at’- + 4 (5.44). and Equation 5.43 can be rewritten T 2 ' I t + B 2 + " 2 > l # ) ^5.451 Equation 5.45 replaces the criteria 5.35 for determining which of the two flow regimes a particular point may be in. Still needed before the consitutive equations 5.39 to 5.40 can be applied are the values of the stresses f and O q . These are the shear and normal stresses necessary in the original coordinate system to just satisfy equation 5.35. They can be found from Tq , by rotating the stress-state T 1 = T , a . = o, and O ' = o back to the original ■ xy o x' y coordinate system. Carrying this out yields (refer to Figure 20) • 9 9x. (I* t P - ) '9y 9xtmax (5.46) and 9u ; 9v T (5.471 max . jl x,y V (3- maxmax R sin© GIVEN R cos© STRESS STATE 9y 9x 9x'max GIVEN \VELOCITY \ STATE 9v . 2 2 % — (J) = TAN Figure 20: MORHR1S CIRCLE FOR THE ROTATIONAL TRANSFORMATION OF THE DEVERATORIC STRESSES AND THE VELOCITY GRADIENTS. 68 To derive the governing.velocity equations requires the substitu­ tion of either equations 5.36 - 5.38, or equations 5.39 - 5.41, de­ pending upon a test of equation 5.45,into Cauchy’s equations of motion. If the first set of equations are used, the Navier-Stokes equations are rederived with the. viscosity coefficient being equal to In the second case a very complex equation results, since the function ( + TT-) may have spatial derivatives. To avoid the cumbersome3y dx max task of finite differing this equation, the stresses themselves are calculated and substituted into the finite difference form of the equations of motion (5.4 and 5.5) BVSMAC This procedure was implemented to produce the fortran computer code BVSMAC. In general, the form and organization of the SMAC program were retained in BVSMAC. However, the portion of SMAC that calculated the tentative velocity field now required rewriting. For this section of the code, equations 5.36 - 5.47 were used to calculate the stresses in each computational cell. These values were substituted into the finite difference form of the equations of motion (5.4 and 5.5) .' For the most part, simple central differencing for spatial, de­ rivatives was employed in the necessary equations... This was facilitated by defining the normal stresses at the cell centers and the sheas stresses at the cell corners. The stress gradients in the equations of motion were thus easily written. For example, in the first equation 69 CS.4)., which is used to compute the u ,velocity, the finite difference approximation for the first stress gradient is. 3y Ay And the second stress gradient term becomes ■Txy ' j+1/2, j-fl/2 • Txy ' i+1/2, j-1/2 ' C-5*48) ^ 0 X gX i+l/j x i, j (5.49)3x Ax Both these terms are centered at i + 1/2, j where the u velocity com­ ponent is defined. A similar arrangement of terms also occurs in the y component equation (5.5). . A slight differencing problem was encountered in calculating the shear stresses at the cell corners. It is seen from equation 5.39, that for this calcuation the term t is needed, and must therefore be found at the cell corners. This requires calculation of thereoy dx max (Equation 5.47). Employing equation 5.44, the term ( ^ + jpr) easily found at the cell corners, but the (^) terjn requires a more elaborate treatment. Since the u velocity components are defined only at the center of the vertical sides of each computational cell, centered dif­ ferencing of yields the results at the cell centers. In order to find the value of .the cell corners, a simple arithmetic average of the result's at the four cell centers ■ surrounding the. corner can be used. This results in the use of a total of six velocity components in 70 calculating , four of which are over one cell dimension away from the location of the calculation. There are inherent numerical stability problems with this type of finite differencing. (See for example Chapter. One of Ames, 19771 • A somewhat less accurate but more stable approximation involves using only the value of at the center of the cell in which the velocity is being calculated. . For example, when calculating * i + 1 / 2 r y and EU j+1/2 the value of used is the value found at the location (i,j). A similar methodology is employed to find values of (^- + ) at the cell centers when calculating the 'U Value required in the normal stress calculation. Programming this methodology into the BVSMAC computer code yielded a program for modeling biviscous material* A copy of this code is reproduced in Appendix A. In order to confirm the operation of the . BVSMAC program,, a simple one-dimensional problem was chosen to simulate. This problem consisted of the gravitational flow of a fluid between parallel infinite plates. A true Bingham material in this configuration will exhibit a central core region moving as a rigid body. Time de­ pendent analytical solutions to even this simple problem are extremely difficult to determine. As a result, only the steady-state configura­ tion will be considered quantitatively, while the time.,dependent motion may be examined qualitatively. In Figure 2 1 , the flow, geometry for this problem is illustrated. The spacing between.the.plates, being in the y direction was designated 71 FLOWING MATERIAL STRESS DISTRIBUTION VELOCITY GRADIENT Figure 21: FLOW BETWEEN PARALLEL PLATES; GEOMETRY, STRESS, AND VELOCITY GRADIENT. 72 d, and the body force per unit mass was assumed to be in the x direction with magnitude g. Simple equilibrium calculations for this plane stress problem permit the steady state shear stress per unit mass density Txy/P' to be foun^ • Measuring the distance in the y. direction from the midpoint between the plates,.this relationship is T (5.50) The value of y where the shear stress increases to T /p is then given by . gp (5.51) Finally, in the region where deformation is occurring, the velocity gradient can be found from the constitutive equation. du dy T --T xy o (5.52) Figure 21 contains plots of the steady-state shear stress and velocity gradient as functions of position. ■ This problem was set up and run using the BVSMAC computer code. Values of the problem parameters were chosen for simplicity, 'and to give reasonable results without excessive calculations. A U quantities are given in dimensionless units. The spacing between plates was chosen as 1.0 unit and was divided computationally into square, cells 0.10 units on a side. The body force g was given a value of 2.0 s and the principle kinematic viscosity, U/P* was initially set to 73 0.2,0. Using these values the steady-state shear stress from equation 5.50 is seen to vary from 0.0 to 1.0,from the midpoint to .the surface of the plates. A TQ/p equal to 0.30 then gives a value of 0.15 for yQ so that the central plug of material should have a. width of 0.30, or 30% of the total width. The steady-state velocity gradient in the . deforming region can be written IOy - 1.5 for y > .15 . (5.53) The sequence of plots shown, in Figure 22 is the output of a plotting program that shows the position of the marker particles calculated by BVSMAC. In these plots the x direction is horizontal and y direction is vertical so that g is pointing to the right. In order to reduce the calculations, the amount of material in the x direction was only given a length of 0.60 . The dense particles at the leading edge of the flow are special marker particles used to give the free surface better definition. This allows the free surface boundary conditions to be applied more accurately (Nichols and Hirt, 1971). The material was allowed to deform in model time for 1.0 second. During this time, essentially steady state flow conditions were achieved. The region of no deformation, calculated earlier to be 0.30 units in width is indicated by the two parallel lines drawn on the t =1.0 s plot in Figure 22. Figure 23 then compares the analyti­ cally calculated velocity profile with the computer generated velocity 74 Q_o 0 . 4 0 0 . 8 0 1 . 20 1 . 6 0 P O S I T I O N (M) T = 0 . 2 0 0 Q_o 0 . 4 0 0 . 8 0 1 . 2 0 P O S I T I O N (M) T= 0.000 Figure 22: BVSMAC PARTICLE PLOT SEQUENCE FOR THE SOLUTION TO THE FLOW OF A BINGHAM MATERIAL BETWEEN PARALLEL PLATES: T . = 0.30, V = 0.20, V 1 = 2.0. o/p 75 CL O 0 . 8 0 I . 20 P O S I T I O N (M) 2.00 0 . 6 0 1 Q-o 0 . 4 0 0 . 8 0 I . 20 P O S I T I O N (M) 2 . 00 0 . 4 0 1 Figure 22: (CONTINUED) 76 Q_o 0 . 4 0 0 . 8 0 1 . 2 0 P O S I T I O N (M) 2.00 I . OOl Q-O 0 . 4 0 0 . 8 0 1 . 20 P O S I T I O N (M) 1 . 6 0 2.00 T = 0 . 8 0 1 Figure 22: (CONTINUED) 77 profile and .copredpondence;is,seen to be good. For this particular.computer run the value of the secondary kinematic viscosity V', was chosen to be 2.0. Xt is this parameter that determines the deformation that takes place in the central'plug regioni The slight curvature of the . leading edge in this region is due to this parameter. A second test was run in which. V ? was doubled to 4.0. This test result differed from the first only in reduced de­ formation in the locked-up region. The particle plot at t = 1.0 s for this test is shown in Figure 24 and it can be seen that the curva­ ture of the leading edge of the flow in the plug region has reduced. Additional tests were run using different T and V parameters.O When EU was doubled to 0.60 the size of the central plug region also doubled. The marker particle plot for this run is shown at t = 1.0 s in Figure 25. It should be noted.that in addition to increasing the size of the locked-up region in the flow, the overall deformation is also reduced. This is due to greater stresses and smaller velocity gradients in the deforming region (see Equation 5.52). Another way of visualizing this result is to plot the shear stress versus velocity gradient curve (the constitutive relation). This is done in Figure 26, and it can be seen that by increasing T , the V viscosity curve is Q displaced upward by the amount T - is increased.. This produces larger stresses for given rates of deformation, and hence more retardation of the flow. One more computer test was run in which the parameter 78 . . .2 . . .6 ------ ANALYTICALLY CALCULATED VELOCITY PROFILE X X X BVSMAC CALCULATED VELOCITY POINTS Figure 23: STEADY STATE VELOCITY PROFILE, COMPARISON BETWEEN ANALYTICAL AND NUMERICAL SOLUTIONS. DE PT H (M ) DE PT H (M ) ,•0 0 O- UO 0. 80 JP .O O 0. 40 0. 80 79 2.000 . 8 0 1 . 2 0 P O S I T I O N (M) 1 . 6 00 . 4 0 T = 1 . 0 0 0 0.00 1 . 6 0 2 . 00 Figure 24: aySMAC PARTICLE PLOT, FLOW BETWEEN PARALLEL PLATES AT t = 1.0 S :TQ/p = Q,3Q, v = 0.20, V ' = 4.0. 80 0 . 4 0 T= 1 . 0 0 1 0 . 8 0 1 . 2 0 P O S I T I O N (M) 1 . 6 0 2.00 CL O 0 . 8 00 . 4 0 1. 20 P O S I T I O N (M) T = 0 . 6 0 1 Figure 25: BVSMAC PARTICLE PLOTS, t = 1.0 S : To/p = 0.60, FLOW BETWEEN PARALLEL PLATES AT V = 0.20, V' = 2.0. 8 1 (Ol) P M U P —> C> Figure 26: CONSTITUTIVE RELATIONS FOR THE BIVISCOUS MODEL SIMULATIONS OF THE FLOW BETWEEN PARALLEL PLATES. 82 values of the previous test, were retained except V was reduced from 0.20 to 0.10. As seen by Figure-26 this should.result in increased overall deformation over the previous test. Also it should be noted that the steady-state velocity should be slightly larger than the first test. This was confirmed by the computer run as can be seen from the particle plot at t =s 1.0 s which is shown in Figure 27. For this one-dimensional example, quantitatively and qualitatively, the biviscous model gives good results. The two-dimensional aspects of the BVSMAC model were only checked qualitatively using symmetry tests, and then indirectly in the modeling of the snow flow tests. Simple two-dimensional analytic solutions for Bingham materials do not exist. Stability and convergence of the BVSMAC model were checked by changing the grid size and varying the integration time step. It was found that doubling the dimensions of the computing grid decreased the resolution of the velocity gradients but did not seriously degrade the shape or magnitude of the velocity profile. Similarly, halving the grid spacing also had little effect on the results. For this problem, varying the time step also produced little variation in the results. However, the.maximum time step that could be used:was .limited to .0015 second .by ■ the viscous diffusion stability condition . (Hirt, 1968)., <$'t. < 4V (5.54). For the velocities in the test problem, this time step resulted in a 83 c L o o 0 . 4 0 T = 1 . 0 0 1 0 . 8 0 1 . 2 0 P O S I T I O N (M) 2.00 CL o 0 . 4 0 0 . 8 0 1 . 2 C P O S I T I O N (M) 0 . 6 0 1 Figure 27: BVSMAC PARTICLE PLOTS, FLOW BETWEEN PARALLEL PLATES AT t = 1.0 S : TQ/p= 0.60, V ’ = 0.10, V = 2.O. 84 maximum, movement of material of less than. 0.01 of a cell, dimension. From past experience with the linear viscous $MA,C computer code, convergence of the solution was adequate as long as the maximum distance moved by any material was less than one-fourth of a cell dimension. Biyiscous Modeling Results The biviscous model was now used to analyze the.snow flow tests described in Chapter Four. The particular test chosen to be simulated was the test designated 2-2-23-80, illustrated in Figures 8, 9 and 10. Before the modeling could commence however, initial inflow conditions had to be determined. Since the flow entered the runout area from the essentially friction-free polyethylene surface, it was allowed that the initial configuration would be a chunk of material moving at constant speed on a horizontal friction free surface. The initial velocity of this material was taken as 17.0 m/s as derived from the initial slope of the position versus time curve shown in Figure 9. The spatial dimensions of this material was difficult to determine from film footage taken of the test. From the dump dimensions given 3 in Figure 8 a total snow mass of 2.2 m was involved in the test. The width of the flowing snow., at the entrance to, the runout zone was about 1.5.m. Assuming.that.the snow was flowing approximately as a rectangular solid, the.longitudinal cross sectional area of..the flow 2 then must have been about 1.5 m . This number agrees with, the longitu­ dinal cross sectional area of the debris that was measured and plotted in Figure 10. The two-dimensional area of the incoming - flow, mass for • 2the computer simulations was. therefore, chosen■to Be 1.5 m . From the 16 mm motion pictures of the snow test the depth, of the incoming snow was estimated to.be about-30 cm. For.the modeling, this was reduced to 25 .cm, for two reasons: first, in,the numerical modeling it was noticed that when the snow initially flowed onto the friction surface, representing the transition from the polyethylene to the show surface, the depth of the flow increased, an initial depth of 25 cm resulted in a flow depth of 30 cm at this point, corresponding to the depth observed on the motion picture film. Second, as the snow accelerated down the polyethylene, it was observed to spread longitudinally to an estimated length of 7-10 m. In order for the initial snow material to obtain 2 this dimension, and still retain a cross sectional area of 1.5 m , the average depth must be less than 30 cm. In fact, the average depth must even be.less than 25 cm. To stay within these constraints it was reasonable to assume that the flow depth tapered from the front to the rear. The profile chosen started with a depth of 25 cm and was reduced to 20 cm after a total length of 8.0 m. This initial, show profile is shown in Figure 28. The numerical■modeling then commenced with,the flow of this material off the frictionless surface onto a surface employing, a no- . slip boundary condition. The computational grid consisted of ah area '28.0 m long and 50. cm high.. The horizontal dimension was divided into 85 86 Z= Q_ 0 . 0 0 2 . 0 0 4 . 0 0 6 . 0 0 UJ P O S I T I O N (M) Q T = 0 . 0 0 0 8.00 O in 2.00 6.00 8.00 ' P O S I T I O N (M) T = 0 . 0 0 0 Figure 28: INPUT FLOW CONFIGURATION FOR BVSMAC SIMULATION OF SNOW TEST NO. 2-2-23-80 (TOP FIGURE IS PLOTTED TO SCALE). 87 140 cells of .0.20 m length each, and the vertical dimension divided into 10 cells each of 0.05 m height. This proved to be about the minimum grid size that was economically feasible.to employ. A smaller cell size was tried.for an abbreviated run and the results showed little overall variance from results.of a similar test on the .0.20 m x 0.05m grid. A larger celled grid was, however, deemed.inappropriate ■ since the boundary layer at the bottom of the flow was on the order of 5 cm. Cells with vertical dimensions larger than 5 cm would be unable to resolve this layer. The horizontal dimension was then chosen to provide reasonable resolution in that direction, and to maintain an aspect ratio between the cell dimensions of no more than five to.one. The three program modeling parameters, T , V, and V , were then adjusted so that the computed flow conformed to the observed motion of test No. 2-2-23-80. It was found that the parameters EU and V were principally responsible for the leading edge motion and total runnout distance. This is in accord with the results and analysis of the parallel plate flow problem. However, many different combinations of EU and V produced the same runout. Flow velocities were not large enough to provide definite distinctions between different combinations of these parameters. Figure 29 shows several one-dimensional equivalent constitutive relations involving combinations, of EU and V that gave good results’ for leading.edge motion. The velocity profile measured in the window test provided another 88 300 du ,dv dy dx Figure 29: CONSTITUTIVE RELATIONS THAT ADEQUATELY MODELED THE LEADING EDGE MOTION OF SNOW TEST 2-2-23-80. 89. criterion for the numerical simulation to satisfy. Jt was found that the computer generated velocity profile was also principally a function of the two parameters EU and V . As' was increased.and V decreased to maintain the same leading edge characteristics, the velocity pro­ file became sharper, with larger gradients, near the surface and ■ smaller gradients above. Conversely combinations of.small:.Tq. and large V produced gradients more closely resembling the parabolic shape expected for pure viscous fluids. Matching the shape of the velocity gradient shown in Figure 16, provided the necessary information to uniquely define the two parameters, T and V. These two values wereO found to be: for T , expressed in units of stress per unit mass, 2.20 m / s ; and for V, the kinematic viscosity, 0.002 m /s. It was also noted that these values provided the best leading-edge-versus-time comparisons with the experimental snow test. This comparison is shown in Figure 30. Figure 31 then shows examples of the velocity gradient calculated by BVSMAC at the point in the flow corresponding to the location of the window in the snow tests. Also shown is the profile found from the window tests. This profile is plotted on a velocity scale twice that of the other plots. This was done because the measured maximum, velocities at the window were about-, half ..those mea­ sured. and calculated for the center of the flow. Jt is believed that this is due primarily to. the boundary drag exerted on the edge of the flowing snow. The window-measured velocities in the upper reaches of Z 9 0 18.0 14.0 10.0. SNOW TEST NO. 2-2-23-80 X X X X BVSMAC I 0.10 TIME (s) Figure 30: LEADING EDGE POSITION VERSUS TIME; COMPARISON BETWEEN EXPERIMENT AND BVSMAC. 91 DEPTH v(m/s) PROTOTYPE (WINDOW TEST) 10.0 v(m/s) .002 m /s V'1 = 0.10 m2/S v(m/s) T O P 3.0 m2/s2, V = .001 m2/s, V1 = .1 m2/s TO P 2 2 2 0.5 m /s , v = .017 m /s, V 1 = 0.10 m2/s Figure 31: VELOCITY PROFILE COMPARISON BETWEEN SNOW TEST AND BVSMAC CALCULATIONS FOR VARIOUS COMBINATIONS OF BVSMAC PARAMETERS. 92 the flow were about 7.0 m/s, even for the. leading edge. Meanwhile at the center of the flow the leading edge was found to be moving near 16 m/s. The magnitude of the third BVSMAC. parameter,.V1, was found:to have very little effect on the leading edge motion. The velocity profile however, was.affected by,/this parameter, though small.adjustments of EU and V could be made to compensate. . It was also found that V 1 had a'.pronounced effect, along with T and V, ,on the final debris distribution. As V 1 increased, deformation in the upper regions of the flow decreased. This resulted in less total overall deformation of the 2initial flow configuration. A value of V 1 equal to 0.10 m /s combined with the previously specified values of andw, provided the best comparison of final debris depth profiles. This result is plotted in 2 Figure 32, along with the results of a simulation with V 1 = 0.20 m /s. Appendix B shows the full time series particle plots for this simulation. In these plots, the friction-free surface extends from the left boundary to the 8.00 m mark. From there onward the surface is no-slip. The vertical dimension (labeled depth) is plotted on a scale exaggerated by a factor of 4 over the horizontal, scale. Also contained in ajppendix I B, to be used for comparison, are additional.time sequenced particle plots for two other sets of computational flow.parameters. The 'vertical scale exaggeration for these plots is a factor of eight. I i As can be seen from examining Figures 30, .31 and 32, and Appendix i m93 CD in = 2.2 m2/s2 U C KQaOO n n } ri^m'irr' 8.004.00 T = 2 . 0 0 4 .002 m /s V Zr1I1 I ■'■ 11* rrV( 1' = 0.10 m /s 12.00 16.00 20.00 P O S I T I O N (M) 24.00 28.00 Oin Q C b. oo 4.00 I* 0.000 8.00 12.00 16.00 P O S I T I O N (M) 20.00 24.00 28.00 Figure 32. FINAL DEBRIS DEPTH PROFILE; COMPARISON BETWEEN EXPERIMENT AND BVSMAC MODEL. 94 B the BVSMAC. results with T./p = 2.2 m ' / s 2 , V - .002 m 2/s, 9,n<3 V = 2 0.10'm /3 closely model the snow experiment results-. Moreover these parameters form a unique set in which variation in one parameter will degrade the modeling results, whatever adjustments may be made in the other two coefficients. Additional validity to the values of these parameters is obtained from.Figure 2, the measure of friction force as a function of speed (Lang et. al. 1981), and the work of the Japanese on fluidized snow. Recall that these experiments produced measurements 2 of the kinematic flow viscosity on the order of .001 m /s. Also, in 2 the first case, a T stress was calculated tb be 540 N/m for snow of 2 2 nominal 40 cm depth. The T /p value of 2.2 m /s corresponds to an 2 3 initiation stress, Tq , of 640 N/m for snow of density 300 KG/m , well within reason considering the difference in speeds involved. The values of the viscosity coefficients used in this modeling. 2 ranged from .001 to 0.20 m /s. At the low end of this scale the viscosities are about that of heavy lubricating oil (SAB .70) at room temperature. The linear viscous motion of material in this viscosity range is demonstrated in Appendix C, where a rectangular blocjc of material 1.8 m by 0.30 m is released from rest and allowed to deform under the influence of earth's gravity. At the upper end of.the viscosity range, the values used for the slow deformation region, the viscosities correspond roughly to that of honey that has been re­ frigerated. Also shown in Appendix C are calculated gravitational 95 viscous, deformations for this material. A last observation is that the tangential boundary condition used in BVSMAC at the bottom boundary was the no-slip condition. The quality of the modeling results .lends credence to the hypothesis that this boundary condition.is appropriate, for flowing, snow. In carrying out the computer, modeling the timer-step between calculational cycles was chosen such that the maximum distance traveled by any part of the fluid remained less than oner-tenth of a cell dimension Using the cell dimensions previously described and this time step cri­ terion, no numerical instabilities were encountered for the range of parameters involved in this modeling. To generate the results exhibited in Appendix B, each modeling run required about 1000 calculational cycles, taking, for the 1400 cell computational system, about 30 minutes of. CPU time on a DEC system VAX 11/780. CHAPTER. VI SUMMARY AND CONCLUSIONS For snow avalanche motion in the speed range below. 20. m/s the biviscous model has provided satisfactory results.. The overall motion of the snow as depicted in the leading.edge motion and final debris depth distribution were well simulated... In addition, details noted in the snow tests were reproduced by the computer model. Quantitatively, the velocity as a function of depth was accurately modeled. Qualita­ tively, the formation of the boundary layer can be seen in the time sequence particle plots in Appendix B. The particles near the front of the flow at the bottom are retarded as the upper part of the flow proceeds over them. These particles are seen to be left in a layer along the bottom boundary, just as the dye placed originally at the front of the flow was seen to be distributed as a layer over the. entire runout area. Examination of the motion of the marker particles in the upper portions of the flow shows that little deformation is taking place in this region. This motion is confirmed by observations and by the small relative displacements of the dye concentrations in the snow tests for dye placed on the surface. Another, aspect ,of the flow seen in both the snow tests and.the computer,modeling is surging motion of the.leading edge. Although it doesn’t show .clearly in the particle.plots, the front of the flowing mass was continually breaking , ' . over the slower moving flow near the surface. This.motion showed up 97 most strikingly by monitoring.the velocity at the leading edge. Xt was found that this velocity was not a smooth, function of. time but exhibited large variations around the average decaying velocity. At the end of each of the., particle plots, in Appendix B are two plots generated by the computer,at the time of execution. They show the leading edge position versus time and the leading edge velocity versus time. The velocity plot shows this surging motion clearly. This motion was also seen when finding the leading edge velocity in the snow flow tests from the 16 mm film. It showed up as anomalous leading edge velocity measurements at sporadic times in the flow. It could also be seen viewing the motion picture film, as surging or jetting of the leading edge, much like the motion of water waves shoaling on a beach after breaking. At speeds above 20 m/s, much conjecture still exists as to what the behavior of flowing snow is. Perla (1980), Mellor (1978), and others have speculated that the flowing snow enters a turbulent flow regime at high speeds. This would be.due to individual snow grains attaining enough energy so that inertial effects become dominant. The collisions between these grains, would become the principal ..energy dissipation mechanism. For this motion to take place,.the snow would have to become.fluidized,, and dilatation of the flowing mass would occur, reducing the mass, density,,.. At this point, turbulent effects would be observed. As evidence of this behavior, the motion of the 98 snow near the front of large, high, velocity avalanches.,is cited, it is seen that this material appears to. he highly turbulent.. .However, it must be pointed out that.what is observed is the very low density dust cloud accompanying the avalanche.. The motion of the.material beneath, this cloud cannot be seen. ..Also cited as evidence of turbulent motion is the existence of forces proportional, to the square of the 2 speed. Most analysis of avalanche motion have proposed these v forces and then derived steady-state equations for the terminal velocity (see equations 3.4 and 3.8). These equations exhibit a square root dependence on the sine and cosine of the slope angle. It is stated that this square root dependence is better represented by data than the linear dependence on the slope angle that would result from forces linear in the velocity (Schaerer, 1975) .. This argument originates from steady state assumptions, ignoring time dependent effects, and could be questioned on that basis. At present, the mechanics of high speed avalanches is not known. The flow could well be of a turbulent nature, but.data is not available to confirm or deny.this. However, for speeds below 20 m/s, no evidence of turbulence motion was observed.in ,the snow flow tests- The dye concentrations, placed, in. the snow did not difuse, but remained in,concentrated form, except for the. dispersion.along the lower boundary.. .Also the linear biviscous model did a very credible job in simulating this motion. Additional experimental data is. needed. 99 Jnvestigatiqn of flow velocities and density distributions in.high speed avalanches must be made, to clarify the mechanics involved. Once this data is collected and available, proposed flow mechanisms can be evaluated. A very useful feature of the computer, model de­ scribed in this thesis, is that. it. may easily be adapted to other flow laws. It is a simple matter, to generalize from a biviscous formulation to a tri-viscous formulation or a flow law that involves more viscosities. In this way even a velocity-squared type of constitutive equation could be approximated. It is the author's opinion that a three-viscosity flow law would probably be more than adequate to provide reasonable simulation of full scale avalanche motion. However, until more data is collected these speculations are academic. Theoretical investigation of grain flow mechanics should also be adapted to snow avalanches. Applying and extending the work done in granular theory could provide a theoretical framework from which constitutive relations could be proposed. These relations could then be tested using the. computer model. The computer simulation methodology, particularly the multi­ viscosity approach, has proved to.be a viable.design for solving flowing.snow problems. The excellent internal consistency shown by the biviscous modeling of the snow flow tests inspires a great deal of confidence in the method. 'The adaptability of the model to other /1 0 0 physical systems, such as rock slides, and mud slides, ,also.should not be overlooked. Many of the same1 types.of problems exist .for other geo physical phenomenon that exist for flowing snow. In addition, the ease with which the model, may be generalized to include more complex constitutive laws indicates a very good, prospect for use as additional information about the mechanics of the flowing material is learned. REFERENCES •Ames, w. F., 1977.. Numerical Methods for. Partial Differential Equations. •Academic Press. Amsden, A. A., and Harlow, F . H. 1970. The SMAC Methods A numerical Technique for Calculating Incompressible Fluid Flows. Los Alamos Scientific Laboratory. Report No. LA-4370. Bagnold, R. A., 1954. Experiments on a gravity-free dispersion of large solid spheres in a Newtonian fluid under shear. .proceedings of the Royal Society of London, A, 225, p. 49-63. Bagnold, R. A. 19.56, The flow of Chesionless Grains in Fluids. Proceedings of the Royal Society of London, A, 249, p. 29-297. Bakkehjdi, S., Cheng, T., Domaas, V., Lied, K., Perla, R., and Schieldrop, B., 1980. On the Computation of Parameters that Model Snow Ava­ lanche Motion. Submitted to Canadian Geotechnical Journal. Briukhanov, A. V., Grigorian, S. S., Miagkov, S. M., Plam, M. Ya., Shurova, I. Ya., Eglit, M. E. and Yakimov, Yu. L., 1967. On Some New Approaches to Dynamics of Snow Avalanches. Phys. of Snow Ice, Conf. Proc., Vol. I, Part 2, p. 1223-1242. Brugnot, G. 1980. Recent Progress and New Applications of the Dynamics of Avalanches. Presented on Snow in Motion Symposium, Fort Collins, Colorado. 12-17 August 1979. Bucher, E. and Roch, A., 1946. Reibungs und Packungswidersland bei raschen Schneebeweguagen (Friction and resistance to compaction of snow under rapid motion), Mitteilungen des Eidgenoss. Institut fur Schnee und Lawinenforschung, Davos-Weissfluhjoch, 9 p. Buser, Othmar and Fruitiger, Hans, 1980. Observed Maximum Rpn-out Dis­ tance of Snow Avalanches and the Determination of the Friction Coefficients y and E, J. Glac.26, 94, p. 121-130. Dent, J. D. and Lang, T. E., 1980. Modeling of Snow Flow. J. Glac. 26, 94, p. 131-140. EybertTrBerard,A., Perroud, P,, Brugnot, G,, Mura, R., Rel, L., 1978, ■ Megures Dynamigues Dans L ' Avalanche. Compte.s Rendus, Deuxieme Rencontre Internationale Sur La Neige Et Les Avalanches, Grenoble, France, p . 203-224. 102 Goodman, M. A. and Cowin, S. C., 1971. Two Problems in the Gravity Flow of Granular Materials.. J. . Fluid Mech. 45,. p. 321-339. Goodman, M. A- and Cowin, S. C., 1972. A Continuum Theory, for Granular Materials. Arch. Pat.'Mech. Anal. 44, p. 249-266. Gubler, H., 1980. Personal Communication. Heimgartner, 'M., 1977. On the Flow of Avalanching Snow. J. Glac. 19, 81, p. 347-363. Hirt, C. W., 1968. . Heuristic Stability Theory for Finite-Difference Equations. J. Comp. Phys. 2, p. 339-355. Hirt, C. W. and J. P. Shannon., 1968. Free-Surface Stress.Conditions for for Incompressible Flow Calculations. J. Comp. Phys. 2, p. 403- 411. Hohenemser, K. and W. Prager, 1932. Uber Die Ansatze Der Mechanik Isotroper Kontinua. Zeits, Angew. Math. U. Mech., 12, p. 216-226. Hopfinger, E. J. and Tochon-Danguy, J. C., 1977. A Model Study of Pow- . der Show Avalanches. J. Glac.19, 81, p. 343-356. Inaho, Y., 1941. Sekisetsu No Do-Masatsukaku Ni Tsuite [Angle of Kinetic Friction of Snow]. SEPPYO, 3, p. 303-307. Kanatani, K . , 1979. A Micropolar Continuum Theory for the Flow of Granular Materials. Int. J. Engng. Sci. 17, p. 419-432. Korner, H. J., 1981. The Energy Line Method in the Mechanics of Ava­ lanches. J. Glac. 26, 94 p . 493. LaChappelle, E. R. and Lang, T. E., 1979. A Comparison of Observed and Calculated Avalanche Velocities. J . Glac♦,. 25, 92,. p. 309-314. Lang, T. E. and Dent, J., 1980. Scale Modeling of Snow-Avalanche Impact on Structures.. J. Glac.26, .94, p. 189-196. Lang, T. E.., Dent, J. and Oakberg, R.. G., 1981. Surface Resistance Measurements Cooperative Agreement RM-81-165-CA U. S. D, A* Forest Service, Rocky Mountain Forest and Range Experiment Station (in progress). 103 Leaf, C. F. and Martinelli, M., Jr., 1977. Avalanche Dynamics: Engineering Applications for Land Use Planning, U . Dept, of Agriculture., Forest Service Research Paper RM-183. Maeno, N. and Nishimura, 1979. Fluidization of Snow. .Cold. Regions ' Science and Technology, I, p. 109-120... Maeno, N., Nishimura, K., and Kaneda, Y., 1980. Viscosity and Heat Transfer in Fluidized Snow. J..Glac. 26, 94, p. 263-274. Malvern, L. E., 1969. Introduction to the Mechanics.of. a Continuous Medium. Prentice Hall. Martinelli, M., Jr., Lang, T. E. and Mears, A. I., 1980. Calculations of Avalanche Friction Coefficients from Field.Data. J. Glac. 26, 94> p. 109-120. Mears, A. I., 1976. Guidelines and Methods for Detailed Snow Avalanche Hazard Investigations in Colorado. Colorado. Geological Survey.. Bulletin 38. Mears, A. I., 1980. A Fragment-Flow Model of Dry-Snow Avalanches. J. Glac. 26, 94, p. 152-164. Mellor, M., 1968. Avalanches. Monograph III A3d. U. S . Army Cold Regions Research Engineering Lab., Hanover, New Hampshire. Nedderman, R. M., and Laohakul, C., 1980. The Thickness of the Shear Zone of Flowing Granular Materials. Powder Tech. 25, p. 91-100. Nichols, B. D. and C. W. Hirt, 1971. Improved Free Surface Boundary Conditions for Numerical Incompressible-Flow Calculations. J. Comp. Phys. 8, p. 434-448. Pao, R. H. F., 1961. Fluid Mechanics. John Wiley and Sons. Perla, R. I., Cheng, T. T., and McClung, D. M., 1980. A Two-Parameter Model of Snow-rAvalanche Motion. J. Glac ■ 26, 9.4, p. 197-20,8. Perla, R. I., 1980. Avalanche Release, Motion, and Impact. In ■'Dynamics of Snow and Ice (S. Colbeck, Ed.I Academic Press. Salm, B., 19,66. Contribution to Avalanche Dynamics. Lnt. Assoc. Sci. Hydrol. Publ. 69, p. 199-214. ■104 S^VEtge, g. g,, 1979.. Gravity Flow Of Cohesionleg^ ,Gpanulaar Materials in CIiutes and. Channels. ' J. Fluid Mech. 92, Part. I, p. 53-9.6. Schaerer, P. A., 1973. Observations of Avalanche ,impact. Pressure. ' Proc. Symp. Adv. North Am. Avilanche Technol. p. 51-54. Schaerer, P. A., 1975. Friction Coefficients and. Speed pf Flowing Avalanches. [Union.Geodesigue et Geophysique, Internationale. Association Internationale des Sciences.Hydrologigues. Commission des Neiges et Glaces.] .Symposium. Mecanique de la neige. Actes du collogue de Grindplwald, avril 1974, p. 425-32 (IAHS-AISH Publication No. 114.) Schaerer, P. A. and Salway, A. A., 1980. Seismic and ,Impact-Pressure Monitoring of Flowing Avalanches. J. Glac. 26, 94, p. 139-188. Shen, H. W., and Roper, A. T., 1970. Dynamics of Snow Avalanches (With Estimation for Force on a Bridge). Bulletin of the International Association of Scientific Hydrology, 15, I, p. 7-26. Shimizu, H., Akitaya, E., Huzioka, T., Nakagawa, M., Kawada, K., 1973. Kurobekyokoku kosoku-nadare no kenkyu. II [Study of High-Speed Avalanches in Kurobe Canyon. II]. Teion-kagaku: Low Temperature Science, A, 33, p. 109-16. Shimizu, H., Huzioka, T. [i.e. Fukioka] ,. Akitaya, E.,. Narita, H., Nakagawa, M., Kawada, K., 1977. Kurobe-kyokoku kosoku-nadare no kenkyu. V [Study of High-Speed Avalanches in Kurobe Canyon. V]. Teion-kagaku: Low Temperature Science, 26, 94, p. 141-152. Shimizu, H. E., Huzioka, T., Akitaya, E., Narita, H., Nakagawa, M., Kawada, K.>1981. A Study on High Speed Avalanches in the Kurobe Canyon Japan. J. Glac. 26, 94, p. 141-152. Shoda, M., 1966. An Experimental Study on Dynamics of Avalanches Snow. IASH. Publ. 69, p. 215-229. Sommerhalder, E., 1972. Deflecting structures-. In "Lawinenschultz in der Schweiz" (F. Castelberg, H. R. in der gand., F , Pfigter, B.. Rageth.,. and G. Yavier» edg-L. Beih. No.. 9.. Genoosenschaft der bundnerischen'HoIsproduzenten. [English t r a n s l . ."Avalanche Protection in Switzerland.'! u.. s. Dept, of Agriculture. Forest Service. Gen. Tech. Rep. RM-9(1975].] 105 Tochon^Panguy, J. C. and Hopfinger, E,. J., 1975. Simulation of the Dynamics of Powder Avalanches. TAHSrAISH. Publ. . 114, p. 369-380. Voellmy, A., 1955. TJher die Zerstorungskraft von. Lawinen. Schweizeri- sche Bauzeitung, Jahrg. 73, Ht. 12, p. 159-62; .Ht. 15, p. 212-17; Ht. 17, p. 280-85. [English.Translation: On the Destructive Force of Avalanches. U.. g. Dept, of Agriculture. Forest Service. Alta Avalanche Study.Center Translation No. 2, 1964.] APPENDICES APPENDIX A BVSMAC FORTRAN COMPUTER CODE o n o o n o n o o n n n n n n o o o o n n n o n n n o n n n n n n n n n n n n n n n n n o n o o o o o o o o o n o 107 I 4-May-I 982 I 4-May-I 982 13:34:35 13:31:52 BVSMAC-SIMPL IFIED MARKER AND CELL-CODE TO SOLVE 2-0 NAVIER-STOKES EONS. ADAPTED FROM AMSDEN AND HARLOW, LOS ALAMOS REPORT LA437# ' MODIFICATIONS BY OIM DENT INITIATION & SET-UP ON LOCAL SYSTEM 10/77 I-DIM ARRAYS 7/79 SURFACE MARKERS 11/81 BIVISCOUS BINGHAM FORMULIZATION 3/82 LOGICAL ASSIGNMENTS " F: 8 INPUT F:10 NUMBER OUTPUT F : I 2 TERMINAL OUTPUT F:13 TERMINAL INPUT (BLANK) F: 15 PLOT OUTPUT LPR OUTPUT CONTROL 0 NO OUTPUT 1 F : 12 TERMINAL ONLY 2 F:10 + F:12 TERMINAL AND NUMBER OUTPUT 3 F: 10 .NUMBERS' ONLY INPUT VARIABLES ( F : 8.) ft * * * K * ft ft ft * * ft ft ftftftftftftftft LINE I (2I5,3F8.3,I5,2F8.3,F8.6) IBAR -NO. CELLS IN X-DIR JBAR - Y-DIR DR -X CELL SIZE DZ -Y "• " , TDT -TIME OF INFLOW SLOPE CHANGE TCYCLE-TOTAL NO. OF CYCLES TO BE DONE PC -LENGTH OF NO-SLIP SURFACE:TIME OF INFLOW ALP -RELAXATION PARAMETER (0-1) EPS -CONVERGENCE CRITERION (.0002) LINE 2 <12A4 > IDENTIFYING NAME (48 CHARACTERS > LINE 3 <4F3.0,3F8.3,5F7.3) BCB -BOTTOM TANGENTIAL B:C. I=FREE SLIP BCR -RIGHT -I=NO SLIP BCT -TOP 0=DYNAMIC BC BCL -LEFT A -BOTTOM TAN BC PARAMETER INFLOW VELOCITY REDUCTION AFTER I SEC. B -BOTTOM TAN BC PARAMETER INITIAL INFLOW HEIGHT CHANGE SPEED C -FINAL INFLOW HEIGHT CHANGE SPEED NU -PRIMARY VISCOSITY NUP -BINGHAM VISCOSITY . . GR -X GRAVITY COMPONENT GZ -Y GAM -BINGHAM COEFF. LINE 4 (7F8.3,12) T -START TIME TWPLT -TIME BETWEEN PLOT OUTPUTS TWPRT -TIME BETWEEN NUMBER OUTPUTS TMIN -MINIMUM TIME STEP TMAX -MAXIMUM TIME STEP TFAC -MAX % OF CELL TRANSITION . . TWF IN -FINISH TIME LPR -SEE LPR OUTPUT CONTROL ABOVE LINE 5-7 <1015,2F8.3) o n o o o o o o n o o o o o o o o n o o o n o o o o o n o o o o o o n o o n o o n n o o -108 i 14-May-1 982 .13:34:35 14-May-1982 13:31:52 TYPE -Sf=NO INFLOW OR OBS >l = INFLOW OBS (TYPE=I LAST OBS) L1-L7 -OBS DIMENSIONS BY CELL NO. NXB -NO. PARTICLES/CELL INFLOW IN X-DIRECTION-- NYB - ' Y-DIRECTION I UL -LEFT INFLOW VELOCITY I TYPE=I CARD UR -RIGHT OUTFLOW VELOCITY (ST = CONT OUTFLOW >- EZ LINE 6-7 (2 I 5,6 F 8.3 ) PARTSECARDS NX -NO. PARTICLES/CELL IN X-DIRECTION (Sf=NO MORE PARTS) NY - Y-DIRECTION (Sf=NO MARKERS) . XC -PARTS DIMENSIONS (REAL UNITS) Y C - XD - YD USf -INITIAL PARTS VELOCITY IN X-DIRECTION LINE 7 ( IX,FlBf.5, 15) DSP -1/NO. OF SURFACE PARTICLES PER CELLf.3 > IF DSP=Sf. NO SURFACE MARK ERS-ELIMI NATE NEXT CARD NPTS -NO. OF SURFACE NODE POINTS LINE 8-?(2X,8F8.3 ) XA1YA -X1Y POSITION OF SURFACE NODE POINTSfCELL DIM) LAST LINE (2I5.F8.3) IBAR - JBAR - DR -OR,A2(2500),A3< 2500),A4(2500).AS(2500I1A6(2500), 1 A7(2500).A8(2500),A9(2500),A10(2500),A1ALP,B,BCB1BCL1BCLTf99),BCR 2 ,BCRT(99 J1BCT1BND,C1COF I,COFZ1COLS,CYCLE,DR1DRODZ1DROU,DSP,DT1DZ1 3 DZODR.EMP,EPS1FUL,GAM1GAMP,GR,GRAD02,GZ1Ht300),I H 2500),12(2500), 4 IALL,IBAR1IFK 300).IPl,IP2.IPHM1JBAR1JPl,JP2,LPB1LPR,NAME(12 >, 5 NIN.NP,MS,NU,NUP,OB,PC,RDR1RDRZ,RDRDZ1RDZ1RDZ2, 6 ROS<99>.SUR,T1TCYCLE,TDT1TFAC,TMAX,TMIN,TPLT.TPRT.TWFIN, 7 TWPLT,TWPRT,TYPE1UL1UR1W1XA(30J1XDIS1XL1 B XP(4000),XR1XSf 4000),YA< 30).YB,VDIS1VFIR1YT1IEND EQUIVALENCE (Al,F ) i(A2.U >,(A3,V ),(A4,UTIL >,(AS,VTIL),(A6,PS I>, * (AT1THETA)1(AS1D)KAglP ).( Il1KF) REAL NU1NUP INTEGER CYCLE,TYPE,TCVCLE I I ■109 14-May-19.82 14-May-l982 I FORMAT(2I5,3F8.3,I5,2F8.3,F8.6> LPB = A8S3B BND=I .J0f FUL=2.0 SUR=3.J0f EMP = 4.0 OB=B.0 10 READ <8,1 > I BAR, JBAR1DR1DZ, TDT1TC VCLE1 PC, ALP, EP1S IP2=IBAR+2 JP2=JBAR+2 IF(DR) 30,20,20 20 CALL PROG ( F , U ,V1 UTIL ,VTI-L , PSI ,THETA, D, KF , P ) GOTO Iflr 30 CONTINUE END 13:34:35 13:31:52 14-M^y-1 982. 13:34:35 14-May-1 982 13:31 :52 SUBROUTINE PROG ( F,U,V,UTIL,VTIL,PSI,THETA,D,KF,P)' COMMON/ZC/Al<25M),A2{25M>,A3(250j0r),A4(25flr0),A5(25M),A6(2W>, 1 A7 < 2500), A8( 2500 >, A9( 2500 ).,A10( 25001, A, ALP, B, BOB, BCL,BCLT( 99), BCR 2 ,BCRT(99>,BCT,BND,C,C0F1,COF2,COLS,CYCLE,DR,DROOZ.DROU.DSP,DT1DZ, 3 DZODR ,EMP , EPS ,FUL ,GAM,GAMP ,GR ,GRAD02 ,GZ, H( 3001', 11 ( 25001,12< 25001, 4 IALL,I BAR,IFK 3001,IPl,IP2,IP HM1JBAR,JPl,JP2,LPB,LPR,NAME!12), 5 NIN,NP,NS,NU1NUP,OB,PC,RDR., RDR2,RDRDZ1RDZ,RDZ2, 6 ROS(99 I,SUR1T1TCVCLE.TDT.TFAC,TMAX1TMIN1TPLT,TPRT,.TWFIN, 7 TWPLT,TWPRT,TYPE,UL,UR,W1XA(30),XDIS,XL, 8 XP(4000),XR,XS(4000),YA(30>,YB,VDIS.YFIR,YT1IEND DIMENSION F(IP2.JP2I,U(IPZ1OPZl1VtIP2,OPZl1UTILtIPZ1OPZl1 * VTILtIPZ1OPZ I.PSItIPZ1OPZ I,THETAtIPZ1OPZI1DtIPZ1OPZ I, * KF ) 15 FORMATtF8.3,15 I 16 FORMATt/.,2X,'MORE CYCLES? (15)') 17 FORMATt/,2X,'MORE TIME? (F8.3)') 18 FORMATtIX,F10.5,15) 19 FORMATt ' T='IPE12.5.' CYCLE='15,' ITER='215,' PLE='E10.3, . I ' VLE='E10.3) 20 FORMATt1015,2F8.3) 21 FORMATt ' DSP=‘F10.5,1 NPTS='15,/,' SURFACE :1 PARTICLE GRID POINTS X-Y' ) 22 FORMAT!'0 NO OBST OR 1/0 BOUNDARIES') 23 FORMATt ' TYPE='12,' LI = ' 15,' L2=' 15, ' L3= ' 15,' 14='15, 1 ' LS='15,' LS='15,' L7= ' 15, ' NXB='15,' NYB='15,' UL='F8.3, 2 ' UR='F8.3 ) 24 FORMATtZX,' T='1PE12.5, ' CYCLE '15) 25 FORMATt 2X,I2F8.3 ) 26 FORMATtIX,14) 27 FORMATtIX,12011) 28 FORMATtIHl ) 29 FORMATtZX,' ITER='15.' MAX NO. ITERATIONS EXCEEDED' ) 30 FORMAT!12X,15,' SURFACE MARKER PARTICLES') 31 FORMATt3X,2(2X,14 I1ZtZX1IPElZ.5 I1ZX115) 32 FORMATt2X.8F8.3) Ill 14-May-1982 13:34:35 I 4-May-1 982 13:31 :52 WRITE (IjBftI) IF (OR.EQ.0)GOTO 460 READ (8,2) (NAME(I),1=1,12) READ (8,3 )BCB,BCR,BCT.BCL1A1B,C1NU,NUP1GR,GZ1GAM . IF(PC.EQ.0. )BCL = 1 .0 READ (8,4) T , TWP L T, TWP R.T, TM I N , TMAX , TFAC , TWF IN , L PR WRI TE<15,2 ) (NAME(I),1 = 1,12) WRITE<15,5)IBAR1JBAR1 DR,DZ1TDT1TCYCLE1PC1ALP,EPS WRITE!I 5,6)BCB,BCR,BCT,BCL1 A,B1C1NU1NUP1GR1GZ1 GAM. WRI TE(15,7)T,TWPLT1TWPRT1TMIN1TMAX1TFAC1TWF IN,LPR I F ( L PR . EQ.. 0 )GOTO 34 IF(LPR.EQ.3)G0T0 33. WRITE! 12,2 ) (NAME( I ), I = I , 12 >• WRITE(12,5)IBAR1JBAR1DR1DZ1TDT1TCVCLE,PC1ALP1 EPS WRITE!12,6)BCB,BCR,BCT,BCL,A1B1C1NU1NUP1GR1GZ1GAM WRITE!12,7) T,TWPLT,TWPRT,TMIN,TMAX,TFAC,TWF IN,LPR IF(LPR.EQ.I>GOTO 34 33 WRITE (10,2) (NAME! D 1I = I1IE) WRITE (10,5) I BAR.JBAR,DR.DZ.TDT.TCYCLE,PC,ALP.EPS WRITE (1 0,6 ) BcB1BcR1BCT1BCL1A1B1C1NU1NUP1GR1GZ1GAM WRITE (10.7) T,TWPLT,TWPRT1TMIN1TMAX,TFAC,TWFIN1LPR 34 CONTINUE IPl=IBAR+! JPl=JBAR+! IP3=IP2+1 IALL=IP2*JP2 IP4-IP2+2 . IAAL=IALL-IP3 RDR=I./DR RDR2=RDR*RDR RDZ=I./DZ RDZ2=RDZ*RDZ DRODZ=DR"RDZ DZODR=DZ*RDR RDRDZ=RDRMRDZ C0F1=2.*NUP*RDR COF2=2.*NUP*RDZ GAMP=NUP/!NUP-NU)*GAM GWMNON=GAMP*!I.-NU/NUP) GRAD0=GAMP/NUP GRAD02=GRAD0*GRAD0 W=I./< 2.*!RDR2 + RDZ2 ) ) TPLT=TWPLT • TPRT=TWPRT TWPLT=TPLT+T TWPRT=TPRT+T CYCLE=0 XR=IBAR*DR YT=JBARMDZ I FS = PC*RDR + 2 I F < B.NE.0..AND.GZ.NE.0.>ALF = 2.*NU/< B*GZ*DZ> AL FS=ALF*ALF IF!GZ.NE.0. )GTE = -(GR+GZ*A)/GZ DO 35 I=I1IALL Al!I )=4. : A2(I)=0. A3!I)=0. 112 14-May-1.982 13:34:35 14-May-1982 13:31:52 A4=0. A5(I ) =0. A6(I)=0. A7(I)=0. A9(I>=0. A10( I J=-DR I 2< I ) = I 35 A8(I>=0. ' DO 36 J=I,JP2 F( I,J J= I . 36 F < 'I P2 , J J = I . DO 37 I = I , IP2 Al(I)=I. 37 AK I ALL - I +1 ) = 1 . READ (8,20) TYPE1Ll,L2,L3,L4,L5,L6,L7,NX8,NVB,UL,OR DO 38 J=I ,JP2 BCLT(J)=BCL ROS(J)=1.0 38 BCRT(J)=BCR . K=I NP=0 I F < TYPE.NE.0)COTO 100 IF(LPR.GE.2) WRITE (10,22) I F(LPR.EQ.I.OR.LPR.EQ.2 > WRITE!12,22 > C C GENERATE MARKER PARTICLES 40 READ (8,8) NX1NY1XC1,YC1XD1YD1UO1VO I F(NX . EQ.0 JGOTO 80 WRI TE(15,9)NX,NY1XC1YC1XD,YD,UO1VO IF( LPR.EQ.I.0R.LPR.EQ.2) WRITE!12,9) NX,NV1XC,VC,XD1VD1UO1VQ IF( LPR.GE .(2 ) WRI TE ( 10,9) NX , N Y, XC , YC , XD , YD , UO, VO I F(NY.NE.0 IGO TO 49 I=XC * RDR+ 2.- IE=XD*RDR+2. JE=YD"RDZ+2. 42 J=YC*KDZ+2. . 44 IF(F(I,J >.EQ.OB >G0 TO 45 IF(F( I ,J).EQ.BND)GO TO 45 F(I1J)=FUL IF( F( I + i, J ) .EQ.BND) U(I1J) = Ud , J > »ROS( J ) + U0*( I.-ROS< J > > IF(F<1 + 1,J).NE.BND.AND.F(1 + 1,J I.NE.OB) U(I1J)=UO IF(F( I-I ,J ). NE . BND . AND . F'( I-I ,J J.NE.OB) U(I-I1J) = UO IF( F< I ,J-I ) .NE .BND.AND.F( I , J-I ) .NE.OB ) V d 1J-I >=.V0 IF( Fd ,J+l ) .NE .BND) ' V(I1J)=VO 45 J=J+I IF(J.LE.JE)GO TO 44 1 = 1 + 1 IF( I .LE.IE >G0 TO 42 GO TO 40 49 XC=XC*RDR YC=YC*RDZ XD=XDrtRDR ■ YD=VDrtRDZ XTE=I./NX VTE=I./NV Y=YTE*.5 1X3 14-May-198.2 13:34:35 14-May-1982 13:31 :52 SZ X=XTE*.5 55 IF U(I1O) = UO IF(FII-I1O).NE.BND. AN D.'FU-I1O >.NE. OB). U(I-I1O) = UO IFtF(I1O-I).NE.BND.AND.F(I1O-I).NE.OB) VII1O-I)=V0 IF(FfI1O+!).NE.BND) V(I1O)=VO .70 X=X+XTE I F(X.L T.I BAR) GO TO 55 Y=Y+YTE IF(Y.LT.OBAR)GQTQ50 GOTO40 72 V I = (Y-YC )**2 Xl=(X-XC )**2 IF((Y1+X1).LE.(XD**2))G0T0 60 GOTO 70 C END MARKER PARTICLE GENERATION 74 WRITE!10,29 ) ITER WRITE!12,29)ITER v . IND=-99 WRI TE(I 5,26 ) I ND WRITE<15,26 )INO GO TO 999 75 WRITE (10.11) ' WRITEt12,11> GOTO 999 80 IF(LPR.GE.2 ) WRITE (10,10) NP IF(LPR.EO.I.0R.LPR.EQ.2)WRITE(12,10)NP WRITE!15,10)NP C GENERATE SURFACE MARKERS NS=0 KS=I READ(S1IS)DSP1NPTS IF(DSP.EQ.0. )XS(2) = L2 IF(DSP.EQ.0.>G0 TO 86 DSPSL=DSP*DSP*1.21 DSPSS=DSP*DSP*.25 DSPSS4=DSPSS*4. READ( 8,32 )( XA( D 1YAU J1I = I1NPTS) X=XA(I> Y=YA(I) DO 85 1=2,NPTS XDF=XA(I)-XA(I-I) VDF=YA(I)-VA(I-I) IFtXDF.EQ.0.)G0 TO 82 U U U 114 I 4-May-1 982. 13:34:39 I 4-May-I 982 13:31:52 SLP=YDFZXDF SEE = YAM I-SLPft-XAl I ) DELX = DSPZSQRTl I.+SLPftSLP) 81 XS(KS)=X XS1KS+1>=Y IA=X+2. 0=Y+2. IFlFlIA,0>.EQ.EMPlFlIA,J)=FVL KS=KS+2 I F I KS.GT.4 0 0 0)G0 TO 75 NS = NS+ I IFI X.GE.XAl I ) >GO TO 85 X=X+DELX I F I X.GT.XAl I ) IX = XAl I> Y=SLPftX+SEE GO TO 81 82 SLPP=XDFZYDF SEEP = XAl I>-SLPPftYA)Y = VA(I ) X=SLPPftY+SEEP GO TO 83 84 I F I Y.LE.YAl I ) >GO TO 85 Y=Y-DELY IFlY.LT.YAl I>)Y = YAl I> X=SLPP*Y+SEEP GO TO 83 85 CONTINUE WL K DK GV hAGEAY e sO EO Ww rhKEkd mpb;Y K ,l'V ;uVE' WRITE!10,32)1XAl I),YAl I),I = I1NPTS) WRITE!10,30)NS 87 IFlLPR.LT.I.OR.LPR.GT.2 >GO TO 88 WRITE!12,21)DSP,NPTS WRITE! 12,32 HXAl I ), YAl I ),I = I1NPTS) WRITE! 12,30)NS- 88 WRITE!15,21IDSP1NPTS WRITE!15,32)1XAl I),VAl I),I = I1NPTS) ■ WRITE!15,30)NS GO TO 250 GENERATE OBSTACLES 100 CONTINUE 101 WRITE!15,23 >TYPE1Ll,L2,L3,L4,L5,L6,L7,NXB1NVB1UL,UR IFlLPR.GE.2)WRITEl10,23)TVPE,L I,L2, n o n 1X5 I 4-May-I 982 13:34 .'35 14-May-1982 13:31:52 * L3,L4,L5,L6iL7,NX9,NYB,UL,UR IF< LPR.EQ. I .0R.LPR.EQ.2 )WRITEU2,23 ) * TYPE ,LI ,L2,1.3, L4,L5,L6,L7, NXBtNYB1UL1UR IF(L7.EQ.0)GOTO 105 IL=L5*2 IR=L6+1 UTE-L7+1 DO 104 J=Z1JTE DO 104 I=JL,IR. 104 FII1J>=0B IFI TYPE,EQ.I>G0 TO 105 READ!8,20ITYPE1L I,L2.L3,L4,L5,L6,L7,NXB,NYB,UL1UR GO TO 101 105 CONTINUE XD IS=0.5/AMAX0INXB1I> YDIS=I.0/AMAX0(NYB,l> YFIR=I .5*YDIS) + L1 NIN = NYB*!L2-L1 ) COL S = 0.0 C SET INFLOW AND OUTFLOW VELOCITIES DO 110 J = I ,JP2 IFIJ.GE.ILl + l ).AND.J.LE.IL2+1)) GOTO 107 106 IFIJ.GE.IL3+2).AND.J.LE.IL4+1)) GOTO 109 GOTO 110 107 UI I ,J ) = UL BCLTIJ I = -I.0 GOTO 106 109 UIIPl1J> = UR BCRT(J)=-I. IFI UR.NE.0. IGOTOl10 ROSIJ)=0. BCRT(J)=I. . 110 CONTINUE DO 115 I=IP4,IAAL IF(A1(I).NE.0B)G0T0 115 IFIAlII + IP2>.NE.OB.AND.All 1 + 1J.EQ.OB) AZ(I)=AZI I+ IPZJftBCB IFIAlII-I ).NE.OB.AND.Al I I+ IPZ).EQ.OB> A3II)=A3(I-IJftBCB IFIAlI I +I ).NE.OB.AND.Al II + IPZ).EQ.OB) A3II)=A3(1 + 1)ftBCB 115 CONTINUE GOTO 40 C * C #***Aw*****A*in******A**************aaaaa*aa**#a*#**aaae******* C * C * C * START OF COMPUTATIONAL CYCLE C * C * C **a«:i**a*)'(« 253 K=K+2 NSP=NSP+! X = XN Y = YN IF(NSP.LE.NS)GO TO 254 GO TO 260 254 XN=XS(K) YN = XS(K+l ) 255 USS=(YN-Y)**2+(XN-X)**2 IFtDSS-CT-OSPSL)G0 TO 2S6 IF(DSS.LT.DSPSS)GO TO 258 GO TO 253 C SPACING TOO GREAT-INSERT PARTICLE 256 XN = (XN + X )* . 5 VN = < YN + V )*-5 OO 257 I=NSP,NS IB=NS+NSP-I+I IKS= I 8+IB- I XS = XS(IKS + 3> NS=NS-I . IF(NSP.LT.NS)G0 TO 254 C CHECK SPACING BETWEEN EVERY OTHER PARTICLE 260 K=I ' N S P = 2 261 NSP=NSP+! IF(NSP-LE-NS)GO TO 263 GO TO 271 263 X=XS(K) Y = XS(K+l ) K=K + 2 264 XN = XS(K + 2 ) YN = XS< K+3 ) DSS = < VN-Y >**2 +(XN-X >**2 IF< DSS-GE-DSPSS )G0 TO 261 DO 265 'I =NSP, NS IKS=2*I-3 XS(IKS)=XStIKS+2) 265 XS(IKS+1 ) = XS(IKS + 3 ) NS=NS-I KDSu'VaGEau'esg EO YLf C w***aA**wnn***********dA****aa********=****a*wa**A*a**#a#*AaAa C * 117 14-May-1 982 13:34:35 I 4-May-1 982 13:31 :52 C " REFLAGGING C * C * £ * DETERMINE WHICH CELLS HAVE PARTS, SET INDICATOR 271 GO TO KRET 273 NPT=O K=I 274 IFINPT.GE.NPlGO TO 275 I=XPIK1+2. J=XPIK+!1+2. KFIItOl=K . K = K + 2 NPT=NPT+! GO TO 274 275 NPT=O KS=I 276 IFINPT.GE.NS IGO TO 277 I=XSIKS 1 + 2. J = XSIKS+I 1 + 2. KF(ItJ)=K K = K+2 KS=KS+2 NPT=NPT+! GO TO 276 C * C * PREFLAGGER C * 277 IFI NP.NE.OlGO TO 283 IFIDSP.EO.O.IGO'TO 283 DO 280 1=2, IPl CSTAT=EMP DO 279 J=EtJPl OW=OP I+2-0 • IFIFI I.JWl.EQ.OBIGO TO 279 IFIKFI I,JWl.EQ.OIGO TO 278 FII1JWl=SUR ­'EcERDF. sO EO Ywb 278 FI I,JWl=CSTAT 279 CONTINUE 280 CONTINUE C * C * SUR TO EMP C * 283 DO 286 I=IP4,IAAL IFIAlI Il.EQ.EMPIGO TO 285 IFIAlI I I.NE.SURlGO TO 286 IFI 111 I I .EQ.OIGO TO 285 IFIAlI 1 + 1 I.NE.EMP.OR.AllI + IP2).NE.EMPlGO TO 286 IFIAK I-I l.EQ.EMP IGO TO 284 rtIFA2m=jer. IF(A1(1-1 l.EQ.EMP)A2(I-I)=0. IF(A1( I + IP2).EQ.EMP)A3( I )=j0f. 1F(A1( I-I P2>. EQ.EMP >A3( I-1P2)=J0T. 286 CONTINUE DO 287 I cT I P4 , IAAL C * C * FUL TO SUR & SUR TO FUL C * I F ( IK I ). EQ. 8 !GOTO 287 IFIAK I ). NE . SUR .AND .Al ( I ). NE . FUL !GOTO 287 AK I I = FUL IFIAK I +U. EQ.EMP. OR,AK I-I >,EQ,EMP,OR. AK I + IP2), EQ.EMP * .OR.All I-IP2).EQ.EMP )AK I ) = SUR. 287 CONTINUE C * . C * SATISFY D=^r FOR NEW SUR CELLS C * 288 DO 299 I=IP4,IAAL IFlAll I).NE.SUR!GOTO 299 N=# IFIAK 1 + 1 !.EQ.EMP ) N = N+1 IFIAK I + IP2).EQ.EMP ) N-N + 2 IFt All I-I>.EQ.EMP ) N-N + 4 IFIAK I-IP2).EQ.EMP ) N = N + 8 GOTO! 289,29jer, 29 I ,292,289,293,290,295,296,290,289,297,289,298T, 289 IN 289 A2II ) = -A3(I-IP2)>) GOTO 299 290 A3II )=A3!I-I P2 )-DZODR*IA21 I)-A2=A2(I> 294 A3!I )=A3(I-IP2 > GOTO 299 295 A3II-IP2 )=A3(I )+DZ0DR*IA2(I)-A2II-I)) GO TO 299 296 A2I I )=A2II-I ) GOTO 298 • 297 A2II-1)=A2!I) 298 A3II-IP2 )=A3(I> 299 CONTINUE I FI DSP.EQ.0.>G0 TO 361 Q ****aA**a*&****w****************#*#a*****##**a****a#**#&A*** C * C * DELETE SURFACE PARTS IN FULL CELLS C * NPN=I NPT=I KN = 3 K=3 313 IFINPT.GE.NS )G0 TO 317 O=XStKN+!1+2. IP = I P2*I O-I> I=XSIKN1+2.+IP 119 14-May-1982 13:34:35 14-May-19.82 13:31:52 IF(AHl).EQ.BND.AND.AK I + IP2).EQ.FUDGO TO 315 IF(Al(I).EQ.FUL)G0 TO 315 NPN=NPN+! XS(K)=XS(KN) XS< K+l ) = XS*IP2+KX+2 IT=I+IP2 I B= I - IP2 I R= I + I IL=I-I XM=KX+ .5 YM=KY+. 5 IF(X.LE.XM.AND.XN.GT.XM)GO TO 342 IF(X.GT.XM.AND.XN.LE.XM)GO TO 342 332 IF + X IF(AKI).EQ.FUDGO TO 337 I F(A K I R ).EQ.F UL)G0 TO 334 IF(Al(IL).NE.FUL>G0 TO 349 DU=ABSIXI-XMtl)*DR ID = -I GO TO 335 334 DU=ABS(XM-X I + I )*DR ID=I 335 IA=I 335 I F(DU.GE.A18(IA).AND.A18(IA).NE.8.)G0 TO 349 A18 ( IA ) = D U I 2(IAJ = ID GO TO 349 120 I 4-May-1982 13:34:35 1.4-May-I 982 13:31 :52 337 DU=ABSIXM-XI)*DR IFIAH IR).EQ.SURlGO TO 333 IA=IL ID=I GO TO 336 338 IA=IR ID = -I GO TO 336 C HORIZONTAL CENTER 342 YI = I YN-Y )/( XN-X )*< XM-XHY IFIAl I I).EQ.FUL>G0 TO 347 IFIAlI IB).EQ.FUDGO TO 343 IFIAH IT),NE,FUDGO TQ 332 DU=ABSIYM-Y I*I. >*DZ I D= I P 2 GO TO 344 343 DU=ABSIYI-VM+!.>*DZ ID=-IPZ 344 IA=I 345 IFG0 TO 332 A10IIAl=DU I 21 IAJ = ID GO' TO 332 347 DU=ABSIVM-YI)*DZ IFIAlI IT).EQ.SURlGO TO 348 IA=IB ID=TPZ GO TO 345 348 IA=IT ID=-IPZ GO TO 345 349 CONTINUE C ****%*****«=***************************&*********************= C * C * THETA CALCULATION C * C* HYDROSTATIC CONDITION ON FUL CELLS UNLESS CONT OUTFLOW 361 IFITYPE.EQ.0.OR.UR.EQ.0.) GO TO 363 DO 362 I=IP4,IAAL A7I I )=£f. IFIAlI I>.NE.FUL ) GO TO 362 J = II-I)/ I PZ+I A71 I)=GZ*(J-1.5)*DZ + GR*(I-IJ-I)*IP2-1.5) 362 CONTINUE C* IMPOSE FREE SURFACE NORMAL STRESS 363 DO 369 I=IP4,IAAL IFIAlI I).NE.SUR)GOTO 369 IFIAlWI I).EQ.0. JGOTO 369 A10II)=DZ/A10.EQ.I.OR.IZI I>.EQ.-1)A10(I)=A10tI)»DRODZ N=0 IFIAlI1+1).EQ.EMP) N=N+I IFIAK I + IP2 ) .EQ.EMP > N = N + 2 IFIAlII-I>.EQ.EMP> N = N + 4 IFIAlII-IPZ).EQ.EMP) N=N+8 GOTO!365,367,369,365,369,369,369, * 367,369,369,369,369,369,369,369) N .121 14-May-1982 13:34:35 14-May-1982 13 :31:52 365 A7( I ) »COF I A2 ( D-AZU-I >)*A1#( I ) + { I.-Aiai D )*A7(.I + I2< D ) GO TO 369 367 A7( I >=COF2*(A3( D-A3U-IP2 > ItlAiai I ) + t I.-Aiai I ) >“A7< I +121 I >) 369 CONTINUE ' C ************************************************************** C * C * FREE SURFACE TANGENTIAL STRESS C * 370 DO 379 I = I P4,IAAL IFlAll I).NE.SURlGOTO 379 IFlAll 1 + 1 l.EQ.EMP.AND.AlI I + IP3).EQ.EMP.AND.Al I I + IP2>.LT.EMP} * A311+1>=A3(I> IFl AllI + IP2 I .EQ.EMP.AND.AltI + IP3 >.EQ.EMP.AND.Al.LT.EMP) * GO TO 372 GO TO 373 372 A21I + IP2)=A21 I) IFlAll I-I >.EQ.EMP.AND.All I + IP1 ).EQ.EMP)A2< I + IP1 >=>A21 I -1 > 373 IFlAll I-I ).EQ.EMP.AND.AllI + IP1 !.EQ.EMP.AND.Al I I +1P2).LT.EMP) * A31I-I>=A31I ) IFl AllI-IP2).EQ.EMP.AND.Al I I-IP I!.EQ.EMP.AND.AllI + ll.LT.EMP) * A21I-IP2)=A2=A2*IP2+I I F I Al I JP ).EQ.FUL.OR.Al I OP).EQ.SUR!GO TO 381 HI I> = =A2*GZ*DZ/NU GO TO 383 382 SN=I . IFlAZlIA>.LT.0.!SN=-I. RAD=ALFS-4.*A*H(I)/B-4.*SN*ALF*A2tIA! IF1RAD.LT.0.!GO TO 335 A21 I>=-A2 + SN*(ALF + SORT(RAD)> ■ 383 IFlABSlAZl I)).LT.ABSlA21IA!!!GO TO 388 385 AZl I) = -A2*BCLT(J) 389 Vl IPZ,J!=VlIPl,0)*BCRT=0. 418 IF(A1(IJ.EQ.FUL.OR.AKIJ.EQ.SUR) AS(I)=RDR* * (A2-A2,1 = 1,12),T,CYCLE WRITE(10,13) DO 434 I=I,IALL IF (A2.EQ.0..AND.A3(I).EQ.0.)SOTO 434 NROWS=II-I)/IP2 OP=NROWS+! IP = I-N,R0WS*IP2 WRITE!10,31>IP1OP1Al(I),A2(I),A3(I),A6(I>,A7(I ), * A10 ( I >, A 9 ( I >, I 2 ( I > 434 CONTINUE I 2P = 2 IC=I 439 IF=IPl I F ( I F . LE . I 20’* IC +1 )GOTO 440 IF=I20*IC+I IC=IC+! 440 DO 442 0=2,OPl 0P=0P2-0+l DO 441 I = I2P , IF 441 IFI(I) = FI I,OP) 442 WRITE!10,27) I I F 11 I ), I = 12P,IF> WRITE!10,26) IF!IF.EQ.IPl)GOTO 443 I2P = I 20*(IC-I) + 2 GOTO 439 443 CONTINUE WRITE!10,26)NP I F(N.NE.3 >G0 TO 448 C WRITE TO PLOTTER 444 WRITE!15,12 X NAME!I),1 = 1,12),T1CVCLE WRITE!15,13) DO 446 I=I,IALL IF!A2(I>.EQ.0..AND.A3!I>.EQ.0..AND.AS!I ).EQ.0.. * AND.A7!I ) .EQ.0. )G0 TO 446 NROWS=!I-I)/IP2 OP=NROWS+! IP=I-NR0WS*IP2 WRITE!15,14)IP,OP,Al(I),A2(I>,A3(I>,A6(I),A7(I>, * AS!I ),A9NS WRITE! 15,25) (XS( I >,I = 1,2*N'S> 448 GO TO KRET 450 CONTINUE IN D = - 9 9 NROWS=!I-I)/IP2 OP=NROWS+! IP = I-NROWS* I P2 WRITE!15,26)IND WRITE!15,26)I ND RETURN 460 CONTINUE 470 CVCLE=CVCLE+! TM I NT= 1000. O O O O O 124 I 4-May-I 982 13:34.'35 14-May-1982 13:31:52 DO 480 I-IP4,IAAL I F(AI.GE.EMPlGO TO 4 80 ABA2=ABS ) ABAS=ABS(A3(I)) IF(ABA2.LT.I.E-5)G0 TO 477 TT=DR/ABA2 I F(TT.LT.TM I NT)TMI NT=TT 477 IF(ABA3.LT.I.E-51G0 TO 480 TT=DZ/ABA3 IF(TT.LT.TM I NT)TMINT=TT 480 CONTINUE DT=TFAClkTMINT IF(DT.LT.TMINIDT=TMIN IF < DT.GT.TMAX)DT=TMAX T=T+DT IXS=0 ITERl=0 Q * C « * IV « A U' rt * W Vl Vt « it ft * I'l n ft * W * rt IS Cl Ur W A Vi * « tt * * W tt a « ft ft ft Q * « a tt * ft * * U ft tt « * * ft t!r ft C * I C * TILDE VELOCITY COMPUTATION (N-S EONS) I C * I 500 I F< GAM.EQ.0. )G0 TO 521 C * C * BIVISCOUS VELOCITY CALCULATION C * C * CALC A=A6=DU/DY+DV/DX DO 501 I=I,IAAL 501 A6 = (A2 )fcRDZ+.GE.EMP>G0 TO 519 IFIAH I > .EQ.BND >G0 TO 519 IL=I-I I R= I + I IB=I-IPZ IT= I + IP2 NROWS=IL/IP2 UP=NROWS+! IP = I -NROVZSfc I P 2 CALCULATE STRESSES AIX=AGI I) GN I =2 . * I A2I I >-A2( ID IfcRDR ANGIS=GNIfcGNI ASIS = I(A6II>+A6(IB ) >fc.5 )fcfc2 IND=I C NORMAL STRESS - SI 502 RSS=GNIfcfc2+ASIS I F I RSS.GT.GRAD02)G0 TO 503 SI=NUPfcGNI GO TO 504 503 RS=SQRTIRSSl SI = IGWMNON/RS+ NUIfcGNI C SHEAR STRESS - VIT 504 RVS=AIXfcfc2+ANGIS 125 I F ( RVS . GT. GRADJ02 >60 TO 505 VIT=NUP*AIX GO TO 506 505 RV=SQRTt RVS > VIT=(GWMNON/RV+NU)*A1X 506 GO TOt 50.7,508,515) IND 507 IFt IP.EQ.IPl >G0 TO 510 IFtAK IR).EQ.EMP >G0 TO 509 . C SHEAR STRESS - VIB & NORMAL STRESS - SIR I ND = 2 S = SI VI=VIT AIX=A6-A2(I ) )*RDR GO TO 502 508 SIR=SI SI = S VIB=VIT VIT=VI C HORIZONTAL EQUATION ZIP=A2(I )*(A2(IL >~A2tIR)) UVBR=0.25*(A2) UVTR=0.25*(A2(I )+A2=A2-A6(IB ) >*RDZ C C SHEAR STRESS - VITL & NORMAL STRESS - SIT 510 I F t OP.EQ.OP I )G0 TO 519 IFtAlt IT).EQ.EMP )G0 TO 511 I ND = 3 S = SI VI=VIT AIX=A6(IL > GNI=2.*tA3(IT)-A3(I)>*RDZ ASIS = < >*.5 >**2 GO TO 502 511 IFtAK IB).NE.EMPlGO TO 519 A5(IB)=A3(IB)+DT*GZ GO TO 519 515 SIT=SI SI=S VITL=VIT VIT=VI C VERTICAL EQUATION ZIP=A3+'DT*tRDZ*-A7t IT) >+GZ * +RDR *(UVTL-UVTR)+tVIT-VITL)#RDR * +(S I+S IT)rtRDZ > 519 CONTINUE DO 520 0 = 2,OPl ' 520 UTIL(IPllO) = UTILt I P I,O ) *ROS(O >+ UTI L( IPl-I,0 )#( I.-ROS(O)) 14-May-1982 .13:34:35 14-May-1982 13:31:52 14-May-1982 1.3:34:35 I 4-May-1 982 13:31 :52 GOTO 532 521 DO 526 I-IP4,IAAL IFlAH I >. GE . FMP )GO TO 526 IFlAH I >.EQ.BNDlGO TO 526 NROWS=II-I1/IP2 OP=NROWS+! I P = I-NROWS1''I P2 I R= I + I IL=I-I IT=I+IP2 I B = I -1P2 : ID=IB-I KDmKVAk7AKVmA7hAT 8 Kh KAskAk4Ves7 EO zYf ­ 5OhKIOuEcG ua' ZIP=AZl I>*=A2(I> + DT*(RDR*(ZIP+A7(I )-A7(IR)) * +RDZ*IUVBR-UVTR)+GR +NU*IRDZ2*(A2(IT )+A2-2.*A2 * -RDRDZ*!A31 I R)-A31 I- IP I)-A31I>+A3!IB)))) C VERTICAL N-S 524 IFlJP.EQ.JPI.0R.A1G0 TO 526 ZIP = A31 I)*(A3(!B)-A31IT)) UVTL=Bf. 25* I A21 I L )+A21 I + IP I ) ) *( A31 IL )+A31 I)) UVTR=JBf. 25*1 A2< I )+A2( IT ) )*( A31 IR )+A3( I ) ) A51I)=A3(I) + DT*(RDZ*(ZIP+A71I )-A7(IT) )+GZ * + RDR*l I UVTL-UVTR)-NU*I IRDZ*(A21IT)-A2(I)) *-RDR*(A3(IR)-A31 I ) ))-(RDZ*(A2(I + I P I )-A21I L) ) *-R D R * I A 31 I >-A3< IL)))))) 526 CONTINUE Q A************************************************************* c * C * CALCULATE DIVERGENCE FOR FUL CELLS C * 532 ALPO=ALP 533 DO 540 I = I P 4,IAAL IFtAll I).NE.FULlGOTO 535 ABl I ) = RDR*( A41 I >-A4( I-I ) ) + RDZ*(A5( I >-A5( I-I.P2 > ) CONTINUE 535 A61 I)=0. 540 CONTINUE 546 CONTINUE (2 *****Vl******G0 TO 554 DO 552 I = IP4,IAAL IFlAll I>.NE:SUR)GO TO 552 AFUL = 4 .*< I .+ALPOHl 4.-11 .+ALPO>*< I .-A101 I ) > >-l. IFIAFUL . L/T\ A71 I + I21 I > > )A7< I + I21 I ) I=AFUL 552 CONTINUE 554 DO 555 1=2,IPl .127 14-May-1982 13:34:35 I 4-May-1982 13:31 :52 A6) 562 CONTINUE 563 IF( IND.EQ.0 >GOTO600 I ND=0 I F(ITER.GT.500 >GOTO 575 ITER=ITER+! DO 570 I = I P4,IAAL IF(Al(I).NE.FUDGOTO 570 1PSIR=A6(1+1) PS IB=A5(I-1 P 2 > PSIL=AG(I-I) IF(TYPE.EQ.0)GOTO 565 IF(Al(I+I).EQ.OB) PSIR=AG(I) IF(.A1( I-IP2 ) .EO.OB ) PSIB=AG(I) IF(A1(I-I ).EQ.0B)PSIL=A6(I) 565 X=W*(I.+A7*<-A8< I ) + RDR2*(PSIR+PSIL) + RDZ2*(PSIB + IAG(I+IP2)))-A7(I)*A6(I) I F(X.GT.I.E + 10)G0 TO 575 Y=ABSiX >-ABS(A6+ABS(A6(I)) AG(I)=X I F < Z.EQ.0 >GOTO 570 I F(ABS(Y/Z).GT.EPS)I ND=I 570 CONTINUE GO TO 554 575 ALPO=0.5*ALP0 WRITEt12,29)ITER WRITE(10,29>ITER IF(ALPO.GE.0.I)G0 TO 533 580 ITERl = ITER IXS=IXS+! IF( I XS . GT.2 )G0 TO 74 DT=DT*.5 T=T-DT GO TO 500 600 CONTINUE £ , * j l < v l * i t r t f t f t * * f t * * f t f t A * f t * f t f t f t f t f t f t f t f t * f t * * * f t ' f t f t f t f t f t f t f t f t f t f t f t f t f t f t f t f t f t f t f t f t f t f t f t f t f t f t f t f t f t C * C * CALCULATE COMPLETE VELOCITIES USING POTENTIAL C * FUNCTION AND BOUNDARY CONDITIONS C * DO 620 I = I P 4,IAAL A9=A7(I )+A6/DT IF(Al(I).GE.EMPFGOTO 620 1F(A1< I ).EQ."BND)GOTO620 IF(A1( 1 + 1 KGE.EMP >'GOTO610 1 28 ' 14-May-1982 13:34:35 I 4-May- I 982- 13:31:52 A4(I)=A4(I)-RDR*(A6(1 + 1>-A6(I> ) 610 IFtAK I + IPZl.EQ.EMP.OR.AK I + IP2 > .EQ.BND )GOTO 620 AS(I)=ASd)-RDZ*(A6( I + I P2 )-A6(I)) 620 CONTINUE- C * SET VEL OM SUR CELL EMP FACES (D=0> 650 DO 661 I = I P4-, IAAL I F( Al ( I ). NE .'SUR )GOTO 661 N=0 IF(Al(1 + 1>.EQ.EMP) N = N+I I F(AI(I + IP2).EQ.EMP) N = N + 2 IFtAKI-I J.EQ.EMP) N = N + 4 IF(AK I-IP2 ) .EQ.EMP ) N = N + 8 GOTOf 651.652„'653.654,65 I,655.652.657.658.652.651.659.651 ,652.651 )H 651 A4(I) = (A4(I-I>-DR0DZ*-A5(I-IP2 > >) GOTO 661 652 A5(I>=A5(I-IP2)-DZODRtt(A4(I>-A4(I-I)> GOTO 661 653 A4(I>=A4) . GO TO 661 655 A4(I-I )=A4(I> 656 A5(I>»A5 GOTO 661 657 A5( I-IP2 >»A5( I l + DZODR^ t A4( I >-‘A4( I-I >) GO TO 661 658 A4.NE.SUR)GOTO 664 IFtAKI + l).EQ.EMP. AND.AKI + IP3 ).EQ.EMP.AND.Alt I +1P2 ). LT. EMP ) * A5(I + I)=A5( I> IFtAltI + IP2).EQ.EMP.AND.AltI + IP3).EQ.EMP.AND.Alt 1 + 1J.LT.EMP) * GO TO 662 GO TO 663 662 A4=A4(I ) IFt Al(I-I).EQ.EMP.AND.Alt I-IPl).EQ.EMP >A4(I + IP1)=A4(I-I ) 663 IFtAK I-I).EQ.EMP.AND.AltI + IPI).EQ.EMP.AND.Alt I +1P2).LT.EMP) * AS(I-I)=AS(I) IFt AltI-IP2).EQ.EMP.AND.Alt I-IPl).EQ.EMP.AND.Al(I + 1).LT.EMP) * A4t I-IP2 )=A4.t I ) 664 CONTINUE Q ***AA********A**MK*******a****aa******************#*******A*** C * C * TANGENTIAL BOUNDARY CONDITION C * DO 675 1=2,IBAR ' IA= I + IP2 IT=IALL-I IB=IT-IPZ A4fIT)=A4(IB )*BCT IF(KLTdFS)GO TO 674 IFtASt I ).EQ.-100.)G0 TO 674 129 IF( BCB.EQ.J0f. )G0 TO 670 A4(I )=A4(IAlwBCB . GO TO 675 670 IF(B.NE.0.I GO TO 671 A4(I>=A4(IA)+A*H(I>*GZ»DZ/NU GO TO 672 671 SN=I. I F(A4(IA).LT.0. ISN = -I. RAD=ALFS-4.ttAwFK I>/B-4.*SN*ALF“A4(IAl IF=A4II + IP2 IttBCB' IFI Al I I-I I.NE.OB.AND.Al II + IP2I.EQ.OB I A5II I=ASII-I IttBCB IFI Al I 1 + 1 I.NE.OB.AND.Al II + IP2l.EQ.OB I A5II>=A5 + UK*DT*RDR XP(KN+I )=XP(K+l)+VK*DT*RD2 O = XP ( KN+ I ).+ 2 . I=XP(KN)+2.+IP2*(0-l) IF( XP(KN). GE. IBAR.OR.XP( KN >.LE.£r. )GOTO 730 IF« XP(KN+I >.CE.JBAR.OR.XP< KN+I ) .LE.0. )COTO 7 30 715 KN=KN+2 NPN=NPN+! 720 IF(Al(I).NE.EMP!GOTO 730 Al(I)=SUR IF (Al ( I + I >. EQ. EMP > A2( D = U(IDfOD) IF(A1( I-I > .EQ.EMP ) A2( I-D = U(ID-IfOD) ■ IF(Al(I+IP2).EQ.EMP) AS(I)=V(IDfOD) IF(A1(I-IP2).EQ.EMP) A3(I-1P2)»V(I D,OD-D 730 K=K+2 NPT=NPT+! GOTO 710 735 NP=NPN C MOVE SURFACE MARKERS 741 NPT=I NPN=I KN = 3 K = 3 742 IF(NPT.GE.NS)GOTO 750 ID=XSlK>+2. HPX=ID-D-XS(K) HMX=D-HPX OR = XSt K+l )+ D 5 HPY = OR-0.5-XS< K+l> HMY=I.-HPY UK=0. IF(NPT-EQ-^ )GO TO 743 UK=HPX*HMY*U(ID-1,0R+1)+HMX*HMY * *U(IDfOR+!) + HPX*HPY*U(ID-I,OR> + HMX*HPY*U( IDfOR) 743 IR=XS(K)+!.5 HPX=IR-0.S-XS(K) HMX=I .-HPX OD = XS(K+l ) + 2. HPV = OD-I.-XS(K+l ) HMY=I.-HPY VK = HPX*HMY*V+HPX*HPY * *V(IRfOD-I)+HMX*HPV*V+VK^ DT*RDZ ' 0 = XS(KN+I ) + 2. I=XS(KN)+2.+IP2"(O-D 744 KN=KN+2 NPN=NPN+! ■ . 746 IF(A1( D.NE.EMP )GOTO 748 A K D = SUR 131 14-May-1 982 13:34:35 14-May-1982 13:31:52 I F(Al<1 + 1 ).E0.EMP ) A2(I>-U(ID1JD) I F ( Al < I-I ). EQ. EMP > A2( I-D = IM ID-1, JD ) I F’( Al < I + IP2>,EQ.EMP > A3( I > = V( ID,JD> IF< AK I-IP2 ) .EQ.EMP ) A3< I - I P 2 )=V( I D , JD-1) 748 K=K+2 NPT=NPT-H GOTO 742 750 NS=NPN Q ********************************* *'» C * C * PARTICLE INFLOW C * IF ( T .GT .TDT HFS=I IF ( TYP E.EQ.0 I GOTO 250 I F(L 2.EQ,0)GO TO 250 ULL=UL-A*iS( 2 ) + BP*DT*RDZ IF(XS(2).LT.L1>XS(2)=L1-.001 IF(XS(2).LT.0. )I F S= I IFXS(2)=L2-.001 L2P = XS(2 ) DO 757 J = 2,JP2 U( I,J >=0. I F(A.EQ.-99 . )ULL = U(2,J ) IF( J.GE.L1+2.AND.J.LE.L2P + 2 ) U(I,J I = ULL 757 CONTINUE I F(NP.EQ.0)£0 TO 250 K N = 2 * N P + I NIN = NYB*(XS(2)-L1 ) 760 X = ULL*RDR*< T + DT)-XDIS*< 2.rtCOLS-H . ) I F(X.LT.0. >GOTO 250 IF(NIN.EQ.01GOTO 250 COLS=COLS+I. Y=YFIR NP N=NP + N I N 770 XP(KN)=X X P ( K N +1 ) = Y KN=KN+2 J = Y-H . I=X + 2.+JMP2 NP=NP+! Y=Y+YDIS IF(A1(I>.NE.EMP)GOTO 775 Al(I)=SUR A2( I DULL 775 IF(KN.GT.LPB)GOTO 75. IF(NP.LT.NPN)GOTO 770 GOTO 760 C «*** A RAB R Ir RBBArtRBBytyiAVrArtABaBBftftBBByrBBBBttAftBftyt **«*»**«* Vt = C * C * 999 RETURN END APPENDIX B TIME SEQUENCE PARTICLE PLOTS OF BVSMAC CALCULATION OF FLOWING SNOW TEST 1. TqZP - 2.2 V = .00.2 m^/s, V 1.= 0.10 m^/s 2. TqZP = 3.0 rtZYs^/ V = .001 m^Zs/ V ' == 0.20 m^Z^ 3. T ZP ” 1.0 m^Zs2, V = .010 Hi2Zsr V 1 = 0.20 m2ZsO B I N G H A M C O E F f I C I E N f G A M = 2 . 2 0 0 V I 5 0 0 . S I ~r \/ N U = 0 . 0 0 2 NUP = 0.100, F R I C T I O N . B O B = - 1 . 0 0 0 G R I D- I B A B = U l O J B A R = I O ..DR = O . .200 ' 0 2 = 0 . 0 5 0 E S f g CL- 0 7 UJ - ' Q JU 3 . 5 0 7 . 0 0 T = 0 . 2 0 1 1 0 . 5 0 1 4 . 0 0 P O S I T I O N (M) 1 7 . 5 0 21 . 0 0 TT-^ n=T 2 4 . 5 0 2= S F § ^ 0 O1. OO Q ^Vv-.rrr/tolvM^r^ 3 . 5 0 7 . 0 0 T = 0 . 1 0 1 T f Tn-TV1.-! \aTrV 1 ■ i 'I H 1 0 . 5 0 1 4 . 0 0 1 7 . 5 0 P O S I T I O N (M) I Ijiii I Tl I I I . . I I i i ; |. i-,— -,— , 2 1 . 0 0 2 4 . 5 0 z: Oin o F S LU Q "0.00 3 . 5 0 7 . 0 0 T= 0 . 0 0 0 1 0 . 5 0 1 4 . 0 0 1 7 . 5 0 2 1 . 0 0 2 4 . 5 0 P O S I T I O N (M) 134 0.00 1 7 . 5 0 21 . 0 0 2 4 . 5 0Lu Q 3 . 5 0 7 . 0 0 T = 0 . 5 0 0 1 0 . 5 0 1 4 . 0 0 P O S I T I O N (M) .=n S ^ j0 O-OO CD I-I 'i t I i i|i I I I1- i 3 . 5 0 T = 7 . 0 0 0 . 4 0 0 1 0 . 5 0 1 4 . 0 0 1 7 . 5 0 P O S I T I O N (M) 21 . 0 0 H_CO' f Ul 2 4 . 5 0 T - O , O C:°o.oo CD I I I i | l H 1. ' T i 1I 1 T i 3 . 5 0 7 . 0 0 1 0 . 5 0 T = 0 . 3 0 1 1 4 . 0 0 1 7 . 5 0 P O S I T I O N (M) 21 . 00 2 4 . 5 0 O3 . 5 0 T = 0 . 8 0 0 21 . 00 rrrrr^-rrr',- 2 4 . 5 0 P O S I T I O N (M) s:S O ,=c S ^ c IO-OO Q T T I Ji I I I I i i i i i T i I i I I I I I i I i I 3 . 5 0 7 . 0 0 T = 0 . 7 0 0 ...... I I I I I I T I I T I T l I I I l I | ' I I I I I i I I I . r 1 0 . 5 0 1 4 . 0 0 P O S I T I O N (M) 21.00 2 4 . ' 5 0 ^ in E P Itj0 O - O O Q rTJTTTTTT 3 . 5 0 TTT^ TTTTTTT^ TTi 7 . 0 0 T = 0 . 6 0 1 1 4 . 0 0 P O S I T I O N (M) .1111 *! i' 11i" I EpKPaQ E 2 1 . 0 0 2 4 . 5 0 136 DE PT H (M ) DE PT H (M ) DE PT H (M o0 . O O 0. 50 r-P ,*0 0 0. 50 Ip. 00 0 .5 0 I r". fi'I ). 00 s'. 5 0 TTTTTTT 7 . 0 0 I = 1 . 101 Trv.1 i rrri 'rrn i R~oRo iRoR 1 0 . 5 0 1 4 . 0 0 1 7 . 5 0 P O S I T I O N (M) 2 1 . 0 0 2 4 . 5 0 . 00 3 . 5 0 T = 7 . 0 0 1 . 0 0 0 1 0 . 5 0 1 4 . 0 0 1 7 . 5 0 2 1 . 0 0 P O S I T I O N (M) TTTTTj-TT-H*" 2 4 . 5 0 1.00 3 . 5 0 T = 0 . 9 0 2 T"1? m ^1Vj I Hl T n 1T1F jf^ f f r ' i' T rM 7 . 0 0 1 0 . 5 0 1 4 . 0 0 1 7 . 5 0 P O S I T I O N (M) 21.00 2 4 . 5 0 137 DE PT H (M ) DE PT H (M ) DE PT H (M: n. O O 0. 50 n. 0 0 0. 50 f3 .0 0 0. 50 i f"'Vi 3 . 0 0 3 . 5 0 -M I'i w w r,T tr,' I T r i1liHiliI ..... .. in,.,. T— 7 . 0 0 1 0 . 5 0 1 4 . 0 0 1 7 . 5 0 2 1 . 0 0 2 4 . 5 0 T = 1 . 4 0 2 P O S I T I O N (M) 3 . 0 0 3 . 5 0 T = 1 . 3 0 1 r"T^ rrT*mT*T"rfr?n 7 . 0 0 'TTfrrmTTT I I i i Ti I I I j 1 0 . 5 0 1 4 . 0 0 P O S I T I O N (M) VPl] Frft !"TfrTrrT -rT^ ri ~r I - rPr,T-i^ 111 .| 1 7 . 5 0 21 . 00 2 4 . 5 0 T'1," i' 'i 'i1 3 . 0 0 I H l l ,'I i p T i', f'r7f T T ? - , , I , I I ,|. , ,'. ,-- 3 . 5 0 7 . 0 0 T = 1 . 2 0 2 1 0 . 5 0 1 4 . 0 0 1 7 . 5 0 P O S I T I O N (M) 21.00 2 4 . 5 0 138 3=8 £ g . . . . . , ii 0.00 3.50 Q 7.00 T = 1 . 7 0 3 10.50 14.00 P O S I T I O N (M) 2 1 . 0 0 24.50 ± 8 T = 1 . 6 0 2 1 0 . 5 0 1 4 . 0 0 1 7 . 5 0 P O S I T I O N (M) 2=8 .zzS - cT T o T Q 3 . 5 0 •TTrfrri^ Trm- 7 . 0 0 T= 1 . 5 0 0 !'ITni! !(,' IVfVf i'fiTi-mn-rii-r^ i-n?3 1 0 . 5 0 1 4 . 0 0 1 7 . 5 0 P O S I T I O N (M) 21 . 0 0 139 O/.g E'-oTuo O 3 .50 i T M M irf 7.00 T = 2 . O O U I rfi'i! I1 11 1 rTTTi . ... . 10.50 14.00 17.50 21.00 24.50 P O S I T I O N (M) J oO1- T r Cl 3.50 T = 1 . 9 0 4 T_i i ^ n i'l r Iinr i i ; i 7.00 10.50 14.00 P O S I T I O N (M) 17.50 21.00 24.50 ’5 ” O Eb^b1. OO CD 3.50 7.00 10.50 J = 1 . 8 0 3 P O S I T I O N (M) 14.00 17.50 21.00 TIliTTTfi 24.50 Om 3 = g I i 0 O T o o 1 3.50 ■ 7 f®l' 11 "i'*"11 ^ r i'v 7.00 T = 2 . 3 0 H r PTjHTT j /NaNN rFrr r- riT~i"n Vfj1^-VrrrM 10.50 14.00 P O S I T I O N (M) 1 7 . 5 0 I ■ ' ' I ' ' > I I * I ‘ ' I 1 ■ I I ! . * 21.00 24.50 : S S "0.00 T t'7 TT Pl 'i *i "i™I ■ T1^rI tyrTTOIFrFTf^ 3.50 7.00 T = 2 . 2 0 H ; - I,,, I, Vi i I' T F 10.50 1 4 . 0 0 1 7 . 5 0 2 1 . 0 0 24.50 P O S I T I O N (M) ES ,^=S E p' E g g Cl 3 . 5 0 T= 2 . 1 0 4 ri I1T BQmEKn 7.00 *. i'I1I1H 1I,. , l'i !' 1 0 . 5 0 1 4 . 0 0 P O S I T I O N (M) 1 7 . 5 0 p21 OO 24.50 141 . 0 0 I „ 5 0 2. 0 0 IjME (SEC) L E A D I N G E D G E P O S I T I O N (M) n.oo 3.00 6.00 S.00 12.00 15.00 18.00 21.00 . — -------- 1------------ 1------------ :------------ 1------------ 1____________ I____________ ! 142 I IML I SEC) L E A D I N G E D G E V E L O C I T Y (M/S) .0. OO 3 . 0 0 6.00 9. G O 12.00 1 5 . 0 0 1 8 . 0 0 143 B N G H R H C O E F F I C I E N T GfiM = 3 . 0 0 0 V I S C O S I T Y M D = O - O U l NUP = 0.20 F R I C T I O N P C B = - 1 . 0 0 0 f i = 2 0 . O O O B = - O . 3 0 3 G R I D I B R R = I 4 0 J B f i B = I O D R = O . 2 0 0 Q Z = O - O S O 144 OLO Of. bo 20.00 28.0012.00 IS.C P O S I T I O N (M) 8.00 T= 0 . 2 0 1 O in 24.00 28.0020.0012.00 16.00 P O S I T I O N (M) 8. CO4.00 0 . 1 0 0 o in 28.0016.00 20.0012.00 P O S I T I O N (M) 4.00 8.00 T= 0 . 0 0 0 145 O in 12.00 16.00 P O S I T I O N (M) 28.008.000.00 4.00 T= 0 . 5 0 1 O I/) f frXV.X'V; .. < I j'-iVn i' i".TT T 1-rrrr; fn i . ..., ■ ■ ■ * 8.00 12.00 16.( P O S I T I O N (M) 20.00 24.00 28.004.00 0 . 4 0 1 O in 20.00 24.00 28.0012.00 16.C P O S I T I O N (M) 4.00 8.00 0 . 3 0 1 146 O LD 24.004.00 8.00 12.00 16.00 P O S I T I O N (M) 20.00 28.00 0 . 8 0 2 O un 4.00 8.00 12.00 16.00 P O S I T I O N (M) 20.00 24.00 T= 0 . 7 0 0 o LD T= 0 . 6 0 0 P O S I T I O N (M) 20.00 28.00 147 28.00 o in U C JmaOO 4.00 T = i -yffrTn 8.00 12.00 1 . 101 P O S I T I O N rrp vn'fHVfTIeITrilTi q ; HNN Q Q fH 16.00 20.00 (M) 24.00 o in L U .. 8.00 12 . 00 16.00 P O S I T I O N (M) 20.00 28. 00 o in 8.00 12.00 16.C P O S I T I O N (M) 20.00 28.0024.00 0 . 9 0 1 148 Oi n YT rf f : V =TZYrnT '0.00 8.00 12.00 16.C P O S I T I O N (M) 20.00 24.00 28.00 o LO LU . 0.00 4.00 8.00 12.00 16.C P O S I T I O N (M) 20.00 o m 24.00 149 Oi n LU . 14. OO 8.00 12.00 16.C P O S I T I O N (M) 20.00 24.00 28.00 T= 1.701 O in 4.00 8.00 12.00 16.C P O S I T I O N (M) • 20.00 T= 1.601 0 m ^rfforfrfrrrrrrTrtn m-THl 12.00 16.00 P O S I T I O N (M) 4.00 20.00 28. 00 T= 1.501 150 O in 4.00 8.00 12.00 16.00 P O S I T I O N (M) 20.00 28.00 2.000 O in L U ; T-Ti VrrfTTf^VffTrYi*! 8.00 12.00 IS.C P O S I T I O N (M) 20.00 24.00 28.00 1 . 900 O in 'tr,-) ff TifiVrrlyM * 1 iVfiT,Vf-p , ft 4.00 20.00 28.00 P O S I T I O N (M)1 . 807 151 T I M E L E A D I N G E D G E P O S I T I O N (M) 3 9.00 12.00B. OO 3 . 0 0 1 5 . 0 0 18.00 21. OO I i -I 152 .00 2.50 L E A D I N G E D G E V E L O C I T Y (M/S) 15.00 18.0012.00 21.009. CO5.00 6.00.0.00 I 153 B I NGHFIH C O E F F I C I E N T G A M = 1 . 0 0 0 V I S C O S I T Y N O = 0 . 0 1 0 NUP = 0.10 F R I C T I O N B O B = - 1 . 0 0 0 0 = 0 . 1 0 0 6 = 0 . 0 1 0 G R I D I B A R = I V O J B A R = I O D R = O . 2 0 0 D Z = O - O S O 154 O in 8.00 12 .00I 16.00 2 P O S I T I O N (M) 26.00 32.00 0 . 2 0 1 o in E g # . I... 111 7 E paOO fA77 WA77 mYA77 mLA77 0 . 102 P O S I T I O N 20.00 (M) 24.00 28.00 32.00 o in ^ TPrFT-T, 0 Tl-OO 4.00 8.00 12.00 T= 0 . 0 0 0 16.00 20.00 P O S I T I O N (M) 24.00 28.00 32.00 155 O Ul LU . mill 111 m i n I nr 8.00 20.00 24.0012.00 15.00 28.00 32.00 T= 0 . 5 0 2 P O S I T I O N (M) O LD 4f " i" [ViIl i-Trrn n7^714%) rT 4.00 8.00 12 . 00 16.00 20.00 24.00 28.00 32.00 0 S I T I 0 NT= 0 . 4 0 0 O LD L U • Tfr'n T r " TnTTT^TrTTr 11 m 8.00 12 . 00 16.00 20.00 24.00 28.00 32.00 T= 0 . 3 0 0 P O S I T I O N (M) 156 O LO 20.0012 . 00 24.00 28.00 32.00 P O S I T I O N (M)0 . 8 0 0 O in L U . 4.00 8.00 12.00 • 16.00 2l P O S I T I O N (M) 24.00 28.00 32.00 T= 0 . 7 0 1 O LO 4.00 8.00 I 16.00 2 P O S I T I O N (M) 12 . 00 32.00 0 . 6 0 2 157 O Lh 20.CO 24.00I 16.00 P O S I T I O N (M) 32.0012 . 008.00 oin I 16.00 2 P O S I T I O N (M) 24.0012 . 00 28.008.00 32.00"0.00 I . O O O o Lh 24.00 28.00I 16.00 2 P O S I T I O N (M) 12.00 32.004.00 8.00 0. 9 0 2 158 Oi n oC b'.bo 4.00 8.00 12.00 1 "l 600.... 2 0 ."bo.... 24.00.... 28.00 32.00 T= 1 . 4 0 0 P O S I T I O N (M) O LD 4.00 8.00 12.00I 16.00 2 P O S I T I O N (M) 28.00 32.00 o in 4.00 12 . 00I 16.00 2 P O S I T I O N (M) 32.00 1 . 2 0 0 159 O in 0.00 4.00 8.00 12.00I 16.00 2 P O S I T I O N (M) 24.00 28.00 32.00 I . 707 o in 4.00 8.00 12.00I 16.00 2 P O S I T I O N (M) 24.00 28.00 T= 1 . 6 0 9 O LO T= 1 . 5 0 6 P O S I T I O N (M) 160 O LO I-U ■ 8.004.00 12.00 16. 00 20.00 32.00 T = 2 . 0 0 7 P O S I T I O N (M) O in 8.00 12. OOI 16.00 2 P O S I T I O N (M) 24.00 28.00 32.00 o in 4.00 8.00 I 16.00 2 P O S I T I O N (M) 24.00 28.00 32.00 161 .00 0.50 1.00 1.50 2.00 2.50 TIME (SEC) L E A D I N G E D G E P O S I T I O N (M) .0.00 6.00 9 . 0 0 12 . 00 1 5 . 0 03 . 0 0 1 8 . 0 0 21.00 162 OG ‘ L E A D I N G E D G E V E L O C I T Y (M/S) .0.00 3 . 0 0 6.00 9 . 0 0 12. OO 1 5 . 0 0 1 8 . 0 0 21 . 0 0 163 APPENDIX C VISCOUS DEFORMATION OF A RECTANGULAR BLOCK. 2 1 . V = .001 m /§ 2. V = . 010 m^/s 3. V = .100 m2/s 4. V = .200 m2/s 165 B I N G H A M C O E F F I C I E N T . G A M = 0 . 0 0 0 V I S C O S I T Y H 3 g R 0.001 F R I C T I O N B O B = - 1 . 0 0 0 0 = 2 0 . 0 0 0 G R I D I B A R = I B ' J B A B = I O D R = O . 2 0 0 D Z = O . 0 5 0 166 O in 5 ° Is O C b .00 ' 0.50 1.00 1.50 2.00 " ' ' ! ' ' 1 ■ 2.50 r * " 3 T = 0 . 2 0 O S I T I O N (M) O LO I'?*'": V f : V f V - eY eI y : : / I I I I j \ : \ y y S . i I i I I I I } J Ax-V' j ‘-T* ^ IcjcT). 00 0.50 I . 00 1.50 2.00 T = 0 . 1 0 2 P O S I T I Q N (M) -I— T- I • 2.50 3.00 o m Q C b . 0 0 1 . 0 0 1 . 5 0 2 . 0 0 T = 0 . 0 0 O ^ 0 S I T I 0 N (M) 2 . 5 0 3 . 0 0 167 oin Q C b 0.50 1.00 1.50 T = 0 . S O C f O S I T I O N (M) 2.00 3.00 O LO 0.50 1.00 1.50 T = O . U O C f O S I T I O N (M) 2.00 2.50 Oin Q C b.oo .--V' i— r 0 . 5 0 1 . 0 0 I . 5 0 2 . 0 0 T = 0 . 3 0 0 P 0 S I T I 0 N (M) T • > • 2.50 3.00 ■168 B I N G H A M C O E F F I C I E N T G A M = 0 . 0 0 0 V I S C O S I T Y N O = 0 . 0 1 0 F R I C T I O N B O B = - l . O O O A = 2 0 . O O O ■ G R I D ■ I B A B = I S J B A R = I O D R = O . 2 0 0 ■ D Z = O . 0 5 0 1 6 9 o in L U a c b / ) . r --- - ... — ?.... i.... ' " I " ' 0.50 l'. OO I I . T I T 1.50 2.00 r T * I ' f r ■ ■ 2.50 3. T= O . 200^06 I T I ON (M) O LO r-fTiVilV: H Tl : I : \J : : : : : :. :: : X V - ........ .— — -r-- - — - --- - - o c b.00 0.50 I . OO 1.50 2.00 T= 0.1 0 2P 05 1 T I ON (M) 2.50 3.00 Q C b . o o 0 . 5 0 1 . 0 0 1 . 5 0 2 . 0 0 T = O . O O Q P O S I T I O N (M) 2 . 5 0 3 . 0 0 170 o LO Q C b.oo 0.50 1.00 1.50 T= 0 . 5 0 ( f O SITION (M) 2.50 3.00 Q C b.oo 0.50 1.00 1.50 2.00 T= 0 . 4 0 Op GS I T I O N (M) 3.00 0 . 5 0 1 . 0 0 1 . 5 0 T = 0 . 3 0 0 P G S I T I 0 N 3.00 171 o in -v.'.v L U 1.00 1.50 0 . 7 0 2 P 0 S I T I 0 N (M) 2.000.50 2.50 o c b. oo 0 . 5 0 1 . 0 0 1 . 5 0 T = 0 . B O S f 3 O S I T I O N ( M ) 2.00 2.50 3.00 172 B I N G H A M C O E F F I C I E N T G A M = 0 . 0 0 0 V I S C O S I T Y ' N O = 0 . 1 0 0 F R I C T I O N B O B = - 1 . 0 0 0 . 0 = 2 0 . 0 0 0 G R I D I B A B = I S J B A R = I O , D R = O . 2 0 0 D Z = O - O S O \ 173 Q in o c b.oo 0 . 50 1. 00 I . 50 T = 0 . 2 0 5 P 0 S I T I G N 2.00 3 . 00 Q C b.oo 0 . 50 1 . 00 1. 50 2 . 00 T= 0.1 OOp OS I T I ON (M) 3 . 00 o in (M ) O X ■• • •• I ««".* ■ •• • •• • •• • •• I— P OOLU ' — r • r- I— -1— ~-|— . ' I —T----'-- • -T • T • t -• - T ■ • •- v• ■ — r T- r T 1 -»,^--V--S— r-^ ' ' ''' ' T " "I r—T'»— + ir 22 f r . f (r ff (r .f 3r f f 3 r . f - r f f T = 0 . O O O p O S I T I O N (M) 174 oi n Q C b'.oo 1. 00 1. 50 5 0 5 P 0 S I T I 0 N 3.00 Q C b.oo 0 . 50 1. 00 1. 50 T= 0 . 4 0 5 P 0 S I T I 0 N 2.00 3.00 Oin oc:b.00 0 . 5 0 1 . 0 0 1 . 5 0 T = 0 . 3 0 5 ^ 0 5 1 T I O N ( M ) 2.00 2.50 3.00 175 oi n Q C b.oo 0 . 50 1. 00 1. 50 T = 0 . S O B p G S I T i a N 2 . 50 3 . 00 O in U C K. g g 0 . 50 1 . 00 1. 50 2 . 00 T = 0 . T O B p G S I T I O N (M) 2 . 50 3 . 00 o in QCb.OO 0 . 50 1. 00 1 . 50 2 . 00 T = 0 . O O B p G S I T I O N (M) "T"'' 2 . 50 3. 00 176 oi n O C ^.00 0 . 50 1. 00 1. 50 2 . 00 T = I . 1 0 5 P 0 S I T I 0 N (M) 3 . 00 Q C b.oo 0 . 50 1. 00 1. 50 2 . 00 T= 1 . 0 0 5 P 0 S I T I 0 N (M) 3 . 00 a c b.oo 0 . 5 0 1 . 0 0 1 . 5 0 T = O . 9 O 5 P 0 S I T I G N 2 . 50 B I N G H A M C O E F F I C I E N T ' G R M = 0 . 0 0 0 ' V I S C O S I T Y N U = 0 . 2 0 0 ' F B I C T I O N B O B = - 1 . 0 0 0 A = 2 0 . O O O G R I D . I B A B = I B J B A R = I O D R = O . 2 0 0 ' D Z = O - O B O 178 oi n 0 . 50 1.00 1.50 2.00 2 . 50 3 . 00 T= 0 . 20[f OS ITION (H) g C K2aOO 0 . 50 1 . 00 1. 50 2 . 00 T= 0.1 O 2P 0 5 1 TION (M) 3. 00 O in I---- O(W) •• • •• • • **»*» • •» • •• 3= 1— P 00 L U .. " • • ' T " ' *"f *"*'* T * ' ' I ' I- , . ...... - . T-- T ■ T - -- I - - -UJ '-I— I— r " ' ' '- -' ' - -- -V » • r * ; 1 • • • ° 0 . 00 0 . 50 1 . 00 1. 50 2 . 00 2 . 5 0 3 . 00 T = O . O O Q P O S I T I O N (M) 179 o LO (M ) 0. X Oin . 00 ' f T f f i 1 ’ I ' ' I I 0 . 50 1. 00 1. 50 2 . 00 T= 0 . 5 0 Cf O S I T I 0 N (M) T 1J e I r e * 2 . 50 f—*"* 3. I----- 0:w ) X 6 )I ~ Oin g J i T J T T I I ' I ' 1 I e e T e 0 . 50 1 . 00 1 . 50 2 . 00 T= 0 . q o c f OS I T I ON (M) 1 ' I ■ ' ■ 2 . 50 3. x ° b OC^ OO ."""T mVm T - - T - ---K*,. «■ .»■ I,' f-...,■»...« ", " ■ l"'" , '0 . 50 1 . 00 1. 50 2 . 00 T= 0 . 3 0 0 P 0 S I T I 0 N (M) 2 . 5 0 T"e"e e 3. 180 o i n OCt).00 1. 00 1. 50 0 0 2 ^ 0 S I T I 0 N 2.00 3 . 00 O LD / f C c r n T D : ? ^ o c T)1.00 ' o'. 50 I'.OO 1.50 ' 2.00 T = 0 . 7 0 2 P 0 S I T I 0 N (M) "*~T ' ’ • 2 . 50 3 . 00 L U - .... . T Q C b.oo 0.50 1.00 1.50 T= o . 602p 0SITI0N D E P T H (M ) D E P T H (M ) JD .. O O 0 . 5 0 _p . 00 0 . 5 0 181 ) 0.50 1.00 1.50 T= I . OODf3OS I T I ON 2.00 3.00 J 0.50 1.00 1.50 2.00 T= 0 . 9 0 Dp 0 S I T I 0 N (M) * • I • ' ■ 2.50 3.00 M ONTANA STA TE U N IV E R SIT Y L IB R A R IE S stks D378.D434@Theses A biviscous modified Bingham model of sn 3 1762 00177748 9