Numerical Treatment of a Diffusion-Convection-Evaporation

Model for Droplets

FRITJOF NILSSON

Master’s Degree Project Stockholm, Sweden 2004

TRITA-NA-E04090

Numerisk analys och datalogi Department of Numerical Analysis KTH and Computer Science 100 44 Stockholm Royal Institute of Technology SE-100 44 Stockholm, Sweden

FRITJOF NILSSON

TRITA-NA-E04090

Master’s Thesis in Numerical Analysis (20 credits) at the School of Engineering Physics,

Royal Institute of Technology year 2004 Supervisor at Nada was Michael Hanke

Examiner was Axel Ruhe

Numerical Treatment

of a Diffusion-Convection-Evaporation Model for Droplets

Abstract

A FOI model for the simulation of moving droplets has been enhanced. Thethree consecutive phases considered in the FOI model of a falling droplet is (1)impact, (2) capillary spreading and sorption and (3) evaporation and redistri-bution. Focus of this work is on the third part, although phase one and twohas been integrated in the code. Improved use of numerical methods has de-creased the discretization errors and increased the speed of the one-dimensionalmodel with at least a factor 100. Furthermore, an extension of the model totwo dimensions, heigth and radius, has been developed and compared with windtunnel experiments. The mean square error of the 2D model is half the value ofthe original 1D model. The improved 1D model gives almost as good results.

The Partial Differential Equation (PDE) of phase 3 is a diffusion-convection-reaction equation, which is semi-discretized with finite differences (FDM) incylindrical coordinates. Boundary conditions are based on the use of infinitedepth and radius, cylindrical symmetry and evaporational flux. The resultingordinary differential equation (ODE) is discretized with variable step size numer-ical differential formulas (NDFs). The nonlinear systems of equations are solvedwith Newtons method plus either the LU solver UMFPACK or Krylow typemethods. An userfriendly implementation including the generation of graphs ismade in Matlab.

Sammanfattning

En FOI-modell för simulering av enstaka vätskedroppars rörelse har förbättrats.FOI-modellen behandlar för en nedfallande droppe tre på varandra följande fy-sikaliska faser; (1) kollision, (2) kapillär utspridning och inträngning i underlagetsamt (3) omfördelning och avdunstning. Detta arbete fokuserar på den tredjefasen, även om de föregående också är integrerade i koden. Användande av för-bättrade numeriska metoder minskar diskretiseringsfelen och ökar hastighetenhos den endimensionella modellen med åtminstone en faktor 100. Därtill harmodellen utvidgats till två rumsdimensioner; höjd och radie. Resultaten harjämförts med vindtunnelförsök. Standardavvikelsen för 2D modellen är cirkahälften så stor som för den ursprungliga 1D modellen. Den förbättrade 1D mo-dellen uppvisar nästan lika goda värden.

Fas tre kan beskrivas som en partiell differential ekvation (PDE) av diffusions-konvektionsreaktionstyp. Denna diskretiseras med finita differenser (FDM) icylindriska koordinater. Randvillkoren skapas med utgångspunkt från avdunst-ning, cylindrisk symmetri och oändligt djup och radie. Den resulterande ordinäradifferential ekvationen (ODE) diskretiseras med numeriska differentieringsform-ler (NDF) med variabla steg. De ickelinjära ekvationssystem som uppstår lösesmed Newtons metod i kombination med antingen LU-lösaren UMFPACK el-ler Krylowmetoder. Implementationen i Matlab är användarvänlig och rik påvisualiseringar.

Zusammenfassung

Ein Model des Schwedischen Verteidigungsforschungsinstituts ist verbessert wor-den. Das Model simuliert Tropfen von Flüssigkeit in drei aufeinander folgendenPhasen: (1) Kollision, (2), kapillares Ausbreiten und Eindringen der Boden und(3) Verdunstung und Verteilung in der Erde. Diese Arbeit konzentriert auf demdritten Teil, aber die anderen sind auch implementiert worden. Anwendungenvon modernen numerischen Methoden haben die Diskretisierungsfehler reduziertund die Geschwindigkeit des Models mit einen Faktor 100 erhöht. Außerdem istein Model in zwei Dimensionen, Höhe und Radius, entwickelt.

Die dritte Phase wird durch eine partielle Differentialgleichung von Diffusions-Konvektions-Reaktions-Typ beschrieben. Diese ist mit finiten Differenzen (FDM)in zylindrische Koordinaten diskretisiert. Die Randbedingungen sind durch Ver-dunstung, zylindrische Symmetrie und unbegrenztes Radius und Tiefe geschaf-fen. Die gewöhnliche Differentialgleichung (ODE) die entsteht, ist mit numeri-schen Differenzenformeln (NDF) mit variablen Schritten behandelt werden. Dieresultierenden nichtlinearen Gleichungssysteme sind mit der Newton Methodein Kombination mit direkten oder iterativen Methoden gelöst werden.

Contents

1 Introduction 3

2 Model and Methods 42.1 Physical Description . . . . . . . . . . . . . . . . . . . . . . . . . 4

2.1.1 Phase 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.1.2 Phase 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.1.3 Phase 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.2 PDE in Cylindrical Coordinates . . . . . . . . . . . . . . . . . . . 62.3 Spatial Discretization . . . . . . . . . . . . . . . . . . . . . . . . 7

2.3.1 Choice of Method . . . . . . . . . . . . . . . . . . . . . . 72.3.2 The Discrete Product Rule . . . . . . . . . . . . . . . . . 82.3.3 Spatial Discretization Performed . . . . . . . . . . . . . . 8

2.4 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . 112.4.1 Boundary Conditions in Two Dimensions . . . . . . . . . 112.4.2 Boundary Conditions in One dimension . . . . . . . . . . 15

2.5 Matrix Generation . . . . . . . . . . . . . . . . . . . . . . . . . . 152.5.1 Lexicographical Ordering . . . . . . . . . . . . . . . . . . 162.5.2 Matrix Generation in Two Dimensions . . . . . . . . . . . 162.5.3 Matrix Generation in One Dimension . . . . . . . . . . . 162.5.4 Compressed Column Storage (CCS) Format . . . . . . . . 17

2.6 Initial Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.6.1 Initial Conditions Two Dimensions . . . . . . . . . . . . . 172.6.2 Initial Conditions One Dimension . . . . . . . . . . . . . . 19

2.7 Remaining Fraction . . . . . . . . . . . . . . . . . . . . . . . . . . 192.8 Evaporated Fraction . . . . . . . . . . . . . . . . . . . . . . . . . 192.9 Jacobian of f . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202.10 Time-Discretization . . . . . . . . . . . . . . . . . . . . . . . . . 24

2.10.1 The trapetzoid formula (Cranck-Nicholson) . . . . . . . . 252.10.2 Backward Differentiation Formulas (BDFs) . . . . . . . . 252.10.3 Numerical Differentiation Formulas (NDFs) . . . . . . . . 27

2.11 Solving Systems of Equations . . . . . . . . . . . . . . . . . . . . 282.11.1 Direct Solvers (GAUSS, LU, CHOL, QR ) . . . . . . . . . 282.11.2 Iterative Solvers . . . . . . . . . . . . . . . . . . . . . . . 302.11.3 Stationary Iterative Methods . . . . . . . . . . . . . . . . 302.11.4 Nonstationary Iterative Methods . . . . . . . . . . . . . . 312.11.5 Preconditioning . . . . . . . . . . . . . . . . . . . . . . . 37

1

3 Implementation 393.1 Choice of Programming Language . . . . . . . . . . . . . . . . . 393.2 General Description of the Implementation . . . . . . . . . . . . 403.3 How to Use the Program . . . . . . . . . . . . . . . . . . . . . . . 403.4 Structure of the Files . . . . . . . . . . . . . . . . . . . . . . . . . 413.5 Description of the Subprograms . . . . . . . . . . . . . . . . . . . 42

4 Results and comparison with experimental data 454.1 Results from the Examination of Matrices . . . . . . . . . . . . . 454.2 Results from the Using of Conservation Laws . . . . . . . . . . . 454.3 Results from the Comparing of Different Solvers . . . . . . . . . 464.4 Results from the Comparison with Experimental Data . . . . . . 49

4.4.1 Wind Tunnel Tests . . . . . . . . . . . . . . . . . . . . . . 494.4.2 Actual Comparisons . . . . . . . . . . . . . . . . . . . . . 50

4.5 General Results of the Project . . . . . . . . . . . . . . . . . . . 50

5 Discussion 54

6 Nomenclature 56

2

Chapter 1

Introduction

The aim of this graduation thesis is to enhance the numerical treatment of anexisting FOI model by Winter et al [27], which describes the performance ofdroplets hitting the ground. Models of this kind can for instance be usefulfor taking preventive measures against pollution, chemical contamination andchemical weapons.

The model is divided into three consecutive steps. The first describes the spread-out caused by kinetic energy and the second describes the short phase when aspherical cap of liquid is still on the ground and the soil is more than saturated.The third phase describes the situation when all liquid has merged into theground and evaporation must be considered. This work focus heavily on thethird phase, even if an Matlab implementation of the first two phases also ismade. Experimental evaporation curves from Prince Mauritz Laboratory, Oese-burg [18], are used as reference data.

Several models have been suggested for describing the motion of liquid in theground, but none has previously been able to describe all three phases properly.Many uses evaporation models based on data for lakes on droplets even thoughother athmospheric processes dominates on the micro scale. The three step FOImodel is promising but has had poor numerical treatment of the equations res-ulting in long execution times and large discretization errors.

The Swedish Defence Research FOI has requested this report in order to improvetheir model, which later on will be incorporated in a larger structure togetherwith models for deposition and atmospherical spreading. FOIs requirementswere: Improve the speed and accuracy of the existing one-dimensional Fortranimplementation, construct eventually a 2D model, design an intuitive user in-terface for the code, generate figures of the results and document the work inan exhaustive report. The goals are achieved completely.

3

Chapter 2

Model and Methods

2.1 Physical DescriptionTransport of non-solid materials in porous medias are caused by physical andchemical phenomena, see for instance Sadhal [23] or Cussler [5]. Liquid chemic-als can be transported by the forces of gravitation, wind and pressure or by cap-pilary forces (diffusion) in porous medias. Gaseous medias can be transportedby diffusion, (movement towards lower concentrations), or convection (termo-dynamic movement caused by differences in temperature or pressure). The FOImodel is divided into three phases. Immediately after impact, when all liquid ison the surface (phase 1), the droplet is affected by wind and gravitation. Whenthe droplet has started to penetrate the surface (phase 2), capilary diffusiondominates. When all liquid has penetrated the surface (phase 3), convection,evaporation, liquid and gaseous diffusion and chemical reactions dominates.

2.1.1 Phase 1When a droplet strikes the ground, an immediate deformation of the droplettakes place due to the impact. The phase is over when all kinetic energy hastransformed into other forms of energy. If v0 is impact velocity (m s−1), γ issurface tension (Nm−1), µ viscosity (kg m−1 s−1), ρ density (kg m−1), re radiusafter impact(m) and r0 radius before impact, the radius is written:

✖✕✗✔

� ✏ ☛✡

✟✠ ✝ ✆

Phase 0 Phase 1 Phase 2 Phase 3

Figure 2.1: A survey of the phases of the falling droplet.

4

re = r0 ·(

re,min

r0

)6

+ 0.1457

√ρ3cr

30v

40

µ2γ

1/6

(2.1)

The radius re,min ≈ 1.61r0 (m) is the smallest possible value of re for v0 = 0,which corresponds to a contact angle of 60 degrees.

2.1.2 Phase 2Immediately after the end of phase 1, a spherical cap of liquid, which has notbeen absorbed by the soil, remains on the ground. Successively it will spreadout on the surface and penetrate the ground until everything is absorbed. Thecombination of physical formulas, approximations and the use of mass conser-vation results in an ordinary differential equation (ODE) of the form r = f(r, t)which reads:

dr

dt=

σ2

2

r −(r10e + 0.7(V0 − 0.86fsβσ

√t · πr2)3 γt

µ

)1/10(2.2)

The term V0 is the volume of the original droplet, fs ≈ φ− θw0 is the saturatedconcentration, φ is total soil porosity, θw0 initial soil water content, t time andzf = σ

√t penetration depth at the symmetry axis of the droplet. The part of

the droplet above the ground is treated as a half-sphere and the absorbed partas a cylinder with height βzf = 0.8zf and concentration 0.86fs. The phase lastsuntil no liquid remains above the surface. This condition can be formulated as:

r ≥√

4πr30

3· 10.86πfsβσ

√t

2.1.3 Phase 3The third phase, on which this project is focused, starts when no liquid remainsabove the ground. When all kinds of relevant transport phenomena in theground are considered, the evolution of the concentration c is governed by thefollowing parabolic Partial Differental Equation (PDE) in which �F (kg m−2 s−1)is the mass flux. Neither A, B nor F can be written in closed form.

c = −∇ · �F − εc = ∇ · (A∇c) + ∇ · ( �Bc) − εc (2.3)

A (m2 s−1) is a (dominating) effective diffusion coefficient, �B (ms−1) is aneffective convection velocity and ε a reaction coefficient. Ds

g,Dsl ,D

slc are effective

soil-gas, soil-dissolved chemical and soil-liquid diffusion coefficients. R−1η , R−1

g

and R−1L are partitioning factors while ρb is the density of the bulk.

A = DsgR

−1g + ρcD

slcR

−1η + Ds

l R−1l

�B = Dsg∇R−1

g + ρcDslc∇R−1

η + Dsl∇R−1

l − �vxyzR−1l (2.4)

5

The numeric solution of this PDE is the main goal of this project. Boundaryconditions and initial conditions are handled in section 2.4 and 2.6.

The derivation of the factors A and B can be found in Winter et al [27] butthe final algorithm is stated. Observe the fixed point iteration with respect ofthe liquid chemical concentration η. The factor a(η) is air filled porosity andfη(η) and R−1

η (η) are partitioning functions. KH is a Henry law constant, KD alinear absorption constant, ρc the density of chemical, εit a breaking condition,cgsat the saturated gas concentration, θw0 soil water content and φ total soilporosity. Dair

