4SIGHT Manual: A Computer Program 

Modelling Degradation of Underground Low L

Waste Concrete Vaults

Kenneth A. Snyder

James R. Clifton

June 1995

Building and Fire Research Laboratory

National Institute of Standards and Technology

Gaithersburg, MD 20899

U.S. Depiutment of Commeme

Ronald H. Brown, ,Secrctqy

Technology Administmtion

Mary L. Good, Under Secrekny for Technology

National Institute of Standards and Technology

Arati Prabhakar, Director
is used to propagate ions by precipitation/dissolution of available salts. The pre

tion/dissolution of salts, in turn, changes the transport properties, which changes th

of ion transport. The result is a model which incorporates the synergism of multiple 

dation mechanisms.

This Report is self-contained. It includes the installation instructions, user m

technical details, and source code. The program was written using a literate program

tool and the “pretty-printing” output of the source code is attached at the end of this 

Keywords: building technology; computer modelling; concrete; corrosion of reinforce

degradation; leaching; radioactive waste; service life; sulfate attack

...

111
3

4

5

6

7

8

2.1 Distributed Files . . . . . . . . . ,

2.2 Hardware Requirements . . . . . .

Program Execution

3.1 Interactive Mode . . . . . . . . . .

3.2 Batch Mode . . . . . . . . . . . . .

3.3 Plotting Results . . . . . . . . . . .

Input Specifications

4.1 Material Properties . . . . . . . . .

4.2 Ion/Fluid Initial Conditions . . . .

4.3 Geometry . . . . . . . . . . . . . . .

4.4 Failure/Termination Limits . . . .

4.5 Sulfate Attack Parameters

4.6 Output Parameters . . . .

Example

5.1 Internal Parameters. . . .

5.2 Initial State . . . . . . . .

5.3 Depth vs. Time . . . . . .

5.4 Failure Data . . . . . . . .

5.5 Final System State . . . .

System

Ion

7.1

7.2

7.3

7.4

7.5

7.6

Transport

Advection-Diffusion . . . .

Darcy Flow . . . . . . . .

Continuity . . . . . . . . .

Dimensionless Variables .

Atkinson-Hearne Model .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

Effects of Cracks and Joints . . . .

Chemical Equilibrium

.

.

.

.

.

.

iv

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . .

. . .

. . .

. . .

. . .

. . .

. . .

. . .

. . .

. . .

. . .

. . .

. . .

. . .

. . .

. . .

. . .

. . .

. . .

. . .

. . .

. . .

. .

,.

. .

. .

. .

. .

. .

. .

. .

. .

. .

. .

. .

. .

. .

,.

. .

. .

. .

. .

. .

. .

. ..’

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

.

.

.

.

.,.

,..

. . .

. . .

. . .

. . .

. . .

. . .

. . .

. . .

. . .

. . .

. . .

. . .

. . .

. . .

.,.

. . .

. . .

. . .

. . .

. . .

.

.

.

.

.

.

.

.

. 

. 

. 

. 

. 

. 

. 

. 

. 

. 

. 

,..

,..

,..

. 

,..

,..

,.,

. 

,..

,..

,..
10 Calculated Material Properties

10.1 Hydration . . . .

10.2 Leaching . . . . .

11 Program Outline

12 Assumptions

13 Acknowledgements

14 References

A

B

Input Parsing

A.1 LEX Specification

A.2 YACCSpecification

CWEB- Source Code

. . . .

. . . .

,.. .

,.. .

.

.

.

.

. . . .

. . . .

,.. .

. . . .

. . . .

. . . .

,.. .

,.. .

v

. . . . .

. . . . .

. . . . .

. . . . .

. . . .

. . . .

. . . .

. . . .

. .

. .

. .

. .

. . . . . . . 

. . . . . . . 

. . . . . . . 

. . . . . . . 
([])ppp

by 4SIGHT for the leaching formation factor. . . . . . . . . . . . . . . . . 

3 4SIGHTflow chart . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ..
of the concrete vault roof given the material properties of the concrete and the extern

concentrations.

4SIGHT incorporates multiple degradation mechanisms by using a single transport 

tion for ions. Individual degradation mechanisms are typically controlled by the conc

tion of a single ion species. However, as various ion species diffuse into the concrete

transport properties change, which changes the rate of ion diffusion. As the pore so

pH changes, available salts either precipitate or go into solution. The precipitation 

dissolution processes changes the porosity, which in turn changes the transport prope

It is the act of maintaining chemical equilibrium that affects the synergism of the di

degradation mechanisms.

Previous efforts [1][2] have yielded mathematical models for various degradation 

anisms. From these equations, the service life could be calculated for any one degrad

mechanism. However, combining the effects of multiple mechanisms was not feasible 

lytically. The advantage of the method used by 4SIGHT is the ability to incorporate 

synergism of multiple degradation processes. Rather than calculate rates of attack, 4

uses a single transport equation to propagate ions through the concrete. Keeping the s

in chemical equilibrium allows interaction of degradation mechanisms.

4SIGHT is a command-driven program that awaits input parameters from the user, 

default values when necessary. Once all known parameters are entered, 4SIGHT reit

the values of all the input parameters used in the calculation with a note as to whethe

value was specified by the user or whether a default value was used. The computations 

by 4SIGHT propagating ions in discrete time increments. At regular intervals 4SIGHT 

to the screen the bulk hydraulic conductivity and diffusivity, and the depth of sulfat

chloride attack. The calculation continues until one of the following occurs: corrosion 

steel reinforcement, sufficient sulfate attack that the roof fails structurally, failure of a 

or a user specified time limit has been exceeded. Upon termination, 4SIGHT prints the 

of the system, composed of the concentration of user selected ions and moles of user se

salts, as a function of depth into the roof.

What follows in this Report is a complete user’s guide to 4SIGHT, including install

information and source code. The goal is an all-inclusive document, acting as a single 

of work, that facilitates review of any aspect of the program. To facilitate source 

review, 4SIGHT was written using a literate programming tool which incorporates bo

programming code and typeset comments. The typeset output does not conform 

structure of a NISTIR. Therefore, the typeset output, as it normally appears, is attach

the end of this report.

1
README Basic installation instructions

4SIGHT . EXE Executable program

ION. DB Ion database

MANUAL. PS 4SIGHT manual in PostScript format

MANUAL.HP 4SIGHT manual in HP-PCL format

MANUAL. TXT 4SIGHT manual in ASCII format

EXAMPLE.DAT Example input file

To install 4SIGHT, simply copy these files to a directory on the hard drive. This manu

been included in three different useful formats for printing.

2.2 Hardware Requirements

To execute the program 4SIGHT, the user must be using a personal computer equipped 

either an 80386 or 80486 microprocessor with a math coprocessor. Lacking either o

requirements, the program will fail to run.

3 Program Execution

4SIGHT was written to be executed in either interactive or batch mode. The inter

mode simply waits for the user to type in the various input parameters at run-time. 

executed in batch mode, the user first edits an ASCII file containing the input param

as if typed during interactive mode. This ASCII file is then included at the comman

when running 4SIGHT.

3.1 Interactive Mode

When run interactively, 4SIGHT simply waits for the user to enter input parameters. 

user specifies the end of the list by entering either quit or exit. For example, assumin

executable program is in the directory C: \4SIGHT, interactive mode is initiated by 

typing the program name:

C:\4SIGHT\4SIGHT

The program responds with:

2
DIFF = 2. OE

12

Wc = 0.45

TIME = 100000

quit

After which, 4SIGHT begins the calculation, printing intermediate results to the scree

3.2 Batch Mode

When run in batch mode the user specifies an input file at the command line. This inp

simply contains any number of input parameters, as in the interactive mode. For ex

assume the following text is saved into a file called input. dat:

DIFF = 5. OE-12

THICKNESS = 1.0

Wc = 0.50

EXTERNALCl = 0.150

EXTERNALNa = 0.150

DEPTH = 0.25

This file can be created by any ASCII editor such as the DOS edit program or the 

vi program. When using a commercial word processing program to create input file

sure that the output is in plain ASCII. As a check, give the C: \type input. dat com

at the DOS prompt. The text should appear as typed. Note that neither quit nor ex

required in the input file. To use input files with 4S IGHT, simply include that file na

the DOS command line:

c: \4s1GHT\4S1GHT input. dat

The output from 4SIGHT goes to the screen. This is useful for initial calculations an

runs. However, if the user desires to save the output to an ASCII file, simply use redir

of the output:

c: \4SIGHT\4S1GHT input. dat >Output. dat

This example uses the commands in the file input. dat as input parameters and sa

output in the file output. dat.

3
mediate results, or the ion concentration as a function of depth from the final state 

system.

4 Input Specifications

The complete list of input parameters is given here to introduce the parameters used 

example to follow. Although many of the parameters are self-evident, the complete m

of all the parameters will not be understood until the user reads the explanatory s

that follow the example.

The input commands to

format of the commands is

Even though every possible

4SIGHT are straightforward.

parameter = value

With very few exceptions

input parameter is listed below, there are sufficient 

parameters for the user to perform a calculation without entering any parameters.

The syntax of the list of parameters below is as follows: Text to be typed verba

shown in typewriter font (e.g., DIFF = ), These words can be typed in any mixture of 

and lower case characters with any number of spaces or tabs between the parameters 

the ‘=’. Words in italics represent user selected input values (e.g., ezpr). There are tw

types: numbers and ions, which are represented by expr and ion, respectively. Nu

can be given in integer, fixed point, or scientific notation. Examples of valid values f

include:

120 120.0

The units for each input value appear under ezpr

the following ions (the valence has been omitted):

H Ca Na K OH

The ion values MUST appear as shown since ion

1. 2E+02

in the definition. 4SIGHT also reco

cl S04 C03

variables are case sensitive. 4S IG

also recognize salts. A salt is specified by concatenating a cation and an anion, sep

by a space. Stoichiometric ratios are not used. For example, the following specificatio

Na Cl Ca OH Na S04

4
4.1 Material Properties

DIFF = ezpr

m2/sec

expr:

default:

PERM = expr

m2

expr:

default:

WC = expr

expr:

default:

The ultimate Cl- diffusivity of the concrete. This value represents 

diffusivity of the undamaged concrete. The default value given below 

suitable for a w/c=O.45 ordinary portland concrete.

DIFF = 5.7E-12

The ultimate permeability of the concrete. This value represents 

permeability of the undamaged concrete.

PERM = 2.5E-18

Concrete water: cement weight ratio during mixing.

Wc = 0.45

4.2 Ion/Fluid Initial Conditions

EXTERNALion = expr

mol/L

ion: The ion of interest.

expr: The external concentration of ion above the roof. This command 

be repeated for each ion external to the roof slab.

example: EXTERNALCl = 0.15

INTERNALion = expr

mol/L
example:

HEAD = expr

m

expr:

default:

INTERNAL Na = 0.10

The equivalent height of ground water (hydraulic head) above th

surface of the roof.

HEAD = 5.00

4.3 Geometry.

THICKNESS = ezpr

m

expr: The vertical thickness of the roof.

default: THICKNESS = 1.00

CRACK= ezprl AT ezprl? DEPTH expr3

m

exprl:

expr2:

expr3:

default:

JOINT PERM =

expr:

example:

m m

The width of the crack.

The spacing between cracks.

The penetration of the crack, relative to the bottom of the slab

CRACK= 0.000100 AT 1.00 DEPTH 0.30

expr

.

m~

Thepermeability of thejoint filling compound between roof slab

value should not reflect the geometry of joint. Rather, it is a 

property of the joint compound.

JOINT PERM = 3. 2E-15

JOINT = exprl AT ezpr2 UNTIL expr3

m m years

6
4.4 Failure/Termination Limits

REBAR = expr

m

expr: The depth of the reinforcement bars (rebars) from the topof the 

When a critical concentration of C/- ions has penetrated to a 

expr the structure fails.

example: REBAR = 0.90

DEPTH = expr

m

expr: Thecritical penetration depth for sulfate attack. Once the sulfate 

has penetrated down to a depth expr from the top, the roof is una

default:

sustain its load and it fails.

DEPTH = 0.20

TIME = expr

day

expr: The maximum number ofdays to continue the calculation. Not

default:

no comma should be used.

TIME = 100000

4.5 Sulfate Attack Parameters

YOUNGS= expr

GPa

expr: The Youngs modulus of the concrete.

default: YOUNGS= 20.0

BETA = expr

expr: Thelinear strain duetoonemoleof Na+ perm30f concrete.

default: BETA = 1.8E-06

7
ROUGHNESS=

ezpr:

default:

GAMMA= expr

expr

The fracture surface roughness factor.

ROUGHNESS= 1.0

J/m2

expr: The fracture surface energy of the concrete.

dejault: GAMMA= 10.0

POISSON = expr

expr: The Poisson ratio of the concrete.

default: POISSON = 0.2

4.6 Output Parameters

OUTPUT= i07?

ion: Output the pore solution concentration (mol/L) of ion in the 

describing the final state of the system. This command lets th

examine the concentrationof specified ions asa functionof depth

command can be repeated for different ions.

example: OUTPUT= C03

OUTPUT= ionl hzl?

Output the quantity (moles) of salt in the volume Vsemple. The 

comprised of the cation ionl and the anion ion2.

