% PDC_Splash_All % 1D model for dilute pyroclastic density current density, temperature, and % runout function [RO,RHOM,H,U,MASS,TM,X,VRVT,PVEL,i,PSAVE]=PDC_1D(type,TR,MP1,RHOP,MA1,CA,CR,R,P,M,E,G,Fr,DT,TA,RHOA,D,PVEL,prescribe) %Mass of particles in different size fractions % MP = MP1; % initial mass of particles (kg) VP = MP./RHOP; VT = sum(VP); % Total particle volume (m^3) MR = sum(MP); % Total particle mass RHOS = RHOP(1); % Substrate density (make same as small particles) mew= 0.5; %tangent to the friction angle nu = 2*10^-5; % dynamic viscosity of air Pa s alpha = 0.007;% Energy loss coefficient % Initial Height of the current % MA(1)=MA1;%Initial mass of air in the current A=2.5*10^-3;% Thermal expansion coefficient of air (1/k) %%% Initial conditions %%% TE(1)=CR*TR*MR+CA*TA*MA(1);%thermal energy TM(1) = TE(1)/(CR*MR+CA*MA(1)); %temperature of the mixture MASS(1)=MR+MA(1);%Total Mass AIRVOLUME(1)=MA(1)/(P/(M*R*TM(1)));%volume of air in the current RHOM(1) = (MASS(1))/(VT+AIRVOLUME(1)); % density of the rock, air mixture H(1)=VT+AIRVOLUME(1); VRVT(1) = VT/H(1);% volume fraction of rock in the current (Volume of rock/total volume) GPRIME = (RHOM(1)-RHOA)/RHOA*G; % Gprime U(1)=Fr*sqrt(GPRIME*H(1)); % Current velocity with FR=1 VENT(1)=E*U(1); % Entrainment velocity X(1)=0;%initial current location Test=RHOM(1)/RHOA; % Density ratio i=2; Height=H; % Now have the current advance (no sedimentation) while Test>1 %while GPRIME>0.01 DX = DT*U(i-1); %Runout distance of the current X(i)=X(i-1)+DX; MENT = DT*VENT(i-1)*RHOA; %Mass of added air MA(i) = MA(i-1)+MENT; %Mass of air in the current %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Calculate settling velocities of particles in the current % u0 = 0; %initialize sedimentation speeds m/s, large particles u2 = 0; %initialize sedimentation speeds m/s, small particles if type ~=1 %Type = 1 corresponds to a current that does not sediment particles % Calculate settling velocities via Equations 4 - 11 in Dufek et al. % % Iterate to find terminal settling velocities of two particle types % preliminaries % d2 = D(1); % small particle diameter d0 = D(2); % large particle diameter rhop2 = RHOP(1); % small particle density rhop0 = RHOP(2); % large particle density rhog = MA(i)/AIRVOLUME(i-1); % density of air in the current alpha2 = VRVT(i-1); % particle volume fraction in the current du = 1; % initialize du value dt = 10^-3; % time step for each iteration (s) % Now iteratively solve for the terminal settling speed of small particles while abs(du) > dt Rep = u2*d2*rhog/nu; % particle Reynolds number f = 1+0.15*Rep^(0.687)+0.0175/(1+42500*Rep^-1.16); % Dufek et al. (2009) eq. 6 tp = (rhop2-rhog)*d2^2/(18*nu); % response timescale Fg = f/tp*u2*(2*rhop2/(2*rhop2+rhog)); % force of the gas on the particle Fp = G*((2*rhop2-2*rhog)/(2*rhop2+rhog)); % force of the particle du = (-Fg+Fp)*dt; % determine if gravitational and drag forces balance u2 = u2 + du; % If now, iterate until the correct terminal velocity is found end % Now iteratively solve for the terminal settling speed of large particles du = 1; while abs(du) > dt Rep = u0*d0*rhog/nu; % particle Reynolds number f = 1+0.15*Rep^(0.687)+0.0175/(1+42500*Rep^-1.16); % Dufek et al. (2009) eq. 6 tp = (rhop0-rhog)*d0^2/(18*nu); % response timescale Fg = f/tp*u0*(2*rhop0/(2*rhop0+rhog)); % force of the gas on the particle Fp = G*((2*rhop0-2*rhog)/(2*rhop0+rhog)); % force of the particle % Now calculate force force from particle - particle impacts % %alpha2 = .01; %Small particle volume fraction, this should come from code input theta = .1*u0; %Granular temperature e = 0.51; % restitution coefficient go = 1; % ~1 for dilute contitions %%% Force from particle - particle impacts Fs = ((6*alpha2*sqrt(theta)*go*(d2+d0)^2)/(sqrt(pi)*d2^3))*((rhop2*d2^3*(1+e)*(u2-u0))/(rhop2*d2^3+rhop0*d0^3)); du = (-Fg+Fs+Fp)*dt; % determine if gravitational and drag forces balance u0 = u0 + du; % If now, iterate until the correct terminal velocity is found end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Assign settling velocity values if you do not prescribe them if prescribe == 0 PVEL(1)=u2;%u2 PVEL(2)=u0; end if i==2 PSAVE(1)=PVEL(1); PSAVE(2)=PVEL(2);%save initial values end %mass of particles sedimented from the current in each size class (MP must be adjusted for the mass lost from each size class) PS(1)=PVEL(1).*MP(i-1,1)./H(i-1)*DT; PS(2)=PVEL(2).*MP(i-1,2)./H(i-1)*DT; %make sure that more particles aren't removed than exist in the current% if PS(1) > MP(i-1,1) PS(1) = MP(i-1,1); end if PS(2) > MP(i-1,2) PS(2) = MP(i-1,2); end %Subtract the particles that were sedimented from the current from the mass of the current by size class MP(i,1)=MP(i-1,1)-PS(1); MP(i,2)=MP(i-1,2)-PS(2); % Mass of particles suspended via a Splash function, only calculate for % entraining current types (i.e., type ~=1 and type ~=2 if type~=1 && type~=2 DC = 0.14*mew^-1.*(RHOP./(RHOS)).^0.5.*D.^(2/3).*(PVEL./(2*G)).^(1/3); %crater depth from particle impacts (m) SPS=alpha.*PS.*PVEL.^2./(2*G.*DC); else SPS = 0; end % Calculate changes in current mass, volume, temperature, and densities MR(i)= MR(i-1)-sum(PS)+sum(SPS); % sum of the particle mass in the current VT(i)= VT(i-1)-sum(PS./RHOP)+sum(SPS)./RHOS; % total volume of particles % Determine the temperature of the entrained particles if type == 3 % cold entrained particles Tentrain = TA; elseif type == 4 % hot entrained particles Tentrain = TM(1); else % currents that don't entrain particles Tentrain = 0; end TE(i)=TE(i-1)+CA*TA*MENT-sum(PS)*TM(i-1)*CR+sum(SPS)*Tentrain*CR; %Thermal energy of the mixture TE(i)=TE(i-1)+CA*TA*MENT; TM(i)=TE(i)/(CR*MR(i)+CA*MA(i)); %Temperature of the mixture % Calculate changes in current height and density due to air % entrainment and thermal expansion MASS(i)=MR(i)+MA(i);%Total Mass AIRVOLUME(i)=MA(i)/(P/(M*R*TM(i)));%volume of air in the current %H(i)=H(i-1)+(VENT(i-1)+((TM(i)-TM(i-1))/DT)*A*(H(i-1)+DT*VENT(i-1))-PVEL(1))*DT;%change in current height due to thermal expansion of air in the current and the change in current height %H(i)=VT(i)+AIRVOLUME(i)-PVEL(1)*DT; %old H(i) H(i)=VT(i)+AIRVOLUME(i); VRVT(i) = VT(i)/H(i);% volume fraction of rock in the current (Volume of rock/total volume) RHOM(i) = (MASS(i))/H(i); % density of the rock, air mixture GPRIME = (RHOM(i)-RHOA)/RHOA*G; % Gprime, a bug was found in this line and corrected on May 10,2017 U(i)=Fr*sqrt(GPRIME*H(i)); VENT(i)=E*U(i); Height=H(i); Test=RHOM(i)/RHOA; if X(i)>50^4 % end if the current runs out too far break end i=i+1;%counter end RO=X(end); % save runout distance end