Next: Example input datafiles
Up: Computer programs for
Previous: Listing for hydreal3d.c
/************************************************************************/
/* */
/* Program disreal3d.c to hydrate three-dimensional cement */
/* particles in a 3-D box with periodic boundaries. */
/* Use cellular automata techniques and preserves correct */
/* hydration volume stoichiometry. */
/* Programmer: Dale P. Bentz */
/* Building and Fire Research Laboratory */
/* NIST */
/* Building 226 Room B-350 */
/* Gaithersburg, MD 20899 USA */
/* (301) 975-5865 FAX: 301-990-6891 */
/* E-mail: dale.bentz@nist.gov */
/* */
/************************************************************************/
/* Three-dimensional version created 7/94 */
/* General C version created 8/94 */
/* Modified to have pseudo-continuous dissolution 11/94 */
/* Difference between disreal3.c and this program, disreal3d.c, */
/* is that in this program, dissolved silicates are located close */
/* to dissolution source within a 17*17*17 box centered at dissolution source */
/* Heat of formation data added: 12/92 */
/* Hydration under self-desiccating conditions included 9/95 */
/* Additions for adiabatic hydration conditions included 9/96 */
#include <stdio.h>
#include <math.h>
#define CUBEMAX 7 /* Maximum cube size for checking pore size */
#define CUBEMIN 3 /* Minimum cube size for checking pore size */
#define SYSIZE 100 /* System size in pixels */
#define DISBIAS 20.0 /* Dissolution bias- to change all dissolution rates */
/* All dissolution rates are divided by this value */
/* To examine early setting, value should be */
/* increased to perhaps 100 or so */
#define DETTRMAX 1200 /* Maximum allowed # of ettringite diffusing species */
#define CHCRIT 50.0 /* Scale parameter to adjust CH dissolution probability */
/* based on number of CH diffusing species */
#define C3AH6CRIT 10.0 /* Scale par. to adjust C3AH6 dissolution prob. */
/* based on potential sulfate (gypsum) in solution */
#define C3AH6GROW 0.01 /* Probability for C3AH6 growth */
#define ETTRGROW 0.002 /* Probability for ettringite growth */
#define C3AETTR 0.001 /* Probability for reaction of diffusing C3A
with ettringite */
/* define IDs for each phase used in model */
#define POROSITY 0
#define C3S 1
#define C2S 2
#define CH 3
#define CSH 4
#define C3A 5
#define GYPSUM 6
#define C4AF 7
#define C3AH6 8
#define ETTR 9
#define AFM 10
#define FH3 11
#define POZZ 12
#define INERT 13
#define INERTAGG 14
#define ABSGYP 15
#define DIFFCSH 16
#define DIFFCH 17
#define DIFFGYP 18
#define DIFFC3A 19
#define DIFFFH3 20
#define DIFFETTR 21
#define EMPTYP 22 /* Empty porosity due to self desiccation */
#define OFFSET 30 /* Offset for highlighted potentially soluble pixel */
/* define heat capacities for all components in J/g/C */
/* including free and bound water */
#define Cp_cement 0.75
#define Cp_pozz 0.75
#define Cp_agg 0.84
#define Cp_h2o 4.18 /* Cp for free water */
#define Cp_bh2o 2.2 /* Cp for bound water */
#define WN 0.23 /* water bound per gram of cement during hydration */
#define WCHSH 0.06 /* water imbibed per gram of cement during chemical
shrinkage (estimate) */
/* data structure for diffusing species - to be dynamically allocated */
/* Use of a doubly linked list to allow for easy maintenance */
/* (i.e. insertion and deletion) */
/* Added 11/94 */
struct ants{
char x,y,z,id;
struct ants *nextant;
struct ants *prevant;
};
/* data structure for elements to remove to simulate self-desiccation */
/* once again a doubly linked list */
struct togo{
int x,y,z,npore;
struct togo *nexttogo;
struct togo *prevtogo;
};
/* Global variables */
/* Microstructure stored in array mic of type char to minimize storage */
/* Initial particle IDs stored in array micpart (for assessing set point) */
static char mic [SYSIZE] [SYSIZE] [SYSIZE];
static short int micpart [SYSIZE] [SYSIZE] [SYSIZE];
/* counts for dissolved and solid species */
long int discount[EMPTYP+1], count[EMPTYP+1];
/* Counts for pozzolan reacted, initial pozzolan, gypsum, ettringite,
and initial porosity */
long int npr,nfill,ncsbar,netbar,porinit;
/* Initial clinker phase counts */
long int c3sinit,c2sinit,c3ainit,c4afinit;
long int nmade,ngoing,gypready,poregone,poretodo;
int cyccnt,cubesize,sealed;
/* Parameters for kinetic modelling ---- maturity approach */
float ind_time,temp_0,temp_cur,time_cur,E_act,beta,heat_cf,w_to_c,s_to_c;
float alpha_cur,heat_old,heat_new,cemmass,mass_agg,mass_water,mass_fill,Cp_now;
/* Arrays for dissolution probabilities for each phase */
float disprob[EMPTYP+1],disbase[EMPTYP+1],gypabsprob;
/* Arrays for specific gravities, molar volumes, heats of formation, and */
/* molar water consumption for each phase */
float specgrav[EMPTYP+1],molarv[EMPTYP+1],heatf[EMPTYP+1],waterc[EMPTYP+1];
/* Solubility flags and diffusing species created for each phase */
int soluble[EMPTYP+1], creates[EMPTYP+1];
int *seed; /* Random number seed */
struct ants *headant,*tailant;
/* Supplementary programs */
#include "ran1.c"
#include "burn3d.c"
#include "burnset.c"
#include "hydreal3d.c"
/* routine to initialize values for solubilities, molar volumes, etc. */
/* Called by main program */
/* Calls no other routines */
void init()
{
int i;
for(i=0;i<=EMPTYP;i++){
creates[i]=0;
soluble[i]=0;
disprob[i]=0.0;
disbase[i]=0.0;
}
/* soluble [x] - flag indicating if phase x is soluble */
/* disprob [x] - probability of dissolution (relative diss. rate) */
gypabsprob=0.01; /* One sulfate absorbed per 100 CSH units */
/* Note that for the first cycle, only the aluminates */
/* and gypsum are soluble */
soluble[C4AF]=1;
disprob[C4AF]=0.1/DISBIAS;
disbase[C4AF]=0.1/DISBIAS;
creates[C4AF]=POROSITY;
soluble[C3S]=0;
disprob[C3S]=0.95/DISBIAS;
disbase[C3S]=0.95/DISBIAS;
creates[C3S]=DIFFCSH;
soluble[C2S]=0;
disprob[C2S]=0.15/DISBIAS;
disbase[C2S]=0.15/DISBIAS;
creates[C2S]=DIFFCSH;
soluble[C3A]=1;
disprob[C3A]=0.8/DISBIAS;
disbase[C3A]=0.8/DISBIAS;
creates[C3A]=POROSITY;
soluble[GYPSUM]=1;
disprob[GYPSUM]=0.1;
disbase[GYPSUM]=0.1;
creates[GYPSUM]=POROSITY;
soluble[CH]=0;
disprob[CH]=0.1/DISBIAS;
disbase[CH]=0.1/DISBIAS;
creates[CH]=DIFFCH;
soluble[C3AH6]=1;
disprob[C3AH6]=0.5/DISBIAS;
disbase[C3AH6]=0.5/DISBIAS;
creates[C3AH6]=POROSITY;
/* Ettringite is initially insoluble */
soluble[ETTR]=0;
disbase[ETTR]=0.10/DISBIAS;
disprob[ETTR]=0.10/DISBIAS;
creates[ETTR]=DIFFETTR;
/* establish molar volumes and heats of formation */
/* molar volumes are in cm^3/mole */
/* heats of formation are in kJ/mole */
/* See paper by Fukuhara et al., Cem. & Conc. Res., 11, 407-14, 1981. */
/* values for Porosity are those of water */
molarv[POROSITY]=18.0;
heatf[POROSITY]=(-285.83);
waterc[POROSITY]=1.0;
specgrav[POROSITY]=1.0;
molarv[C3S]=71.0;
heatf[C3S]=(-2927.82);
waterc[C3S]=0.0;
specgrav[C3S]=3.21;
/* For improvement in chemical shrinkage correspondence */
/* Changed molar volume of C-S-H to 108 5/24/95 */
molarv[CSH]=108.;
heatf[CSH]=(-3283.0);
waterc[CSH]=4.0;
specgrav[CSH]=2.12;
molarv[CH]=33.1;
heatf[CH]=(-986.1);
waterc[CH]=1.0;
specgrav[CH]=2.24;
molarv[C2S]=52.;
heatf[C2S]=(-2311.6);
waterc[C2S]=0.0;
specgrav[C2S]=3.28;
molarv[C3A]=89.1;
heatf[C3A]=(-3587.8);
waterc[C3A]=0.0;
specgrav[C3A]=3.03;
molarv[GYPSUM]=74.2;
heatf[GYPSUM]=(-2022.6);
waterc[GYPSUM]=0.0;
specgrav[GYPSUM]=2.32;
molarv[C4AF]=128.;
heatf[C4AF]=(-5090.3);
waterc[C4AF]=0.0;
specgrav[C4AF]=3.73;
molarv[C3AH6]=150.;
heatf[C3AH6]=(-5548.);
waterc[C3AH6]=6.0;
specgrav[C3AH6]=2.52;
/* Changed molar volume of FH3 to 69.8 (specific gravity of 3.0) 5/23/95 */
molarv[FH3]=69.8;
heatf[FH3]=(-823.9);
waterc[FH3]=3.0;
specgrav[FH3]=3.0;
/* Changed molar volue of ettringite to 735 (spec. gr.=1.7) 5/24/95 */
molarv[ETTR]=735;
heatf[ETTR]=(-17539.0);
/* Each mole of ETTR which forms requires 32 moles of water, */
/* six of which are supplied by 3 moles of gypsum in forming ETTR */
/* leaving 26 moles to be incorporated from free water */
waterc[ETTR]=26.0;
specgrav[ETTR]=1.7;
molarv[AFM]=313;
heatf[AFM]=(-8778.0);
/* Each mole of AFM which forms requires 12 moles of water, */
/* two of which are supplied by gypsum in forming ETTR */
/* leaving 10 moles to be incorporated from free water */
waterc[AFM]=10.0;
specgrav[AFM]=1.99;
molarv[POZZ]=27.0;
/* Use heat of formation of SiO2 (quartz) for unreacted pozzolan */
/* Source- Handbook of Chemistry and Physics */
heatf[POZZ]=-907.5;
waterc[POZZ]=0.0;
specgrav[POZZ]=2.2;
/* Assume inert filler has same specific gravity and molar volume as SiO2 */
molarv[INERT]=27.0;
heatf[INERT]=0.0;
waterc[INERT]=0.0;
specgrav[INERT]=2.2;
molarv[ABSGYP]=74.2;
heatf[ABSGYP]=(-2022.6);
waterc[ABSGYP]=0.0;
specgrav[ABSGYP]=2.32;
molarv[EMPTYP]=18.0;
heatf[EMPTYP]=(-285.83);
waterc[EMPTYP]=0.0;
specgrav[EMPTYP]=1.0;
}
/* routine to check if a pixel located at (xck.yck,zck) is on an edge
/* (in contact with pore space) in 3-D system */
/* Called by passone */
/* Calls no other routines */
int chckedge(xck,yck,zck)
int xck,yck,zck;
{
int edgeback,x2,y2,z2;
int ip;
edgeback=0;
/* Check all six neighboring pixels */
for(ip=1;((ip<=6)&&(edgeback==0));ip++){
x2=xck;
y2=yck;
z2=zck;
switch (ip){
case 1:
x2+=1;
if(x2>=SYSIZE){x2=0;}
break;
case 2:
x2-=1;
if(x2<0){x2=SYSIZE-1;}
break;
case 3:
y2+=1;
if(y2>=SYSIZE){y2=0;}
break;
case 4:
y2-=1;
if(y2<0){y2=SYSIZE-1;}
break;
case 5:
z2+=1;
if(z2>=SYSIZE){z2=0;}
break;
case 6:
z2-=1;
if(z2<0){z2=SYSIZE-1;}
break;
default:
break;
}
if(mic [x2] [y2] [z2]==POROSITY){
edgeback=1;
}
}
return(edgeback);
}
/* routine for first pass through microstructure during dissolution */
/* low and high indicate phase ID range to check for surface sites */
/* Called by dissolve */
/* Calls chckedge */
void passone(low,high,cycid)
int low,high,cycid;
{
int i,xid,yid,zid,phid,iph,edgef;
/* gypready used to determine if any soluble gypsum remains */
if((low<=GYPSUM)&&(GYPSUM<=high)){
gypready=0;
}
/* Scan the entire 3-D microstructure */
for(xid=0;xid<SYSIZE;xid++){
for(yid=0;yid<SYSIZE;yid++){
for(zid=0;zid<SYSIZE;zid++){
/* Identify phase and update count */
phid=50;
for(i=low;((i<=high)&&(phid==50));i++){
if(mic [xid][yid][zid]==i){
phid=i;
/* Update count for this phase */
count[i]+=1;
if((low<=GYPSUM)&&(GYPSUM<=high)&&(i==GYPSUM)){
gypready+=1;
}
/* If first cycle, then accumulate initial counts */
if(cycid==1){
if(i==POROSITY){porinit+=1;}
else if(i==C3S){c3sinit+=1;}
else if(i==C2S){c2sinit+=1;}
else if(i==C3A){c3ainit+=1;}
else if(i==C4AF){c4afinit+=1;}
else if(i==GYPSUM){ncsbar+=1;}
else if(i==POZZ){nfill+=1;}
else if(i==ETTR){netbar+=1;}
}
}
}
if(phid!=50){
/* If phase is soluble, see if it is in contact with porosity */
if((cycid!=0)&&(soluble[phid]==1)){
edgef=chckedge(xid,yid,zid);
if(edgef==1){
/* Surface eligible species has an ID OFFSET greater than its original value */
mic [xid][yid][zid]+=OFFSET;
}
}
}
} /* end of zid */
} /* end of yid */
} /* end of xid */
}
/* routine to locate a diffusing csh species near dissolution source */
/* at (xcur,ycur,zcur) */
/* Called by dissolve */
/* Calls no other routines */
int loccsh(xcur,ycur,zcur)
int xcur,ycur,zcur;
{
int xpmax,ypmax,effort,tries,xmod,ymod,zmod;
struct ants *antnew;
effort=0; /* effort indicate appropriate location found */
tries=0;
/* 500 tries in immediate vicinity */
while((effort==0)&&(tries<500)){
tries+=1;
xmod=(-8)+(int)(17.*ran1(seed));
ymod=(-8)+(int)(17.*ran1(seed));
zmod=(-8)+(int)(17.*ran1(seed));
if(xmod>8){xmod=8;}
if(ymod>8){ymod=8;}
if(zmod>8){ymod=8;}
xmod+=xcur;
ymod+=ycur;
zmod+=zcur;
/* Periodic boundaries */
if(xmod>=SYSIZE){xmod-=SYSIZE;}
else if(xmod<0){xmod+=SYSIZE;}
if(zmod>=SYSIZE){zmod-=SYSIZE;}
else if(zmod<0){zmod+=SYSIZE;}
if(ymod<0){ymod+=SYSIZE;}
else if(ymod>=SYSIZE){ymod-=SYSIZE;}
if(mic[xmod][ymod][zmod]==POROSITY){
effort=1;
mic[xmod][ymod][zmod]=DIFFCSH;
nmade+=1;
ngoing+=1;
/* Add this diffusing species to the linked list */
antnew=(struct ants *)malloc(sizeof(struct ants));
antnew->x=xmod;
antnew->y=ymod;
antnew->z=zmod;
antnew->id=DIFFCSH;
/* Now connect this ant structure to end of linked list */
antnew->prevant=tailant;
tailant->nextant=antnew;
antnew->nextant=NULL;
tailant=antnew;
}
}
return(effort);
}
/* routine to count number of pore pixels in a cube of size boxsize */
/* centered at (qx,qy,qz) */
/* Called by makeinert */
/* Calls no other routines */
int countbox(boxsize,qx,qy,qz)
int boxsize,qx,qy,qz;
{
int nfound,ix,iy,iz;
int hx,hy,hz,boxhalf;
boxhalf=boxsize/2;
nfound=0;
/* Count the number of requisite pixels in the 3-D cube box */
/* using periodic boundaries */
for(ix=(qx-boxhalf);ix<=(qx+boxhalf);ix++){
hx=ix;
if(hx<0){hx+=SYSIZE;}
else if(hx>=SYSIZE){hx-=SYSIZE;}
for(iy=(qy-boxhalf);iy<=(qy+boxhalf);iy++){
hy=iy;
if(hy<0){hy+=SYSIZE;}
else if(hy>=SYSIZE){hy-=SYSIZE;}
for(iz=(qz-boxhalf);iz<=(qz+boxhalf);iz++){
hz=iz;
if(hz<0){hz+=SYSIZE;}
else if(hz>=SYSIZE){hz-=SYSIZE;}
/* Count if porosity, diffusing species, or empty porosity */
if((mic [hx] [hy] [hz]<C3S)||(mic[hx] [hy] [hz]>ABSGYP)){
nfound+=1;
}
}
}
}
return(nfound);
}
/* routine to create ndesire pixels of empty pore space to simulate */
/* self-desiccation */
/* Called by dissolve */
/* Calls countbox */
void makeinert(ndesire)
long int ndesire;
{
struct togo *headtogo,*tailtogo,*newtogo,*lasttogo,*onetogo;
long int idesire;
int px,py,pz,placed,cntpore,cntmax;
/* First allocate the linked list */
headtogo=(struct togo *)malloc(sizeof(struct togo));
headtogo->x=headtogo->y=headtogo->z=(-1);
headtogo->npore=0;
headtogo->nexttogo=NULL;
headtogo->prevtogo=NULL;
tailtogo=headtogo;
cntmax=0;
printf("In makeinert with %ld needed elements \n",ndesire);
fflush(stdout);
/* Add needed number of elements to the end of the list */
for(idesire=2;idesire<=ndesire;idesire++){
newtogo=(struct togo *)malloc(sizeof(struct togo));
newtogo->npore=0;
newtogo->x=newtogo->y=newtogo->z=(-1);
tailtogo->nexttogo=newtogo;
newtogo->prevtogo=tailtogo;
tailtogo=newtogo;
}
/* Now scan the microstructure and rank the sites */
for(px=0;px<SYSIZE;px++){
for(py=0;py<SYSIZE;py++){
for(pz=0;pz<SYSIZE;pz++){
if(mic[px][py][pz]==POROSITY){
cntpore=countbox(cubesize,px,py,pz);
if(cntpore>cntmax){cntmax=cntpore;}
/* Store this site value at appropriate place in */
/* sorted linked list */
if(cntpore>(tailtogo->npore)){
placed=0;
lasttogo=tailtogo;
while(placed==0){
newtogo=lasttogo->prevtogo;
if(newtogo==NULL){
placed=2;
}
else{
if(cntpore<=(newtogo->npore)){
placed=1;
}
}
if(placed==0){
lasttogo=newtogo;
}
}
onetogo=(struct togo *)malloc(sizeof(struct togo));
onetogo->x=px;
onetogo->y=py;
onetogo->z=pz;
onetogo->npore=cntpore;
/* Insertion at the head of the list */
if(placed==2){
onetogo->prevtogo=NULL;
onetogo->nexttogo=headtogo;
headtogo->prevtogo=onetogo;
headtogo=onetogo;
}
if(placed==1){
onetogo->nexttogo=lasttogo;
onetogo->prevtogo=newtogo;
lasttogo->prevtogo=onetogo;
newtogo->nexttogo=onetogo;
}
/* Eliminate the last element */
lasttogo=tailtogo;
tailtogo=tailtogo->prevtogo;
tailtogo->nexttogo=NULL;
free(lasttogo);
}
}
}
}
}
/* Now remove the sites */
/* starting at the head of the list */
/* and deallocate all of the used memory */
for(idesire=1;idesire<=ndesire;idesire++){
px=headtogo->x;
py=headtogo->y;
pz=headtogo->z;
if(px!=(-1)){
mic [px] [py] [pz]=EMPTYP;
count[POROSITY]-=1;
count[EMPTYP]+=1;
}
lasttogo=headtogo;
headtogo=headtogo->nexttogo;
free(lasttogo);
}
/* If only small cubes of porosity were found, then adjust */
/* cubesize to have a more efficient search in the future */
if(cubesize>CUBEMIN){
if((2*cntmax)<(cubesize*cubesize*cubesize)){
cubesize-=2;
}
}
}
/* routine to implement a cycle of dissolution */
/* Called by main program */
/* Calls passone, loccsh, and makeinert */
void dissolve(cycle)
int cycle;
{
int nc3aext,ncshext,nchext,nfh3ext,ngypext,plok,edgef;
int xpmax,ypmax,phid,phnew,plnew,cread;
int i,xloop,yloop,zloop,ngood,ix1,iy1,xc,yc,valid,xc1,yc1;
int iz1,zc,zc1;
long int water_left,ctest;
int placed,cshrand,ntrycsh,maxsulfate;
long int nsurf,suminit;
long int xext,nhgd;
float pdis,plfh3,fchext,fc3aext,mass_now,tot_mass,heatfill;
float heatsum,molesh2o,molesdh2o,h2oinit,alpha,heat2,heat3,heat4;
FILE *heatfile,*phfile;
struct ants *antadd;
/* Initialize variables */
nmade=0;
cshrand=0; /* counter for number of csh diffusing species to */
/* be located at random locations in microstructure */
heat_old=heat_new; /* new and old values for heat released */
/* Initialize dissolution and phase counters */
nsurf=0;
for(i=0;i<=EMPTYP;i++){
discount[i]=0;
count[i]=0;
}
/* Pass one- highlight all edge points which are soluble */
soluble[C3AH6]=0;
soluble[CH]=0;
passone(0,EMPTYP,cycle);
/* If first cycle, then determine all mixture proportions based */
/* on user input and original microstructure */
if(cycle==1){
/* Mass of cement in system */
cemmass=(specgrav[C3S]*(float)count[C3S]+specgrav[C2S]*
(float)count[C2S]+specgrav[C3A]*(float)count[C3A]+
specgrav[C4AF]*(float)count[C4AF]+specgrav[GYPSUM]*
(float)count[GYPSUM]);
/* Total mass in system neglecting single aggregate */
tot_mass=cemmass+(float)count[POROSITY]+specgrav[INERT]*
(float)count[INERT]+specgrav[POZZ]*(float)count[POZZ];
/* water-to-cement ratio */
w_to_c=(float)count[POROSITY]/cemmass;
mass_water=(1.-mass_agg)*(float)count[POROSITY]/tot_mass;
/* pozzolan-to-cement ratio */
s_to_c=(float)(count[INERT]*specgrav[INERT]+
count[POZZ]*specgrav[POZZ])/cemmass;
/* Conversion factor to kJ/kg for heat produced */
heatfill=(float)(count[INERT]+count[POZZ])/cemmass;
heat_cf=0.001*(0.3125+w_to_c+heatfill);
mass_fill=(1.-mass_agg)*(float)(count[INERT]*specgrav[INERT]+
count[POZZ]*specgrav[POZZ])/tot_mass;
printf("Calculated w/c is %.4f\n",w_to_c);
printf("Calculated s/c is %.4f \n",s_to_c);
printf("Calculated heat conversion factor is %f \n",heat_cf);
printf("Calculated mass fractions of water and filler are %.4f and %.4f \n",
mass_water,mass_fill);
}
heatsum=molesh2o=molesdh2o=0.0;
alpha=0.0; /* degree of hydration */
/* heat2 is heat of hydration based on heats of hydration of pure */
/* compounds (c3s,c2s,c3a,c4af) from Neville and Brooks */
/* page 14, Table 2.4 */
/* Originally from Lerch and Bogue, */
/* See book (The Chemistry of Cement) by Lea for complete ref. */
/* heat3 is heat of hydration based on heats of hydration of pure */
/* compounds (c3s,c2s,c3a,c4af) from Cement Chemistry by Taylor */
/* page 232, Table 7.5 */
/* heat4 contains measured heat release for C4AF hydration from */
/* Fukuhara et al., Cem. and Conc. Res. article */
heat2=heat3=heat4=0.0;
mass_now=0.0; /* total cement mass corrected for hydration */
suminit=c3sinit+c2sinit+c3ainit+c4afinit+ncsbar;
/* ctest is number of gypsum likely to form ettringite */
/* 1 unit of C3A can react with 2.5 units of Gypsum */
ctest=count[DIFFGYP];
if((float)ctest>(2.5*(float)count[DIFFC3A])){
ctest=2.5*(float)count[DIFFC3A];
}
for(i=0;i<=EMPTYP;i++){
if((i!=0)&&(i<=ABSGYP)&&(i!=INERTAGG)){
heatsum+=(float)count[i]*heatf[i]/molarv[i];
/* Tabulate moles of H2O consumed by reactions so far */
molesh2o+=(float)count[i]*waterc[i]/molarv[i];
}
/* assume that all C3A which can, does form ettringite */
if(i==DIFFC3A){
heatsum+=((float)count[i]-(float)ctest/2.5)*heatf[C3A]/molarv[C3A];
}
/* assume all gypsum which can, does form ettringite */
/* rest will remain as gypsum */
if(i==DIFFGYP){
heatsum+=(float)(count[i]-ctest)*heatf[GYPSUM]/molarv[GYPSUM];
heatsum+=(float)ctest*3.30*heatf[ETTR]/molarv[ETTR];
molesdh2o+=(float)ctest*3.30*waterc[ETTR]/molarv[ETTR];
}
else if(i==DIFFCH){
heatsum+=(float)count[i]*heatf[CH]/molarv[CH];
molesdh2o+=(float)count[i]*waterc[CH]/molarv[CH];
}
else if(i==DIFFCSH){
heatsum+=(float)count[i]*heatf[CSH]/molarv[CSH];
molesdh2o+=(float)count[i]*waterc[CSH]/molarv[CSH];
}
else if(i==DIFFETTR){
heatsum+=(float)count[i]*heatf[ETTR]/molarv[ETTR];
molesdh2o+=(float)count[i]*waterc[ETTR]/molarv[ETTR];
}
else if(i==C3S){
alpha+=(float)(c3sinit-count[i]);
mass_now+=specgrav[C3S]*(float)count[C3S];
heat2+=.502*(float)(c3sinit-count[i])*specgrav[C3S];
heat3+=.517*(float)(c3sinit-count[i])*specgrav[C3S];
heat4+=.517*(float)(c3sinit-count[i])*specgrav[C3S];
}
else if(i==C2S){
alpha+=(float)(c2sinit-count[i]);
mass_now+=specgrav[C2S]*(float)count[C2S];
heat2+=.260*(float)(c2sinit-count[i])*specgrav[C2S];
heat3+=.262*(float)(c2sinit-count[i])*specgrav[C2S];
heat4+=.262*(float)(c2sinit-count[i])*specgrav[C2S];
}
else if(i==C3A){
alpha+=(float)(c3ainit-count[i]);
mass_now+=specgrav[C3A]*(float)count[C3A];
heat2+=.867*(float)(c3ainit-count[i])*specgrav[C3A];
heat3+=1.144*(float)(c3ainit-count[i])*specgrav[C3A];
heat4+=1.144*(float)(c3ainit-count[i])*specgrav[C3A];
}
else if(i==C4AF){
alpha+=(float)(c4afinit-count[i]);
mass_now+=specgrav[C4AF]*(float)count[C4AF];
heat2+=.419*(float)(c4afinit-count[i])*specgrav[C4AF];
heat3+=.418*(float)(c4afinit-count[i])*specgrav[C4AF];
heat4+=.725*(float)(c4afinit-count[i])*specgrav[C4AF];
}
else if(i==GYPSUM){
alpha+=(float)(ncsbar-count[i]);
mass_now+=specgrav[GYPSUM]*(float)count[GYPSUM];
}
}
alpha=alpha/(float)suminit;
/* Current degree of hydration on a mass basis */
alpha_cur=1.0-(mass_now/cemmass);
h2oinit=(float)porinit/molarv[POROSITY];
/* Update water consumed to reflect that used in pozzolanic reaction */
/* 2.1 moles of water per mole of pozzolanic C-S-H */
/* assuming that molar volume of pozzolanic C-S-H is 80.87 */
molesh2o+=2.1*((float)(count[POZZ]-nfill)+(float)(npr)/1.35)/80.87;
/* Update heatsum for pozzolanic conversion */
heatsum+=((float)count[POZZ]-(float)nfill+(float)npr/1.35)*
(heatf[CSH]/80.87-heatf[POZZ]/molarv[POZZ]);
/* Assume 780 J/g S for pozzolanic reaction */
heat2+=0.78*((float)npr/1.35)*specgrav[POZZ];
heat3+=0.78*((float)npr/1.35)*specgrav[POZZ];
heat4+=0.78*((float)npr/1.35)*specgrav[POZZ];
/* Adjust heat sum for water left in system */
water_left=(long int)((h2oinit-molesh2o)*molarv[POROSITY]+0.5);
heatsum+=(h2oinit-molesh2o-molesdh2o)*heatf[POROSITY];
heatfile=fopen("heat.out","a");
if(cyccnt==0){
fprintf(heatfile,"Cycle alpha_vol alpha_mass heat1 heat2 heat3 heat4\n");
}
fprintf(heatfile,"%d %f %f %f %f %f %f\n",
cyccnt,alpha,alpha_cur,heatsum,heat2,heat3,heat4);
heat_new=heat4; /* use heat4 for all adiabatic calculations */
/* due to best agreement with calorimetry data */
fclose(heatfile);
if(molesh2o>h2oinit){
printf("All water consumed at cycle %d \n",cyccnt);
fflush(stdout);
exit();
}
/* Attempt to create empty porosity for account for self-desiccation */
if((sealed==1)&&((count[POROSITY]-water_left)>0)){
poretodo=count[POROSITY]-water_left;
makeinert(poretodo);
poregone+=poretodo;
}
/* Output phase counts */
phfile=fopen("phases.out","a");
fprintf(phfile,"%d ",cyccnt);
for(i=0;i<=EMPTYP;i++){
fprintf(phfile,"%ld ",count[i]);
printf("%ld ",count[i]);
}
printf("\n");
fprintf(phfile,"%ld\n",water_left);
fclose(phfile);
cyccnt+=1;
if(cycle==0){
return;
}
/* See if ettringite is soluble yet */
if(ncsbar!=0.0){
if((soluble[ETTR]==0)&&((count[AFM]!=0)||(((float)count[GYPSUM]/
((float)ncsbar+((float)netbar/3.21)))<0.10))){
soluble[ETTR]=1;
printf("Ettringite is soluble beginning at cycle %d \n",cycle);
/* identify all new soluble ettringite */
passone(ETTR,ETTR,2);
}
} /* end of soluble ettringite test */
/* Adjust ettringite solubility */
if(count[DIFFETTR]>DETTRMAX){
disprob[ETTR]=0.0;
}
else{
disprob[ETTR]=disbase[ETTR];
}
/* Adjust solubility of CH */
/* based on amount of CH currently diffusing */
/* Note that CH is always soluble to allow some */
/* Ostwald ripening of the CH crystals */
soluble[CH]=1;
passone(CH,CH,2);
if((float)count[DIFFCH]>=CHCRIT){
disprob[CH]=disbase[CH]*CHCRIT/(float)count[DIFFCH];
}
else{
disprob[CH]=disbase[CH];
}
/* Address solubility of C3AH6 */
if((count[GYPSUM]>(int)((float)ncsbar*0.05))||(count[ETTR]>500)){
soluble[C3AH6]=1;
passone(C3AH6,C3AH6,2);
/* Base C3AH6 solubility on maximum sulfate in solution */
/* from gypsum or ettringite available for dissolution */
/* The more the sulfate, the higher this solubility should be */
maxsulfate=count[DIFFGYP];
if((maxsulfate<count[DIFFETTR])&&(soluble[ETTR]==1)){
maxsulfate=count[DIFFETTR];
}
/* Adjust C3AH6 solubility based on potential gypsum which will dissolve */
if(maxsulfate<(int)((float)gypready*disprob[GYPSUM]*
(float)count[POROSITY]/(float)(SYSIZE*SYSIZE*SYSIZE))){
maxsulfate=(int)((float)gypready*disprob[GYPSUM]*
(float)count[POROSITY]/(float)(SYSIZE*SYSIZE*SYSIZE));
}
if(maxsulfate>0){
disprob[C3AH6]=disbase[C3AH6]*(float)maxsulfate/C3AH6CRIT;
if(disprob[C3AH6]>0.5){disprob[C3AH6]=0.5;}
}
else{
disprob[C3AH6]=disbase[C3AH6];
}
}
else{
soluble[C3AH6]=0;
}
/* See if silicates are soluble yet */
if((soluble[C3S]==0)&&((cycle>1)||(count[ETTR]>0)||(count[AFM]>0))){
soluble[C2S]=1;
soluble[C3S]=1;
/* identify all new soluble silicates */
passone(C3S,C2S,2);
} /* end of soluble silicate test */
/* Pass two- perform the dissolution of species */
nhgd=0;
/* Once again, scan all pixels in microstructure */
for(xloop=0;xloop<SYSIZE;xloop++){
for(yloop=0;yloop<SYSIZE;yloop++){
for(zloop=0;zloop<SYSIZE;zloop++){
if(mic[xloop][yloop][zloop]>OFFSET){
phid=mic[xloop][yloop][zloop]-OFFSET;
/* attempt a one-step random walk to dissolve */
plnew=(int)(6.*ran1(seed));
if((plnew<0)||(plnew>5)){ plnew=5;}
xc=xloop;
yc=yloop;
zc=zloop;
switch(plnew){
case 0:
xc-=1;
if(xc<0){xc=(SYSIZE-1);}
break;
case 1:
yc-=1;
if(yc<0){yc=(SYSIZE-1);}
break;
case 2:
xc+=1;
if(xc>=SYSIZE){xc=0;}
break;
case 3:
yc+=1;
if(yc>=SYSIZE){yc=0;}
break;
case 4:
zc-=1;
if(zc<0){zc=(SYSIZE-1);}
break;
case 5:
zc+=1;
if(zc>=SYSIZE){zc=0;}
break;
default:
break;
}
/* Generate probability for dissolution */
pdis=ran1(seed);
if((pdis<=disprob[phid])&&(mic[xc][yc][zc]==POROSITY)){
discount[phid]+=1;
cread=creates[phid];
mic[xloop][yloop][zloop]=POROSITY;
if(phid==C3AH6){nhgd+=1;}
/* Special dissolution for C4AF */
if(phid==C4AF){
plfh3=ran1(seed);
if((plfh3<0.0)||(plfh3>1.0)){
plfh3=1.0;
}
/* For every C4AF that dissolves, 0.5453 */
/* diffusing FH3 species should be created */
if(plfh3<=0.5453){
cread=DIFFFH3;
}
}
if(cread!=POROSITY){
nmade+=1;
ngoing+=1;
phnew=cread;
count[phnew]+=1;
mic[xc][yc][zc]=phnew;
antadd=(struct ants *)malloc(sizeof(struct ants));
antadd->x=xc;
antadd->y=yc;
antadd->z=zc;
antadd->id=phnew;
/* Now connect this ant structure to end of linked list */
antadd->prevant=tailant;
tailant->nextant=antadd;
antadd->nextant=NULL;
tailant=antadd;
}
if((phid==C3S)||(phid==C2S)){
plfh3=ran1(seed);
if((phid==C2S)||(plfh3<=0.521)){
placed=loccsh(xc,yc,zc);
if(placed!=0){
count[DIFFCSH]+=1;
}
else{
cshrand+=1;
}
}
}
if(phid==C2S){
plfh3=ran1(seed);
if(plfh3<=0.077){
placed=loccsh(xc,yc,zc);
if(placed!=0){
count[DIFFCSH]+=1;
}
else{
cshrand+=1;
}
}
}
}
else{
mic[xloop][yloop][zloop]-=OFFSET;
}
} /* end of if edge loop */
} /* end of zloop */
} /* end of yloop */
} /* end of xloop */
/* Now add in the extra diffusing species for dissolution */
/* Expansion factors from Young and Hansen and */
/* Mindess and Young (Concrete) */
ncshext=cshrand;
if(cshrand!=0){
printf("cshrand is %d \n",cshrand);
}
/* CH, Gypsum, and diffusing C3A are added at totally random */
/* locations as opposed to at the dissolution site */
fchext=0.61*(float)discount[C3S]+0.191*(float)discount[C2S]+
0.2584*(float)discount[C4AF];
nchext=fchext;
if(fchext>(float)nchext){
pdis=ran1(seed);
if((fchext-(float)nchext)>pdis){
nchext+=1;
}
}
fc3aext=discount[C3A]+0.5917*(float)discount[C3AH6]+
0.696*(float)discount[C4AF];
nc3aext=fc3aext;
if(fc3aext>(float)nc3aext){
pdis=ran1(seed);
if((fc3aext-(float)nc3aext)>pdis){
nc3aext+=1;
}
}
ngypext=discount[GYPSUM];
count[DIFFGYP]+=ngypext;
count[DIFFCH]+=nchext;
count[DIFFCSH]+=ncshext;
count[DIFFC3A]+=nc3aext;
for(xext=1;xext<=(nc3aext+ncshext+nchext+ngypext);xext++){
plok=0;
do{
xc=(int)((float)SYSIZE*ran1(seed));
yc=(int)((float)SYSIZE*ran1(seed));
zc=(int)((float)SYSIZE*ran1(seed));
if(xc>=SYSIZE){xc=0;}
if(yc>=SYSIZE){yc=0;}
if(zc>=SYSIZE){zc=0;}
if(mic[xc][yc][zc]==POROSITY){
plok=1;
phid=DIFFCH;
if(xext>nchext){phid=DIFFCSH;}
if(xext>(nchext+ncshext)){phid=DIFFC3A;}
if(xext>(nchext+ncshext+nc3aext)){phid=DIFFGYP;}
mic[xc][yc][zc]=phid;
nmade+=1;
ngoing+=1;
antadd=(struct ants *)malloc(sizeof(struct ants));
antadd->x=xc;
antadd->y=yc;
antadd->z=zc;
antadd->id=phid;
/* Now connect this ant structure to end of linked list */
antadd->prevant=tailant;
tailant->nextant=antadd;
antadd->nextant=NULL;
tailant=antadd;
}
} while (plok==0);
} /* end of xext for extra species generation */
printf("Dissolved- %ld %ld %ld %ld %ld %ld \n",count[DIFFCSH],
count[DIFFCH],count[DIFFGYP],count[DIFFC3A],count[DIFFFH3],
count[DIFFETTR]);
printf("C3AH6 dissolved- %ld with prob. of %f \n",nhgd,disprob[C3AH6]);
fflush(stdout);
}
/* routine to add nneed one pixel elements of phase randid at random */
/* locations in microstructure */
/* Called by main program */
/* Calls no other routines */
void addrand(randid,nneed)
int randid;
long int nneed;
{
int ix,iy,iz;
long int ic;
int success;
/* Add number of requested phase at random pore locations in system */
for(ic=1;ic<=nneed;ic++){
success=0;
while(success==0){
ix=(int)((float)SYSIZE*ran1(seed));
iy=(int)((float)SYSIZE*ran1(seed));
iz=(int)((float)SYSIZE*ran1(seed));
if(ix==SYSIZE){ix=0;}
if(iy==SYSIZE){iy=0;}
if(iz==SYSIZE){iz=0;}
if(mic[ix][iy][iz]==POROSITY){
mic[ix][iy][iz]=randid;
success=1;
}
}
}
}
/* Calls dissolve and addrand */
main()
{
int ncyc,ntimes,icyc,valin;
int cycflag,ix,iy,iz,phtodo,setflag;
int iseed,burnfreq,setfreq,oflag,adiaflag;
long int nadd;
int xpl,xph,ypl,yph,fidc3s,fidc2s,fidc3a,fidc4af,fidgyp,fidagg;
float pnucch,pscalech,pnuchg,pscalehg,pnucfh3,pscalefh3,krate;
float mass_cement,mass_cem_now,mass_cur;
FILE *infile,*outfile,*adiafile;
char filei[80],fileo[80];
ngoing=0;
cycflag=0;
heat_old=heat_new=0.0;
time_cur=0.0; /* Elapsed time according to maturity principles */
cubesize=CUBEMAX;
poregone=poretodo=0;
/* Get random number seed */
printf("Enter random number seed \n");
scanf("%d",&iseed);
printf("%d\n",iseed);
seed=(&iseed);
printf("Dissoluton bias is set at %f \n",DISBIAS);
printf("Do you wish to output final microstructure to file (0-No 1-Yes)\n");
scanf("%d",&oflag);
printf("%d\n",oflag);
if(oflag==1){
printf("Enter name of file to create \n");
scanf("%s",fileo);
printf("%s\n",fileo);
}
/* Open file and read in original cement particle microstructure */
printf("Enter name of file to read initial microstructure from \n");
scanf("%s",filei);
printf("%s\n",filei);
/* Get phase assignments for original microstructure */
/* to transform to needed ID values */
printf("Enter IDs in file for C3S, C2S, C3A, C4AF, Gypsum, Aggregate\n");
scanf("%d %d %d %d %d %d",&fidc3s,&fidc2s,&fidc3a,&fidc4af,&fidgyp,&fidagg);
printf("%d %d %d %d %d %d\n",fidc3s,fidc2s,fidc3a,fidc4af,fidgyp,fidagg);
fflush(stdout);
infile=fopen(filei,"r");
for(ix=0;ix<SYSIZE;ix++){
for(iy=0;iy<SYSIZE;iy++){
for(iz=0;iz<SYSIZE;iz++){
fscanf(infile,"%d",&valin);
mic[ix][iy][iz]=POROSITY;
if(valin==fidc3s){
mic[ix][iy][iz]=C3S;
}
else if(valin==fidc2s){
mic[ix][iy][iz]=C2S;
}
else if(valin==fidc3a){
mic[ix][iy][iz]=C3A;
}
else if(valin==fidc4af){
mic[ix][iy][iz]=C4AF;
}
else if(valin==fidgyp){
mic[ix][iy][iz]=GYPSUM;
}
else if(valin==fidagg){
mic[ix][iy][iz]=INERTAGG;
}
}
}
}
fclose(infile);
fflush(stdout);
/* Now read in particle IDs from file */
printf("Enter name of file to read particle IDs from \n");
scanf("%s",filei);
printf("%s\n",filei);
infile=fopen(filei,"r");
for(ix=0;ix<SYSIZE;ix++){
for(iy=0;iy<SYSIZE;iy++){
for(iz=0;iz<SYSIZE;iz++){
fscanf(infile,"%d",&valin);
micpart[ix][iy][iz]=valin;
}
}
}
fclose(infile);
fflush(stdout);
/* Initialize counters, etc. */
npr=0;
nfill=0;
ncsbar=0;
netbar=0;
porinit=0;
cyccnt=0;
setflag=0;
c3sinit=c2sinit=c3ainit=c4afinit=0;
/* Initialize structure for ants */
headant=(struct ants *)malloc(sizeof(struct ants));
headant->prevant=NULL;
headant->nextant=NULL;
headant->x=0;
headant->y=0;
headant->z=0;
headant->id=100; /* special ID indicating first ant in list */
tailant=headant;
/* Allow user to iteratively add one pixel particles of various phases */
/* Typical application would be for addition of silica fume */
printf("Enter number of one pixel particles to add (0 to quit) \n");
scanf("%ld",&nadd);
printf("%ld\n",nadd);
while(nadd>0){
phtodo=25;
while((phtodo<0)||(phtodo>INERT)){
printf("Enter phase to add \n");
printf(" C3S 1 \n");
printf(" C2S 2 \n");
printf(" CH 3 \n");
printf(" CSH 4 \n");
printf(" C3A 5 \n");
printf(" GYPSUM 6 \n");
printf(" C4AF 7 \n");
printf(" C3AH6 8 \n");
printf(" ETTR 9 \n");
printf(" AFM 10 \n");
printf(" FH3 11 \n");
printf(" POZZ 12 \n");
printf(" INERT 13 \n");
scanf("%d",&phtodo);
printf("%d \n",phtodo);
}
addrand(phtodo,nadd);
printf("Enter number of one pixel particles to add (0 to quit) \n");
scanf("%ld",&nadd);
printf("%ld\n",nadd);
}
fflush(stdout);
init();
printf("After init routine \n");
printf("Enter number of cycles to execute \n");
scanf("%d",&ncyc);
printf("%d \n",ncyc);
printf("Do you wish hydration under 0) saturated or 1) sealed conditions \n");
scanf("%d",&sealed);
printf("%d \n",sealed);
printf("Enter max. # of diffusion steps per cycle (500) \n");
scanf("%d",&ntimes);
printf("%d \n",ntimes);
printf("Enter nuc. prob. and scale factor for CH nucleation \n");
scanf("%f %f",&pnucch,&pscalech);
printf("%f %f \n",pnucch,pscalech);
printf("Enter nuc. prob. and scale factor for C3AH6 nucleation \n");
scanf("%f %f",&pnuchg,&pscalehg);
printf("%f %f \n",pnuchg,pscalehg);
printf("Enter nuc. prob. and scale factor for FH3 nucleation \n");
scanf("%f %f",&pnucfh3,&pscalefh3);
printf("%f %f \n",pnucfh3,pscalefh3);
printf("Enter cycle frequency for checking pore space percolation \n");
scanf("%d",&burnfreq);
printf("%d\n",burnfreq);
printf("Enter cycle frequency for checking percolation of solids (set) \n");
scanf("%d",&setfreq);
printf("%d\n",setfreq);
/* Parameters for adiabatic temperature rise calculation */
printf("Enter the induction time in hours \n");
scanf("%f",&ind_time);
printf("%f \n",ind_time);
time_cur+=ind_time;
printf("Enter the initial temperature in degrees Celsius \n");
scanf("%f",&temp_0);
printf("%f \n",temp_0);
temp_cur=temp_0;
printf("Enter apparent activation energy for hydration in kJ/mole \n");
scanf("%f",&E_act);
printf("%f \n",E_act);
printf("Enter kinetic factor to convert cycles to time for 25 C \n");
scanf("%f",&beta);
printf("%f \n",beta);
printf("Enter mass fraction of aggregate in concrete \n");
scanf("%f",&mass_agg);
printf("%f \n",mass_agg);
printf("Hydration under 0) isothermal or 1) adiabatic conditions \n");
scanf("%d",&adiaflag);
printf("%d \n",adiaflag);
fflush(stdout);
adiafile=fopen("adiabatic.out","w");
fprintf(adiafile,"Time Temperature Alpha Krate Cp_now Mass_cem\n");
for(icyc=1;icyc<=ncyc;icyc++){
if(icyc==ncyc){cycflag=1;}
dissolve(icyc);
printf("Number dissolved this pass- %ld total diffusing- %ld \n",nmade,ngoing);
fflush(stdout);
if(icyc==1){
printf("ncsbar is %ld netbar is %ld \n",ncsbar,netbar);
}
hydrate(cycflag,ntimes,pnucch,pscalech,pnuchg,pscalehg,pnucfh3,pscalefh3);
temp_0=temp_cur;
/* Handle adiabatic case first */
/* Cement + aggregate +water + filler=1; that's all there is */
mass_cement=1.-mass_agg-mass_fill-mass_water;
if(adiaflag==1){
/* determine heat capacity of current mixture, */
/* accounting for imbibed water if necessary */
if(sealed==1){
Cp_now=mass_agg*Cp_agg;
Cp_now+=Cp_pozz*mass_fill;
Cp_now+=Cp_cement*mass_cement;
Cp_now+=(Cp_h2o*mass_water-alpha_cur*WN*mass_cement*(Cp_h2o-Cp_bh2o));
mass_cem_now=mass_cement;
}
/* Else need to account for extra capillary water drawn in */
/* Basis is WCHSH(0.06) g H2O per gram cement for chemical shrinkage */
/* Need to adjust mass basis to account for extra imbibed H2O */
else{
mass_cur=1.+WCHSH*mass_cement*alpha_cur;
Cp_now=mass_agg*Cp_agg/mass_cur;
Cp_now+=Cp_pozz*mass_fill/mass_cur;
Cp_now+=Cp_cement*mass_cement/mass_cur;
Cp_now+=(Cp_h2o*mass_water-alpha_cur*WN*mass_cement*(Cp_h2o-Cp_bh2o));
Cp_now+=(WCHSH*Cp_h2o*alpha_cur*mass_cement);
mass_cem_now=mass_cement/mass_cur;
}
/* Determine rate constant based on Arrhenius expression */
/* Recall that temp_cur is in degrees Celsius */
krate=exp(-(1000.*E_act/8.314)*((1./(temp_cur+273.15))-(1./298.15)));
/* Update temperature based on heat generated and current Cp */
temp_cur=temp_0+mass_cem_now*heat_cf*(heat_new-heat_old)/Cp_now;
}
/* Else isothermal conditions */
else{
krate=1.0;
mass_cem_now=mass_cement;
}
/* Update time based on simple numerical integration */
/* simulating maturity approach */
/* with parabolic kinetics (Knudsen model) */
if(cyccnt>1){time_cur+=(2.*(float)(cyccnt-1)-1.0)*beta/krate;}
fprintf(adiafile,"%f %f %f %f %f %f\n",time_cur,temp_cur,
alpha_cur,krate,Cp_now,mass_cem_now);
fflush(adiafile);
/* Check percolation of pore space */
if((icyc%burnfreq)==0){
burn3d(0,1,0,0);
burn3d(0,0,1,0);
burn3d(0,0,0,1);
}
/* Check percolation of solids (set point) */
if((icyc%setfreq)==0){
setflag=burnset(1,0,0);
setflag+=burnset(0,1,0);
setflag+=burnset(0,0,1);
}
}
fclose(adiafile);
dissolve(0);
/* Output final microstructure if desired */
if(oflag==1){
outfile=fopen(fileo,"w");
for(ix=0;ix<SYSIZE;ix++){
for(iy=0;iy<SYSIZE;iy++){
for(iz=0;iz<SYSIZE;iz++){
fprintf(outfile,"%d\n",(int)mic[ix][iy][iz]);
}
}
}
fclose(outfile);
}
}