g and Dwaterl are diffusion coefficients for the chemical in air and

water.

Algorithm DCR (Diffusion-Convection-Reaction)ηnew = εit, p = false

while p == false

ηold = ηnew

a = φ− θw0 − ηnew

fη = cgsat/(max(ηnew, ηcrit))

R−1η = 1/(ρc + fη(a + θw0 +

ρbKD

KH))

ηold = ηnew

ηnew = c ·R−1η

ηdiff = |(ηnew − ηold)/ηold|if (ηdiff ≥ εit), p = true, endif

endwhile

Dsg =

a10/3 ·Dairg

φ2

Dsl =

θ10/3w ·Dwater

l

φ2

Dslc = 1.07 · 10−3 · σ

2

18·(e

8.1η(φ−θ) − δk

)R−1

g = R−1ηnew

· fηR−1

L = R−1ηnew

· fη/KH

Discretize and solve the PDE

O

2.2 PDE in Cylindrical CoordinatesSince the model problem has cylindrical symmetry, cylindrical coordinates willbe used on the PDE. If the definition ∂x ≡ ∂

∂x is made, according to the formulacollection Beta [22], the formulas for differentiation in cylindrical coordinatesare:

6

∇u = ∂r(u)er +1r∂r(u)eθ + ∂z(u)ez

∇ · �u =1r∂r(rur) +

1r∂θ(uθ) + ∂z(uz)

These identities applied to the PDE leads to the following relationships:

�B = (Dsg∂r(R−1

g ) + ρcDslc∂r(R−1

η ) + Dsl ∂r(R−1

l ) − vrR−1l )er +

1r(Ds

g∂θ(R−1g ) + ρcD

slc∂θ(R−1

η ) + Dsl ∂θ(R−1

l ) − vθRl)eθ +

(Dsg∂z(R−1

g ) + ρcDslc∂z(R−1

η ) + Dsl ∂z(Rl) − vzR

−1l )ez =

Brer +1rBθeθ + Bzez

c = −∇ · �F − εc = ∇ · (A∇(c) + �Bc) − εc =

∇ · (A∂r(c)er +A

r∂θ(c)eθ + A∂z(c)ez +

(Brer +1rBθeθ + Bzez)c) − εc =

1r∂r(rA∂r(c)) +

1r2

∂θ(A∂θ(c)) + ∂z(A∂z(c)) +

1r∂r(rcBrer) +

1r2

∂θ(cBθeθ) + ∂z(cBzez) − εc

Symmetry along the z-axis gives that: ∂θu = ∂θ · �u = 0, which leads to the 2Dformulation of the PDE in cylindrical coordinates:

c =1r∂r(rA∂r(c)) + ∂z(A∂z(c)) +

1r∂r(rcBr) + ∂z(cBz) − εc (2.5)

If the concentration is assumed to be independent of the radius, the 1D formu-lation is derived.

c = ∂z(A∂z(c)) + ∂z(cBz) − εc

2.3 Spatial Discretization

2.3.1 Choice of MethodThe discretization of the spatial-coordinate dependent factors can either be donewith Finite Element Methods (FEM), pseudo-spectral methods, compact finitedifferences or classical Finite Difference Methods (FDM). While the computa-tional grid is rectangular and not smooth or polygonal, the standard classicalfinite differences (FDM) were chosen in agreement with Carlsund [4].

7

The computational domain in 3D is a cylinder (r, z, θ), in 2D a rectangle (r,z) and in 1D a line z. The continous domain is split with the folllowing FDMdiscretizaton with N equally spaced gridpoins in the vertical direction and M inthe horizontal.

zi = (i− 1) · ∆z, i = 1..N, ∆z = −|∆z| = |zmax|/(N − 1) = (zi+1 − zi)

ri = (i− 1) · ∆r, j = 1..M, ∆r =rmax

(M − 1)= (rj+1 − rj)

Centered discretization with second order accuracy are used for the approxim-ation of derivatives. ∂x(xi) ≈ (xi+1/2−xi−1/2)

∆x . Average values are denoted byxi = (xi+1/2 + xi−1/2)/2 and xi±1/2 = (xi + xi±1)/2.

2.3.2 The Discrete Product RuleAccording to Weinan E [26], the discrete product rule reads:

∂x(fg) = f∂x(g) + g∂x(f)

With the centered difference scheme the discrete product rule becomes:

∂x(f · g)i = f∂x(g)i + g∂x(g)i =

(fi+1/2 + fi−1/2

2)(gi+1/2 − gi−1/2

∆x) + (

gi+1/2 + gi−1/2

2)(fi+1/2 − fi−1/2

∆x) =

(2fi+1/2 · gi+1/2 − 2fi−1/2 · gi−1/2)2∆x

+

(fi+1/2 · gi−1/2 − fi+1/2 · gi−1/2)2∆x

+(fi−1/2 · gi+1/2 − fi−1/2 · gi+1/2)

2∆x=

(fi+1/2 · gi+1/2 − fi−1/2 · gi−1/2)∆x

, i = 1..n

The three term extension is:

∂x(f · g · w)i =(fi+1/2 · gi+1/2 · wi+1/2 − fi−1/2 · gi−1/2 · wi−1/2)

∆x

2.3.3 Spatial Discretization PerformedThe last few results are enough for determinating the shapes of the sought FDMstencils.

Nota bene that in the 2D model, A, B, Dsg, R−1

g , Dslc, R

−1η , Ds

l , R−1l and c

are dependent of z(i) and r(j). In order to keep the equations as simple as pos-sible, at this stage only the index, which is currently worked with, is denoted.For example: Bi±1 ≡ Bi±1,j and ∂rcj ≡ ∂rci,j . If the radial parts are ignoredthe formulas for the 1D model are obtained.

8

Initially, the discretizations of A and B with respect to z(i) are done. Toachieve the corresponding two radial equations, just exchange every z with rand every i with j.

Ai = Dsg(i)Rg(i) + ρcD

slc(i)Rη(i) + Ds

l (i)Rl(i) (2.6)

Bi±1/2 = ± 12∆z

· ((Dsg(i± 1) + Ds

g(i)) · (R−1g (i± 1) −R−1

g (i)) +

ρc · (Dslc(i± 1) + Ds

lc(i)) · (R−1η (i± 1) −R−1

η (i)) +

(Dsl (i± 1) + Ds

l (i)) · (R−1l (i± 1) −R−1

l (i))) − vzR−1l (i± 1)(2.7)

Applying the relations above on the z-dependent terms in the model-problemwill result in a three-point stencil. From the section of cylindrical coordinates:

c =1r∂r(rA∂r(c)) + ∂z(A∂z(c)) +

1r∂r(rcBr) + ∂z(cBz) − εc

The spatial differentiation of the diffusion, convection and degradation terms inthe z-direction are deduced below.

∂z(Ai · ∂z(ci)) ≈∂z(Ai · (ci+1/2 − ci−1/2))

∆z=

∂z(Ai · ci+1/2)∆z

− ∂z(Ai · ci−1/2)∆z

≈(Ai+1/2 · ci+1 −Ai−1/2 · ci)

(∆z)2− (Ai+1/2 · ci −Ai−1/2 · ci−1)

(∆z)2=

((Ai+1 + Ai) · ci+1 − (Ai+1 + Ai) · ci)2(∆z)2

+((Ai + Ai−1) · ci−1 − (Ai + Ai−1) · ci)

2(∆z)2=

12(∆z)2

· ((Ai+1 + Ai) · ci+1 + (Ai + Ai−1) · ci−1 − (Ai+1 + 2Ai + Ai−1) · ci =

Ka1ci−1 + Ka2ci + Ka3ci+1

∂z(ci ·Bi) ≈ (Bi+1/2 · ci+1/2 −Bi−1/2 · ci−1/2)∆z

≈(Bi+1/2 · ci+1 + Bi+1/2 · ci

2∆z− (Bi−1/2 · ci + Bi−1/2 · ci−1

2∆z=

(Bi+1/2 · ci+1 + (Bi+1/2 −Bi−1/2) · ci −Bi−1/2 · ci−1)2∆z

=

Kb1ci−1 + Kb2ci + Kb3ci+1

−εc = −εci = Kc2ci

If K1 = Ka1 + Kb1, K2 = Ka2 + Kb2 + Kc2 and K3 = Ka3 + Kb3 the vertical

9

componenents of the 2D five-point stencil or the 1D three-point stencil can bedenoted: K1ci−1 + K2ci + K3ci+1. Analog treatment of the radial parts gives:

1r∂r(rjAj∂r(cj))) ≈

∂r(rjAj · (cj+1/2 − cj−1/2))rj∆r

=∂r(rjAj · cj+1/2)

rj∆r− ∂r(rjAj · cj−1/2)

rj∆r≈

(Aj+1/2rj+1/2 · cj+1 −Aj−1/2rj−1/2 · cj)rj(∆r)2

− (Aj+1/2rj+1/2 · cj −Aj−1/2rj−1/2 · cj−1)rj(∆r)2

≈1

4rj(∆r)2· ((Aj+1 + Aj)(rj+1 + rj) · cj+1 + (Aj−1 + Aj)(rj−1 + rj) · cj−1

− ((Aj−1 + Aj)(rj−1 + rj) + (Aj+1 + Aj)(rj+1 + rj)) · cj) =Ha1cj−1 + Ha2cj + Ha3cj+1

1rj

∂r(rjBjcj)) ≈(rj+1/2 ·Bj+1/2 · cj+1/2

rj(∆r)− rj−1/2 ·Bj−1/2 · cj−1/2)

rj(∆r)≈

14rj(∆r)

· ((rj+1 + rj) ·Bj+1/2 · cj+1 − (rj−1 + rj) ·Bj−1/2 · cj−1 +

((rj+1 + rj) ·Bj+1/2 − (rj−1 + rj) · Bj−1/2) · cj) =Hb1cj−1 + Hb2ci + Hb3ci+1

If H1 = Ha1 + Hb1, H2 = Ha2 + Hb2 and H3 = Ha3 + Hb3, the five-point sten-cil of the 2D model can now be merged together. H1ci,j−1 + (H2 + K2)ci,j +H3ci,j+1 + K1ci−1,j + K3ci+1,j . The factors of the expression can be simpli-fied if the following notation is used: rj±1/2 = (rj + rj±1)/2 = r±, Aj±1/2 =(Aj + Aj±1)/2 = A±

j , Ai±1/2 = A±i .

The Bj-terms contains a factor 1/∆r and Bi a factor 1/∆z. With the definitionsB±

j = Bj±1/2 ·∆r/2 and B±i = Bi±1/2 ·∆z/2, the factors are moved outside the

terms. Finally the 2D and 1D stencils can be stated. The 1D stencil consist ofthe vertical components K1, K2 and K3.

K1 =(A−

i −B−i )

(∆z)2

K2 = K2+ + K2− =((−A+

i + B+i ) − (A−

i + B−i ))

(∆z)2

K3 =(A+

i + B+i )

(∆z)2

H1 =r−(A−

j −B−j )

rj(∆r)2

H2 = H2+ + H2− =(r+(−A+

j + B+j ) − (r−(A−

j + B−j )))

rj(∆r)2

H3 =r+(A+

j + B+j )

rj(∆r)2(2.8)

10

∂rc = 0 c = 0

c = 0

evaporation

inner regions

❄

✲r

z

Figure 2.2: Principle sketch of the boundary conditions of the two-dimensionalmodel. The upper edge corresponds to the surface, the left edge to the symmetryaxis and the right and lower edges to infinite radius and depth.

The spatial semidiscretization of the PDE into ODE form is almost completed,but the boundary conditions are still to be considered. That is a crucial step inthe model.

2.4 Boundary Conditions

2.4.1 Boundary Conditions in Two DimensionsIn the two-dimensional model boundary conditions occurs at the surface, alongthe symmetry axis, at z = ∞ and at r = ∞. Infinite depth and radius areapproximated with depth and radius which are large compared to the expectedmaximum values of the droplet. Thus, the boundaries are shaped like a hollowcylinder which due to symmetry along the z-axis can be translated into a rect-angle.

The right edge symbolizes the infinite radius, the left edge the center of thecylinder, the upper edge the surface and the lower edge the infinite depth. Theconcentration at infinite depth or radius is zero, ∂r(c) in the center is zero dueto symmetry and the flux on the surface must equal the flux of the evaporation.

Inner Regions

In the inner regions of the rectangle the whole five-point stencil from the previouschapter is used.

S1 : cni,j = (H2 + K2)cni,j + K1cni−1,j + K3c

ni+1,j + H1c

ni,j−1 + H3c

ni,j+1

Right Boundary Conditions

On the right edge (infinite radius) the concentration is always zero. This meansthat the time derivative on the edge is also always zero. This corresponds tothe identity K1 = K3 = H1 = H3 = 0. (K2 + H2) does not influence the result

11

and can for example be set to zero.

S2 : cn+1i,M = cni,M = 0 =⇒ cn+1

i,M − cni,M∆t

= cni,M = 0

Lower Boundary Conditions

Similary, the concentration at infinite depth is also always zero. Thus, the timederivative at the bottom is zero. K1 = K3 = H1 = H3 = 0. K2H2 = 0.

S3 : cn+1N,j = cnN,j = 0 =⇒ cn+1

N,j − cnN,j

∆t= cnN,j = 0

Left Boundary Conditions

The boundary conditions on left edge corresponds to the fact that the concen-tration in the cylindrical body is independent of the angle, which for instancemeans that: cθ=0 = cθ=π. This means that the concentration (and other radius-dependent factors) one (imaginary) discretization point step left of the centerr1 = 0 would equal those one step right of the center. As a result there willbe no net flow from one side of the cylinder to the other; the derivative of theradius is zero.

∂rcni,1 = 0 =⇒ (cni,2 − cni,0)

2∆r= 0 =⇒ (cni,2 = cni,0)

Both the radial diffusion term and the radial convection term has the radius rjin the denominator. The singularity along the symmetry axis can be avoidedby using the symmetry. At r1 = (r0 + r2)/2 = 0 one has A0 = A2, v0.5 = −v1.5,B0.5 = −B1.5 and r0 = −r2. Note that:

r0 + 2r1 + r2r1

=2(r0 + r2)(r0 + r2)/2

= 4

Thus, on the left edge the following holds.

1r∂r(rjAj∂r(cj))) ≈ 1