NOTE: The two ion specifications are separated bya space and stoichiom

8
using Clions gave 5.0 x 1012 rn2/.sec (DIFF5. OE12). No permeability measure

are available. The roof slab is 1.0 meter thick (THICKNESS=1. 000). Regularly spaced 

have been observed on the bottom of the slab. The cracks are approximately 100 pm 

spaced 2 meters apart, and are assumed to extend upwards to the neutral axis w

25 centimeters from the bottom of the roof slab (CRACK= 0.000100 AT 2.0 DEPTH 

Also, the roof is buried, giving an effective pressure head of 2 meters (HEAD = 2. O)

Soil analysis indicates the presence of SO:- at a concentration of 1.0 moles p

(EXTERNALS04 = 1. O). Engineering analysis indicates that if the sulfate degradation 

etrates down 20 centimeters from the top surface of the roof then the vault will c

(DEPTH = O. 20). Additionally, chloride ions are present in the soil at a concentration 

moles per liter (EXTERNALCl = O. 40) and engineering drawings indicate that the re

ment bars are located 60 centimeters from the top of the roof (REBAR = O. 60). Th

soil analysis indicates that sodium ions are also present at a concentration of 2.40 mo

liter (EXTERNALNa = 2.40), giving a nearly neutral soil pH.

Internal to the concrete the pH is approximately 13. Therefore, internal potassium 

sodium concentrations are approximately 0.1 (INTERNAL K = O. 1) and 0.05 (INTER

= O. 05) moles per liter, respectively.

To monitor the ingress of the external ions, calcium (OUTPUTCa) and chloride (

Cl) ions will be included in the output of the final state of the system. Also, to m

leaching, the solid calcium hydroxide content will also be included in the output (OUT

OH).

The input file for this example, example. dat is included in the distribution disket

reiterated here:

DIFF = 5. OE-12

Wc = .40

THICKNESS = 1.000

EXTERNALNa = 2.40

EXTERNALCl = 0.400

EXTERNALS04 = 1.00

INTERNAL K = 0.1000

INTERNAL Na= 0.0500

OUTPUTCa

OUTPUTCl

OUTPUTCa OH

REBAR = .8000
To use this input file, simply include example.dat at the DOS command line. To 

results, redirect the output to afile:

C:\4SIGHT\4SIGHT EXAMPLE.DAT>EXAMPLE.OUT

In this example, the output is stored in EXAMPLE.OUT,an ASCII file which can be im

into any spreadsheet program.

NOTE: This calculation required 90 seconds to complete on a personal computer eq

with a80486 microprocessor operating at66MHz.

5.1 Internal Parameters.

4SIGHToutputs all of theparameters theuser has specified (“USER’’),along with defa

calculated values (“DEFAULT”).

Ii This is 4SIGHT (Version 1.0)”

“THICKNESS “ 1.00000 “(m) “ “USER”

“DIFF tl 5. Oe-12 “(mA2/see)” “USER”

“PERM II 9.8e-13 “(m/see) “ “DEFAULT”

“WC II 0.40000 “

II I,USERII

“HEAD II 2.00000 “(m) “ “USER”

Sulfate Attack Parameters:

“YOUNGS “ 2.Oe+10 “(N/m-2) “ “DEFAULT”

“BETA 81 1.8e-06 “ “ “DEFAULT”

“CE “ 350.00000 “(Mel/m-3)” “DEFAULT”

“ROUGHNESS“ 1.00000 “ “ “DEFAULT”

“GAMMA “ 10.00000 “(J/m-2) “ “DEFAULT”

“POISSON “ 0.20000 “ “ “DEFAULT”

“DEPTH II 0.20000 “(m) “ “USER”

“REBAR II 0.80000 “(m) ,, I,USERII

“TIME II 100000 “(day) “ “USER”

“CRACK = 0.00010 AT 2.00000 DEPTH 0.25000”

“Chloride failure (yr)” 1743

“Sulfate failure (yr)” 447

10
y

“ ION “ “EXTERNAL” “INTERNAL”

tl H: “ 0.00000 0.00000

“ Ca:” 0.00000 0.00035

“ Na:” 2.40000 0.05000

*I K:” 0.00000 0.10000

“ OH:” 0.00000 0.15070

“ cl:” 0.40000 0.00000

“ S04:” 1.00000 0.00000

“ pH:” 7.00000 13.17810

These results indicate that the pH of the environment is 7, while thepH of th

solution isinitially 13. To change the pHof the environment the user can simply chan

concentration of EXTERNALanions or cations.

5.3 Depth vs. Time

As 4SIGHTis calculating ion transport it regularly prints the current status of the 

dation. The sulfate and chloride penetration depths as afunctionof time are given 

columns labelled S04and Cl, respectively. Listhe remaining thickness of the slab, K

hydraulic conductivity, D is the diffusivity, Flux is the flux of pore solution out the b

of the slab, and pHis thepH of the flux.

11
15055

20023

25031

30018

35011

40048

45104

50017

55045

60015

65003

70002

75013

80069

85037

90050

95069

100017

0.982 1.3e-12 5.Oe-12 0.018 0.129 0.027

0.976 1.3e-12 5.Oe-12 0.025 0.150 0.028

0.969 1.3e-12 5.Oe-12 0.031 0.173 0.028

0.963 1.3e-12 5.Oe-12 0.037 0.192 0.028

0.957 1.3e-125.Oe-12 0.043 0.211 0.028

0.951 1.3e-12 5.Oe-12 0.049 0.229 0.029

0.945 1.3e-12 5.Oe-12 0.055 0.246 0.029

0.939 1.3e-12 5.Oe-12 0.061 0.263 0.029

0.933 1.3e-12 5.Oe-12 0.067 0.280 0.029

0.927 1.3e-12 5.Oe-12 0.074 0.295 0.030

0.920 1.3e-12 5.Oe-12 0.080 0.311 0.030

0.914 1.3e-12 5.Oe-12 0.086 0.327 0.030

0.908 1.3e-12 5.Oe-12 0.092 0.342 0.031

0.902 1.4e-12 5.Oe-12 0.098 0.357 0.031

0.896 1.4e-12 5.Oe-12 0.104 0.372 0.031

0.890 1.4e-12 5.Oe-12 0.110 0.387 0.031

0.884 1.4e-12 5.Oe-12 0.116 0.401 0.031

0.878 1.4e-12 5.Oe-12 0.122 0.416 0.032

13.2

13.2

13.2

13.2

13.2

13.2

13.2

13.2

13.2

13.2

13.2

13.2

13.2

13.2

13.2

13.2

13.2

13.2

Asper the linear model for sulfate attack, the sulfate front increases linearly with 

The chloride depth has a tl/2 because of the low Peclet number (the ratio of D

‘diffusive’ flow [2]), At sufficiently high Peclet numbers (greater permeability or hy

head) the chloride depth approaches a linear relationship totime.

5.4 Failure Data

Afterreporting the time dependent behavior, 4SIGHTreports the reason fortermination 

thestatus of the sulfate and chloride penetration:

Exceeded TIME limit.

II II T 17.283

“Day” 100017

“S04 (m)” 0.122

“Cl (m)” 0.416

12
5.5 Final System State

Upon termination of the calculation, the final state of the system for the slab is printed. 

L(m) is measured from the top of the slab. Psi is the hydraulic pressure, VDis the D

velocity, xi is the inverse of the formation factor, phi is the porosity, fc is the estim

compressive strength using ACI 211.

Final System state:

“L(m)” “Psi” “vD” “xi” “phi” “pH” “fc” “Ca” “cl” llca(JHll

0.0000 0.392 0.000 0.0024 0.1634 7.000 5052 0.0000 0.4000 0.000

0.0500 0.392 0.000 0.0024 0.1631 7.000 5057 0.0000 0.4000 0.000

0.1000 0.392 0.031 0.0024 0.1658 7.000 5011 0.0000 0.4000 0.000

0.1500 0.374 0.031 0.0024 0.1663 13.277 5003 0.0001 0.3740 31.599

0.2000 0.342 0.031 0.0025 0.1671 13.351 4989 0.0001 0.3253 31.860

0.2500 0.311 0.031 0.0025 0.1672 13.3924988 0.0001 0.277331.886

0.3000 0.279 0.031 0.0025 0.1672 13.4054987 0.0001 0.231431.899

0.3500 0.248 0.031 0.0025 0.1673 13.3944986 0.0001 0.1889 31.912

0.4000 0.217 0.031 0.0025 0.1673 13.365 4986 0.0001 0.1509 31.918

0,4500 0.186 0.031 0.0025 0.1673 13.3264986 0.0002 0.117831.919

0.5000 0.155 0.031 0.0025 0.1673 13.286 4986 0.0002 0.0899 31.919

0.5500 0.124 0.031 0.0025 0.1673 13.250 4986 0.0002 0.0671 31.919

0.6000 0.093 0.031 0.0025 0.1673 13.220 4986 0.0003 0.0488 31.919

0.6500 0.062 0.031 0.0025 0.1673 13.199 4986 0.0003 0.034831.919

0.7000 0.031 0.031 0.0025 0.1673 13.1854986 0.0003 0.024231.919

0.7500 0.000 0.031 0.0025 0.1673 13.1764986 0.0004 0.016431.919

0.8000 0.000 0.031 0.0025 0.1673 13.1734986 0.0004 0.0110 31.919

0.8500 0.000 0.031 0.0025 0.1673 13.1724986 0.0004 0.007331.919

0.9000 0.000 0.031 0.0025 0.1673 13.172 4986 0.0004 0.0050 31.919

0.9500 0.000 0.031 0.0025 0.1673 13.172 4986 0.0004 0.0039 31.919

1.0000 0.000 0.031 0.0025 0.1673 13.172 4986 0.0004 0.0039 31.919

Note that theexternal quantities areduplicated toadepthof O.1 meters. As the su

penetrates the concrete, the concrete fails and the external boundary conditions move 

the concrete. Also note that the pressure (P) is zero from 0.75 to 1.00 meters. Sinc

concrete was cracked, the permeability of the cracked portion of the concrete overwhel

13
Although each LLW facility may be unique, most underground facilities can be repre

schematically y as in Figure 1. Below the ground surface there will be an engineered 

to deflect surface water away from the concrete vault below. The entire facility is site

geologically suitable location.

Figure 1: Schematic of underground LLW facility

From the schematic in Figure 1 it appears as though the roof is the most critical e

because it is most likely to have a moist environment, especially upon the failure 

engineered barrier. Because of this, the analysis of the entire concrete vault can conce

on the roof. If the roof is treated as a simple slab, the analysis simplifies further. Sin

flow through the roof slab will be approximately uniform over the surface of the slab, 

dimensional analysis of transport vertically through the slab should serve as a suffi

accurate model for transport through the roof of the vault.

7 Ion !llansport

A single transport equation is developed to propagate ions through the slab. This eq

can be converted into a finite difference equation so that it can be implemented in a com

program. A result of this approach is that time will advance in discrete intervals. After 

time interval each computational element is put in chemical equilibrium using vo

products and a charge balance. This step is achieved through dissolution/precipitation 

available salts. Any change in the quantity of solid salts in the pore space will chan

14
At the core of 4SIGHT is the advection-diffusion equation which establishes the tra

of ions through the slab. The advection-diffusion equation is simply the diffusion eq

with an extra term to account for Darcy flow. The development of the equation starts 

simply from the relation bet ween flux, j, and concentration, C:

j = –llVc + cu 

The parameter D is the effective diffusion coefficient, and the quantity u is the averag

fluid velocity. The rate of change in concentration is the negative divergence of eqn. 

(9C —= V” DVC– U”VC

at

after neglecting the divergence of the volume averaged velocity, since the fluid is vi

incompressible and the rate change in porosity is small.

7.2 Darcy Flow

The average pore fluid velocity can be calculated from the Darcy velocity, which 

volume-averaged pore fluid velocity. Given a porous media with permeability y k an

fluid viscosity p, the Darcy velocity, VD, is proportional to the pressure gradient, V

the density of the fluid [3]:

VD = -~ (Vp-pg)

P

The Darcy velocity VD can be related to the average pore velocity, u, from a geom

argument. Let v(x) represent the velocity of the pore fluid at some point x in the sla

sume that the pore fluid completely fills the available pore space. Representing the p

by ~, u is defined as:

1

/() u= P4vxdV 

The Darcy velocity is the average over the entire volume V:

1

/() ‘D=vv

vxdV

Since v(x) = O outside of the porosity, eqn. 5 can be limited to the pore space:

1

vD=— /() Vfb

vxdV 

15
Once the Darcy equation is incorporated into eqn. 2, the continuity equation is 

in order to update the pressures. Upon execution of the computer program, the 

ing sequence is continuously reiterated: transport ions, attain chemical equilibrium, 

transport properties, update boundary conditions. As the porosity of the concrete ch

the transport properties also change. The change in the transport properties will 

the hydraulic pressure distribution in the concrete since the transport properties w

by uniform throughout the concrete. The pressure distribution will be updated usi

continuity equation for porous media.

The continuity equation for a fluid is

ap

~+v”pv=o 

To develop a continuity equation for porous media, find the average integral value of 

over a representative volume, V:

Rearranging the order of integration gives

Finally, assuming p is constant and using the definition of Darcy velocity gives

This is the result obtained by Slattery [4] for porous media. Substituting for VD from 

gives

3$—=

at

V.qvp-pg) 

P

7.4 Dimensionless Variables

Eqn. 2, 3, and 11 can be combined to form a system of equations to propagate ions th

a porous media, However, the system of equations can be condensed by a transform

into dimensionless variables. Consider the following definitions:

t

X =—; T=—

kO/DO

P= p

pDO/kO

16
ac

“ VP+ & .Vc F“v Do o

This equation will be simplified further later using material properties related to por

7.5 Atkinson-Hearne Model

As yet, a dependable mechanistic model does not exist for sulfate attack. Therefo

incorporate sulfate attack, the depth of sulfate degradation is calculated using the At

and Hearne [5] model:

~ = E@ CocED

W)’(1 – v)

E Youngs modulus

co External sulfate concentration

CE Concentration of sulfate as

D Sulfate diffusion coefficient

a Roughness factor

‘Y Fracture surface energy

v Poisson ratio

Since the Atkinson-Hearne model

ettringite

gives the location of the sulfate front, 4SIGHT pre

that all of the concrete behind the sulfate front has been completely disintegrated, g

the properties of the surrounding soil. Since the transport coefficients of soil are much 

than concrete, the external boundary conditions are advanced to the sulfate front, c

a moving boundary condition.

The Youngs modulus, roughness factor, fracture surface energy, and Poisson ratio m

determined from experimental measurements. The quantity CE can be either calculated 

the cement composition, or more accurately from experiment. To experimentally dete

CE, the sulfate reacted per unit mass hydrated cement is plotted against the logarit

time. This data is fit to the equation [5]

tc

()

‘m = ?-nologlo —t,ck

where m is the moles of sulfate reacted in the cement, m. is the free parameter 

regression, t is time, c is the concentration of sulfate in liquid, t, is the characteristic 

for reaction, and ck is the concentration in kinetic experiments. The maximum valu

17
complete consumption of the C3A. The time to spalling is [5]

where

x:pa,,cE

t.pa[l=

2Dc0

X;pa,l =

2CY7(1– v)

E(/3C’E)2

In the event that t,pall< -tm,CE must be calculated in a self consistent manner using 

16,17, and 18. More complete details can be found in [5]. For ordinary portland ceme

is likely that tspall> tw. Therefore, since it is assumed that the quantity of external 

k Sufficiently great tO aCt aS an infinite reSerVOir (eXternal COnCentratiOn iS COnStan

can be calculated from the C3A content of the cement. Since each mole of ettringite r

one mole of A12~3 [6], the molar concentration of A12~3 will be the molar concentrati

.

ettringite. Given a concrete mix design having z~~~ kilograms of cement per cubic m

concrete, and the cement having a weight fraction 4A1203 of aluminum oxide, the m

ettringite formed per cubic meter of concrete is

CE =

xcem4Alz03

0.10196

with the A12~3 molar weight of 0.10196 kg per mole.

7.6 Effects of Cracks and Joints

Although not a degradation mechanism per se, effects due to the presence of crack

joints can be incorporated into the models for transport of the ions. Since the roof sl

likely be a supported member such that the bottom of the roof slab will be in tens

is assumed that cracks will appear on the bottom of the slab and extend upwards 

neutral axis of the slab. In the case of joints, 4SIGHT assumes that the joint extends th

the entire depth of the roof slab. The joint is filled with a joint compound with a 

permeability and service life.

The permeability of a cracked slab is calculated assuming the crack walls are s

and parallel. Given a square slab with width L and depth D“ having cracks with w

penetrating the full depth D of the slab, the permeability of the slab is a weighted 

18
L+w12 L+wko

Since w3/12 istypically far greater than LkO, the permeability of the slab can beapp

mated by the permeability due to the crack. Further, if each of the cracks of width 

spaced a distance a apart, the permeability y of the slab, k,, is

W3

k~=~ 

Joints can be handled in a similar manner as cracks. However, joints will typica

very much wider than cracks Since joints will presumably extend the entire thickness 

slab, once the joint fails, the flow through the joint would overwhelm the transport 

through the central portion of the slab. In fact, the transport coefficients could be a

as, or greater t ban, those of the soil. Therefore, upon failure of the joint, 4S I GHTa

that the roof fails to impede the flow of water into the vault, the transport propert

the concrete should be approximated by the transport properties of soil, and the calcu

ceases.

8 Chemical Equilibrium

After each time step, each computational element is brought to chemical equilibriu

satisfying two conditions:

1.

2.

If a salt exists as solid, the constituent ion concentration product equals the vo

product.

The sum of the free charges from all available ions equals zero, insuring local 

neutrality.

Given the salt A@C., composed of anion A“- and cation C’@+,condition 1 above 

that if AOCa exists as a solid then

[Aa-]@[CP+]@ = K,, 

where K~P is the ion volubility product. Condition 2 implies that with m anions 

cations present in the pore solution:

~ QA;- = ~/?C;+ 

i=l j=l

19
eqn. 24

These two steps are repeated, as necessary, until both equations are satisfied.

9 Failure Criteria

Once the rules for propagating ions through the slab have been defined, criteria for ter

ing the calculations are needed for each degradation mechanism. Termination occurs 

suitable failure criteria have been met.

9.1 Reinforcement Corrosion

Corrosion of the rebar is due to the presence of chloride ions and proceeds in two 

initiation and corrosion. The initiation stage is the time during which chloride io

diffusing through the concrete. During this time there is an insufficient concentrati

chloride ions at the rebar for corrosion to occur. As the concentration of chloride ions 

rebar increases, the pH decreases to insure charge neutrality. When a sufficient concent

of chloride ions reaches the rebar, corrosion begins. The corrosion stage begins at the 

of corrosion and lasts until failure of the steel reinforcement, which is typically only 

years hence. Since the initiation stage may last for hundreds of years, the duration 

corrosion stage is insignificant to the total lifetime of the reinforcement. Therefo

accurate estimate of the service life of the concrete can be approximated from the du

of the initiation period.

9.2 Sulfate Attack

As the sulfate front proceeds into the concrete, the effective thickness of the concre

creases. At some point in time the roof has an insufficient thickness to support its lo

this time, structural failure of the roof will occur. This critical depth of sulfate p

tion can only be determined through detailed structural analysis. Therefore, it is the 

responsibility to provide this information.

20
pj

10 Calculated Material Properties

All of the concrete physical parameters (e.g., diffusivity, permeability, etc.) are user sp

inputs to 4SIGHT. However, in cases where not all properties are available, missing qua

must be approximated using existing correlations. The physical properties must be 

lished due to both hydration and leaching. The physical properties due to hydration 

initial conditions. However, as the porosity changes due to leaching, corrected values 

physical parameters are needed.

10.1 Hydration

The hydration of cement can be approximated by a reaction between tri-calcium 

(C3S) and water, forming a calcium silicate hydrate (CSH). A more elaborate model 

porating multiple mineral phases would require chemical analysis of the cement and 

relatively little additional information concerning degree of hydration. The weight r

water to cement, ~, is the oft reported quantity to characterize the concrete mix. 

some period of hydration, the fraction of the initial C3S which has hydrated is the 

of hydration, a, The relation between these two properties and porosity can be deter

stoichiometrically [7],

1 + 1.31a

@=l-1+32y 

c

and is valid for ~ values used in practice.

The diffusivity can be related to either ~ or ~. After the first 100 days of hydratio

transport properties of most cement pastes are near their asymptotic values. Althoug

values are still changing after 100 days, these changes are small compared to the ac

with which these transport measurements can be made. Due to this apparent steady 

an empirical relation between Do, the chloride diffusivity in m2/s, and ~ was develop

cement paste by Atkinson, Nickerson, and Valentine[8] and Walton, Plansky, and Sm

loglo DO = 6.0; – 9.84 

for values of ~ in the range (0.2–0.6).

21
y, y 

and l?(x) is the Heaviside function. The quantity 0 is a constant for all ions.

The above equations relate the diffusivity of cement paste to water: cement ratio o

ity. A relationship is now needed between cement paste diffusivity and concrete diff

Experimental results of Luping and Nilsson[l 1] for cement paste and mortar sugge

their diffusivities are approximately equal. Additionally, results from numerical exper

by Garboczi, Schwartz, and Bentz[12] investigating the effects of aggregate-paste int

zone diffusivity upon bulk diffusivity suggest that concrete diffusivity is approximately 

to the paste diffusivity.

Given ~, the permeability can be approximated from the data in Hearn, et al.[13

for ~ in

10.2

Once D

k = 10500~ X 10-21 rn2 

the range (0.35,0.80).

Leaching

and k have been established, changes due to leaching can be calculated from 

where +0 is the initial porosity due to hydration. As the Ca(OH)2 is leached, the p

increases. Let the value of porosity after leaching be 4’, and the diffusivit y be D’. 

tunately, the ratio D’/D~ is not simply 0(#’) because as calcium hydroxide is leache

the paste the ratio lY/@ does not retrace eqn. 27. Rather, the ratio D’/Dj is greate

0(~’), as demonstrated by the NIST microstructural model.

The NIST cement microstructural model[14] was used to determine the relation for 

upon leaching. Results for a ~ =0.35 paste are shown in Figure 2. The formation 

decreases with decreasing porosity due to hydration, denoted by circles, lJpon leac

the calcium hydroxide, the formation factor follows the curve denoted by squares. Gi

following definitions:

00= O(+.) ti’ = 0(($) 

An empirical relation was developed to relate the leached pore structure to the unda

pore structure.

( = ~o + 2“0($ – ~o) 

and is shown by the solid curve in Figure 2. Therefore, the ratio of the leached v

diffusivity, D’, to the initial value D. is

D’ t

Z=K

22
0

1

J

o e,

I I I

o 0.05 0.1 0.15 0.2 0.25 0.3

)

Figure 2: Formation factor, n/@, due to hydration (circles) and to leaching (squa

w/c==O.35 (after [14]). The solid line represents the approximation used by 4SIGHT 

leaching formation factor.

The relative change in permeability can be calculated using ~. The Katz-Tho

equation relates permeability to d [15]:

where dC is the diameter of the largest sphere which can pass through the pore spac

sample. Also, given that dC cc 0 [16], the relative change in permeability y becomes

The final dimensionless advection-diffusion equation is a combination of eqns. 14

33:

aci

—=V

()

~i 3

~T

.$vp+ ~ vD “VG 

o 0

There is one such equation for chemical species i.

11 Program Outline

The flow of the program is summarized by the pseudo-code program shown in Figure 

program naturally divides itself into three parts: initialization, ion propagation, a

state output.

23
I

ion Irmsporl

I

Y

NO

*

GE)

Figure3: 4SIGHTfiow chart.

24
1.

2.

3.

4.

..

5.

6.

13

The vault system is, and remains, water saturated. This should be valid for all 

most arid locations.

Darcy flow is valid, even at the bottom of the slab which might not be in direct 

with liquid.

There is sufficient oxygen present for corrosion. This a conservative approach 

rosion.

The concrete vault begins service approximately 100 days after casting. Therefo

transport properties have reached about 9070 of their asymptotic values.

The external ion concentration and hydraulic head are constant over the lif

vault.

As the sulfate front advances, the concrete in its wake can be treated like soil 

external conditions can be advanced to this point.

Acknowled~ements u

The authors wish to acknowledge the financial support of the U.S. Nuclear 13egulatory 

mission (NRC!). Mr. Jacob Philip was the NRC Program Manager and his advice was 

appreciated. The authors would also like to express their thanks to Dr. James Pomme

of Bucknell University and Mr. Dale Bentz of NIST for their discussion and commen

25
[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]

[12]

[13]

J.M. Pommersheim, and J.R. Clifton, “Models of Transport Processes in Con

IVIS7’ll? 4405,U.S. Department of Commerce, (1990).

R.B. Bird, W.E. Stewart, and E.N. Lightfoot, Transport Phenomenon, John W

Sons, New York, (1960).

J.C. Slattery, Mass Transport, McGraw-Hill Book Company, New Yorlc, (1972).

A. Atkinson, and J.A. Hearne, “Mechanistic Models For The Durability Of C

Barriers Exposed To Sulfate-Bearing Groundwaters,” JZl?S’Proc., 176, 149–156, 

H.F.W. Taylor, Cement Chemistry, Academic Press, New York, (1990).

Young and Hansen, “Volume Relationships for C-S-H Formation Based on Hy

Stoichiometries,” Proc. MRS, Eds. Struble and Brown, 85, 313-322, (1987).

A. Atkinson, A. I<. Nickerson, and T.M. Valentine, “The Mechanism of Leaching 

Some Cement-Based Nuclear Wasteforms,” Radioactive Waste Management 

Nuclear Fuel Cycle, 4 (4), 357-378, (1984).

J. Walton, L. Plansky, and R, Smith, “Models for Estimation of Service Life 

crete Barriers in Low-Level Radioactive Waste Disposal,” NUREG/CR-5542, 

2597, (1990).

E.J. Garboczi and D.P. Bentz, “Computer simulation of the diffusivity of cement

materials,” J. Materials Science, 27,2083-2092, (1992).

T. Luping, and L. Nilsson, “Rapid Determination of the Chloride Diffusivity in C

by Applying an Electric Field,” ACI Materials Journal, 89 (l), 49-53, (1992).

E.J. Garboczi, L.M. Schwartz, and D.P. Bentz, “Modelling the influence of the int

zone on the D. C. electrical conductivity y of mortar,” (Submitted to) JournaZ 

ACBM.

N. Hearn, R.D. Hooton, and R.H. Mills, “Pore Structure and Permeability,” Signif

of Tests and Properties of Concrete and Concrete-Making

Lamond, ASTM, 240-262, (1994).

26

Materials, Elds. Klieg
Rock,” Phys. Rev. B, 34 (11), 8179-8181, (1986).

[16] S. Li, and D.M. Roy, “Investigation Of Relations Between Porosity, Pore St

and Cl- Diffusivity Of Fly Ash and Blended Cement Pastes,” Cement and C

Research, 16 (5), 749, (1986).

[17] T. Mason and D. Brown, Lex @ Yacc, O’Reilly & Associates, Inc., Sebastopo

(1991).

27
pg

The input parsing was accomplished using public domain versions of LEX and YACC

for the PC. These tools simplify parsing user input, facilitating modification and ext

in program capabilities. The user unfamiliar with LEX and YACCshould read the b

Mason and Brown listed in the references. To obtain a copy of the public domain 

access the following directory through anonymous ftp:

machine: omnigate. clarkson. edu

directory: pub/msdos/dj gpp

LEX and YACCare programs that take as input user specifications for a particular 

and create a C file which performs the functions specified in the user LEX and YA

The resulting C! file is then simply ttinclude-ed into the main C program, Of cour

user can code the parsing directly in C, but LEX and YACCmake the coding less err

and more extensible.

A.1
integer

dreal

ereal

real

cliff

perm

Wc

thick

neuax

youngs

beta

ce

rough

gamma

poisson

ext rn

int rn

time

depth

rebar

output

head

crack

joint

at

until

assign

wspace

kill

nl

%%

{wspace}

OH

cl

S04

[0-9] +

( [0-9]+”. “ [O-9]*) [ ([O-9]*” .“ [0-9]+)

({dreal} I {integer}) [eE] [+-1? [o-91+

{dreal} I{ereal}

[IM] [11] [Ff] [Ff]

[Ppl [Ee] [Rr] [Ml

[Wwl[cc]

[Tt] [Hh] [Ii] [Ccl [Kk] [Nn] [Ee] [Ss1 [Ss1

[Nn] [As] [XX] [Ii] [ss]

[Yyl [001 [Uul [Nn] [Gg] [Ss1

[Bb] [Ee] [Tt] [As]

[Ccl [Ee]

[Rr] COO][Uul [Gg] [Hh] [Nn] [Ee] [Ss] [Ss]

[Gg] [As] [Mm][h] [As]

[Ppl [Ool [Ii] [Ss] [ss1 [Ool [Nn]

[Ee] [XX] [Tt] [Ee] [Rr] [Nn] [As] [Lll

[Ii] [Nn] [Tt] [Ee] [Rr] [Nn] [As] [Lll

[Tt] [11] [Mm] [Ee]

[Dal] [Ee] [Ppl [RI [Hh]

[Rr] [Ee] [Bb] [As] [Rr]

[00] [Uu] [Tt] [PP] [Uu] [Tt]

[Hh] [Ee] [As] [Dal]

[cc] [Rr] [As] [Ccl [Kk]

[Jj] [001 [Ii] [Nn] [Tt]

[As] [Tt]

[UUI [Nn] [Tt] [Ii] [Lll

It=,l I !, :,,

[ \tl +

(( [Qq] ( [Uul [Ii] [Tt] )7) I ( [Ee] ( [XX] [Ii] [Tt])?))

\n

.

~ return HYDROXIDE;}

{ return CHLORINE;}

{ return SULFATE;}

A.2
{cliff}

{perm}

{WC]

{thick}

{neuax}

{youngs}

{beta}

{cc>

{rough}

{gamma>

{poisson}

{extrn}

{intrn}

{time}

{depth}

<rebar}

{output}

{head}

{crack}

{joint}

{at>

{until}

{integer}

{real}

{assign}

{nl}

{kill}

i%

{ etuOSSU}

{ return DIFF;}

{ return PERM;}

{ return WC;}

{ return THICKNESS;}

{ return NEUAX;}

< return YOUNGS;}

{ return BETA;}

{ return Cc;}

{ return ROUGHNESS;}

{ return GAMMA;}

{ return POISSON;}

{ return EXTRN;}

{ return INTRN;}

{ return TIME;}

{ return DEPTH;}

{ return REBAR;}

{ return OUTPUT;}

{ return HEAD;}

{ return CRACK;}

{ return JOINT;}

{ return AT;}

{ return UNTIL;}

{ sscanf(yytext, “~ld”, &yylval.lnt32);

return INT32;}

{ sscanf(yytext, “%lf”, &yylval.real);

return REAL;}

{ return EQUALS;}

{ return NEWLINE;}

{ return KILL;}

{ return yytext[Ol;}

A.3
void yyerror(s)

char *s;

{

/* fprintf(stderr,”~s \t see line %d\n’’,s,input-file-line); 

void syntax-error(char *);

void syntax-error(char *s)

{

fprintf(stderr,”~s \t see line %d\n’’,s,lnput-flle-llne);

3

%union {

double real; /* real value */

long int32; /* DoubleWord Integer 

int int16; /* SingleWord Integer 

3

%-token <real> REAL

%token <int32> INT32

%token KILL

)!token EQUALS

%token NEWLINE

%token DIFF

%token PERN

%token WC

%token THICKNESS

%token NEUAX

%token YOUNGS

%token BETA

%token Ce

%token ROUGHNESS

%token GAMMA

A.4
%token CALCIUM

xtoken SODIUM

%token POTASSIUM

%token HYDROXIDE

%token CHLORINE

%token SULFATE

%token CARBONATE

%token TIME

%token DEPTH

%token REBAR

%token OUTPUT

%token HEAD

%token CRACK

%token JOINT

%token AT

%token UNTIL

%type <real> expr

~type <int16> anion

~type <int16> cation

~left ‘+’ ‘-’

Xleft ‘*’ ‘/’

lines: /* NOTHING*/

I lines line

.

*

line: NEWLINE {input-file.line++;}

I terminate { printf(’’%d lines parsed.\n’’,lnput-fil~:-llne);

return 0;}

I assign

I recvr NEWLINE

A.5
1

.

*

expr:

I

.

*

assign:

I

I

I

I

I

I

I

I

I

I

I

I

I

recvr EQUALS expr { syntax-error(’’UnknoWn Command’’);}

INT32 {$$ = (real)$I; }

REAL {$$ = $1;}

DIFF EQUALS expr {

PERM EQUALS expr {

WC EQUALS expr {

THICKNESS EQUALS expr {

NEUAX EQUALS expr {

YOUNGS EQUALS expr {

BETA EQUALS expr <

Ce EQUALS expr {

ROUGHNESS EQUALS expr {

GAMMA EQUALS expr {

POISSON EQUALS expr {

EXTRN cation EQUALS expr

EXTRN anion EQUALS expr

INTRN cation EQUALS expr

Dlnfty.value = $3;

Dinfty.is_default = FALSE;}

kinfty.value = $3;

kinfty.is-default = FALSE;}

wc value = $3;

wc.is.default = FALSE;}

thickness value = $3;

thickness.is-default = FALSE;

sample-length = thickness.value;}

neutral.axis-depth.value = $3;

neutral-axis-depth .is-default=FALS

Youngs.value = $3*1.0E+09;

Youngs.is_default = FALSE;}

beta. value = $3;

beta.is-default = FALSE;}

CE. value = $3;

CE.is-default = FALSE;}

roughness.value = $3;

roughness.ls-default = FALSE;}

gamma.value = $3;

gamma.is-default = FALSE;}

nu value = $3;

nu.is-default = FALSE;}

{ catlon[$21 .c[O1 = $4;}

{ anlon[$21 .c[O1 = $4;}

{ catlon[$21.c[ll =$4;}

A.6
I

I

I

I

I

I

I

I

I

.

>

cation:

I

I

I

.

$

DEPTH EQUALS expr {

REBAR EQUALS expr {

OUTPUT cation anion {

OUTPUT cation

OUTPUT anion

HEAD EQUALS expr

CRACK EQUALS expr

JOINT EQUALS expr

{

{

{

MaxDay.value = MIN(MaxDay.value,

MaxDay.is-default = FALSE;}

termination-type = STRUClr-LIMIT

MaxDepth = $3;}

rebar-depth=$3;

termination.type = STRUCT-LIMIT

sol-array[$2] [$3] .output..flag=TRUE

cation[$2] .output-flag = TRUE;}

anion[$2] .output-flag = TRUE;}

head. value = $3;

head.is-default = FALSE;)P

AT expr DEPTH expr

{ crack-width.value = $3;

crack-width. is-default = FALS

crack-spacing = $5;

crack-depth = $7;}

AT expr UNTIL expr

{ joint-width= $3;

joint-spacing = $5;

joint-lifetime = $7*365;

if(MaxDay. is-default == TRUE)

MaxDay.value = joint-lifetime;

else

MaxDay.value = MIN(MaxDay.value,joint-lifeti

MaxDay.is-default = FALSE;

joint-is-specified = TRUE;}

JOINT PERM EQUALS expr { joint-permeability= $4;}

HYDROGEN { $$ = H;}

CALCIUM { $$ = Ca; >

SODIUM { $$ = Na; }

POTASSIUM <$$ = K;}

A.7

of two programs: cweave and ct angle. The user writes a CWEBfile using any ASCII 

This file contains both the C code and the documentation. ct angle takes this inpu

and extracts the C code. Similarly, cweave takes the input file and creates a TEXfile 

can be ~-ed, creating a formatted version of the complete documentation, includin

source code. Therefore, to use these tools, the user must have both CWEBand TEXf

PC. Fortunately, both of these tools are public domain, and can be accessed by anony

ftp:

~ machine: f tp. nj it. edu

directory: pub/msdos/emt ex

CWEB machine: labrea. stanford. edu

directory: pub/cweb

The only modification made to CWEBfor the development of 4SIGHT was the exte

to L TEX. To use the extensions of L TEXwith CWEB,one only needs a L TEXstyle sheet. 

too is available by anonymous ftp:

machine: ftp.th-darmstadt. de

directory: pub/programming/literate-programming/c. c++

The pages that follow are verbatim output from cweave, via L TEX:

B.1
1.

2.

14.

18.

22.

27.

31.

39.

53.

54.

56.

61.

63.

65.

68.

69.

73.

INTRODUCTION

MAIN

User defined data structures..................................

Global variables.........................................

Interrupt handlers........................................

INITIALIZATION

Initializesolubilityarray....................................

Initializeparameters......................................

ADVECTION-DIFFUSION

Ion transport...................................... ...”

Chemicalequilibrium......................................

Advancesulfatefront.....................................

Chloridepenetration.....................................

Update Pressures.......................................

Adjust physical parameters..................................

Proper time.................................. ..””.”.””

OUTPUT
objectiveisacheivedby consideringprominentdeterioration mechanismsto predicthlei

permeability, effluent, and workingthicknessofthevaultslab.

The conditionofthevaultslabisassesssedby modellingtheingressofdeleterious ion

by bothdiffusionad convectiveflowdue tohydraulicpressuregradients.Additionally, ch

physicalpropertiesdue to the leaching/deposition of salts are incorporated into the program

propagation of ions are calculated using an advection-diffusion equation. After each discret

step, each computational cell is brought to chemical equilibrium through precipitation/dissoluti

of ionic salts.

The program incorporates a one-dimensional model of the vault since the critical deterio

will most likely occur at the ceiling of the vault. The ceiling is the most likely element to exp

moisture, stress, and external ions such aa Cl- and S04. Additionally, the majority of the tra

will occur at the center of the roof slab since the edges of the slab will be supported. There

simple one-dimension model of the flow at the center of the slab will give a conservative estim

the transport for the entire slab.

The program proceeds in a strainghtforward manner; The program gets input parameters 

the user; initializea the ion information and scaling parameters; iterates diffusion and con

calculations, bringing each computational element to chemical equilibrium between iteration

termines whether failure has occured; and prints pertenent information about the per fcxma

the vault.

The output of the program can be printed to either the screen or to an ASCII file. The 

file can then be used as input to a spreadsheet program for further analysis.
below.

#define banner “\”\tTh%isU4SIGHTU( Versional. O)\’’\n\n”

#define NUM.CELLS 20

#define NUM.SURFACES(NUM-CELLS+ 1)

#define DLIMIT 0.0005

(Include files s);

( Preprocessor definitions 11);

( User defined data types 14);

( Global variables 4);

( Function declarations 23);

(Interrupt handlers 24);

(Input parsing routine 13);

int main (int argc, char **argv )

{

int i, j, k;

( Systemsetup 5);

( Parse input 6);

( Parameter initialization 7);

( Print header 75);

( Print intermediate results 76);

do {

if (kbhii( )) {

i= geich( );

cmiroLbreak( );

}

(Ion transport 8);

( Adjust parameters g);

( Print intermediate results 76);

(Assess termination 10);

} while (miermination);

( Print termination information 77);

( Print desired output 78);

return 1;

}

3. The following are *include files needed for some of the intrinsic functions such as ,pr

scanf. The math functions pow and log and the math error routine matherr require 

Additional include files are declared when needed. The routine kbhit ( ) requires the he

conio. h.

(Include files 3) s

#include <stdio .h>

#include <math .h>

#include <conio .h>

#include <bios .h>

#include <time. h>

See also sections 22, 28, and 32.

This code is used in section 2.
See also sections 12, 19, 20, 21, 29,31, 40, 63, and 69.

This code is used in section 2.

5. System Setup.

#define CTRL.BRK.ON 1

#define GET.TIXE O

(Systemsetup5) S

printf(banner);

elapsed-time = biostime (GET-TIME, O~); /* get BIOS timer value */

( Determine STDIN 30);

i = setcbrk(CTRL-BRK_ ON); /* check ctrl-brk every system call */

ctrlbrk(control-brea k);

sigraal(SIGFPE, div-O);

iniiialize_ion-data ( );

This code is used in section 2.

6. Parse Input.

( Parse input 6) s

i = yyparse( );

This code is used in section 2.

7. Parameter Initialization.

( Parameter initialization 7) E

initialize-parameters( );

update-pressures ();

aci-guidlines( );

chemical-equilibrium (FALSE);

chemical-eguilibrium (FALSE); *

chemicaLequilibrium (FALSE);

( Print initial system state 74);

This code is used in section 2.

8. Ion Transport.

(Ion transport 8) s

ion-diflusion [ );

chemical_eguilibrium(TRUE);

advance_sulfate-front ();

This code is used in section 2.

9. Adjust Parameters.

( Adjust parameters 9) s

( Advance global clock 70);

(Adjust physical parameters 68);

update-pressures ();

( Update the porosities 67);

This code is used in section 2.
This code is used in section 2.

11. The following are useful macro definitions which save a lot of typing.

( Preprocessor definitions 11)=

#define SqR(a) ((a)* (a))

#define CUB(a) ((a) * (a)* (a))

#define MAX(a, b) ((a > b) ? a : b)

#define MIN(a, b) ((a <b)? a : b)

See also section 46.

This code is used in section 2.

12. Declare yyparse which is written in YACC and LEX. The function yyparse takes no arg

and simply parses the input, assigning values to global variables. The routine does ncjt “

until the entire input stream has been parsed.

( Global variables 4)+-

int input. file.line = 1; /* line number of input file (lex/yacc) */

13. (Input parsing routine R) s

int yyparse (void);

#include “yyparse. c“

This code is used in section 2.
Using real lets the programmer easily change between float and double data types for real n

The boolean types are defined {FALSE, TRUE} because Borland assigns the first entry the 

O. Variables defined as boolean can be used in a conditional statement as if(myboolvw) 

direct comparison to FALSE or TRUE.

( User defined data types 14) =

t ypedef double real;

typedef enum {

FALSE,TRUE

} boolean;

See also sections 15, 16, end 17.

This code is used in section 2.

15. The ION data type contains all important information about a single ion. This info

includes:

co: concentration at each computational cell surface.

q: the ratio Dji/DO foreachion.

Dj: freeiondiffusivity.

valence: the valence of the ion.

output -f lag: a flag todenotewhetherthefinalionconcentration, asafunctionofdistancei

insystem.out.

nanaeo: the ASCII name of the ion.

( User defined data types 14)+=

typedef struct {

}

16.

real c[MJM.SURFACES]; /* concentration at each surface */

real moles [WJM-SURFACES]; /* moles ion in solution */

real q; /* Dfi/Do */

real Dj; /* ion free diffusivity */

real valence; /* ion valence */

boolean ouiput-j?ag; /* print to system .out? */

char name[5]; /* ASCII name of ion */

ION ;

The SOLARRAY data type contains the volubility products and stiochiometric ratios 

anion-cation combination. The salt is C~A~, where C represents the cation, and A the ani

quantity of solid is defined on the cell surfaces.

( User defined data types 14)+=

typedef struct {

1

real sIIWILSURFACES]; /* moles solid salt in element */

real ksp; /* salt volubility constant */

real m; /* stiochiometric cation factor */

real n; /* stiochiometric anion factor */

real molar-density; 1* mn3/mole */

boolean output-jlag; /* print salt to system. out? */

SOLARRAY;

real @JMCELLS]; /

changing formation factor /

real <’[NUM-SURFACES]; /* changing formation factor at the surface */

real WINIJILSURFACES]; /* Pressure */

real & [NUN_CELLS]; /* current porosity */

real #’n [NUN-SURFACES]; /* porosity at the surface */

real #n-’[NUM_SURFACES]; /* previous porosity at concentration */

real litre [NUH.URFACES]; /* litres of solution (surface area)*/

real strength [NUM-SURFACES]; /* strength in psi */

int FIRST_CELL = O; /* first active cell number */

20. User specified parameters. These parameters determine the ultimate properties of the 

Although default values are given, user specified values will override these values.

( Global variablea A) +-

MATERIAL thickness = {1.0, TRUE}; /* roof thickness */

MATERIAL neutral.axis-depth = {0.75, TRUE}; /* depth of neutral axis */

MATERIAL km = {2.5. 10-18, TRUE}; /* ultimate permeability */

MATERIAL Dm = {5.7 o10-12, TRUE}; /* ultimate diffusivity */

MATERIAL head = {5. O,TRUE}; /* external pressure head */

MATERIAL WC = {0.45, TRUE}; /* water: cement */

MATERIAL Youngs = {20.0. 10+0’, TRUE}; /* Young’s modulus */

MATERIAL ~ = {1.8 010-06, TRUE}; /* linear strain coefficient */

MATERIAL CE = {350., TRUE}; /* cone. sulfate aa ettringite */

MATERIAL roughness = {1.0, TRUE}; /* roughness factor */

MATERIAL y = {10.0,TRUE}; /* fracturesurfaceenergy*/

MATERIAL v = {0.2, TRUE}; /* Poissonration*/

21. Scaling parameters.

( Global variables 4)+=

real L = 0.010; /* length of ELEMENTS so that AX = 1 */

real DO = 1.0. 10-ll; /* scaling diffusivity */

real kO = 1.0. 10-18; /* scaling permeability */

real AT= 1.0; /* dimensionless time increment */

real AX[NUH-SURFACES]; /* dimensionless length between each cell */
The

routine signal requires signal. h.

(Include files 3)+=

#include <dos. h>

#include <signal .h>

23. ( Function declarations 23) s

int control. break(void);

int matherr(struct ezception *);

void div-O (void);

See also sections 33, 39, 43, 51, 54, 56, 59, 60, 61, and 65.

This code is used in section 2.

24. control-break handles interruptions due to the use striking CTRL-BRK during execution.

(Interrupt handlers 24) s

int control-break (void) {

int i, j, k;

printf(’’\n\n\tUSERuCTRL-BRK ! . \n\n” ); ( Print desired output 78)

return O; }

See also sections 25 and 26.

This code is used in section 2.

25. The routine matherr handles exceptions from the math coprocessor.

( Interrupt handlers 24)+-

int matherr(struct exception *e)

{

print~(’’\n\n”);

prini~(”\tmatherr: u“);

switch (e-type ) {

case DOMAIN: printj(’’DOMAIN”);

break;

case SING: pm”ntf(’’SINGULARITY”);

break;

case OVERFLOW:printj ( “OVERFLOW”);

break;

case UNDERFLOW:printf(’’UNDERFLOW”);

break;

case TLOSS: printf(’’TLOSS” );

break;

default: print$(’’UNKNOWNUERRORUTYPE!”);

break;

}

printf(’’\n”);

print~(’’\n\n\t “);

print~(”f unction: \t%s (%lf ,%lf ) \n”, e-name, e-argl, e-arg2);

printf(’’\n”);

ezii (0);

}

#include <errno .h>

29. If redirection is not used, set input.jile_handle to the file specified at the commandl lin

( Global variables 4 ) +=

int input.jiie-handie; /* file containing input parameters */

30. If an input file was given without redirection, use dup2 ( ) to copy the input file handle 

to the existing stdin file handle.

( Determine STDIN 30) s

if (argc > 1) { /* redirect from input file */

fprintf(stderr, ‘“Calculat ing. . . “);

if (( inputJile-handle = open (argv [1], OJDONLY)) = –1) {

switch (ermo) {

case ENOENT: fprinif(siderr, “Nousuchuf ile: &\n”, ar’gv[1]);

break;

case EMFILE: jpn’ntf (stderr, “Tooumanyuopenuf iles. \n” );

break;

case EACCES: fpra”ntf(stderr, “Permissionudenied. \n” );

break;

case EINVACC: fprintf (stderr, “Invaliduaccess ucode. \n” );

break;

default: fprintf(stderr, “Error: JJnknownuerrorucode. \n”);

break;

}

ezit(0);

}

dup2 (input-file-handle, STDIN);

}

else /* interactive mode */

printf (“EnterUcommands: \n”);

This code is used in section 5.
initializes thetwo ION arrays:anion and cation

( Globalvariables4)+=

int num. cations;

int num_anions;

SOLARRAY **sol_ array;

ION *cation;

ION *anion;

32. The routine malloc( ) requires stdlib. h.

( Include files 3)+-

#include <stdlib. h>

33. ( Function declarations 23) +-

void iniiialize_ion_ data (void);

34. Initializing ion and salt array. The data for the ions and salts is in the file ion. db. M

determine the number of anions and cations, including OH- and H+, and then allocate 

for anion, caiion, and sol-array.

void initialize-ion-data ( )

{

int i, j, k;

FILE *ion_ database; /* ion and salt data base */

if ((ion-database = ~open(’’ion. db”, “rt”)) s A) {

print~(”\n\tERROR: ~uUnableutouopenuion. db !\n\n”);

ezit (0);

}

fscanf (ion_database, “%du%d”, &num-cations, &num-anions);

( Allocate ion and salt arrays 35)

( Input ion data 36)

( Input salt data 37)

fclose (ion-database);

( Zero-out ion and salt data 38 )

}

35. Allocate memory for all of the arrays. The total number of cations/anions includes OH

H+.

( Allocate ion and salt arrays M) s

anion = (ION *) malloc(num-anions * sizeof(ION));

cation = (ION *) ma/ioc(num-cations * sizeof(ION)); /* allocate 2-D array */

sol-array = (SOLARRAY w) malloc(num_cations * sizeof(SOLARRAY *));

for (i = O; i < num-cations; i++) {

sol-array [i] = (SOL ARRAY T.) malloc(num-anions * sizeof (S OLARRAY));

}

This code is used in section 34.
for (j= O; j < num.anions; ji+) {

fscanf (ion-database, “%s”, anion [j].name);

fscanf (ion-database, “%lf”, &anion ~].valence);

fscanf (ion-database, “!!le”, &anion ~].ll~);

anion ~]. output-flag = FALSE;

}

This code is used in section 34.

37. Read salt data in cation major order.

( Input salt data 37) s

for (i= O; i < num-cations; ii+)

for (j = O; j < num-anions; j++) {

fscanf (ion-database, “%le’’,&soi-away[i]~].ksp);

fscanf (ion-database, “%lf”, &soLarray[i]~].m);

fscanf (ion-database, “?!lf’’,&soLarra3f[i]~].n);

fscanf (ion-database, ‘Wf”, &so/-array [i]~].molar-density);

sol-array [i]~]. ordput$ag = FALSE;

}

This code is used in section 34.

38. Initialize the ion and salt data to zeros.

( Zero-out ion and salt data 38) s

for (k= O; k < EUM.SURFACES; k++) {

for (i= O; i < num-cations; ii-t) {

cation [i].c[k] = 0.0;

caiion[i].moles[k] = 0.0;

for (j= O; j < num-anions; j-l+) {

anion ~].c[k] = 0.0;

anion ~] moles [k] = 0.0;

sol-array [i]~].s[k] = 0.0;

}

}-

}

This code is used in section 34.
#define OH O

#define Cl 1

#define S04 2

#define C03 3

#define H O

#define Ca 1

#define Na 2

#define K 3

#define TIllE_LIMIT O

#define STRUCT_LIMIT 1

#define PRINTOUTS 20.O

( Function declarations 23)+-

void initialize.pararn eters ();

40. ( Global variables 4) +-

real sample.length = 1.0; /* length of sample in meters */

real sulfate_ failure.year;

real chlom’de.failure_y ear;

real R[NUM_CELLS]; /* perm. coefficient due to cracking */

real sulfate-rate; /* dimensionless rate */

real sulfate-rate.ms; /* m/see */

real sulfate-depth = 0.0; /* depthofsulfatepenetration*/

real chloride-depth = 0.0; /* depthofchloridepenetration*/

real rebar-depih = 1.0;

boolean termination = FALSE; /* simulation termination */

int termination-type = TIME-LIMIT; /* limiting condition for termination */

MATERIAL MazDay = {10000.0, TRUE}; /* time limit (day) */

real MaxDepth = 100.0; /* penetration depth limit (m) */

MATERIAL crack-width = {0.000100, TRUE};

real crack-spacing = 100.0;

real crack-depth = 0.0;

real crack_permeabilit y;

real joint-width = 0.0;

real joint-spacing = 100.0;

real joint-lifetime = 20000.0;

boolean joint.change_fiag = FALSE;

boolean joint-is_specijied = FALSE;

real joint_permeability = 0.0;

real Pout; /* dimensionless ext. pressure */

real dP; /* dimensionless rho g h pressure */

real vINI.KCELLS]; /* Darcy velocity */

real p = 0.001; /* pore fluid viscosity */

real g = 9.8; /* gravitational constant(rn/sec/see) */

real p = 1000; /* density of pore fluid (kg/m3) */

real kB; /* bulk permeability */

real DB; /* bulk diffusivity */
real

real

real

real

real

real

real

real

real

real

0; /

degree of hydration /

f; /* initial solids fraction */

Vsample; /* volume of element */

pCH = 33.1; /* density of CH */

molesCH = 13.0; /* moles CH/ltr soln. */

19rJ= 0.00281; /* initial l/F */

19; /* current l/F */

print.day; /* next day to printinterm.results*/

print. day-interval; /* days between interm. results */

earliest-failure-day; /* earliest structural failure */

41. initialize-parameiemo. This routine initializes the scaling parameters and physical 

ters.

void initialize-parameterso

{

int i, j, k;

int n, m;

int iter;

real moles-cation, moles-anion, charge, x, dx, z, criiical.C’l-concentration, T-crit

concen-ratio, darcy-velocity, perm-iime, chloride-time, perm-factor,

/* modify DT to account for high permeability */

perrmdepih, A, crack-factor; /* modify DT to account for cracks */

for (k = FIRST.CELL; k < NUM.SURFACES; k++)

&x[k]= 1.0;

( Establish material parameters of concrete 42)

aci-211 ( );

( Propagate ion concentration in concrete 45)

( Print parameters to sidout 47)

( Initialize ion ‘eta’ parameter 48)

L = sampleJength/NUM_ CELLS; /* the universal length scale */

MaxDepth = MaxDepth /L; /* dimensionless */

( Calculate crack and joint adjustment to permeability 49)

cm”tical-Cl-concentration = 0.000400 * Vsample * (1.0 – 4) * 2.5/35.4;

if (anion [cl] .c[FIRST_CELL] > critical-cl-concentration ) {

/* calculate T for [C.3-crit]/[Cl(x=O)] = erfc(X/2 sqrt(DT)) */

concen_ratio = inv-erfc( criticaLCILconcentration /anion [Cl] .c[OJ);

T-crit = SQR((rebar-depth /L)/(2.O * concen-ratio));

chloride-failure-year = 3.171. 10-Os * T-crit * SQR(L)/D.;

if (head value > 0.0) { /* correct diffusion estimate to include darcy penetratio

darcy-velocity = (kO/p) * (p * g * head ,value)/sample_length;

A = 1.0/(4.0* Do * SQR(concen-raiio));

perm-depth = rebar-depth + 1.0/( 2.0 * A * darcy_velocity) – sqrt(rebar_depth /(A 

darcy-veiocity ) + 1.0/(4.0 * SQR(A * darcy-velocity)));

perndime = perm-depth /darcy-velocity;

chlon”de-failure-y ear = 3.171.10-08 * perm-time;

}

}

else chloride-failure_y ear = 1.00 10+lO;
}

42.

() {

crack.factor= 1.0;

}

else {

if (sampleJength – crack-depth >0.05*sample.length)

crack. factor =sampleJength/( sample_ lengih - crack-depth);

else crack.factor= 10.0;

}

AT= 1.0/(64.0 * perm-factor * crack-factor * anion [Cl].q * t90);

sulfa te-raie.ms = Youngs value * SQR(~. value) * anion [S04] .cIFIRST-CELL] * CE.vah6e 

Do/(roughness value * y.value *(1 – v.value)); /* m/see */

sulfate_rate = suifaie-rate-ms * (L/DO); /e AX/AT */

if (sulfate-rate > 0.0) suljate-failure-y ear = 3.171. 10-Os * (MaxDepth * L)/sulfat e-rate

else sulfaie-faihre-year = 1.0. 10+lO;

printf (“\’’Sulf at euf ailureu(yr) \“\t\t ~X8. Olf \n”, sulfate-failure-year);

/* calculate number of days until termination */

earliest-failure-day = 365. * UIN(su)fate-failure_y ear, chloride_ failure_year);

/* determine intervals between print outs */

print-day-interval = MIN( earliest-failure-day, MaxDay. value )/PRINT_• UTS;

print-day = print-day-interval;

Pout = (kO/(p * DO)) *p* g * head value;

dP = (kO/(p * DO)) * p * g $ (sampleJength/NUM-CELLS); /* change per cell */

for (k= O; k < NUMJ3URFACES; ki+) {

<’[k] = 190;

~’n[k] = 4;

qWl[k] = (p;

V[k] = POd * (real) (NUM_SURFACES– k – 1)/(NUM-SURFACES – 1);

}

( Initial estimate of pH 50)

The material parameters of the concrete must be determined based upon the information 

by the user. Missing pieces of information must be calculated from estab~shed relationships. 

is-default portion of the material data determines whether the user specified the value.

( Establish material parameters of concrete 42) s

if (Dm .is_default # TRUE) { /* calculate porosity */

if (WC.is-default s TRUE) wc.value = (loglO(10000. * Dw.value) + 9.84)/6.0;

if (km is-default s TRUE)

/* kinfty.value = 1.0E18*(0.8904+ .002525 *exp(15.07*wc.value)); */

k~.value = pow(10.0, (5.0* wc.value – 21.0));

}

if (km. is-default # TRUE) {

if (&is-default - TRUE) {

D~.value = 1.0010-04* POW(10.O,6.O* wc.value - 9.84);

if (we. is-default s TRUE) wc. value = (loglO(10000. * D~.value) + 9.84)/6.0;

}

}

if (wc.is-default # TRUE) { /* calculate diffusivity and porosity */
if (km, is-default s TRUE)

/* kinfty.value = 1.0E-18*(0.8904+.002525 *exp(15.07*wc.value)); */

km. value = pow(10.0, (5.0* wc. value – 21.0));

}

OiI= Dm.value/anion [Cl].Df;

if (t90 c DLIMIT) L90 = DLIMIT + 0.0001;

if (00 < DLIMIT + 0.07 * SQR(.18)) # = sqrt ((00 – DLIMIT)/.O7);

else ~ = .17326+ sqd(.03002 – (.05832 + DL1~lT – ~0)/1.8’7);

a = (3.2* wc.vaiue)/l.36-@ * (1+ 3.2* wc.value)/l.36;

~;;l~(l + 3.2* wc.value);

= 1000.0/4;

molesCH = cr * (0.61) * f * vSU7??ple/f@fI;

F = l.()/tie;

DO = Da.value;

k. = km.value;

litre[O] = 0.001* Vsample *+;

for (k= 1; k < NUM.SURFACES; k++) {

soLarray [Ca] [OH].s[k] = molesCH;

Jitre[k] = 0.001* Vsample *~;

}

for (k= O; k < NUM_CELLS; l?++) {

<[k] = t90;

@n[k] = ~;

}

This code is used in section 41.

43. The routine aci.211 estimates the strength, in psi, of the concrete in the individual 

putational elements. The estimate is based upon Table 5.3.4(a) of ACI 211.1-81 for air-entr

concrete. The strength can be adjusted to account for changes in porosity using:

( F’Unction declarations 23)+-

void aci-21J (void);

44. void aci-211 ( )

{

int i;

real @;

$ = pow(ezp(wc.value – 3.575), –2.6817);

for (i= O; i < NUMJSURFACES; i++) strength [i] = v;

}

45. The function yyp.rse sets the concentration in the concrete by putting the user-specified 

into the . c[l] cell. The internal concentration is established by copying this value into the rem

cells. However, this procedure k not performed for OH and H.
cation [~. moles [0] = cation [zj.c[O] 

litre[O];

cation [~. moles [l] = cation [i].c[l] * We[l];

for (k= 2; k < NUM.SURFACES; k++) {

caiion [i],c[k] = cation [i].c[l];

cation [i]. moles [k] = cation [i].c[k] * Iitre [k];

}

}-

for (j= 1; j < num.anions; jH) {

anion ~]. moles [0] = anion ~] .c[O] * litre [0];

anion ~]. moles [l] = anion ~].c[l] * Iiire[l];

for (k= 2; k < NUM.SURFACES; kt+) {

anion ~] .c[k] = anion ~] .c[l];

anion ~].moles[k] = anion [j].c[k] * litre[k];

}

}

This code is used in section 41.

46. Report material parameters to stdout.

( Preprocessor definitions 11)+=

#define GET_STAT(a) ((a. is.defatdt s TRUE) ? “DEFAULT” : “USER”)

47. ( Print parameters to stdout 47) s

printf (“\” THICKNESSu\’’\tu%8. 51f \t \“ (m) uuuuuu\”\t\’’~s\ n”\n”, sample-length,

GET-STAT(thickn ess));

pm”ntf(“\’’DIFFuuuuuu\’’\ tu%8. lle\t\” (m-2/s ee)\’’\t\’’!!s\ n”,n”, Do, GET_STAT(&));

pn”ntf (“\”PERMuuuuuu\ ’’\tu%8. lle\t\’’(m/sec) uu\’’\t%s\%\n”, \n”, 9.80 10+06 * Ito,

GET_STAT(kw));

pm”ntf(’’\’’WCuuuuuuuu\tu%8u%8 .51f\t\’’uuuuuuuuu \!$\t\~~~~\ll\ntC, wc.value, GET_STAT(W

printf(’’\’’HEADuuuuuu\tu%8u%8 .51f\t\’’(m)uuuuuu\’’\t%s\%\n”’\n” ,head.value,GET_STAT(he

prinif(’’SulfateuAttackuParameters :\n”);

printf(’’\’’YOUNGSuuuu\tu%8u%8. lle\t\’’(N/m”2)uu\’’\t%s\%\n”,\n”, Youngs.vahte,

GET_STAT(Youngs ));

prinif(’’\*’BETAuuuuuu\’’\tu%8 .lle\t\’’uuuuuuuuu \ll\t\l~~~\~$\n”, ~.value, GET_STAT(P))

prhtf(’’\”CE uuuuuuuu\”\tuX8. 51f \t\” (Mel/m-3)\’’\t\’’Xs\n”,n”, CE.vahe, GET_STAT(C

printf(’’\’’ROUGHNESSu\tu%8u%8 .51f\t\” uuuuuuuuu\’’\t\’’%s\n”\n” ,roughness.value,

GET-STAT(roughness));

pn”ntf(’’\’’GAMMAuuuuu\tu%8u%8 ,51f\tY’’(J/m-2 )uu\’’\t\’’%s\n”\n” ,7.value, GET_STAT(7))

printf(’’\’’POISSONuuu\tu%8u%8. 51f\t\’’uuuuuuuuu\’’\tXs\X\n”,\n”, ~.value,GET-STAT(v))i

printf(’’\n”);

if(MazDepth <sample-length)

printf(’’\’’DEPTHuuuuu\tuX8uX8 .51f\t\’’(m)uuuuuu\’’\tk\’\n”’\n” ,MazDepth, ’’USER”

if (rebar_depth <sampleJength)

printf(’’\’’REBARuuuuu\tu%8u%8 .51f\t\’’(m)uuuuuu\’’\t%s\%\n”’\n” ,rebar-depth,’’USER” 

prinif(’’\’’TIMEuuuuuu\tu~8u~8 .Olf\t\’’(day)uuuu\’’\t%s\%\n”’\n” ,MaxDay.value,

GETJ3TAT(MazDay ));

if (crack-widih.is-de fault ~TRUEAneuiral-axis_depth .is-default -TRUE) {

crack-depth =0.25 *sampleJength;
if (neutral_ axis. depth .isdefauli FALSE)

crack-depth = sample-length – neutral_ azis_depth. value;

if (crack.width.value > 0.0) printf(’’\’’CRACKu=u%8 .51fuATu%8.SlfuDEPTHu%8 .51f\’’\n”,

crack-widih.value , crack-spacing, crack-depth);

if (joini-permeabili@ >O) {

prinif (“\”JOINTuPERMu=uX8. ile\’’\n”, joini.permeability);

pm’ntf(’’\’’JOINT=X8 .51fuATuX8.51euUNTILu%8 .51f\’’\n’’,joinwidthth ,joint-spacing,

joini_lifetime/365);

}

This code is used in section 41.

48. Once the parameter DO has been established, thevalue ofqforeach ionmust reinitializ

( Initialize ion ‘eta’ parameter 48) -

for(i=O; i<num_cations; ii-t) caiion[i].q =cation[i].D~/DO;

for(j=O; j<num-anions; j+t) anion ~].q=anion~].D~ /DO;

This code is used in section 41.

49. The presence of cracks is reflected in K which a multiplicative adjustment for permea

In the absence of cracks x[i]=l. The quantities crack_width and crack-spacing are specified 

user.

Theequation forthepermeability ofcracks ofwidthw spaced a distancea apart is

(Calculatecrackandjointadjustmentopermeability49 ) =

crack-permeability = (SQR(crack-widih .vahw)/12);

X=o; /* zstarts frominside surface */

dx = sample-length/NUM_ CELLS;

for(k=NUM.CELLS–l; k~O; k--) {

z+=dz;

if (ZS crack-depth) {

K[k] = k.* (1. – (joint-width/joint-spacing ) – (crack-width value/crack-spacing ));

/*weighted length *I

K[k] +=joini.permeability * (joint-width/joint-spacing);

~[k] += crack-pemeability *(crack-width .value/crack-spacing);

}

else {

if (z–dz <crack-depth) { /* crack ends inthis cell*/

K[k] = k. * (1. – (joint-width/joint-spacing ) – (crack-width .vahe/crack-spacing));

/* weighted length *I

K[k] +=joint_permeability *(joint_width /joint.spacing);

tc[k] = (z – crack-depth )/~[k];

K[k]+= (dz-(z- crack-depth ))/(crack-pemeability *(crack-width .value/crack-spacin

@]/= dx;

K[k]= l,o/lc[k];

}

else {
1

;[k]=/c[k]/ko;

}

This code is used in section 41.

50. Need to make an initial estimate of the equilibrium concentrations and the pH of eac

( Initial estimate of pH 50) ~

/* make a rough guess for concentrations in cases of solids present */

for (k= O; k < NUM.SURFACES; k+t) {

for (i= O; i < num-cations; ii-t) {

for (j= O; j < num-anions; ji+) {

if (soLarray[i]~]. s[k] > 0.0) {

m = sol-array [i]~] .m;

n = sol.array[i]~].n;

if(m SO VnSO){

if (m S O) printf(’’m=Ou\n”);

if (n S O) printf(’’n=Ou\n”);

printf (“\n\tProblemuwithu%sXs\n” , cation [~. name, anion ~].name);

ezii(0);

} ““

moles-cation = pow (sol-array [il~] .ksp,

l./(m + n)) * pow(anion ~].valence/cation [i].valence, n/(n + m));

moles-anion = moles-cation * (cation [i]. valence/anion ~]. valence);

cation [i].c[k] = moles-cation;

anion ~] .c[k] = moles-anion;

cation [i]. moles [k] = cation [i].c[k] * litre [k];

anion ~]. moies[k] = anion [~.c[k] * Iitre[k];

}

}

}

} /* calculated [OH] and [H] concentrations */

for (k= O; k < NUM.SURFACES; k++) {

charge = 0.0;

for (j= 1; j < num-anions; j+t) charge += anion ~].valence * anion ~].c[k];

for (i = 1; i < num-cations; ii-t) charge –= cation [i] valence * cation [i].c[k];

anion [OH].c[k] = 0.5* (–charge + sqrt(SQR(charge) + 4.O * soi-array[lf][OH] .ksp));

cation [H] .c[k] = sol-array [H] [OH].ksp/ anion [OH].c[k];

anion [OH]moles [k] = anion [OH].c[k] * litre [k];

cation [H] .moles[k] = cation [H] .c[k] * iitre [k];

)

This code is used in section 41.

51. Inverse complementary error function. Find z given C such that C’= erfc(z).

(Function declarations 23) +=

real inv-erjc(real);
, , z[] {, ,,,, ,,,, , , 

1.60,1.80, 2.00}, y-table[] = {1.00, .888,.777,.671,.572, .480,.396,.322,.258, .203,.1

.090,.048,.024,.011, .005};

int i;

i = 16;

while (i >0 A y.table [--i] < C) ; /* find bounds */

if (z < 15) {

frac = (y-table [i] - C)/(y-tabie[il - y-table[i + l]);

z = z-tabie[i] + frac * (z-table[i + 1] – z-iable[i]);

}

else z = z-i!able[15]; /* if c >2.0 simply use last c value */

return x;

}
and Darcy flow of the pore solution due to hydrostatic head. The flux of ions due to both gr

in the ion concentration and to a volume average flow of pore solution is

j= –DVC + CU

where j is the ion flux, D is the diffusivity, c is the ion concentration, and u is the volume-ave

velocity of the pore solution. The time dependent change in concentration is the negative div

of the flux:

ac—= V. DVC– U. VC– CV” U

&

Given a hydrostatic pressure head on a vertical column of porous media, the pore volume-ave

flow VD is

vD = ‘:(vp–~g)

For the hydrostatic heads considered here, the body force term, pg is non-negligible. This e

can be cast into the more familiar form, assuming a constant density pore fluid, using a m

pressure potential:

~=p–pgz

This gives the more familiar Darcy equation

The volume-averaged velocity u can be related to the Darcy flow velocity:

vD = @_l

where ~ is the porosity.

Finally, the above equations can be combined to give

ac

% . Vc-1-Cv E=V”DVC% .-.$W

This equation gives the spatial and temporal behavior of the concentrations.

calculations a means is needed to update the hydrostatic pressure potential, +.

Continuity Equation

The temporal behavior of+ is calculated using the continuity equation:

ap—=–v. pv

at

where v is the intrinsic velocity of the pore solution. After averaging over

continuity equation becomes

ad

Z=–V” VD

To comple

the microstructure

This can be related back to pressure using the Darcy equation once again:
void ion.difltkm (void);

55. ion-diffusion.

void ion-diffusion ( )

{

int i, j, k;

real cation.jluz [10][NUM_SURFACES],anion$uz [1O][NUM-SURFACES],upper-limit = 5.0. 

lowerJimit = 5.0. 10+03, max_jiux-ratio;

real tmp-cation [10][NUM_SURFACES], tmp-anion [10][NUM_SURFACES];

AXIFIRST-CELL] = (real) (FIRST_CELL + 1) – sulfate-depth;

max$uz-ratio = 0.0;

for (i = O; i < num-cations; ii-t) {

for (k= FIRST_CELL + 1; k < NUM-SURFACES– 1; ki+) {

caiion-jiwr [~[k] = cation [i].q * (f[k] * (cation [i].c[k + 1] – cation [i].c[k])/AX[k] – f[

(cation [i].c[k] - cation [i].c[k - 1])/AX[k - 1])/(0.5 * (AX[k] + AX[k - l]));

caiion_fZuz [~[k] –= (l./#n[k]) * v[k] * ((cation [i].c[k + 1] – cation [i].c[k – 1])/(AX

AX[k - l]));

if (cation [i].c[k] > 0.0)

maz.jlwr.ratio = MAX(maz-fluz-ratio, AT * cationJ7uz [i][k]/cation [i].c[k]);

}

}

for (j = O; j < num-anions; j++) {

for (k= FIRST-CELL +1; k < NUM_SURFACES– 1; kH) {

anionJiuz ~][k] = anion ~].q * (~[k] * (anion fj].c[k + 1] – anion ~].c[k])/AX[k] – <

(anion ~].c[k] - anion ~].c[k - 1])/AX[k - 1])/(0.5 * (AX[k] + AX[k - l]));

anion_f7uz ~][k] -= (l./qYn[k]) * v[k] * ((anion ~].c[k + 1] - anion ~].c[k - 1])/(AX

AX[k - l]));

if (anion ~].c[k] > 0.0)

ma@uz_ratio = MAX(mazJluz-ratio, AT xcanion-fluz ~][k]/anion ~].c[k]);

}

}

if (mazJluz-ratio > upper-limit) AT *= upperJimit/maz.j7uz_ ratio;

if (max-jluz-raiio < lowerJimit ) AT *= upper_limit /maxJ7uz-ratio;

for (k = FIRST_CELL + 1; k < NUM_SURFACES– 1; k++) {

for (i = O; i < num-cations; ii+ ) cation [i],c[k] += AT* cation-flux [i][k];

for (j= O; j < num-anions; j+t) anion ~].c[k] += AT * anion-jluz ~][k];

}

k = NUM-SURFACES– 1;

for (i= O; i < num-cations; ii+) caiion [i].c[k] = cation [i].c[k – 1];

for (j = O; j < num-anions; j++) anion ~].c[k] = anion ~].c[k – 1];

for (i = O; i < num-cations; i+t)

for (k= FIRST-CELL; k < NUM_SURFACES;k-i+)

cation [~moles [k] = cation [i].c[k] * Iitre [k] * AX[k];

for (j = O; j < num-anions; j+t)

for (k = FIRST_CELL + 1; k < NUM-SURFACES; k++)

anion ~] moles [k] =

}

anion ~] .c[k] * Iitre [k]* AX[k];
go in to/out of solution. Determine the pH and adjust H and OH, accordingly.

( Function declarations 23)+=

void chemical-equilibrium (boolean);

57. Minimizing function. This is the function to minimize for the chemical. equilibrium routine.

The objective of chemical_equilibrium is to determinehow many molesofsaltshouldbeleached/pr

Therefore,minfunc must adjust for the pore volume.

The equation minfunc is minimizing the square of

[c+~]m[A+~]”-kv=o

where z is moles salt, m and n are stoichiometric ratios, and V is the pore solution volume.

real minfunc(real, int, int, int );

real minfunc(real z, int i, int j, int k)

{

real C, A;

int m, n;

m = sol-array[i]~].m;

n = soi-array[fl~].n;

C = cation [i].c[k];

A = anion ~].c[k];

return SQR(pow(C + m * c/(litre[k] * AX[k]), m) * POV(A + n * z/(iitre[~] * AX[~]),

n) – soLarrag [i][j].l.xp);

}

58. chemical-equilibriurn( ). Cycle throughsol-array and determineifanyionsshouldgo into/ou

ofsolutionbasedupon theconcentration, thevolubility constant,and thepresenceofsolidsalt.

#define TOL 1.0010-06

void chemicaLequilibn”um (boolean Chandelorosity)

{

int i, j, k, n, m, iterations;

real w-maz, 6;

real xa, zb, xc, fa, fb, fc, moles, imp;

real dsolid;

real Alit re; /* change in pore volume */

real old.cOH;

real charge, temp;

boolean NEG_FLAG;

for (k = FIRST.CELL + 1; k < NUM_SURFACES; k+t ) {

iterations = O;

do {

w_max = 0.0;

for (j = O; j < num-anions; jt+) {

for (i= O; i < num-cations; ;*) {

if(m(i=OAj~())) {
p p([][], () )p(][], () );

if ((sol-array [i]~].s[k] > 0.0) V (temp > soLarray[i]~].ksp)) {

za = –MIN(caiion[i].c[k], anion~].c[k])/10.O;

xb = —xa;

xc = 0.0;

mnbrak(&xa, &xb, &xc, &fa, &fb, &fc, minfunc, i, j, k);

imp = brent (xa, xb, xc, minfunc, TOL, &moles, i, j, k);

if (fabs(moles) > 0.05)

moles *= 0.10; /* no drastic changes */

if (sol-array [i]~].s[k] – moles < 0.0)

moles = sol_array [~~].s[k]; /* insufficient salt */

if (cation [i]. moles [k] + m * moles < 0.0)

moles = – cation [i]. moles [k]/m + 0.000001;

if (anion [j]. moles [k] + n * moles < 0.0)

moles = –anion [j]. moles [k]/n + 0.000001;

Alitre = 0.001 * moles* sol-array [i]~].molar.density/AX[ k];

if (Iitre[k] + Alitre < 0.0) {

Alitre = –iitre[k]; /* insufficient pore space */

moles = Alit re * AX[kJ * 1000 .O/sol-array [~~] molar-density;

}

cation [~. moles [k] += m * moles;

anion ~].moles[k] += n * moles;

if (Change-Porosity s TRUE) {

Jitre[k] += Alitre;

soLarray[i]~]. s[k] –= moles;

}

~ = fabs(moles);

w-maz = MAX(6,w-maz);

}

}

}-

} /* readjust the concentrations */

~’n[k] = 1000.0* litre[k] * AX[k]/ Vsample;

iitre[k] = MAX(litre[k],0.0001); /* avoid zero volume problems */

for (i = O; i < num-cations; i+t)

cation [i].c[k] = cation [i].moJes[kJ/(liire [k] * AX[k]);

for (j = O; j < num-anions; j++)

anion ~J.c[k] = anion ~].moles [k]/(litre[k] * AX[k]);

/* Determine [OH] and [ll]concentration */ /* store old value */

old-cOH = anion [OH].c[k];

charge = 0.0;

for (j = 1;j < num.anions; j++) charge += anion ~].valence * anion ~].c[k];

for (i = 1; i < num-cations; iH) charge -= caiion[i].oalence * cation [i].c[k];

anion [OH].c[k] = 0.5* (– charge + sqrt(SqR(charge) + 4.0* soLarray[H][OH].ksp))

anion [OH].moles[k] = anion [OH].c[k] * litre[k] * AX[k];

if (fabs((anion [OH].c[k] – old-cOH)/oid-cOH) > 1.00) {

w.max = 1.0;

anion [OH].c[k] = 0.5* (anion [OH].c[k] + oid-cOH);

anion [OH].moies[k] = anion [OH].c[k] * Iitre[k] * AX[k];
#n[FIRsT-cELL + 1]);

#n[FIRST_CELL] = MAX(gYn[FIRST_CELL], 0.0); /* extrapolate #n[NUM_CELL] value 

g5’n[NUM-CELLS]= 2.0 * #n[NUM_CELLS – 1] – #n[NUM_CELLS – 2];

/* extrapolate htre[FIRST_CELL] value */

litre [FIRsT_CELL] = 0.001* Vsample * #n[FIRST_CELL];

}

59. Minimization. The routine mnbrak brackets a minimum.

( Function declarations 23)+=

void mnbrak (real *, real *, real *, real *, real *, real *, real(*~unc)(real, int, int, int 

int, int);

#include %nbrak. c“

60. Parabolic approximation. The routine brent uses parabolic approximation to determ

minimum given three points bracketing the minimum.

( Function declarations 23)+=

real brent (real, real, real, real(xfvnc)(real, int, int, int ), real, real *, int, int, int );

#include “brent. c“
62. void advance-srdfate-front( )

{

int i, j, k, new$rst.cell;

sulfate-depth += sulfate.rate * AT;

if (sulfate-depth > (real) (FIRST_CELL + 1) v AXIFIRST_CELL] < 0.10) {

newJirst-cell = (int) (sulfate-depth + 0.4); /* penetration beyond next cell */

/* copy info from old first cell into new first cell */

for (k= FIRST_CELL + 1; k ~ newJirst-cell; ki+) {

for (i= O; i < num-cations; i~) cation[i].c[k] = catzon[i].c[FIRST_ CELL];

for (j = O; j < num-anions; j++) anion ~].c[k] = anion ~]. CIFIRST_CELL];

for (i= O; i < num-caiions; ii-t)

for (j = O; j < num-anions; j+t) sol.array[i]~]. s[k]= soi_array[i]~]. s[FIRST_C

}

W[new-first-cell] = WIFIRST_CELL];

FIRST-CELL = new-first-cell;

AXIFIRST_CELL] = (real) (FIRST-CELL+ 1) – sulfate.depth;

}

}
mass of chlorine, ~ is porosity, and 2.5 is the density of concrete.

( Global variables 4 ) +-

real C1.ratio [NUM_SURFACES];

64. ( Assess chloride penetration 64) ~

for (k = FIRST-CELL; k < NUM_SURFACES;k+t)

CLratio[k] = 35.5* anion [Ci].moles[k]/(2.5 * (1.0 – #’n[k]) * Vsample);

k = NUM.SURFACES– 1;

while (k ~ OA CLratio[k] < 0.0004) k--;

if (k < O) chloride-depth = 0.0; /* k at lower limit */

else {

if (k s FIRST_ CELL) { /* account for distance equals AX */

chloride_depth = AX IFIRST_CELL] * ( CLratio [k] - 0.0004)/( CLraiio [k] – Cl_ratio [k + 

chloride_depth += (real) FIRST_ CELL;

1

~lse chloride-depth = (real) k + ( CILratio [k] – 0.0004)/( Cl-ratio [k] – CILratio [k + l]);

}

This code is used in section 10.
The pressures are updated using the continuity equation:

ap

l%=-v”pv

where v is the intrinsic pore fluid velocity and p is the pore fluid density. An equation for 

material can be obtained from a volume average over a represet ative volume V:

1

/

OQ_ 1

7“& v“/

v .p’ d3x

For these equations, assume p is a constant since water is virtually incompressible. Rearrangin

simplifying the above equation gives:

——al

atvv J

pd3x=–V. ~

v“J

/)’ d3x

Since p is zero outside the pore volume, if VPrepresents only the pore volume then the above e

simplifies to

84 4

%=–v”~ I“,vd3x

which finally gives

*=_v ~u=_v ‘D

at”

where ‘D is the Darcy velocity. Substituting for the Darcy velocity gives

—=

a~ V.;VP

at

( Function declarations 23)+=

void updat e_pressures (void);

66. update_pressures.

void update-pressures( )

{

int k, iteration;

real perm [NUM_CELLS], tmp [NUM-SURFACES], tol= 0.0005, c;

for (k = FIRST_CELL; k < NUM_CELLS;kt+) perm [k] = K[k] * CUB(C[k]/Oo);

WIFIRST.CELL] = Pod + (NUM_CELLS– 1 + AX IFIRST.CELL]) * dP;

W[NUM-SURFACES – 1] = O;

iteration = O;

do {

for (k= FIRST_CELL + 1; k < NUM.SURFACES– 1; k+t) {

tmp[k] = perm[k] * Q[k + 1]/AX[k] + perm[k – 1]* w[k – 1]/AX[k – 1];

/* **** tmp[k] -= 0.5*( DX[k-l]+DX[k])* (PHIn[k]-PHIn[k-l] )/DT; *** */

tmp[k] *= l. O/(perm[k]/AX[k] + perm[k – 1]/AX[k – 11);

}

c = 0,0;

for (k= FIRST_CELL + 1; k < NUM-SURFACES – 1; k+t) {

if (w[k] > 0.0) c = MAX(c,fabs(~[k] – tmp[k])/!l[k]);
v[k] = –/c[k] * CUB(~[k]/Oo) * (W[k + 1] – W[k])/AX[k];

}

67. Copy the new porosities into the n – 1 values.

( Update the porosities 67) s

for (k= FIRST-CELL; k < NUW3URFACES; k~) #’n-l [k] = qY”[k];

This code is used in section 9.
+n[k] = 0.5 * (#’n[k] + fb’n[k+ l]); /* interpolate phi’ values */

d = DLIMIT + 0.07 * SQR(&[k]);

if (@n[k] > 0,180) d += 1,8* SQR(qV[k] – .180);

if (@n[k] > ~) /* leaching */

<[k] = d, + 5.0 * (19- t9,);

else /* precipitation */

f[k] = O;

}-

<’[FIRST_ CELL] = <[FIRST_ CELL]– (AXIFIRST_CELL] /(1. O+ AXIFIRST-CELL])) * (t[FIRST

1] – ~[FIRST_CELL]);

for (k= FIRST-CELL +1; k < NUM_SURFACES– 1; k-t+) (’[k]= 0.5* (f[~] + ([k – 1]);

<’[NuM_suRFAcEs – 1] = 1.5* &[NUM_CELLS– 1] – 0.5 * &[NUM_CELLS- 2];

This code is used in section 9.
( Global variables 4 ) +=

real Time = O; /* dimensionless cumulative time */

real Day = O; /* cumulative time (day) */

70. ( Advance global clock TO) =

Time += AT;

Day = Time * SQR(~)/(DO * 86400);

This code is used in section 9.

71. Check time dependent variables.

( Check time dependencies 71) ❑

if (Day > joint.lifetime A joint. is-specijied ❑ TRUE) {

joint-change-flag = TRUE;

}

This code is used in section 10.

72. (Assess simulationtermination72) s

if (Day > MaxDay value) termination = TRUE;

if (sulfate-depth > MaxDepth ) termination = TRUE;

if (( chlom”de-depth * L) > rebar-depth ) termination = TRUE;

This code is used hasection 10.
74.

g p

Print the state of the system for the first 4 concentration surfaces.

( Print initial system state 74) ~

printf (“\n\nInitialustateUo fusystem:\n\n”);

printf(’’\’’uIONu\t\\EXTERNAL\N\t\”’\t\” INTERNAL\’’\n”);

for(i=O; i<num.caticms; ii+) {

if (cation [i].CIFIRST_CELL]# Ov cation [i].C[FIRST_CELL + 1] # O) {

printf (“\”%4s: \“\t”,cation [~. name);

for (k= FIRST_ CELL; k < FIRST_CELL + 2; ktl-) print} (”%8. 51f \t”, cafion[i].c[k]);

printf(’’\n”);

}

}-

for (j = O; j < num.anions; ji+) {

if (anion ~]. CIFIRST_CELL] # Ov anion ~]. CIFIRST_CELL + 1] # O) {

printf (“\”%4s: \“\t”,anion ~].name);

for (k= FIRST-CELL; k < FIRST-CELL+ 2; kH) prinif (“%8. Slf\t”, anion ~].c[k]);

printf(’’\n”);

}

}

pm”ntf (“\”UupH: \“\t”);

for (k = FIRST_ CELL; k < FIRST-CELL+ 2; ki+) prinif (“%8. 51f \t”, –logl O(cation [H].c[

pn”ntf(“\n”);

printf (“\n”);

This code is used in section 7.

75. Print header for permeability, sulfate penetration, and chloride penetration depths.

( Print header 75)=

printf(’’\’’Day \t\\L\’’L\’’\t\ ’’Kuu\''\t\''Duuuu\'' \t\''S04\''\t\'' Cl\'' \t\''Flux\''\tuu\ ''pH

pn"nif(''\''uuu\ ''\t\''m\''\t\' 'm/s\ ''\t\''m-2/s\'' \t\''%u\''\t\'' mu\'' \t\''ti/dy/m2\' '\n'');

printf (“\n”);

This code is used in section 2.

76. Print intermediate results for permeability, sulfate and chloride penetration depths. The 

values of permeability and diffusivity must be calculated here.

(Print intermediate results 76) s

if (Day S 0.0 V Day > print-day) {

print-day += print-day-interval;

printf (“%ldu\t%6. 31f”, (long) Day, (AX IFIRST-CELL] + NUM_CELLS– FIRST_CELL – 1

kB = AX IFIRST_CELL] * CUB(OO/~[FIRST_CELL] )/tCIFIRST_CELL];

11~ = AXIFIRST_CELL] * (t10/<[FIRsT_cELL]);

for (k= FIRST_CELL + 1; k c NUM_CELLS; k+t) {

kB += CUB(Oo/~[k])/~[k];

DB += ticI/~[k];

}

kB = (AXIFIRST_CELL] + NuM-cELL5 – FIRST-CELL – 1) * kJk~;
}

This

77.

p

prinij(’’\t%6. ile”, DB);

prinif (“\t%6. 31f \t%6. 31 f”, sulfate.depih * L, chloride.depih * L);

printf(”\t%6. 31f”, 1000. * 86.4* v[NUM_CELLS– 1] * Do * Vsa?nple/SQR(L));

printf (“\t%6.Ilf\n”,–Ioglo(cation [H]. C[NUM_SURFACES– l]));

code is used in section 2.

Print termination information to include the day of termination and the depth of chlori

sulfate penetration.

( Print termination information 77) S

prinif (“\n\n”);

if (Day > MaxDay value A joinLchange_j7ag s FALSE) printf (“ExceededuTIMEulimit. \

if (Day > MazDay. value A joinLchangeJ7ag s TRUE)

printf (“ExceededuJOINTuLIFETIME. \n\n”);

if (sulfate-depth > MaxDepth ) pn”ntf (“Sulf at eUFailure. \n\n” );

if (chloride. depth * L > rebar-depth ) printf (“Chlori.deuFailure. \n\n” );

printf (“\”T\’’\t%8. 31f \n”, Time);

prinif (“\”Day\’’\t%8. Olf \n”, Day);

printf (“\”S04u(m)\’’\t %8. 31f\n”, sulfate-depth * L);

prinif (“\”Clu(m)\’’\t%8. 31f \n”, chioride-depih * L);

prinif(’’\n\n”);

This code is used in section 2.

78. Print final system state for any desired ion or salt.

( Print desired output 78) z

printf (“FinalUSyst emustat e: \n\n” );

prinif(’’\’’L(m) u\t”);”);

prin~f(’’\’’Psiu\ t”);”);

pn”nif(’’\’’vDu\ t”);”);

