Up: Main Previous: Bibliography
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define STEFAN 5.669e-8 /* Stefan-Boltzmann constant */
#define EMISS 0.90 /* Emissivity of concrete surface */
#define SOLARABS 0.65 /* Solar absorptivity of concrete */
/************************************************************************/
/* */
/* Program: Tempcalc.c */
/* Program to calculate concrete surface temperature during a */
/* yearly weather cycle (based on a 1-D finite difference solution */
/* */
/* Inputs: File inputs.par with material properties of concrete */
/* and soil. */
/* Weather data file XXXXX.tm2 from NREL database */
/* Outputs: Files with */
/* 1) surface temperature history (temphist.out) */
/* 2) time-of-wetness history (wettime.out) */
/* 3) time-of-freezing history (freezetime.out) */
/* */
/* Programmer: Dale P. Bentz */
/* Last revision: July 19, 2000 */
/* Effort funded by the Federal Highway Administration */
/* */
/************************************************************************/
/* Global variables */
double tempnew[22],tempold[22],t_out=0.0,sum_t,tk_sky;
double conc_heat_cap,conc_thermal_k,conc_density;
double h_coeff=0.0,alpha_conc,t_cum=0.0;
double soil_heat_cap,soil_thermal_k,soil_density,alpha_soil;
/* Routine to read in input property parameters for concrete and soil */
/* Called by: main program */
/* Calls: none */
read_inputs(){
FILE *parfile;
parfile=fopen("inputs.par","r");
/* Note that all inputs are in SI units */
fscanf(parfile,"%lf",&conc_heat_cap);
printf("Concrete heat capacity is %lf J/kg/C\n",conc_heat_cap);
fscanf(parfile,"%lf",&conc_thermal_k);
printf("Concrete thermal conductivity is %lf W/m/C\n",conc_thermal_k);
fscanf(parfile,"%lf",&conc_density);
printf("Concrete density is %lf kg/m/m/m\n",conc_density);
alpha_conc=conc_thermal_k/conc_heat_cap/conc_density;
printf("Calculated concrete alpha value is %.12lf m*m/s\n",alpha_conc);
fscanf(parfile,"%lf",&soil_heat_cap);
printf("Soil heat capacity is %lf J/kg/C\n",soil_heat_cap);
fscanf(parfile,"%lf",&soil_thermal_k);
printf("Soil thermal conductivity is %lf W/m/C\n",soil_thermal_k);
fscanf(parfile,"%lf",&soil_density);
printf("Soil density is %lf kg/m/m/m\n",soil_density);
alpha_soil=soil_thermal_k/soil_heat_cap/soil_density;
printf("Calculated soil alpha value is %.12lf m*m/s\n",alpha_soil);
fclose(parfile);
}
/* Routine to update all temperatures in the finite difference grid */
/* Called by: main program */
/* Calls: none */
upd_temp(q_in,t_amb,t_dew,cloud_cover,delta_x,delta_t)
double q_in,t_amb,t_dew,cloud_cover,delta_x,delta_t;
/* q_in is in W/m*m
t_amb and t_dew are in degrees Celsius
cloud_cover is in tenths (0.1-1.0)
delta_x is in meters
delta_t is in seconds */
{
int i;
FILE *tempfile;
double q_sun,q_sky,q_conv,q_cond,q_stor,q1,q2;
double tk_amb,tk_dew,e_s,cloud_fact;
tk_amb=t_amb+273.15; /* Conversion to Kelvin */
tk_dew=t_dew+273.15; /* Conversion to Kelvin */
/* Sky temperature computed based on algorithm of Walton et al. (TARP) */
cloud_fact=1.0+0.024*cloud_cover-0.0035*cloud_cover*cloud_cover;
cloud_fact+=0.00028*cloud_cover*cloud_cover*cloud_cover;
/* Compute the sky emissivity e_s */
e_s=(0.787+0.764*log(tk_dew/273.))*cloud_fact;
/* Now estimate the sky temperature tk_sky */
tk_sky=sqrt(sqrt(e_s))*tk_amb;
/* Compute the four heat flow contributions at the top surface */
q_sun=q_in*SOLARABS;
q_sky=STEFAN*EMISS*(pow((tempold[0]+273.15),4.0)-pow(tk_sky,4.0));
q_conv=h_coeff*(tempold[0]-t_amb);
q_cond=(conc_thermal_k/delta_x)*(tempold[0]-tempold[1]);
q_stor=q_sun-q_sky-q_conv-q_cond;
/* Update the surface temperature */
tempnew[0]=tempold[0]+q_stor*2.*delta_t/delta_x/conc_density/conc_heat_cap;
/* First, the top concrete layer */
/* Equations based on those presented in Heat Transfer by Holman */
for(i=1;i<10;i++){
tempnew[i]=((alpha_conc*delta_t)/(delta_x*delta_x))*
(tempold[i+1]+tempold[i-1])+(1.0-2.*alpha_conc*delta_t/
(delta_x*delta_x))*tempold[i];
}
/* Second, the interface between concrete and soil */
q1=(conc_thermal_k/delta_x)*(tempold[10]-tempold[9]);
q2=(soil_thermal_k/delta_x)*(tempold[10]-tempold[11]);
tempnew[10]=tempold[10]-(q1+q2)*2.*delta_t/delta_x/
(conc_density*conc_heat_cap+soil_density*soil_heat_cap);
/* Last, the soil subbase */
for(i=11;i<20;i++){
tempnew[i]=((alpha_soil*delta_t)/(delta_x*delta_x))*
(tempold[i+1]+tempold[i-1])+(1.0-2.*alpha_soil*delta_t/
(delta_x*delta_x))*tempold[i];
}
/* Last soil node is constant at a temperature of 13 C */
tempnew[20]=13.0;
/* Copy the new temperatures over the old ones to prepare for the next update */
for(i=0;i<=20;i++){
tempold[i]=tempnew[i];
}
/* Output temperature results once per hour */
if((t_cum-t_out)>=3600.){
tempfile=fopen("temphist.out","a");
fprintf(tempfile,"%.2lf %.2lf %.2lf %.2lf ",
t_cum/3600.,t_amb,t_dew,tk_sky-273.15);
t_out=t_cum;
/* Output surface temperature */
for(i=0;i<=1;i++){
fprintf(tempfile,"%lf ",tempnew[i]);
}
/* And temperature for every fifth node into concrete and soil */
for(i=5;i<=20;i+=5){
fprintf(tempfile,"%lf ",tempnew[i]);
}
fprintf(tempfile,"\n");
fclose(tempfile);
}
}
/* Main program to read in weather data file and call routines to */
/* set material properties and update temperatures */
/* Calls: read_inputs and upd_temp */
main(){
FILE *infile,*wetfile,*frfile;
char wba[80],city[80],state[80],latd[10],longd[10];
char filein[80];
register int i,j;
int tzone,lat1,lat2,long1,long2,elev;
int wetflag=0,duration,numcum=0,freezeflag=0;
long int timeorg,timefin;
long int timeforg,timeffin;
int year,month,day,hour,exhr,exdnr,glhr,glhrflu,dnr,dnrflu;
int dfhr,dfhrflu,glhi,glhiflu,dhi,dhiflu,dni,dniflu;
int zi,ziflu,tskycov,tskycovflu,opqskycov,opqskycovflu,tempdry,tdflu;
int tempwet,twflu,rh,rhflu,atmp,atmpflu,wdir,wdirflu,wspd,wspdflu;
int vis,visflu,ceilflu,days[13];
long int ceil,durcum;
double qsun,temp_amb,temp_dew,del_x,del_t,del_t_soil,tfact;
double temp_amb_old,wspd_old,temp_dew_old,wspd_new;
double temp_amb_cur,wspd_cur,temp_dew_cur,tfmin=0.0;
float rhsum=0.0,rhave,rhinit,twinit=0.0;
int pw0,pw1,pw2,pw3,pw4,pw5,pw6,pw7,pw8,pw9,prewat,prewatflu;
char glhrfls,dnrfls,dfhrfls,glhifls,dhifls,dnifls,rest[240];
char zifls,tskycovfls,opqskycovfls,tdfls,twfls,rhfls,atmpfls,wdirfls;
char wspdfls,visfls,ceilfls,prewatfls;
read_inputs();
del_x=0.02; /* 2 cm spacing in nodes for 40 cm total */
durcum=0;
/* Compute cumulative days per month */
days[1]=0;
days[2]=31;
days[3]=days[2]+29;
days[4]=days[3]+31;
days[5]=days[4]+30;
days[6]=days[5]+31;
days[7]=days[6]+30;
days[8]=days[7]+31;
days[9]=days[8]+31;
days[10]=days[9]+30;
days[11]=days[10]+31;
days[12]=days[11]+30;
/* Get name of weather data file to use for this simulation */
printf("Enter name of file to read \n");
scanf("%s",filein);
printf("%s\n",filein);
infile=fopen(filein,"r");
/* Open log files for recording time-of-wetness and time-of-freezing */
wetfile=fopen("wettime.out","w");
frfile=fopen("freezetime.out","w");
/* Read in the header for this weather data file and echo values */
fscanf(infile,"%s %s %s %d %s %d %d %s %d %d %d",
wba,city,state,&tzone,latd,&lat1,&lat2,longd,&long1,&long2,&elev);
printf("Number is %s\n",wba);
printf("City is %s\n",city);
printf("State is %s\n",state);
printf("Time zone is %d\n",tzone);
/* Set time step in finite difference routine */
/* based on thermal properties of concrete and soil */
del_t=del_x*del_x/4./alpha_conc;
del_t_soil=del_x*del_x/4./alpha_soil;
if(del_t_soil<del_t){del_t=del_t_soil;}
/* Scale time step to a factor of 3600 (seconds in an hour) */
if(del_t<50.0){
del_t=10.;
}
else if((del_t>50.0)&&(del_t<100.0)){
del_t=50.0;
}
else if((del_t>100.0)&&(del_t<200.0)){
del_t=100.0;
}
else if((del_t>200.0)&&(del_t<300.0)){
del_t=200.0;
}
else if((del_t>300.0)&&(del_t<400.0)){
del_t=300.0;
}
else if(del_t>400.0){
del_t=400.0;
}
printf("del_t is %lf \n",del_t);
fflush(stdout);
/* Read in the first record */
/* Format statement is directly from TMY2DATA (NREL) Manual */
fscanf(infile,"%2d%2d%2d%2d%4d%4d%4d%1s%1d%4d%1s%1d%4d%1s%1d%4d%1s%1d%4d
%1s%1d%4d%1s%1d%4d%1s%1d%2d%1s%1d%2d%1s%1d%4d%1s%1d%4d%1s%1d%3d%1s%1d%4d%1s%1d%3
d%1s%1d%3d%1s%1d%4d%1s%1d%5ld%1s%1d%1d%1d%1d%1d%1d%1d%1d%1d%1d%1d%3d%1s%1d%s",
&year,&month,&day,&hour,&exhr,&exdnr,&glhr,&glhrfls,&glhrflu,&dnr,&dnrfls,
&dnrflu,&dfhr,&dfhrfls,&dfhrflu,&glhi,&glhifls,&glhiflu,&dni,&dnifls,&dniflu,
&dhi,&dhifls,&dhiflu,&zi,&zifls,&ziflu,&tskycov,&tskycovfls,&tskycovflu,
&opqskycov,&opqskycovfls,&opqskycovflu,&tempdry,&tdfls,&tdflu,&tempwet,
&twfls,&twflu,&rh,&rhfls,&rhflu,&atmp,&atmpfls,&atmpflu,&wdir,&wdirfls,
&wdirflu,&wspd,&wspdfls,&wspdflu,&vis,&visfls,&visflu,&ceil,&ceilfls,
&ceilflu,&pw0,&pw1,&pw2,&pw3,&pw4,&pw5,&pw6,&pw7,&pw8,&pw9,&prewat,
&prewatfls,&prewatflu,rest);
/* Extract the relevant values */
temp_amb_old=(double)tempdry/10.;
wspd_old=(double)wspd/10.;
temp_dew_old=(double)tempwet/10.;
/* Initialize temperature data to ambient temperature value */
for(j=0;j<=20;j++){
tempold[j]=temp_amb_old;
}
tempold[20]=13.0;
/* Now process each weather record in turn */
for(i=1;i<8760;i++){
fscanf(infile,"%2d%2d%2d%2d%4d%4d%4d%1s%1d%4d%1s%1d%4d%1s%1d%4d%1s%1d%4d
%1s%1d%4d%1s%1d%4d%1s%1d%2d%1s%1d%2d%1s%1d%4d%1s%1d%4d%1s%1d%3d%1s%1d%4d%1s%1d%3
d%1s%1d%3d%1s%1d%4d%1s%1d%5ld%1s%1d%1d%1d%1d%1d%1d%1d%1d%1d%1d%1d%3d%1s%1d%s",
&year,&month,&day,&hour,&exhr,&exdnr,&glhr,&glhrfls,&glhrflu,&dnr,&dnrfls,
&dnrflu,&dfhr,&dfhrfls,&dfhrflu,&glhi,&glhifls,&glhiflu,&dni,&dnifls,
&dniflu,&dhi,&dhifls,&dhiflu,&zi,&zifls,&ziflu,&tskycov,&tskycovfls,
&tskycovflu,&opqskycov,&opqskycovfls,&opqskycovflu,&tempdry,&tdfls,&tdflu,
&tempwet,&twfls,&twflu,&rh,&rhfls,&rhflu,&atmp,&atmpfls,&atmpflu,&wdir,
&wdirfls,&wdirflu,&wspd,&wspdfls,&wspdflu,&vis,&visfls,&visflu,&ceil,
&ceilfls,&ceilflu,&pw0,&pw1,&pw2,&pw3,&pw4,&pw5,&pw6,&pw7,&pw8,&pw9,&prewat,
&prewatfls,&prewatflu,rest);
/* Extract the relevant values */
temp_amb=(double)tempdry/10.;
temp_dew=(double)tempwet/10.;
wspd_new=(double)wspd/10.;
qsun=glhr;
sum_t=0.0;
/* Now process one hours worth of data */
while(sum_t<3600.){
/* Update the hourly and cumulative times */
sum_t+=del_t;
t_cum+=del_t;
tfact=(sum_t)/3600.;
/* Linearly interpolate the temperatures and wind speeds */
/* within any 1 h period */
temp_amb_cur=temp_amb_old+tfact*(temp_amb-temp_amb_old);
temp_dew_cur=temp_dew_old+tfact*(temp_dew-temp_dew_old);
wspd_cur=wspd_old+tfact*(wspd_new-wspd_old);
/* Convection coefficient based on wind speed from FEMMASSE */
/* of Intron in The Netherlands (3/00) */
if(wspd_cur<5.0){
h_coeff=5.6+4.0*wspd_cur;
}
else{
h_coeff=7.2*pow(wspd_cur,0.78);
}
/* Call the routine to update all temperatures */
upd_temp(qsun,temp_amb_cur,temp_dew_cur,(double)tskycov/10.,del_x,del_t);
}
/* Now look for freezing and wetting events */
temp_amb_old=temp_amb;
temp_dew_old=temp_dew;
wspd_old=wspd_new;
/* Freeze events */
if((freezeflag==0)&&(tempold[0]<(0.0))){
freezeflag=1;
timeforg=(day*24+days[month]*24+hour);
tfmin=tempold[0];
}
/* Log minimum temperature achieved during freezing */
if((freezeflag==1)&&(tfmin>tempold[0])){tfmin=tempold[0];}
if((freezeflag==1)&&(tempold[0]>=(0.0))){
freezeflag=0;
timeffin=(day*24+days[month]*24+hour);
duration=timeffin-timeforg;
fprintf(frfile,"%ld %d %lf\n",timeforg,duration,tfmin);
}
/* Wet if rain (precipitation) or dew (condensation) event */
if(((wetflag==0)&&(pw0==0)&&((pw2!=9)||(pw6!=9)||(pw3!=9)))||
((wetflag==0)&&((tempold[0]<temp_dew)&&(tempold[0]>0.0)))){
wetflag=1;
rhsum+=(float)rh;
rhinit=(float)rh;
twinit=tempold[0];
printf("Wet from %d:00 %d/%d/%d with T= %lf ",hour,
month,day,year,tempnew[0]);
timeorg=(day*24+days[month]*24+hour);
}
if((wetflag==1)&&(pw0==0)&&((pw2==9)&&(pw6==9)&&(pw3==9))&&
(tempold[0]>=temp_dew)){
wetflag=0;
printf("until %d:00 %d/%d/%d ",hour,month,day,year);
timefin=(day*24+days[month]*24+hour);
duration=timefin-timeorg;
durcum+=duration;
numcum+=1;
printf("of duration %d hours\n",duration);
fprintf(wetfile,"%ld %d %f %f \n",timeorg,duration,
(twinit+tempold[0])/2.,rhinit);
}
}
printf("\n");
printf("total cumulative time of wetness is %ld in %d separate events\n",
durcum,numcum);
rhave=rhsum/(float)numcum;
printf("average rh before start of rainfall is %f \n",rhave);
fclose(infile);
fclose(wetfile);
fclose(frfile);
}
Up: Main Previous: Bibliography