4r1(∆r)2· ((A2 + A1)(r2 + r1)c2 + (A0 + A1)(r0 + r1)c0

− ((A0 + A1)(r0 + r1) + (A2 + A1)(r2 + r1))c1) =1

2r1(∆r)2· ((A2 + A1)(r0 + r2)c2 − (A2 + A1)(r0 + r2)c1) =

1(∆r)2

· (−(A2 + A1)c1 + (A2 + A1)c2)

1rj

∂r(rjBjcj)) ≈ 14r1(∆r)

· ((r2 + r1)B3/2c2 − (r0 + r1)B1/2c0 +

((r2 + r1)B3/2 − (r0 + r1)B1/2)c1) =

12

12r1(∆r)

((r0 + r2)B3/2c2 + (r0 + r2)B3/2)c1) =

1(∆r)

· (B3/2c1 + B3/2c2)

This can be written:

H1 = 0

H2 =2(−A+

1 + B+1 )

(∆r)2

H3 =2(A+

1 + B+j1)

(∆r)2(2.9)

The result follows.

S4 : cni,1 = (H2 + K2)cni,j + K1cni−1,j + K3c

ni+1,j + H3c

ni,j+1

Upper Boundary Conditions

On the upper edge S1 the evaporation on the surface must be considered. Ac-cording to Baines and James [1], the vertical flux (kg m−2s−1) from an in-dividual small wetted spot with radius rend is described by the the followingequation, in which bar indicates the averaged value:

F1 = 0.66716(R−1g,1c

n1 − cb) · (

(Dairg,1 · v∗)2rendν

)1/3 = k1cn1 − k2cb (2.10)

In the two dimensional extension all factors with bars are replaced with theactual values at the points.

F1,j = 0.66716(R−1g,1,jc

n1,j − cb) · (

(Dairg,1,j · v∗)2rendν

)1/3 = k1j cn1,j − k2jcb

The factor rend is the approximation of the radius of the wetted circle at thesurface. In this work it is defined as the radius where the concentration at thesurface is 100·rtol per cent of the concentration in the middle of the surface. Thefactor influences the rate of evaporation and consequently the fraction remainingof the substance. A standard value of rtol is 0.01.

rend = r : (cn1,rcn1,1

= rtol)

An approximation is made by stepping outwards from the center comparingconsecutive quotes of and picking the last one which is larger than rtol. Theindex of that object is called rjstop, which initially is equal to the number rjend.

The fact that the flux is known at the surface and immediately below the sur-face can be used for deriving an expression for c at the surface. The equation

13

below is obtained for the vertical components with index j. The horizontal com-ponents are calculated as usual. The background concentration directly abovethe ground can be considered to be proportional to the concentration on theground. Therefore the approximation c0/R

−1g ≈ αc1 is made. The background

concentration is usually very small and can for most applications be set to zero.

c1,j + ∇r(F1) = −∇z(F1) − εc1 = −F3/2 − F1

∆z/2− εc1 =

(A3/2∂z(c3/2) + c3/2B3/2 + k1jc1 − k2c0)∆z/2

− εc1 =

2A3/2(c2 − c1)(∆z)2

+2B3/2(c2 + c1)

2∆z+

2k1j(1 − α)c1∆z

− εc1 =

2A+1 (c2 − c1)(∆z)2

+2B+

1 (c2 + c1)(∆z)2

− 2k1j(1 − α)c1|∆z| − εc1

This can be stated in the following stencil.

K1 = 0

K2 =2(B+

1 −A+1 )

(∆z)2− 2k1j(1 − α)

|∆z| − ε

K3 =2(B+

1 + A+1 )

(∆z)2(2.11)

Remember that:

k1j = 0.66716R−1g,1j · (

(Dairg,1j · v∗)2rendν

)1/3

Two facts that must be considered are that an decrease in concentration res-ults in an increase in evaporation divided by concentration and that the fluxfrom Baines [1] is based on experiments with static radius. When the liquidspreads out horizontally, the flux expression is no longer completely valid. Thespreadout leads to lower concentrations on a larger area which implies too highevaporation rates. Deriving an precise formula for the damping factor is not aneasy task. A simple but useful approximation is setting the damping functionproportional to the surface concentration. It can then simply be merged intothe α-term.

The approximation does not lack drawbacks. The contribution to the flux fromthe edge of the droplet becomes too large compared to the local flux rate in thecenter of the surface. For certain problems this can result in an initial unphys-ical increase of evaporation up to a proper rate. Anyhow, despite its drawbacksthe approximation gives for example surprisingly nice evaporation curves.

At last, the upper boundary condition can be written with the following fourpoint stencil:

14

S5 : cn1,j = (H2 + K2)cn1,j + K3cn2,j + H1c

n1,j−1 + H3c

n1,j+1

The whole three-dimensional system including the boundaries is now trans-formed into an ordinary differential equation (ODE) on the form c = f(c).

2.4.2 Boundary Conditions in One dimensionIn the one-dimensional model there will be just one upper and one lower bound-ary condition. The standard stencil has three points. In practice, the stencilsequals the horizontal parts of the 2D model. So in the central parts there is anrelationship as below with the same K values as in the 2D part.

S1 : cni = K1cni−1 + K2c

ni + K3c

ni+1

The concentration at the lower boundary (infinite depth) is set to zero. In 1Dthis equals that K1 = K2 = K3 = 0 or:

S2 : cn+1N = cnN = 0 ⇒ cnN = 0

In the one-dimensional model the radius Rend is considered to be constant witha value equal to the final radius re from the spreading-sorption submodel. Useof the same idea as in the three-dimensional model gives the upper boundarycondition:

S3 : cn1 = K2cn1 + K3c

n2

The K value equals those from the 2D model except from the facts that theyare independent of the radius and that no damping factor is included in theα-term.

K1 = 0

K2 =2(B+

1 −A+1 )

(∆z)2− 2k1(1 − α)

|∆z| − ε

K3 =2(B+

1 + A+1 )

(∆z)2

k1 = 0.66716R−1g,1 · (

(Dairg,1 · v∗)2rendν

)1/3

At this stage when an equation ˙ci,j = f(ci,j) is stated for each point on thegrid, the obvious next step is gathering all equations into a system of equations.This must be done if the equations are to be solved.

2.5 Matrix GenerationA necessity before gathering all equations into a system of equations is that theconcentrations can be written as a vector. If the dimension of the grid is one,this is already done. Otherwise lexicographical ordering may be used.

15

2.5.1 Lexicographical OrderingThis kind of ordering is the one that occurres in dictionarys. The system is: Firstsort due to the first position, then sort due to the second position inside eachgroup that has the first value in common and so on. The following expressionis an example of lexicographical order:

sort(xi,j : i = 1..2, j = 1..2) ⇒ {x1,1, x1,2, x2,1, x2,2}

2.5.2 Matrix Generation in Two DimensionsIn the two-dimensional case with a grid consisting of N rows and M columns wehave N times M equations. With the lexicografical ordering of c, the equationcan be written cn = f(cn) = W (cn) · cn. The matrix W is a NM*NM pentadi-agonal matrix. The value on the main diagonal on row (i-1)M+j correspondsto the H2 +K2 value of the gridpoint (i,j). The subdiagonal corresponds to H1

and the super diagonal to H3. The subdiagonal number M corresponds to K1

and the superdiagonal M to K3. A pedagogical sketch of a system with N=3and M=4 is shown below.

c1,1 c1,2 c1,3 c1,4c2,1 c2,2 c2,3 c2,4c3,1 c3,2 c3,3 c3,4

=⇒

∗ ∗ 0 0 ∗ 0 0 0 0 0 0 0∗ ∗ ∗ 0 0 ∗ 0 0 0 0 0 00 ∗ ∗ ∗ 0 0 ∗ 0 0 0 0 00 0 ∗ ∗ ∗ 0 0 ∗ 0 0 0 0∗ 0 0 ∗ ∗ ∗ 0 0 ∗ 0 0 00 ∗ 0 0 ∗ ∗ ∗ 0 0 ∗ 0 00 0 K1 0 0 H1 H2 + K2 H3 0 0 K3 00 0 0 ∗ 0 0 ∗ ∗ ∗ 0 0 ∗0 0 0 0 ∗ 0 0 ∗ ∗ ∗ 0 00 0 0 0 0 ∗ 0 0 ∗ ∗ ∗ 00 0 0 0 0 0 ∗ 0 0 ∗ ∗ ∗0 0 0 0 0 0 0 ∗ 0 0 ∗ ∗

·

c1,1c1,2c1,3c1,4c2,1c2,2c2,3c2,4c3,1c3,2c3,3c3,4

=

f1,1

f1,2

f1,3

f1,4

f2,1

f2,2

f2,3

f2,4

f3,1

f3,2

f3,3

f3,4

2.5.3 Matrix Generation in One DimensionIn one dimension the gathering of the equations is very simple. Just put oneequation on each row to obtain a three-diagonal N*N system of equations onthe form W · c = f with W as the N*N matrix and c and b as N*1 vectors. Forexample, if N=6:

16

W =

∗ ∗ 0 0 0 0∗ ∗ ∗ 0 0 00 H1 H2 H3 0 00 0 ∗ ∗ ∗ 00 0 0 ∗ ∗ ∗0 0 0 0 ∗ ∗

Note that the nonzero elements in both cases are gathered into diagonal stripes.This is not a random coincident. The finite differential method always resultsin matrices with diagonal bands and is therefore also called the band matrixmethod.

2.5.4 Compressed Column Storage (CCS) FormatFor large matrices with relatively few non-zero elements it is extremely ineffi-cient to store all elements as in a usual full matrix. Fortunately a large amountof sparse matrix formats are at hand, see for instance Robert [21] and Don-garra [10]. Matlab uses Compressed Column Storage (CCS) format for sparsematrices. This format stores one vector V with the values of the matrix takencolumnwise from the upper left corner, one with the column indices of the ele-ments and one with the element numbers of V for the first element in eachcolumn.

With the compressed column (or row) storage format no zeros at all are storedand the matrix does not need a specific shape. The drawback is that the ad-dressing of the elements are not straightforward and extra calculations must bemade when the matrix is used.

2.6 Initial ConditionsMany scientists, like Youngs [28], Reichardt [20] and Gardner [13], have ex-amined the shape of the profile of concentrations which occurs when a dropletis put on the top of a vertical column and is absorbed into the column. Their res-ults coincides. The FOI phase 3 submodel starts to calculate when the heigthof the spherical cap of the droplet has become infinitely small, which meansthat all of the liquid has penetrated the surface. As a result the shape of theexperimental curves at t = 0 can be used as initial data in the model.

2.6.1 Initial Conditions Two DimensionsFrom the experiments of Gardner [13] a graph of the curve c(z), c = 0..1, z = 0..1is given for t = 0. In order to bring the curve into a function on closed forma polynomial was interpolated from a couple of points from the experimentalcurve translated from z = −1..0. to z = 0..1.

c(z) = −1.515 · z1/2 + 3.4514 · z1/4 − 0.9005 · z1/40, z = 0..1 → c = 0..1

17

0 0.2 0.4 0.6 0.8 1 1.2 1.4−1

−0.9

−0.8

−0.7

−0.6

−0.5

−0.4

−0.3

−0.2

−0.1

0Interpolated Curve

Concentration

Depth

interpolatedexperimental

Figure 2.3: Interpolation of a Inital Profile from Experimental Data from Gard-ner

In the two-dimensional case the initial values are set to a constant value insidea rotational body symmetric around the z-axis. In the first test the body wassliced with horizontal cuts and the normed area was set equal to the normedconcentration at the same height. Noting that c ∝ r2, use of a three term leastmean square interpolation could be used do deduce a pleasant formula. Thez-values are still translated one step.

r(z) = −1.4561 · z1/2 + 2.8716 · z1/4 − 0.4094 · z1/40, z = 0..1 → c = 0..1

Surprisingly the speed and accuracy of the program was improved when theinitial curve was shaped like a plug instead of a relatively smooth function.Because the plug works better, it was chosen. This means that all gridpointsinside the cylinder πr2

end · zend have the same nonzero concentration, while allother gridpoints have ct=0 = 0. When the curve z is integrated from zero to onea factor occurs. 1 − ∫ 1

0c(z)dz ≈ 0.84. The maximum initial concentration is

given by (φ− θ)ρc Winter,[27]. Inside the plug the average initial concentrationought to be 0.84(φ− θ)ρc, but while the factor 0.86 is previously used by FOI,the final expression for average initial concentration reads:

c1i,j ={

0.86(φ− θ)ρckcomp 1 ≤ i ≤ ziend, 1 ≤ j ≤ rjend0 (all other i and j)

A compensation factor kcomp is introduced in order to simplify comparisonswith the 1D model. In 2D, if c(rend) = c and c(rend + ∆r) = 0, the averageconcentration between the points are c/2. In 1D it is zero. The difference comesthus from the dicretization and not the physics.

kcomp =r2jend

r2jend + (r2

jend+1 − r2jend)/2 + 1.5

18

2.6.2 Initial Conditions One DimensionThe intuitive choice for the 1D initial conditions is taking the previously gener-ated function c(z), evaluating it at a properly choosen amount of equally spacedgridpoints from 0 to 1 and simply multiplying with Zend and (φ− θ)ρc. Thus,practical tests show that a plug works better in 1D too. The result is as follows.

c1i ={

0.86(φ− θ)ρc 1 ≤ i ≤ ziend0 i > ziend

If smoother initial values are needed anyway, one could run the program a fewseconds without diffusion to spread out the droplet a bit.

2.7 Remaining FractionThe remaining fraction of substance in the ground at time tn, where n is thecurrent timestep number, is determined by the relationship:

Fr(tn) = Fnr =

∫VcndV∫

Vc1dV

=mn

m1

In three dimensions with cylindrical coordinates it can be approximated by:

Fnr =

π∑N−1

i=1

∑M−1j=1 (r2

j+1 − r2j )(c

ni+1,j + cni,j + cni,j+1 + cni+1,j+1)∆zi/4

π∑N−1

i=1

∑M−1j=1 (r2

j+1 − r2j )(c

1i+1,j + c1i,j + c1i,j+1 + c1i+1,j+1)∆zi/4

=

∑N−1i=1

∑M−1j=1 (r2

j+1 − r2j )(c

ni+1,j + cni,j + cni,j+1 + cni+1,j+1)∑N−1

i=1

∑M−1j=1 (r2

j+1 − r2j )(c

1i+1,j + c1i,j + c1i,j+1 + c1i+1,j+1)

In the one dimensional case it can be approximated by:

Fnr =

πR2end

∑N−1i=1 (cni+1 + cni )∆zi/2

πR2end

∑N−1i=1 (c1i+1 + c1i )∆zi/2

=∑N−1

i=1 (cni+1 + cni )∑N−1i=1 (c1i+1 + c1i )

Note that the remaining fraction is not explicitely dependent on the time andcan thus be calculated outside the timeloops in order to save calculation time.

If the times are stored in the vector T, n = 1..nmax, changes in the frac-tion remaining per timestep can be written like:

∂tFnr =

(Fn+1r − Fn

r )(Tn+1 − Tn)

2.8 Evaporated FractionIf the degeneration term of the equation is set to zero, the remaining part of thesubstance is only dependent of the accumulated evaporation. A way to keep an

19

eye on errors of integration and discretization is to use a conservation law. Theevaporated fraction plus the remaining fraction must equal one. If the fractionevaporated is called Fev then 1−Fev ≈ Fr. The index rjstop is abbreviated withjs and in 3D the following approximation can be used.

∂tFnev = −π(1 − α)

M1(

js−1∑j=1

(r2j+1 − r2

j )(kn1,j+1 + kn

1,j)(cnj+1 + cnj )/4) +

(r2js+1 − r2

js)(kn1,js)(c

njs+1 + cnjs)/2) )

F 1ev = 1

F 2ev = 1 + ∂tF

1ev · T (2)

Fn+1ev = 1 + ∂tF

nev ·

(Tn+1 − Tn−1)2

The corresponding 1D relationships are even simpler. The α1 term in 1D isusually zero while it does not contain any damping.

∂tFnev = −πr2

end(1 − α1)kn1 c

n1

M1

F 1ev = 1

F 2ev = 1 + ∂tF

1ev · T (2)

Fn+1ev = 1 + ∂tF

nev ·

(Tn+1 − Tn−1)2

2.9 Jacobian of fThe Jacobian of f, J = ∇c(f(c), is of major importance for many ODE solvers.While the Jacobian in this case not can be derived on closed form a numericapproach will be used.

The current ODE can be written as: cn = f(cn) = W (cn)cn. Note that fdoes not depend explicit on the time; just implicit through c. This kind ofequation is called an autonomus ordinary differential equation. If the numberof rows in the small original matrix is called N = Zimax and the number ofcolumns M = Rjmax we have K = M ∗N . The lengths of f and c are both Kand the size of J , V and W are K ∗K. The vertical index of W is called i andthe horizontal j. The Jacobian can be stated like:

Jn = ∇c(f(cn)) =

∂c1(f1) . . . ∂cK(f1)

.... . .

...∂c1(fK) . . . ∂cK

(fK)

20

Jik = ∂ck(fi) = ∂ck

N∑

j=1

wij(c)cj

=

N∑j=1

∂ck(wij(c)cj) =

wik(c) +N∑

j=1

∂ck(wij(c))cj = wik(c) + vik(c)

The Jacobian can thus be written as: Jn = Wn + V n. From the precedingderivation of the stencils it can be noted that the dependence of c for wij in 3Dis as follows:

wi,i−M (ci, ci−M )wi,i−1(ci, ci−1)wi,i(ci, ci−1, ci+1, ci−M , ci+M )wi,i+1(ci, ci+1)wi,i+M (ci, ci+M )wi,j = 0 (j �= i±M, j �= i± 1, j �= i)

It is possible to state the formulas for vik more explicit. The following holds forrow i of the matrix V.

vi,i−M = (∂ci−M(wi,i−M ))ci−M + (∂ci−M

(wi,i))civi,i−1 = (∂ci−1(wi,i−1))ci−1 + (∂ci−1(wi,i))civi,i = (∂ci

(wi,i−M ))ci−M + (∂ci(wi,i+M ))ci+M +

(∂ci(wi,i−1))ci−1 + (∂ci

(wi,i+1))ci+1 + (∂ci(wi,i))ci

vi,i+1 = (∂ci+1(wi,i+1))ci+1 + (∂ci+1(wi,i))civi,i+M = (∂ci+M

(wi,i+M ))ci+M + (∂ci+M(wi,i))ci

vi,k = 0 (All other k)

Note that the pentadiagonal structure is preserved and that the number of non-zero elements does not increase.

Derivatives of v are approximated with forward differences in order to reducethe amount of calculations needed; the second of the two terms in the expressionis already evaluated.

∂ck(f(c)) =

f(ck + h, cj �=k) − f(c)h

An example of the use of the formula is given.

∂ci−1(wi,i−1(ci, ci−1)) =(wi,i−1(ci, ci−1 + h) − wi,i−1(ci, ci−1))

h

The factor hmust be chosen carefully to minimize numerical errors. The approx-imation of the derivative becomes better with decreasing h until the cancellation

21

error due to the limited computer precision starts to dominate.

The concentration c has a smooth behaviour, which implies that the magnitudeof c will be approximately the same in the whole region. If it is about one,the smallest possible h before the cancellation errors starts to dominate is thesquare-root of the machine precision. If the magnitude of c differs from one, theexpression should be multiplied with the norm of c. The maximum norm willbe used because it is more suited for computers than the 2-norm. Eventually asmall extra term could be added to eliminate the risk of dividing by zero if theconcentration is zero everywhere. The condition

√εm · εk > εm ⇒ εk >

√εm is

used. Finally:

h =√εm · (||c||max + εk)

Remember that for example wi,j−1 = H1. The terms K1, K3, H1, H3, K2+,K2−, H2+ andK2− are previously evaluated. If the concentration ci is disturbedby h the corresponding factors are denoted Kh

1 , Hh3 ,etc. On the contrary, if ci

remains constant but ci+1, ci−1, ci+M or ci−M are perturbed by h, the factorsare called Hh1

3 etc. The elements on row i in the V matrix can therefore bestated as in the system below.

vi,i−M =(Kh1

1 −K1

h

)ci−M +

(Kh1

2− −K2−h

)ci

vi,i−1 =(Hh1

1 −H1

h

)ci−1 +

(Hh1

2− −H2−h

)ci

vi,i =(Kh

1 −K1

h

)ci−M +

(Hh

1 −H1

h

)ci−1 +

(Hh

3 −H3

h

)ci+1 +(

Kh3 −K3

h

)ci+M +

(Hh

2 −H2

h

)ci +

(Kh

2 −K2

h

)ci

vi,i+1 =(Hh1

3 −H3

h

)ci−1 +

(Hh1

2+ −H2+

h

)ci

vi,i+M =(Kh1

3 −K3

h

)ci−M +

(Kh1

2+ −K2+

h

)ci

vi,k = 0 (All other k)

A couple of additional calculations are needed to obtain the values of Hh11 , Hh

3

and the other unknown factors. First the matricesDhsg,Dhs

lc,Dhsg, Rhg, Rhη, Rhl

are constructed by recalculating Dsg etc with cnew = c + h . Thereafter values

for among others Ah±r , Ah±

z and Bh1±r , of which Hh1

1 and the other unknownsconsists, are searched for on the original grid 1 ≤ j ≤ M , 1 ≤ i ≤ N . As usualthe radial components can be extracted by simply replacing all i with j for Aand B.

Ahi = Dhs

g(i)Rhg(i) + ρcDhslc(i)Rhη(i) + Dhs

l (i)Rhl(i) =⇒

Ah±i =

Ahi + Ai±1

2

22

Ah1±i =

Ai + Ahi±1

2

Bh±i = ±1

4· ((Ds

g(i± 1) + Dhsg(i)) · (Rg(i± 1) − Rhg(i)) +

ρc · (Dslc(i± 1) + Dhs

lc(i)) · (Rη(i± 1) − Rhη(i)) +

(Dsl (i± 1) + Dhs

l (i)) · (Rl(i± 1) − Rhl(i))) − ∆z

2VzR

sl (i± 1)

Bh1±i = ±1

4· ((Dhs

g(i± 1) + Dsg(i)) · (Rg(i± 1) − Rg(i)) +

ρc · (Dhslc(i± 1) + Ds

lc(i)) · (Rhη(i± 1) − Rη(i)) +

(Dhsl (i± 1) + Ds

l (i)) · (Rhl(i± 1) − Rl(i))) − ∆z

2VzRhs

l (i± 1)

The 3D stencils can then be written:

Kh1 =

(Ah−i −Bh−

i )(∆z)2

Kh2 = Kh

2− + Kh2+ =

((−Ah+i + Bh+

i ) − (Ah−i + Bh−

i ))(∆z)2

Kh3 =

(Ah+i + Bh+

i )(∆z)2

Hh1 =

r−(Ah−j −Bh−

j )rj(∆r)2

Hh2 = Hh

2− + Hh2+ =

(r+(−Ah+j + Bh+

j ) − (r−(Ah−j + Bh−

j )))rj(∆r)2

Hh3 =

r+(Ah+j + Bh+

j )rj(∆r)2

At the upper boundary we have:

Kh1 = 0

Kh2 = 0 + Kh

2+ =2(Bh+

1 −Ah+1 )

(∆z)2− 2k1(1 − α)

|∆z| − ε

Kh3 =

2(Bh+1 + Ah+

1 )(∆z)2

And at the left boundary:

Hh1 = 0

Hh2 = 0 + Hh

2+ =2(−Ah+

1 + Bh+1 )

(∆r)2

Hh3 =

2(Ah+1 + Bh+

1 )(∆r)2

23

At the right boundary:

Hh1 = Hh

2− = Hh2+ = Hh

3 = 0

And at the lower boundary:

Kh1 = Kh

2− = Kh2+ = Kh

3 = 0

To obtain the remaining stencils, just replace ”h” with ”h1” in the formulasabove. The whole 3D Jacobian Jn = Wn + V n can then be gathered. The 1DJacobian is obtained by setting M = 1 and simply ignoring every factor H andthe terms which exclusively depend on H.

The amount of work can be reduced by storing values of w(c = 0) and v(c = 0)and using the precalculated values instead of recalculating them. Rememberthat changes in h still gives new values of the v-terms.

Another way of reducing the calculations needed is done by noting that insidethe grid, where no boundary conditions are present, the following holds:

Ah1+i = Ah−

i+1

Ah1−i = Ah+

i−1

Bh1+i = Bh−

i+1

Bh1−i = Bh+

i−1

=⇒

Kh11 (i) = −Kh

2+(i− 1)Kh1

2−(i) = −Kh3 (i− 1)

Kh12+(i) = −Kh

1 (i + 1))Kh1

2+(i) = −Kh2−(i + 1)

The principle is the same for the radial terms, but the index of the explicitoccurances of r does not change even though the rest does.

2.10 Time-DiscretizationThe choice of method for the discretization of the time derivatives an OrdinaryDifferential Equation (ODE) is heavily dependent of the properties of the Jac-obian. If the largest negative eigenvalues of the Jacobian is large in comparisonwith the reciprocal of the time span, the system is stiff. [3]. The present Jac-obians have typically very large γ and very small 1/T , which confirms stiffness.

Explicite methods, which only uses previously known information in each timestep,are not efficient on stiff problems while the time steps must be very small toassure stability. Instead implicit methods, in which a non-linear system of equa-tions must be solved in each time step, ought to be used. Each implicit methodcan be written on the following form with the values of k1, k2 and k3 known.

Φ(cn+1) = (I + k1)cn+1 − k2fn+1 − k3 = 0 (2.12)

Non-linear systems of equations can preferably be solved with a Newton iter-ation, in which linear systems are solved in each iteration step. The searchdirection in Newton iteration number k is denoted pk.

24

∇cn+1Φ(cn+1(k) )pk = −Φ(cn+1

(k) )

cn+1(k+1) = cn+1

(k) + pk (2.13)

The following time-discretization schemes are described: The Trapetzoid method,Backward Differential Formulas (BDF) and Numerical Differential Formulas(NDF).

2.10.1 The trapetzoid formula (Cranck-Nicholson)The trapezoid formula is a semi-implicit one-step method for ODE solving. Eventhough it might not be the best choice for stiff problem, it will still be describedwhile it can be used to initialize multistep methods. It has a truncation errorproportional to the square of the timestep. Applied on c = f the trapezoidformula can be written:

∂tcn =

(cn+1 − cn)∆t

=fn + fn+1

2(2.14)

On the specific diffusion-convection problem with f(cn) = W (cn)cn this rela-tionship yields that:

I · (cn+1 − cn)∆t

=(W (cn+1)cn+1 + W (cn)cn)

2⇐⇒

Φ(cn+1) = (I − ∆t

2W (cn+1)) · cn+1 − (I +

∆t ·W (cn)2

)cn = 0 =⇒

∇cn+1Φ(cn+1) = (I − ∆t

2Jn+1)

The original FOI model used a explicit-implicit variant of Crank-Nicholson. Itused the trapezoid rule with central differences in time (cn+1 − cn−1)/(2∆t) =(fn+fn+1)/2 on the diffusion terms and explicit Forward Euler (cn+1−cn)/∆t =fn on the convection and reaction terms. The approximationW (cn+1) ≈ W (cn)was made. The model could be described as a quasi two-step method becausetwo previous timesteps must be known to run the method. Observe that thismethod is not used in the present models.

2.10.2 Backward Differentiation Formulas (BDFs)The backward differential formulas (BDFs) are linear multistep method espe-cially useful for stiff problems. BDF formulas can be of any order, but only thefirst six are zero stable, which means that the methods have small sensitivityto small changes in indata even when ∆t → 0. For higher orders of BDFs, thestability region decreases.