prin~f(’’\’’xiu\ t”);”);

prinif(’*\’’phi\ ’’\t”);

printf (“\’’pH\’’u\t “);

pn”ntf(“\”f c\ ’’u\t “);

for (i = O; i < num_cations; i+i-)

if (cation [i]. output-jlag s TRUE)

for (j= O; j < num_anions; j+)

if (anion ~]. output-flag z TRUE)

for (i= O; i < num-caiions; i-l+)

printf(’’\’’%su\t”,t”, cation [i].name);

prinif (“\’’%s\’’u\t”, anion ~].name);

for (j = O; j < num-anions; j++)

if (soLarray [~ ~]. output-flag = TRUE)

printf(’’\’’%s\su\t”,t”, cation [i].name, anion ~].name);

printf(’’\n”);

for (k= (); k < NUM.SURFACES;kt+) {

printf (“%7.41f\t*’,L * k);

printf (0’%7.31f\t”,W’[k]);

if (k < WM.CELLS) {

printf(’%$. 31f u\t”, v[k]);
p(\, []);

printf (“%6. 31f u\t”, 14.0 + loglo(anion [o].c[k]));

printf(”%6. Olf ~\t”, strengih[k] * (log (#n[k])/lo9(@)));

for (i = O; i < num.cations; i+t)

if (cation [i].ouiput.flag s TRUE) printf(”%6. 41f ~\t”, cation [i].c[k]);

for (j = O; j < num-anions; j++)

if (anion ~]. outptitJiag s TRUE) printf(”%6. 41f ~\t”, anion ~].c[k]);

for (i = O; i < num-caiions; ii+)

for (j= O; j < num-anions; jt+)

if (sol-array [i]~]. output-flag s TRUE) printf(”%6. 31f u\t”, sol-array [i]~].s[k]);

printf(’’\n”);

}

elapse d.time = biostime(CET_TIME, OL] – elapsed-time;

printf (“\’’Elapsedutimeu( see) : \“\t~6. If \n”, elapsed.time/CLK_TCK);

This code is used in sections 2, 24, and 26.