Alpine3D

Alpine3D Svn Source Tree

Root/trunk/alpine3d/snowdrift/Sublimation.cc

  • Property svn:executable set to *
1/***********************************************************************************/
2/* Copyright 2009-2015 WSL Institute for Snow and Avalanche Research SLF-DAVOS */
3/***********************************************************************************/
4/* This file is part of Alpine3D.
5 Alpine3D is free software: you can redistribute it and/or modify
6 it under the terms of the GNU Lesser General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
9
10 Alpine3D is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU Lesser General Public License for more details.
14
15 You should have received a copy of the GNU Lesser General Public License
16 along with Alpine3D. If not, see <http://www.gnu.org/licenses/>.
17*/
18/*************************************************************************/
19/* calculate sublimation drifting snow */
20/*************************************************************************/
21/*------------------------------------------------------------------------+
22 | |
23 +-----------------------------------------------------------------------*/
24/*************************************************************************/
25#include <alpine3d/snowdrift/SnowDrift.h>
26
27using namespace mio;
28
29/**
30 * @brief Sublimation
31 * Calculates drifting snow sublimation at every node of the concentration field
32 */
33void SnowDriftA3D::Sublimation()
34{
35nodes_Subl.grid3D = 0.;
36
37for (unsigned int ix=0; ix<nx;ix++) {
38for (unsigned int iy=0; iy<ny; iy++) {
39for (unsigned int iz=1; iz<nz; iz++){ //start at second level since below > BC saturation
40//calculate magnitude windspeed instead of u,v,w
41const double WindVelocity = sqrt(mio::Optim::pow2(nodes_u.grid3D(ix,iy,iz))+mio::Optim::pow2(nodes_v.grid3D(ix,iy,iz))+mio::Optim::pow2(nodes_w.grid3D(ix,iy,iz)));
42if ( nodes_c.grid3D(ix,iy,iz)>1e-9 && nodes_RH.grid3D(ix,iy,iz)<1.0){
43const double repRadius = 6.25e-5; //representative radius for sublimation of a particle size distribution
44const double dmdt_ini = calcSubldM(repRadius, nodes_Tair.grid3D(ix,iy,iz), nodes_RH.grid3D(ix,iy,iz), WindVelocity, nodes_z.grid3D(ix,iy,iz));
45nodes_Subl.grid3D(ix,iy,iz) = calcS(nodes_c.grid3D(ix,iy,iz),repRadius,dmdt_ini);
46}
47}
48}
49}
50}
51
52/**
53 * @brief Calculate sublimation mass change
54 * This function calculates the mass change of a single ice particle for the given conditions
55 * @param Radius radius of the ice particle
56 * @param AirTemperature
57 * @param RH relative humidity
58 * @param WindSpeed magnitude of the wind
59 */
60double SnowDriftA3D::calcSubldM(const double Radius,const double AirTemperature,const double RH, const double WindSpeed, const double altitude)
61{
62const double sigma = RH - 1.; //Undersaturation of the surrounding air
63
64//Nusselt number Nu
65double Nu;
66const double Re = reynoldsNumberforFallingParticles(Radius, WindSpeed, AirTemperature, RH, altitude); //Reynolds number
67if (Re<=10.) {
68Nu = 1.79+0.606*sqrt(Re);
69} else {
70Nu = 1.88+0.580*sqrt(Re);
71}
72
73//Sherwood number
74const double Sh = Nu;
75
76//diffusivity of water vapor in the atmosphere
77const double diffusivity = 2.11e-5 * pow((AirTemperature/273.15), 1.75);
78
79//water vapor density in the atmosphere
80double const SatVaporDensity = Atmosphere::waterVaporDensity(AirTemperature,Atmosphere::vaporSaturationPressure(AirTemperature));
81
82//Calculate dm/dt according to Thorpe and Mason (we split it in 3 parts)
83const double temp0 = (Constants::lh_sublimation*molecularWeightofWater) / (Constants::gas_constant_mol*AirTemperature) - 1.;
84const double temp1 = 2.*Constants::pi*Radius*sigma;
85const double temp2 = Constants::lh_sublimation / (thermalConductivityofAtm*AirTemperature*Nu) * temp0 +
86 1. / (diffusivity * SatVaporDensity * Sh) ;
87
88//And return the value for dm/dt for 1 particle
89const double dmdt = temp1/temp2;
90return dmdt;
91}
92
93/**
94 * @brief Calculate sublimation
95 * The mass change of a single particle is now extended to all particles
96 * @param concentration snow concentration
97 * @param sublradius radius of a single particle
98 * @param dmdt mass change of a single particle due to sublimation
99 */
100double SnowDriftA3D::calcS(const double concentration,const double sublradius, const double dmdt)
101{
102const double subllossrate = (concentration / (4./3.*Constants::pi*pow(sublradius,3)*Constants::density_ice)) * dmdt;
103return subllossrate;
104}
105
106/**
107 * @brief Reynolds number
108 * Calculate the Reynolds number for a falling particle
109 * @param Radius radius of the particle
110 * @param Windspeed magnitude of the wind
111 * @param AirTemperature
112 * @param RH relative humidity of the air
113 */
114double SnowDriftA3D::reynoldsNumberforFallingParticles(const double Radius, const double Windspeed,const double AirTemperature,const double RH, const double altitude)
115{
116const double tmp_venVelocity = ventilationVelocity(Radius, Windspeed, AirTemperature, RH, altitude);
117const double ReynoldsNumber = (2.*Radius*tmp_venVelocity) / kinematicViscosityAir;
118return ReynoldsNumber;
119}
120
121/**
122 * @brief Ventilation velocity
123 * @param Radius radius of the particle
124 * @param Windspeed magnitude of the wind
125 * @param AirTemperature
126 * @param RH relative humidity of the air
127 * @param altitude Altitude above sea level
128 */
129double SnowDriftA3D::ventilationVelocity(const double Radius, const double Windspeed, const double AirTemperature, const double RH, const double altitude)
130{
131const double MeanTerminalFallVelocity = terminalFallVelocity(Radius, AirTemperature, RH, altitude);
132const double FluctuatingComponent = 7.5e-3*sqrt(2)*pow(Windspeed, 1.36);
133const double VentilationVelocity = MeanTerminalFallVelocity + FluctuatingComponent;
134
135return VentilationVelocity;
136}
137
138/**
139 * @brief Terminal fall velocity
140 * NB assumed to be 0.5
141 * @param Radius radius of the particle
142 * @param Temperature air temperature (K)
143 * @param RH relative humidity of the air (between 0 and 1)
144 * @param altitude Altitude above sea level
145 */
146double SnowDriftA3D::terminalFallVelocity(const double /*Radius*/, const double /*Temperature*/, const double /*RH*/ , const double /*altitude*/)
147{
148//problems for radius larger than 0.000054>> try constant velocity 0.5
149return 0.5;//settling velocity assumed 0.5 m/s
150
151/*const double tmp1 = (6.203*kinematicViscosityAir) / Radius;
152const double Psat = Atmosphere::vaporSaturationPressure(Temperature);
153const double sigma = Constants::density_ice / ( Atmosphere::stdDryAirDensity(altitude, Temperature) +
154 Atmosphere::waterVaporDensity(Temperature, Psat*RH) );
155const double vr = .5 * (-tmp1 + sqrt((tmp1*tmp1)-4.*(-1.379)*mio::Cst::gravity*sigma*Radius) );
156return vr;*/
157}
158
159/**
160 * @brief Relative humidity
161 * Gives the relative humidity for a given air temperature and specific humidity
162 * @param AirTemp Air temperature
163 * @param qi specific humidity
164 * @param altitude Altitude above sea level
165 */
166double SnowDriftA3D::RH_from_q(const double AirTemp, const double qi, const double altitude)
167{
168const double SatVaporDensity = Atmosphere::waterVaporDensity(AirTemp,Atmosphere::vaporSaturationPressure(AirTemp));
169const double tmp = (qi/(1.-qi)) * Atmosphere::stdDryAirDensity(altitude, AirTemp);
170
171const double RH = std::min(tmp/SatVaporDensity , 1.5);
172return RH;
173}
174

Archive Download this file

Revision: HEAD