All BDF methods used with backward differences applied on cn+1 = f(cn+1)can be written on the following form: (See Quarteroni [19], Sjö [25] or Deuflard[9])

25

p∑m=1

(∆t)m−1

m∇m

t cn+1 =cn+1 +

∑p−1j=0 ajc

n−j

β∆t= f(cn+1) ⇐⇒

Φ = cn+1 − β∆t · f(cn+1) +p−1∑j=0

ajcn−j = 0 ⇒ {f = Wc} ⇒

Φ = (I − β∆tW (cn+1))cn+1 +p−1∑j=0

ajcn−j = 0

∇cnΦ = (I − β∆tJn+1) (2.15)

The coefficients for all zero-stable fixed step BDFs are stated below.p a0 a1 a2 a3 a4 a5 β1 -1 0 0 0 0 0 12 - 43

13 0 0 0 0 2

33 - 1811

911 - 2

11 0 0 0 611

4 - 48253625 - 1625

325 0 0 12

255 - 300137

300137 - 200137

75137 - 12

137 0 60137

6 - 360147450147 - 400147

225147 - 72

14710147

60147

A good choice for starting a Newton iteration is setting:

cn+1(0) =

p∑m=0

(∆t)m∇mt cn

An approximation of the truncation error of a BDF of order p is given by:

∆tp+1

p + 1cn+1(p+1) ≈

1p + 1

(∆t)2p+1∇p+1t cn+1

As an example, the second order BDF2 is deduced. The truncation error ofBDF2 is proportional to the square of the timestep. Two preceeding steps mustbe evaluated to initialize it.

2∑m=1

(∆t)m−1

m∇m

t cn+1 =((cn+1 − cn) + (cn+1 − 2cn + cn−1)/2)

∆t=

(cn+1 − 43c

n + 13c

n−1)23∆t

= f(cn+1) =⇒

Φ(cn+1) = (I − 23∆tWn+1)cn+1 − 4

3cn +

13cn−1 = 0

∇cn+1Φ(cn+1) = (I − 2∆t

3Jn+1)

cn+1(0) = cn + (cn − cn−1) +

(cn − 2cn−1 + cn−2)2

=5cn − 4cn−1 + cn−2

2

26

For nonstiff problems, BDFs are relatively inefficient compared to the best ex-plicit methods. Anyhow, the present problem is very stiff. Another problemfor all multistep algorithms is that it is rather complicated to control the time-discretization errors if the stepsize is not constant. If constant stepsize is used,a lot of unnecessary operations will have to be done. For this reason, Matlabsbuild-in adaptive ODE solvers are used in the main program. The change ofstep size in these are based on the fact that a certain interpolation polynomialought to give the same value for different ∆t. In the alternative home-madeversion the step size can only be multiplied or divided by 2.

For convection domianated diffusion convection equations, a semiimplicite BDF2is suggested by Hundsdorfer [16], but the diffusion term dominate in the presentequation, and the method is thus not chosen.

There also exist a quasi first order BDF called TR-BDF2 described in Hosea[15], which combines the use of BDF2 and Cranck-Nicholson in a cleaver way.This could be a good alternative, but practical tests on the matrices showsthat the best was using a kind of fine tuned BDFs called Numerical DifferentialFormulas, NDFs.

2.10.3 Numerical Differentiation Formulas (NDFs)The NDFs presented by Shampine et al [24] equals the BDFs with a slightdifference of an additional term which reduces the truncation error and allowsfor larger timesteps. The factor κ is a scalar.

Φ(cn+1) = cn+1 − β∆t · f(cn+1) +p∑

j=0

ajcn−j −

p∑j=1

κ

j(cn+1 − cn+1

(0) ) = 0 =⇒

Φ(cn+1) = (I − β∆tW (cn+1) − I

p∑j=1

κ

j)cn+1 +

p∑j=0

ajcn−j +

p∑j=1

κ

j(cn+1

(0) ) = 0

∇cn+1Φ = (I − β∆tJn+1 − I

p∑j=1

κ

j)

cn+1(0) =

p∑m=0

(∆t)m∇mt cn (2.16)

The truncation error is approximated with:

εT ≈ (p∑

j=1

κ

j+

1p + 1

)∆tp · cn+1(p+1) (2.17)

A table of κ coefficients and stability angles for the first 5 orders of fixed stepBDFs and NDFs are attached.

27

order(p) NDF koeff κ BDF (degrees) NDF (degrees)1 -1.850 90 902 -1/9 90 903 -0.0823 86 804 -0.0415 73 665 0 51 51

When the time discretization of the ODE is completed, the next step is tohandle the solving of linear and non-linear systems of equations.

2.11 Solving Systems of EquationsThere are two principle ways of solving linear systems of equations: Direct meth-ods and iterative methods. Among the direct ones classical Gauss eliminationand factorization methods as LU, CHOLESKY, and QR can be mentioned. Theiterative ones can be divided into stationary iterative methods like Gauss-Seideland non-stationary as the conjugated gradient method (CG), GMRES (gener-alized minimum residual) and other Krylow type methods. Nonlinear systemscan be solved by Newtons method in an iterative procedure, in which linearsystems must be solved each iteration step.

For each method a problem can be constructed for which another algorithmworks better. Every variant has its strengths and weaknesses. Direct methodsare usually to prefer for dense problems and iterative for sparse, but for matriceswith large condition numbers, especially if they are neither symmetric, positivedefinit or easy to find good preconditioners for, sparse direct solvers can becomparatively. The efficiency of iterative methods is strongly dependent of thechoice of preconditioner, the eigenvalues and the condition number.

2.11.1 Direct Solvers (GAUSS, LU, CHOL, QR )The group of direct solvers compound among others Gauss-elimination and fac-torization methods. Lower and upper triangular matrices are easy to solve withgauss elimination. The basic idea for factorization is to split one difficult equa-tion into two simple ones which can be easily solved with Gauss elimination.QR factorization is used for rectangular (nonquadratic) matrices. For symmet-ric matrices, Cholesky factorisation is most efficient. Anyhow, since the presentmatrices are quadratic and unsymmetric, Lower Upper (LU)-factorisation isused. If L is a lower triangular matris, U a upper triangular and the systemAx = b is to be solved, the basic strategy can be displayed with:

Ax = LUx = Ly = b =⇒{Ux = yLy = b

Most of the work is thus done in the factorization step. This is an advantage ifthe equation shall be reevaluated with constant system matrix but with variableright hand sides.

28

0 200 400 600 800 1000 1200 1400 1600

0

200

400

600

800

1000

1200

1400

1600

nz = 7992

Sparsitypattern J

Figure 2.4: Sparcity pattern of J without reordering. NNZ(J)=7992,NNZ(L)=4824, NNZ(U)=11965.

There are many ways to improve the performance of the basic solvers. Oneof the most commonly used is pivoting, especially partial pivoting. The back-ground is that if the pivot, the frontal term according to which the other termsare normed, in a step of gauss elimination is small, the roundoff errors will growunacceptably much. If the rows (or/and columns) are interchanged in such away that the row with the largest scalar in the current pivot column is inter-changed with the current pivot raw (or vice versa), this will not happend. Ifboth rows and colums are exchanged, the method is called complete pivoting.Else if only one of them is swiched it is called partial pivoting. The reorderingsare stored in a row interchange matrix or/and a column interchange matrix. Incontrast to the partial pivoting, the number of operations for complete pivotingis usually too high for practical applications.

The rows or columns can also be reordered to preserve sparsity of the LUfactors. Some reordering strategies like Reversed Cutchill McKee (RCM), whichworks best for “long and skinny” problems, attempt to gather the non-zero ele-ments of the original matrix around the main diagonal and thereby reducingfill in. Other types are (symmetric or column) minimum degree permutations(SYMMMD/COLMMD) and (symmetric or column) approximate minimum de-gree permutation (SYMAMD/COLAMD). While our matrices are unsymmet-rical only COLRCM, COLMMD and COLAMD are of interest. COLAMD tendsto generate lower amount of fill in than the others. A few examples are given inthe graphs below. Note that the fill in of a non-reordered bandmatrix is stronglydependent of the bandwidth of the band. As a consequence, problems of higherdimension than 2 are rarely well suited for direct methods.

For unsymmetrical matrices, the natural strategy is to use column permuta-tions to preserve sparsity, and row permutations to reduce roundoff errors. Thefactorisation would read PAQ = LU ⇔ A = P−1LUQ−1.

29

0 200 400 600 800 1000 1200 1400 1600

0

200

400

600

800

1000

1200

1400

1600

nz = 7992

J med colamd

Figure 2.5: Sparcity pattern of J with COLAMD reordering. NNZ(J)=7992,NNZ(L)=1730, NNZ(U)=7992.

Efficient algorithms can be applied for the actual numerical factorization. Two ofthe best LU solvers of today, SuperLU [7], [8] and UMFPACK [6], which bothincorporate row and column reordering, use different strategies. UMFPACKuses a multifrontal strategy with small, dense and easily solved frontal matricesused together with a column elimination tree. On the other hand, SuperLUuses a supernodal elimination tree in order to perform calculations on smalldense submatrices in each super column. Both programs are implemented inone Matlab version, one Fortran version and one C version. On the presentmatrices the performance of SuperLU and UMFPACK are approximately equi-valent. UMFPACK is chosen.

2.11.2 Iterative SolversIterative methods for solving linear or non-linear systems of equations Ax = bare all based on the fundamental idea that if x is the correct values, a convergentsequence gives that x = limk→∞ xk. Generally every iterative scheme can bedescribed on the form:

x0 = f0(A, b)xk+1 = fk+1(xk, xk−1, ..., xk−l, A, b), (k ≥ l)

The method is stationary if every fk is independent of k, otherwise it is non-stationary. If fk depends lineary on the x values, the method is lineary, otherwiseit is nonlineary. The order of the method equals the number of previous stepsused in each iteration.

2.11.3 Stationary Iterative MethodsStationary variants are obtained by adding zero in a clever way to the originalsystem. (Carlsund [4]) Schematically the system is obtained by:

30

Ax = b + (Px− Px) ⇐⇒Px = (P −A)x + b ⇐⇒x = (I − P−1A)x + P−1b = Bx + c =⇒x(k+1) = (I − P−1A)x(k) + P−1b

The error is given by:

e(k+1) = Be(k) ⇒ e(k+1) = Bke(0)

Obviously only iff ||B|| < 1 the error will grow smaller. The smaller ||B||, thefaster convergence up to a rate given by the spectral radius, which is denoted:χ(B) = maxi|λi| < 1.

The splitting of a matrix A into a lower triangular, an upper triangular anda diagonal part, gives rise to submatrices A = (L+D +U). These submatricesare used as the matrix P in different stationary methods.

Jacobi : P = D

Gauss− Seidel : P = L + D

Successive−Over −Relaxation : P = L +D

ω, (0 ≤ ω ≤ 2)

The performance can be improved by use of preconditioners which are men-tioned further on. Anyway, in a competition with the nonstationary solvers, thestationary ones will lose severly.

2.11.4 Nonstationary Iterative MethodsIterative methods of nonstationary kind can basically be described with a iterat-ive scheme which reads that the next iteration point is a sum of the value of theinitial point plus the correction terms of all iterations. The correctors are vec-tors, which means they have both magnitude and direction. Even optimizationprocedures like Gauss-Newtons method can be described in this way.

xk+1 = xk + sk · pk = x0 +k∑

i=1

sipi

Krylow methods are based on the use of the Krylow space Uk, see Hanke [14].For a linear system Ax = b with residual r = Ax − b, the Krylow space isdenoted:

U0 = {0}Uk(A, r0) = lin{r0, Ar0, ..., A

k−1r0}The energy inner product and the energy norm are defined by:

< x, y >A:=< x,Ay >=< Ax, y >= yT (Ax)

||x||A =< x, x >1/2A

31

All Krylow formulas use either of the two principles the Orthogonal ResidualCondition (OR) or the Minimal Residual condition (MR). The most famous ORformula is probably the conjugated gradient method (CG), even though it onlyworks for symmetric positive definite matrices. The most useful MR method oftoday is likely the generalized minimal residual method GMRES.

The Minimal Residual condition can be formulatated: Find xk such that:

< rk, Au >= (Au)T (b−Axk) = 0, ∀u ∈ Uk ⇔rk ⊥ AUk ⇐⇒||rk|| = ||b−Axk|| = min

z∈Uk

||b−A(x0 + z)|| = minz∈Uk

||r0 −Az||

The Orthogonal Residual condition can be formulatated: Find xk such that:

< rk, u >=< b−Axk, u >= uT (b−Axk) = 0, ∀u ∈ Uk ⇔rk ⊥ Uk

Note that for the OR case if rk �= 0 the residuals are mutually orthogonal andtogether span Uk+1.

< ri, rj >= δij < ri, ri >

Uk+1 = lin{r0, r1, ..., rk}

Conjugated Gradient method (CG)

The Conjugated Gradient method (CG) is a classic fundamental method onwhich many others are based. Basically in each iteration step the new valuexk+1 equals the previous value plus αk times a directional factor pk equal tothe residual plus a fraction β of the previous direction factor. αk and βk aredetermined by the orthogonal residual condition.

xk+1 = xk + αkpk

pk+1 = rk + βk+1pk

Even though it only works for positive definite matrices and is thus not of in-terest for this specific problem, a description of the algorithm for CG is addedfor the curious. This and all subsequent algorithms is from Hanke [14].

32

Algorithm CGp0 = r0 = b−Ax0

for k = 1 : kmax

αk =rTk rkpTk Apk

xk+1 = xk + αkpk

if OK

BREAK

endif

rk+1 = rk − αkApk

βk+1 =rTk+1rk+1

rTk rkpk+1 = rk + βk+1pk

endfor (2.18)

The number of operations needed are: 1 matrix vector product, 2 scalar productsand 3 axpy operations. One axpy operation is defined a · x + y.

Algorithm CGSu1 = p1 = r0 = b−Ax0

Choose r∗0 : < r0, r∗0 >�= 0

p∗1 = r∗0for k = 1 : kmax

αk =r∗T0 rk−1

r∗T0 Apkqk = uk − αkApk

xk = xk−1 + αk(uk + qk)if OK

BREAK

endif

rk = rk−1 − αkA(uk + qk)

βk+1 =r∗T0 rk

r∗T0 rk−1

uk+1 = rk + βk+1qk

pk+1 = uk+1 + βk+1(qk + βk+1pk)endfor (2.19)

The number of operations needed are: 2 matrix vector products, 2 scalarproducts and 7 axpy operations.

33

Conjugated Gradient Squared method (CGS)

An extension of of the CG approach is the Conjugated Gradient Square method(CGS), which usually is a rapidly converging Krylov type method. It works forlinear systems with quadratic matrices even if they are neither symmetric norpositive definite. One of the nice features of CGS is that it does not demandthe presence of the transpose of the matrix, AT .

GMRES

The Generalized Minimal RESidual method (GMRES) is based on the minimalresidual principle. The matrix Vk consists of vectors v1, v2, · · · , vk which togetherspan Uk. Define:

v1 =r0

||r0|| =r0β

Uk � z =k∑

j=1

yjvj = Vky

The nestled for-loop in the GMRES algorithm is a orthogonalization schemecalled Arnoldi with Modified Gram-Schmidt. It produces a recursion whichgives that:

AVk = Vk+1Hk

The matrix Hk is on Hessenberg form, which means that all non zero elementsare contained in either the upper triangular part or on the first subdiagonal.

Hk =

h11 h12 h13 . . . h1k

h21 h22 h23 . . . h2k

h32 h33 . . . h3k

. . . . . ....

hk,k−1 hkk

hk+1,k

In GMRES, the value of xk does not have to be evaluated before this is reallywanted.

The objective function of the minimal residual condition can be written: Findz = Vky ∈ Uk which minimizes J(y):

J(y) = ||b−Axk|| = ||b−A(x0 + z)|| = ||r0 −Az|| =||βv1 −AVky|| = ||βe1 − Hky||

The system which has to be solved for J is usually very small compared to thedimension of the problem.

34

GMRES stops when Arnoldi does, and Arnoldi terminates iff the computedsolution is exact.

Algorithm GMRESr0 = b−Ax0

β = ||r0||v1 = r0/β

for j = 1 : mwj = Avj

for i = 1 : jhij =< wj , vi >

wj = wj − hijvi

endfor

hj+1,j = ||wj ||if hj+1,j = 0 then m = j goto(∗), endif

vj+1 = wj/hj+1,j

endfor

(∗) Compute ym asminimizer of ||βv1 −AVky||xm = x0 + Vmym (2.20)

The number of operations needed are: 1 matrix vector product, j+1 scalarproducts, j axpy operations and 1 vector scaling.

NonLinear-GMRES

For non-linear systems a modified GMRES method called NonLinear GMRES(NL-GMRES) can be used.

Nonlinear systems of equations on the form F (x) = 0, are (as previously noted)usually solved with Newtons method. An initial guess x0 is assumed.

F ′(xk)zk = −F (xk)xk+1 = xk + λkzk

A linear system must be solved for each step of the iteration. It can either besolved with direct methods, which is the case in the main program, or in com-bination with Krylow type methods, which is the case in the additional part ofthe program.

Preconditioning is an important feature that will be examined in the nextchapter. Anyhow, if GMRES is used on the linear system only right side pre-conditioning (F ′M−1) · (Mz) = −F can be applied, otherwise the residual willbe improperly scaled and totally meaningless for practical use.

35

In the case of Krylow type methods, for a well chosen small σ �= 0 the Jac-obian of the right side of the linear equation can be approximated by:

F ′(xk)z ≈ (F (xk + σz) − F (xk))σ

The nonlinear system of equations can be viewed as an equivivalent minimationproblem. If F is the nonlinear mapping D ⊆ �n → �n and x ∈ D:

F (x) = 0 ⇔ minx∈D

||F (x)||22

= minx∈D

g(x) = 0

A necessity for finding a minimum F (x) is that a descent direction can beobtained. At position x, that is a direction which for sufficiently small α thefollowing holds:

g(x + αp) < g(x) ⇔ g′(x) = ∇xg(x)T p < 0

For the Newton direction we have that:

g′(x) = −(FTF ′)(F−1F ) = −||F ||2 < 0

This means that the Newton direction is always a decent direction. The direc-tion obtained by GMRES is decent if it is close enough to the Newton one orif z0 = 0. Note that the optimal solution are seldom found, while just a fewGMRES steps are taken.

The residual of the GMRES iteration number k is denoted rk = −(F + F ′zk).By construction the orthogonal relation < rk, F

′zk > is present.This gives that:

g′zk = FTF ′zk = −FT (rk + F ) =− < rk,−(rk + F ′zk) > −||F ||2 = ||rk||2 − ||F ||2

Consequently, zk is a decent direction if ||rk|| < ||F ||. For GMRES we havethat ||rk+1|| < ||rk|| = ρk ∀k ≥ 1. If z0 = 0 is the startvalue, decent directionswill always be found. The whole algorithm is stated before preconditioningis considered. A good Matlab NL-GMRES implementation called nsolgm.m ispresented by Kelly [17].

36

Algorithm NL-GMRES(1) Choose x0 ∈ D, ε0 > 0k = 0f0 = F (x0)(2) r0 = −fk

Choose σ = σ(xk, fk)β = ||r0||v1 = r0/β

for j = 1 : mSolve Mvj = vj

wj =(F (xk + σvj) − fk)

σfor i = 1 : j

hij =< wj , vi >

wj = wj − hijvi

endfor

hj+1,j = ||wj ||if hj+1,j = 0 then m = j goto(∗), endif

vj+1 = wj/hj+1,j

Compute residualnorm ρj

if ρj ≤ εk then m = j goto(∗), endif

endfor

(∗) Compute ym asminimizer of ||βv1 −AVky||Solve Mzk = Vmyk

(3) choose λk ∈ (0, 1| : ||F (xk + λkzk)|| � ||F (xk)||xk+1 = xk + λkzk

fk+1 = F (xk+1)(4) if OK BREAK endif

k = k + 1Choose εk > 0goto(2) (2.21)

2.11.5 PreconditioningOne way to improve the convergence rate and efficiency of iterative methodsis to use preconditioning matrices. The system is multiplied with some otherproperly chosen matrix with the aim to receive a new system which is easierto solve because of smaller condition number. In fact, a good preconditioneris often a necessity for iterative methods to work at all. To be of any use, thepreconditioning matrixM must have something in common with A. The rate of

37

similarity between A andM strongly effects the efficiency of the preconditioning.

Preconditioning can be applied from the right and from the left. This cor-responds to the relations below.

Left : M−1Ax = M−1b

Right : AM−1Mx = AM−1y = b

Five types of strategies will be described.

Jacobi (diagonal) preconditioning P = diag(A) demands very few additionalcalculations to be used. It is simple but often useful. It is used in the iterativevariant of the implementation.

Incomplete factorisation preconditioning exploits the fact that the inverse ofa matrix can be computed by using LU- och Cholesky factorisation, but thatthe computional cost of a complete factorization is too high in practice. Anapproximative factorization can be performed if all small elements in the copyof A are set to zero during factorization and only the dominant terms are saved.If the distribution of the magnitudes of the elements range over a large interval,this method is very efficient.

Polynomial preconditioning uses the fact that the applying a few numbers of astationary iterative method xk+1 = Bxk + c to the system can be viewed as aninverse of A. The iterative procedure approximates an expression for x, and sodoes an operation A−1Ax = A−1b ≈ Ix = x. One common choice with acceler-ator terms is A−1 ≈ ℘(q) =

∑kj=0 αkjqj ,

∑kj=0 αkjqj = 1,∀k. Another one is

based on Neumann series. If the spectral radius of the matrix G is less than onethe following holds: (I − G)−1 =

∑∞j=0 Gj . The first terms can be used as an

approximation of A−1. Polynomial preconditioning could be of interest in thiscase.

Domain decomposition preconditioning is useful if the nonzero structure of Acan be divided into submatrices of simple shape like diagonal matrices, raw andcolumn vectors and zeros. The inverse of A can be formed from the inverses ofthe submatrices.

Least-squares preconditioners, which are only useful for symmetric positive def-inite problems, notices that optimaly (I−P−1) = {0} yields. An approximationps ≈ A−1 is obtained by minimizing ϕ(x) = (1 − ps(x) · x) for x.

Next comes a description of the implementation of the model and methods.

38

Chapter 3

Implementation

3.1 Choice of Programming LanguageThe determing whether the computer implementation should be programmed inC, Fortran 90 or Matlab 6.5 has consumed a lot of time. All three have both nicefeatures and drawbacks. C is fast and widely adopted. Fortran is fast, handlesmatrices efficient and has a glorious history as the dominating programminglanguage for scientific calculations, for decades. Matlabs routines for numericallinear algebra are very efficient, it is userfriendly and has powerful routines forgeneration of graphs. Matlab is unfortunately an interpreted language, whichmeans that the code must be interpreted before it is executed, which makes itslower.

At the beginning of the project I intended to use Fortran 90 and spent a lotof time learning Fortran, coding the project in Fortran, solving Fortran specificprogramming problems and searching for pleasant Fortran functions. I foundout that the debugging of Fortran code was a little bit less efficient than theMatlab equivalent. Error information in Matlab is more specific and printingthe values of a symbol on the screen corresponds to removing a “;” and not towriting a long command sequence. Realizing that Matlab or some other graphgenerating program had to be used for visualizing the results I decided to makeboth a Matlab and a Fortran version of the code.

At that time I also compared the speed of different operations in the two lan-guages. The main results was that the speed for numerical linear algebra wasabout the same for all of the languages, loops in Matlab < 6.5 was about 100times slower than in Fortran but only a few (5-10) percent slower than Fortranin Matlab 6.5. Function calls seems to be the large bottle neck left in Matlab.Therefore the number of function calls in the program is as low as possible.

When the Matlab version was finished I decided not to complete the Fortranversion even though a lot of work was spent on it already. The main reasonswas that in matlab, all functions of a certain type, for example iterative solversfor linear systems of equations or ODE solvers, has the same syntax. Thereby,it is much simpler to update a Matlab version than a fortran one. Furthermore,

39

if a Matlab solver is improved, it still keeps its old name, which means thatold programs are automatically updated if they use that structure. Finally, theimprovement rate and speed acceleration of Matlab seems to be higher than forC and Fortran. Neither C nor Fortran will disappear in a long time, but I guessthat Matlab will have a dominating position in the future.

If somebody would like to make a Fortran version for the present type of prob-lems in spite of that, UMFPACK and SUPERLU are nice for LU factorizationand Byrnes et al DVODPK [3] or eventually Browns et al DASPK [2] can beused as iterative (GMRES) ODE solvers.

3.2 General Description of the ImplementationThe program MainDrop.m with subprograms is an implementation of the nu-merical treatment of a model which describes how droplets of liquid behaveon the ground. The programming language used is Matlab 6.5. The mainprogram uses finite differences for the spatial discretization, NDFs in combin-ation with Newtons method for the time-discretization and LU-factorizationwith UMFPACK for solving the linear systems of equations. Future versions ofMatlab will most likely automatically result in the use of NL-GMRES insteadof UMFPACK. Graphs are generated automatically when running the programand some additional features like versions with Krylow type methods are alsoimplemented.

3.3 How to Use the ProgramIn order to use the program, begin by checking that your Matlab version isversion 6.5 or better, otherwise the for-loops will take too long time. Check alsothat all files are in the same catalog. Next, open the file indata.m and add newvalues of the constants. The syntax is on the form:

----------------------------------------------------

function[phi,theta,...,alpha1]=indata(inputnumber)if (inputnumber==1)

phi= X; %Total soil porositytheta= X; %Soil water content...r0= X %Initial radius before phase 1

elseif (inputnumber==2)phi= Z..

end---------------------------------------------------------

40

Next, open the file MainDrop.m. That is the control file and the only file thatshould be opened during normal execution of the code. Dont touch the commandsequences at the bottom of the file. The syntax of the raws of MainDropp.mthat should be considered is:

---------------------------------------------------------% MAIN FILE DROP PROGRAM (c) FRITJOF NILSSON 2004clear,clf %Remove old constants and graphs

Dimension=X; %Dimension of problem, (1) or (3)Tmax= X; %Timespan in seconds

Ziend= X; %Gridpoints in z-direction with initial conc.kz= X; %Ziend*kz=gridpoints in z without initial conc.Rjend= X; %Gridpoints in r-direction with initial conc.kr= X; %Rjend*kr=gridpoints in r without initial conc.

inputNr= X; %Number of current set of constantsreltol= 1e-3; %Relative Tolerance, accuracy criteriaabstol= 1e-6; %Absolute Tolerance, accuracy criteria

standard= X; %Run mainprogram? 1==’yes’, 0==’no’,optional= XX; %1X=BDF,2X=NDF,X1=LU,X2=GMRES,X3=CGS,X4=NL-GMRES,0=no---------------------------------------------------------

When running the program, control that program does not write "Warning,Infinite Radius/Depth Reached". Otherwise the Fr curve and the Fev curvewill show too low values at distant times. The values of Ziend and Rjend arepreferably chosen somewhere between 3 and 20. If the timespan is one day andthe ground is of concrete, kz=14 and kr=3 are good choices. With 24h on sand,try kz=80 and kr=20.

If additional experimental curves are to be added to the graphs, add a columnvector with times and a column vector with fraction values at column “nr”,respective “nr+1” in the MS-EXCEL file Experimental.xls.

3.4 Structure of the FilesThe main program MainDropp.m is the file from which everything is controlled.From there a call to either main1D.m or main2D.m is done. Almost all sub-programs exists in a 1D and a 2D version. Usersupplied indata is fetched fromindata.m. If the indata constants rmax and Zmax are set to zero, they are calcu-lated in the file phase1o2.m. Initial values are calculated in startvalues1D.m (or2D). The optional file optional.m executes some interesting but non-central pro-cesses like generation of sparsity patterns. Then comes the intrinsic ODE solverode15s.m which uses the files jacobian1D.m,cdot1D.m, stencilCalc1D.m and sten-cilCalch1D.m. The output is processed in the file evaporation1D.m and graphsare constructed in graphs1D.m. Alternative ODE-solvers are implemented inthe optional file alternatives1D.m which calls the same files as ODE15s.

41

35,95

MainDrop.m phase1o2.m

main2D.mmain1D.m indata.m

startvalues2D.m

optional2D.m

ode15s.m

alternatives2D.m

jacobian2D.m

cdot2D.m

stencilCalch2D.m

stencilCalc2D.m

evaporation1D.m

graphs2D.m

✲

✲

✲

✲

✲

✘✘✘✘✘✘✘✘✘✘✘✘✘✘✿

�

❆❆❆

❆❆❆✁

✁✁✕

✄✄✄✄✗

✄✄✄✄✄✄✗

✄✄✄✄✄✄✄✄✄✄✗

✄✄✄✄✄✄✄✄✄✄✄✄✗

✄✄✄✄✄✄✄✄✄✄✄✄✄✄✄✗

✄✄✄✄✄✄✄✄✄✄✄✄✄✄✄✄✄✄✗

❄

Figure 3.1: Structure of the files of the program MainDrop.m.

3.5 Description of the SubprogramsThe subprograms of the project are here described briefly. When both a 2Dand a 1D version exists, the 2D version is described but all 2D specific parts areprinted in italics. The functionality is the same. Most of the information fromthe file indata is handled global.

MainDrop.mInput :(-)Output : (inputnumber, Tmax, Ziend, kz, Rjend, kr, standard, reltol, abstol,NDFhm, BDFhm, Trapetshm)Function: Gives the user the opportunity to choose the functions and files whichare to be used. Calls either main3D or main1D.

main2D.mInput : (inputnumber, Tmax, Ziend, kz, Rjend, kr, standard, reltol, abstol, ND-Fhm, BDFhm, Trapetshm)Output : resultsFunction: Intermediate buffer which structures the flow of operations. Callsat most indata.m, startvalues3D.m, optional3D.m, ode15s.m alternatives3D.mand graphs3D.m.

indata.mInput : inputnumber

42

Output : phi, theta,...,alpha1results (see above)Function: Contains all userdefined constants.

phase1o2.mInput : r0Output : rE, Rmaxabs, ZmaxabsFunction: Calculates the steps of phase 1 and 2 of the model. This includesfunction calls to the ODE-solver ode45.m and a function F. Note that the func-tion only is called if Zmaxabs and Rmaxabs from indata are zero and r0 is not.

startvalues2D.mInput : Ziend, kz, Zmaxabs, Rjend, kr, Rmaxabs, phi, theta, densitychem, cgsOutput : c0, Zimax, dz, z, Rjmax, dr, rFunction: Interpolates a curve from experimental data, integrates the curve andgenerates a initial grid with nonzero concentrations in a rectangle from (1,1) to(Rjend,Ziend). Rjmax and Zimax are the total numbers of gridpoints in eitherdirection.

ode15s.mInput : c0,Jpattern, cdot.m, jacobian.m, Tspan, options, Zimax, dz, Rjmax, dr,rOutput : T, YFunction: Intrinsic Matlab ODE-solver with NDFs for stiff problems.Returns avector with times and a vector with concentrations.

alternatives2D.mInput : c0, cdot.m, jacobian.m, Tspan, Zimax, dz, Rjmax, dr, rOutput : T, YFunction: Contains homemade implementations of NDFs, BDFs and CN witheither LU or GMRES. Returns a vector with times and a vector with concen-trations.

jacobian2D.mInput : t, c, Zimax, dz, Rjmax, dr, rOutput : JFunction: Calculates the Jacobian in each timestep. Calls both stencilCalc3D.mvariants.

cdot2D.mInput : t, c, Zimax, dz, Rjmax, dr, rOutput : cdot, WFunction: Approximates the time derivate of the concentration in each timestep.cn = W (cn)cn. Calls stencilCalc3D.m.

stencilCalc2D.mInput : c, Zimax, dz, Rjmax, dr, rOutput : K1,H1,H2K2,H3,K3, K2plus, K2minus, H2plus, H2minus, Dgs, Dlcs,RgammaInvert, RGInvert, RLInvert, AFunction: Calculates the stencils of the spatial discretization of each timestep.

43

stencilCalch2D.mFunction: In principle the same function, input and output as stencilCalc3D.mexcept from a small disturbance h in the concentration.

evaporation1D.mInput : Y, T, dz, Zimax, r, Rjend, RjmaxOutput : Fr, Fev, d_dtFr, d_dtFev, Fluxtotal, Fluxlocal, AKRG, massFunction: Postprocesses the gathered information from the ODE solver for thewhole timeinterval.

graphs2D.mInput : Zimax, Rjmax, Zmaxabs, Rmaxabs, Y, T, Fr, Fev, d_dtFr, d_dtFev,Fluxtotal, Fluxlocal, nr)Output : graphsFunction: Generates graphs of among others concentrations profiles, evapora-tion curves and fluxprofiles.

44

Chapter 4

Results and comparison withexperimental data

The results of this work can be divided into the following subcategories: Resultsfrom the examination of matrices, from the use of conservation laws, from com-parisons of different ODE and equation system solvers and from comparisonswith experimental data. A final category is general results of the project

4.1 Results from the Examination of MatricesThe eigenvalues of the matrices are of great importance for the efficiency of allnumerical methods for solving systems of equations. Remember that systemsare considered stiff if some eigenvalues of the Jacobian is negative and large incomparison with the reciprocal of the timespan of interest. Here comes a tablewith the largest and smallest eigenvalue of the initial Jacobian from one sandand one concrete test. The number of gridpoints are 350 (1D) respective 20000(2D).

Eigenvalue of J concrete 1D sand 1D concrete 2D sand 2Dlargest -0.8554 -3797 -0.3885 -1775smallest -0.2930 -1303 -0.3849 -1721

The magnitudes of the eigenvalues are much larger than the reciprocal of thetimeinterval 24*3600s. Thus the systems are really stiff, especially the sandcase.

4.2 Results from the Using of Conservation LawsSome important test of the implementations can be performed through the useof conservation laws. The first check is to confirm that if no chemical degrada-tion is present, the amount of evaporated chemical equals the loss of chemicalin the ground. A graph of Fr − (1 − Fev) for one concrete example in 2D isshown. The Time interval is 24h and the number of initially filled gridpointsare 5. Note that the differences are small compared to one.

45

0 1 2 3 4 5 6 7 8 9

x 104

−1.5

−1

−0.5

0

0.5

1

1.5

2

2.5

3x 10

−3 Fr−(1−Fev) vs Time

Time (s)

Fr−(

1−Fe

v) (−

)

Figure 4.1: Differences of Fraction Remaining minus Fraction Evaporated during24h on concrete (3D). As in all tests, the differences are a few promille of thetotal amount. 1 ≥ |Fr| > 0.2. Absolute tolerance is 106 and relative is 10−3.

If the differences are larger than a few promille something is wrong. If the Fr

becomes lower than Fev at the end of the interval, the grid is probably too small,which leads to a loss of mass at infinity. If Fev rapidly grows larger than Fr

during the first hour and then the proportions remains approximately constant,then Fev is wrong because of a summation with too large timesteps. In the 2Dsand case it can be tempting to save too few values because lack of memory.

Next test is to see that mass conservation holds during spread out. In order toconfirm this, the evaporation rate and degradation term are set to zero for thewhole timeinterval. As can be seen in the fraction remaining graph of a 24h 2Dconcrete curve, practically no mass dissappears because of discretization errors.

The initial mass of the 1D case should equals the 2D mass. From a concreterun for 1D: M0=1.1871 kg, and for 2D: M0=1.1842 kg .

4.3 Results from the Comparing of Different Solv-ers

Two test are performed in order to compare different solvers. The first com-pares variable size ODE-solvers on concrete with 100 (1D) respective 2000 (2D)gridpoints and on sand with 350 gridpoints (1D). All tests are for 24h.

• ode15s NDF (implicit)• ode15s BDF (implicit)• ode23tf BDF/Trapets (implicit)• ode23s Rosenbrock (implicit)

46

0 1 2 3 4 5 6 7 8

x 104

0.9998

0.9999

0.9999

1

1

1.0001

1.0001

Fraction Remaining vs Time

Time (s)

Frac

tion R

emain

ing (−

)Numerical ModelExperimental ValuesChemical Fraction

Figure 4.2: Fraction remaining without evaporation for concrete 24h (2D). Thediscretization errors are too small to be notized.

• ode45 Runge-Kutta (explicit)

Note that NDFs are fastest on the present matrices and that the explicit methodworks bad on this stiff problem. NDFs have the lowest number of timesteps onconcrete and Rosenbrock on sand. In the figure the execution time is abbrevatedt, the number of timesteps with nr, concrete with c and sand with s.

Method 1D c. t 1D c. nr 3D c. t 3D c. nr 1D s. t 1D s. nrode15s NDF 0.516 184 15.031 231 11.29 1597ode15s BDF 0.656 216 18.797 304 12.50 1874

ode23tf 1.312 245 30.890 273 14.72 1122ode23s 1.266 202 34.953 227 13.31 1060ode45 4.9210 2929 84.041 1757 210.39 23937

A couple of graphs which shows the step length in time for the variable stepNDFs are shown. Note that small steps are taken initially and large are takenat the end.

The second test is for the solving of equations. A couple of solvers havebeen used on a concrete case in 1D and 2D. 100 gridpoints are used in 1D and2000 in 2D. The tested methods and the results are as follows.

• UMFPACK LU+colamd (direct)• GMRES no preconditioning (iterative)• GMRES diagonal preconditioning (iterative)• GMRES incomplete (lu(θ)) preconditioning (iterative)• CGS diagonal preconditioning (iterative)

47

0 20 40 60 80 100 120 140 160 180 20010

−4

10−3

10−2

10−1

100

101

102

103

104

Size of dt vs Timestepnumber k

Timestepnumber (−)

dt

Figure 4.3: Size of dt versus timestepnumber for 24h on concrete (1D). Thesteps are very small in the beginning but becomes soon larger.

0 200 400 600 800 1000 1200 1400 160010

−7

10−6

10−5

10−4

10−3

10−2

10−1

100

101

102

103

Size of dt vs Timestepnumber k

Timestepnumber (−)

dt

Figure 4.4: Size of dt versus timestepnumber for 24h on sand (1D). Note thatthe steps are ultimately smaller than in the concrete case.

48

• NL-GMRES no preconditioning (iterative)

Method 1D(a) 2D(a) 1D(b) 2D(b)UMFPACK 0.844 9.109 0.609 6.875

GMRES no prec. 1.515 8.875 0.640 5.891GMRES diag prec. 1.734 10.047 0.703 6.406NL-GMRES no prec. 2.969 18.532 0.500 4.375

In the (a) runs, NDF2 has stepped 48h with 125 iterations. 16 < ∆t < 2048.Except for NL-GMRES, only one Newton iteration is done each time step. Inthe (b) tests, 100 steps with NDF2 with ∆t = 1 are taken.

Incomplete LU factorization with 10−6 < θ < 10−3 did not work on this par-ticular problem even though it is usually the best choice. The CGS did neverconverge, which is a bit strange while it usually is rather good. For small timesteps, NL-GMRES was fastest but otherwise UMFPACK was fastest. Anyway,since more than one Newton steps may be used in NL-GMRES due to errorcontrol, the comparison is not really fair. The diagonal preconditioning doessuprisingly not improve the speed of GMRES on this problem.

According to Hanke [14], direct solvers are usually best suited for one-dimensionalproblems (bandwidth 3) and iterative solvers for problems with three space di-mensions (bandwidth n2). Iterative methods are usually just slightly betterthan direct methods on two-dimensional problems with bandwidth n. The testdoes not contradict this.

4.4 Results from the Comparison with Experi-mental Data

Three important checks can be done with the present experimental data; Com-parisons with evaporation curves, initial flux rate and general concentrationprofile. Before this is done a short reiteration of how the wind tunnel experi-ment was done is presented.

4.4.1 Wind Tunnel TestsThe wind tunnel test are made at Prince Mauritz Laboratories in the Nether-lands. In each run three droplets of some liquid are dropped on a substrate.The evaporation of each droplet are measured gravimetrically during the wholeprocess and in the end the remaining mass is measured chemically as a control.The soils that are used in the test are concrete, sand, artificial grass and glass.The diffusion into glass is to low to be of interest and the grass is not symmetricenough for the task, thus sand and concrete remains.

Most of the runs are in one way or another erronous; the two evaporation curvesdo not coincide, the measurements are irregular or the curves are contaminatedby evaporated water which implies that the chemically and the gravimetricallydetermined endpoints do not coincide. The runs with least amount of errors are

49

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

x 105

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1Fraction Remaining vs Time

Time (s)

Frac

tion R

emain

ing (−

)Numerical ModelExperimental ValuesChemical Fraction

Figure 4.5: Fraction Remaining for concrete fx310100 48h (1D). The star marksthe more reliable chemically determined remaing fraction of substance.

chosen as reference material for the model. In practice we have chosen the sandtests fx10069 and fx01069 and the concrete tests fx310100 and fx29129. Thedroplets in the test consist of methyl salicylate.

4.4.2 Actual ComparisonsEvaporation curves for two sand experiments and two concrete experiments in1D and 2D will follow. It can be notized that the 2D curves fits the experimentaldata even better than the 1D curves.

The concentration profiles coincides with experimental shapes by Baines andJames [1]. After 48h on concrete the curves below are obtained. The profile inthe center of the 2D graph looks similar.

Finally it can be shown that the initial rate of evaporation in the 1D modelagrees with the experimental value. Because of the damping term, the test isnot really applicable for 2D, but it can be seen that if horizontal movements aretaken away and the damping is zero, the 2D values coincides too.

4.5 General Results of the ProjectAn enhanced 1D program has been constructed and a 2D model has been de-veloped. Both programs can use either direct or iterative methods when solvinglinear system of equations.

An increase of speed has taken place in the 1D model. In spite of the factthat Fortran is generally faster than Matlab, the new Matlab version is in prac-tice more than 100 times faster than the old Fortran version. For 1D curvesnot more than a few seconds are needed for the execusion time and in 2D some-

50

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

x 105

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1Fraction Remaining vs Time

Time (s)

Frac

tion R

emain

ing (−

)

Numerical ModelExperimental ValuesChemical Fraction

Figure 4.6: Fraction Remaining for concrete fx310100 48h (3D).

0 1 2 3 4 5 6 7 8 9

x 104

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1Fraction Remaining vs Time

Time (s)

Frac

tion R

emain

ing (−

)

Numerical ModelExperimental ValuesChemical Fraction

Figure 4.7: Fraction Remaining for sand fx10069 48h (1D)

51

0 1 2 3 4 5 6 7 8 9

x 104

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1Fraction Remaining vs Time

Time (s)

Frac

tion R

emain

ing (−

)

Numerical ModelExperimental ValuesChemical Fraction

Figure 4.8: Fraction Remaining for sand fx10069 48h (3D)

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5−4

−3.5

−3

−2.5

−2

−1.5

−1

−0.5

0x 10

−4 Concentration at t=86400 vs Depth

Concentration (kg/m3)

Depth

(m)

Figure 4.9: Concentration profile after 48h on concrete (1D). It coincides withexperimental data.

52

0 0.5 1 1.5 2 2.5 3 3.5

x 10−3

−4

−3.5

−3

−2.5

−2

−1.5

−1

−0.5

0x 10

−4 Concentration at t=86400(s), cmax

=5.0066

Radius (m)

Depth

(m)

Figure 4.10: Concentration vs depth and radius after 48h on concrete (3D). Theconcentration profile along the symmetry axis coincides with experimental data.

where between a few seconds and half an hour is needed on a 1.8 GHz, 256 MbWindows PC.

Improvements for speed includes for example vectorizing loops, avoiding unnec-cessary recalculations, reducing the number of function calls and using variablestepsize NDFs and efficient equation solvers.

The number of discretization errors has been reduced notably, which allowsthe use of fewer discretization points. This increases the speed further.

The difference of the curves for Fr and 1−Fev are now negliable. Mass conser-vation holds during spread out which can be checked by setting the evaporationrate to zero for the whole process. The initial flux for both the 1D and 2D modelcoincides with the experimental value if horizontal movements are forbidden andα1 is zero.

Evaporation curves for the 1D model has good agreement with experimentaldata. The average distance between the 2D model curves and the experimentalcurves is approximately half the average distance between the 1D model curvesand the experimantal curves. If δk is set to zero, the 1D concrete curves are notaffected but the 1D sand curves becomes much better and the execusion timesare four times faster.

The program is easy to use and generates automatically a large amount ofuseful graphs and warning messages. Furthermore, it is coded in a way that willtake advantage of future improvements of Matlab and build-in Matlab functions.

This report documents the project. Comparisons with the windtunnel experi-ments of Oeseburg et al[18] are never published outside FOI before.

53

Chapter 5

Discussion

In the 2D model a damping function kα occures in the upper boundary condi-tions. A constant kα ≈ 0.7 works fairly well for evaporation curves and similar.Anyhow, it is an simplification that gives too large evaporation rate from theedge of the droplet and too small in the center, even if the total evaporation iscorrect. The background is that the formula of the upper boundary condition isdeveloped from experiments with drops with constant radius and no horizontalmovement. The formula gives that the relative flux (flux/concentration) is lowerfor points on the surface with high concentration than for those with low, evenif the reverse holds for the absolute flux.

In order to find a damping function that describes the two-dimensional fluxmore properly, experiments that measures the evaporation as a function of theradius must be performed. A testplant could consist of a soil above whichseveral concentrical hollow cylinders are placed. In each cylindrical segment adetector could measure the evaporated amount from a ring shaped segment onthe ground and the sum of the measurements are the total evaporated mass.If a droplet is placed in the center of the site, data about the relative amountof evaporation as a function of the normed radius could be received, and thatcould be used for finding a better damping function.

It can be noted that if the physical constant δk is set to zero instead of one,the 1D curves suprisingly agrees much better with the windtunnel test. Thephysical background has not yet been found, but based on the empirical tests,the factor is suggested to be taken away. This also reduces the run times by 75 %.

More comments of the physics of the model can be found in Winter [27].

Neither the avoidance of recalculating identical stenciles nor a graphical in-terface is implemented. This could be done in a future improvement.

The choice of Matlab instead of Fortran as programming language can be dis-cussed, but as Matlab improves so rapidly, it should be an acceptable choice.Moreover, Matlab allows for an easy implementation of a grafical user interface(GUI) if this should become necessary.

54

0 0.2 0.4 0.6 0.8 1 1.2 1.4−7

−6

−5

−4

−3

−2

−1

0

1x 10

−4 Concentration at t=86400 vs Depth

Concentration (kg/m3)

Depth

(m)

Figure 5.1: Concentration profile after 24h in sand (1D, delta=0). The edge ofthe droplet is just about to reach the lower boundary.

0 1 2 3 4 5 6 7 8 9

x 104

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1Fraction Remaining vs Time

Time (s)

Frac

tion R

emain

ing (−

)

Numerical ModelExperimental ValuesChemical Fraction

Figure 5.2: Fraction remaining of the chemical after 24h in the sand (1D, δk =0). Note that the curve is very close to the experimental data.

55

56

Chapter 6

Nomenclature

a Air filled porosity=φ− η − θw m3m−3

cg Concentration of chemical vapour phase in soil air kg m−3

cg0 Vapour concentration at the surface kg m−3

cgsat Volatility=saturated vapour concentration kg m−3

cL Concentration of dissolved chemical in soil water kg m−3

cs Concentration of absorbed chemical kg m−3

c Total concentration of chemical in the ground kg m−3

δk Constant in Dslc. Should be zero! (−)

Dg Diffusion coefficient for chemical vapour in the air m2 s−1

Dsg Effective vapour diffusion coefficient in the soil m2 s−1

Dsl Effective diffusion coefficient for dissolved chemical in the soil m2 s−1

Dwaterl Diffusion coefficient for chemical in water m2 s−1

A Effective diffusion coefficient m2 s−1

B Effective convection coefficient m s−1

Dslc Diffusion coefficient for liquid phase chemical in the soil m2 s−1

fs Effective porosity of the substrate=φ− θw0 m3 m−3

fη Partitioning function for vapour and liquid phases=cg/η kg m−3

F Massflux of chemical kg m−2s−1

F Mean flux of chemical from the surface kg m−2s−1

KD Linear adsorption coefficeient=cs/cL m3 kg−1

KH Non-dimensional Henry law constant (-)m Mass kgr Radius of spreading droplet patch mr0 Radius of spherical droplet mrmaxabs Horizontal radius of the wet spot at the end of sorption mre Droplet radius after impact mre,min Smallest value of re for settling velocity u0 = 0 mrtol Constant for radius determining (about 0.01) (-)rjend Number of iteration points with conc. on surface. (-)Rg Partitioning function for total and vapour concentration=c/cg (−)R−1

l Inverse of Rg (−)Rc Horizontal radius of liquid spherical cap mRl Partitioning function for total and dissolved concentrations= c/cL (−)R−1

l Inverse of Rl (−)Rη Partitioning function for total and liquid concentrations= c/η (−)R−1

l Inverse of Rη (−)t Time stmaxabs Total time phase 3 s57

t1 Time for sorption to be completed sT Temperature Kv Wind velocity m s−1

v0 Settling velocity of falling droplet m s−1

v∗ Friction velocity m s−1

vz Convective vertical velocity of soil water m s−1

V0 Initial volume of the spherical droplet. m3

Vc Volume of liquid spherical cap m3

Vs Liquid volume contained in the wetted zone m3

zf Vertical distance from the soil surface to the wetting front mzf1 Depth of the wetting front at end of sorption mz Vertical position below surface mr Horizontal position mε Time constant for degradation of the chemical in the soil s−1

εT Truncation error (−)εm Precision of the actual computer, machine epsilon (−)εk Small value for avoiding division by zero >

√εm (−)

γ Surface tension of liquid chemical N s−1

η Liquid chemical concentration m3 m−3

ηcrit Limit for liquid concentration, (cg = cgsat) m3 m−3

θw Volumetric content of soil water concentration m3 m−3

θw0 Volumetric soil water content before attachment of chemical m3m−3

κv air velocity shear near the surface s−1

κ NDF koefficient (−)κ∞ Condition number maximum-norm (−)κ2 Condition number 2-norm (−)λ Eigenvalue (−)µ Dynamic viscosity of liquid chemical kg m−1s−1

ν Kinematic viscosity of atmospheric air m2 s−1

ρb Soil bulk density kg m−3

ρc Density of liquid chemical kg m−3

ρw Liquid denisity of water kg m−3

ρs Density of water solution kg m−3

σ Penetrability =zf/√t =konstant m s−1/2

φ Total soil porosity m3 m−3

χ Spectral radius (−)zmaxabs Mean depth of the wet layer at the end of sorption m

58

Bibliography

[1] W.D.Baines and D.F.James, Evaporation of a Droplet on a Surface, Ind.Eng. Chem. Res. 33. 411-416, 1994.

[2] P.N.Brown A.C.Hindmarsh L.R. Petzold, DASPK Package - 1995 Revision,1995.

[3] G.Byrne A.Hindmarsh and P.Brown, DVODPK: Variable Coefficient Or-dinary Differential Equation Solver with the Preconditioned Krylow MethodGMRES for the Solution of linear Systems, Illinois Institute of Technology,2002.

[4] N.Carlsund, Lecture Notes: Finite Element Method, KTH Department ofNumerical Analysis and Computer Science, 2003.

[5] E.L.Cussler, Diffusion-Mass Transfer in Fluid Systems, Cambridge Univer-sity Press, 1997, ISBN 052156477

[6] T.A.Davis, UMFPACK - a Column Preordering Strategy for theUnsymmetric-Pattern Multifrontal Method, Univerity of Florida, Depart-ment of Computer and Information Science and Engineering, 2003.

[7] J.W.Demmel et.al, SUPERLU - A Supernodal Approach to Sparse PartialPivoting, University of Califonia, 1999.

[8] J.W.Demmel et.al, SUPERLU - Users Guide, University of Califonia, 1999.

[9] P.Deuflard, Scientific Computing with ODE, Texts in Applied Mathematics,Springer-Verlag NY, YEAR? ISBN?

[10] J.Dongarra, Survey of Sparse Matrix Formats,http://www.netlib.org/linalg/html_templates/node90.html, last vis-ited 040126, 1995.

[11] G.Eriksson, Kompendium i Tillämpade Numeriska Metoder,KTH Depart-ment of Numerical Analysis and Computer Science, 2002.

[12] G.Eriksson, Numeriska Algoritmer med MATLAB,KTH Department of Nu-merical Analysis and Computer Science, 1998.

[13] W.R.Gardner et.al, Post-Irrigation Movement of Soil Water-Simultaneous Redistribution and Evaporation, University of Wisconsin,Madison,Wisconsin 53706, 1970.

59

[14] M.Hanke, Lecture Notes: Advanced Numerical Methods, Royal Instituteof Technology, Department of Numerical Analysis and Computer Science,2003.

[15] M.E.Hosea et.al, Analysis and Implementation of TR-BDF2, Applied Nu-merical Mathematics 20 (1996) 21-37, Elsevier Science B.V. SSDI 0168-9274(95)00115-8, 1996.

[16] W.Hundsdorfer, Partially Implicit BDF2 Blends for Convection DominatedFlows, Siam J.Numerical Analyses, Vol 38 No.6, pp.1763-1783, Society forindustrial and Applied Mathematics, 2001, ISBN 9162838989.

[17] C.Kelly, Iterative Methods for Linear and Nonlinear Equations, Siam Fron-tiers of Numerical Analalysis, 1995.

[18] F.Oeseburg et.al, Evaporation of CW Agents from Various Sub-strates; Wind Tunnel Experiments and Validation of Evaporation Model,FOI/TNO-Prince Mauritz Laboratorium, 1998.

[19] A.Quarteroni and R.Sacco and F.Salieri, Numerical Mathematics, Texts inApplied Mathematics, Springer-Verlag NY, 2000, ISBN 0-387-98959-5.

[20] K.Reichardt and D.R.Nielsen and J.W.Biggar, Scaling of Horizontal Infilt-ration into Homogeneous Soils, Soil Sci. Soc. American Proc. 36, 1972.

[21] Robbert, Sparse Matrix Compression Formats,http://ce.et.tudelft.nl/robbert/sparse-matrix-compression.html, lastvisited 040126.

[22] L.Råde and B.Westergren, Beta, Mathematics Handbook for Science andEngineering , Studentlitteratur, 1998, ISBN 91-44-00839-2.

[23] S.S.Sadhal et.al, Transport Phenomena with Drops and Bubbles, Mechan-ical Engineering Series, Springer Verlag NY, 1997, ISBN 0387946780

[24] L.F.Shampine and M.W. Reichelt, The Matlab ODE Suite, Siam Journalof Scientific Computing Vol 18 No 1 pp 1-22, Society for Industrial andApplied Mathematics, 1997.

[25] A.Sjö, Analysis of computational Algorithms for Linear Multistep Methods,Lunds University, Centre for Mathematical Sciences Numerical Analyses,1999, ISBN 91-628-3898-9.

[26] Weinan.E, Finite Difference Schemes for incompressible Flows in theVelocity-Impuls Density forumulation, Courant Institute of MathematicalScience. NY10012, www.math.princeton.edu/ weinan/papers/cfd8.pdf, lastvisit 031010

[27] S.Winter et.al, On Modeling of the Evaporation of Chemical Warfare Agentson the Ground, Journal of Hazardous Materials A 63 5-24, 1998.

[28] E.G.Youngs and A.Poulovassilis, The Different Forms of Moisture ProfileDevelopment During the Redistribution of Soil Water after Infiltration. ,American Geophysical Union Vol 12 No 5, 1976.

60