Geiser J. Decomposition Methods for Differential Equations: Theory and Applications

Подождите немного. Документ загружается.

Numerical Experiments 157

6.4.2.1 First Test Example: Wave Equation with Constant

Coeﬃcients in Two Dimensions

In the ﬁrst example, we compare the classical splitting method (ADI method)

and the iterative splitting methods for a nonstiﬀ and a stiﬀ case.

We compare the nonstiﬀ case for the operators A and B.

We b egin with D

=1,D

,Δx =

,Δy =

, two iterations per

time-step, and Ω = [0, 1]

, t ∈ [0, 3 · (1/

√

2)].

We discretize with the temporal and spatial derivatives using second-order

ﬁnite diﬀerence methods.

The results are given in Table 6.31 for the classical method (ADI method)

and in Table 6.32 for the iterative operator-splitting method.

TABLE 6.31: Classical operator-splitting

method; D

=1,D

,Δx =

,Δy =

,and

Ω=[0, 1]

, t ∈ [0, 3 · (1/

√

2)].

η err

err

tsteps 44 45 46

0 0.0235 2.3944e-004 4.8803e-007

0.0100 2.0405e-004 5.3070e-007 7.9358e-008

0.0800 4.1689e-005 4.1885e-005 4.2068e-005

0.0833 4.5958e-005 4.5970e-005 4.5982e-005

0.1000 7.0404e-005 6.9235e-005 6.8149e-005

0.2000 3.2520e-004 3.0765e-004 2.9168e-004

TABLE 6.32: Iterative operator-splitting

method; D

=1,D

,Δx =

,Δy =

,two

iterations per time-step, starting solution with an

explicit time-step, and Ω = [0, 1]

t ∈ [0, 3 · (1/

√

2)].

η err

err

tsteps 44 45 46

0 0.0230 2.3977e-004 4.9793e-007

0.0100 1.4329e-004 1.2758e-006 7.9315e-008

0.0800 4.1602e-005 4.1805e-005 4.1994e-005

0.0833 4.5858e-005 4.5879e-005 4.5898e-005

0.1000 7.0228e-005 6.9074e-005 6.8003e-005

0.2000 3.2375e-004 3.0636e-004 2.9052e-004

158Decomposition Methods for Diﬀerential Equations Theory and Applications

In the next computations we compare the stiﬀ case for the operators A and

We deal with D

, D

100

,Δx =

,Δy =

,andΩ=[0, 1]

×[0, 3 ·

(1/

√

2)]. We discretize with the time and space derivatives with second-order

ﬁnite diﬀerence methods.

The results are given in Table 6.33 for the classical method (ADI method)

and in Table 6.34 for the iterative operator-splitting method.

TABLE 6.33: Classical operator-splitting method; D

100

,Δx =

,Δy =

,andΩ=[0, 1]

, t ∈ [0, 3 · (1/

√

2)].

η err

err

tsteps 3 5 6 7 15 100

0.0 41.0072 0.1205 9.2092e-004 0.0122 0.1131 0.1475

0.01 29.0629 0.0626 5.5714e-004 0.0212 0.1171 0.1476

0.05 6.9535 0.0100 0.0491 0.0764 0.1336 0.1479

0.08 0.9550 0.0984 0.1226 0.1331 0.1465 0.1482

0.09 0.1081 0.1409 0.1515 0.1539 0.1508 0.1483

0.1 0.3441 0.1869 0.1819 0.1755 0.1552 0.1484

0.2 0.7180 0.6259 0.4960 0.4052 0.2005 0.1494

0.3 0.2758 0.7949 0.7126 0.6027 0.2479 0.1504

0.4 0.0491 0.7443 0.7948 0.7316 0.2963 0.1513

0.5 0.0124 0.6005 0.7754 0.7905 0.3448 0.1523

TABLE 6.34: Iterative operator-splitting method; D

100

,Δx =

,Δy =

, starting solution with an explicit

time-step, and Ω = [0, 1]

t ∈ [0, 3 · (1/

√

2)].

η err

err

tsteps 3 5 6 7 15 100

#iter 4 8 10 12 28 198

0.0 11.5666 0.0732 6.1405e-004 0.0088 0.0982 0.1446

0.01 9.0126 0.0379 3.6395e-004 0.0151 0.1017 0.1447

0.05 3.3371 0.0047 0.0315 0.0540 0.1161 0.1450

0.08 1.4134 0.0520 0.0785 0.0935 0.1272 0.1453

0.09 0.9923 0.0756 0.0970 0.1081 0.1310 0.1454

0.1 0.6465 0.1019 0.1166 0.1231 0.1348 0.1455

0.2 0.7139 0.4086 0.3352 0.2877 0.1737 0.1465

0.3 0.7919 0.6541 0.5332 0.4466 0.2141 0.1474

0.4 0.6591 0.7963 0.6792 0.5787 0.2552 0.1484

0.5 0.5231 0.8573 0.7734 0.6791 0.2961 0.1493

Numerical Experiments 159

REMARK 6.13 The iterative splitting method is in both cases more

accurate, and the best results are given for the discretization parameter η ≈

0.01. Two iteration steps are suﬃcient to obtain the second-order results. For

the stiﬀ case, more iteration steps are important to obtain the same results

as for the nonstiﬀ case.

6.4.2.2 Second Test Example: Wave Equation with Constant

Coeﬃcients in Three Dimensions

In the second example, we compare the iterative splitting methods for dif-

ferent starting solutions u

i−1

n+1

We consider the explicit method ﬁrst for the starting solution u

i−1,n+1

,with

i−1,n+1

− 2u

+ u

n−1

= τ

∂

∂x

+ D

∂

∂y

In the second method, we calculate u

i−1,n+1

using the classical operator-

splitting method (ADI method), described in Section 6.4.1.

For our model equation, we apply the three-dimensional wave equation:

∂

∂t

= D

∂

∂x

+ D

∂

∂y

+ D

∂

∂z

, (6.66)

where we have exact initial and boundary conditions given in the equation

(6.47).

We use D

=1,D

and D

100

,Δx =Δy =Δz =

,three

iterations per time-step, and Ω = [0, 1]

, t ∈ [0 , 6·(1/

√

3)]. We discretize with

the time and space derivatives with second-order ﬁnite diﬀerence methods.

The results are given in Table 6.35 for the classical method (ADI method)

and in Table 6.36 for the iterative operator-splitting method.

TABLE 6.35: Iterative operator-splitting method;

=1,D

and D

100

,Δx =Δy =Δz =

three iterations per time-step, starting solution with an

explicit time-step, and Ω = [0, 1]

, t ∈ [0 , 6 ·(1/

√

3)].

η err

err

tsteps 8 9 10 11 12

0 169.4361 8.4256 0.5055 0.2807 0.3530

0.2000 0.0875 0.1315 0.1750 0.2232 0.2631

0.3000 0.1151 0.0431 0.0473 0.0745 0.1084

0.4000 0.3501 0.2055 0.0988 0.0438 0.0454

0.4500 0.4308 0.3002 0.1719 0.0844 0.0402

0.5000 0.4758 0.3792 0.2510 0.1402 0.0704

160Decomposition Methods for Diﬀerential Equations Theory and Applications

TABLE 6.36: Iterative operator-splitting method;

=1,D

and D

100

,Δx =Δy =Δz =

three iterations per time-step, starting solution with

classical splitting method (ADI), and Ω = [0, 1]

t ∈ [0 , 5 ·(1/

√

3)].

η err

err

tsteps 6 7 8 9 10

0 79.0827 10.1338 0.2707 0.1981 0.2791

0.2000 0.1692 0.2658 0.3374 0.3861 0.4181

0.3000 0.0311 0.0589 0.1182 0.1983 0.2692

0.4000 0.1503 0.0400 0.0389 0.0793 0.1322

0.4500 0.2340 0.0832 0.0316 0.0502 0.0930

0.5000 0.3101 0.1401 0.0441 0.0349 0.0645

REMARK 6.14 The starting solution prestep combined with the

ADI method is more accurate and gives the best results for the discretization

parameter η ≈ 0.01. In a three-dimensional case, one more iteration step is

needed compared to the two-dimensional case.

Numerical Experiments 161

6.5 Real-Life Applications

In this section we deal with the real-life applications that were introduced

in Chapter 1. We discuss the eﬃciency and accuracy of our proposed splitting

methods for the multi-physics problems.

6.5.1 Waste Disposal: Transport and Reaction

of Radioactive Contaminants

In the next subsections, we describe the two- and three-dimensional simu-

lations of waste disposal.

6.5.1.1 Two-Dimensional Model of Waste Disposal

We calculate some waste disposal scenarios that help us to reach new con-

clusions about the waste disposals in salt domes.

We consider a model based on a salt dome with a layer of overlying rock,

with a permanent source of groundwater ﬂow that becomes contaminated with

the radioactive waste. Based on our model, we calculate the transport and

the reaction of these contaminants coupled with decay chains, see [76]. The

simulation time is 10000[a], and we calculate the concentration of waste in

the water that ﬂows to the top of the overlying rock. With these dates we

can determine if the waste disposal is safe. We present the two-dimensional

test case with the dates of our last project partner GRS in Braunschweig

(Germany), cf. [74] and [75].

We have a model domain with the size of 6000[m] × 150[m] consisting of

four diﬀerent layers with diﬀerent permeabilities, see [74]. Groundwater is

pooled in the domain from the right to the left boundary. The groundwater

ﬂows faster through the permeable layer than through the impermeable layers.

Therefore, the groundwater ﬂows from the right boundary to the middle half

of the domain. It ﬂows through the permeable layer down to the bottom of

the domain and pools at the top left of the domain to an outﬂow at the left

boundary. The ﬂow ﬁeld with the velocity is calculated with the program

package D

F (Distributed-Density-Driven-Flow software toolbox, see [72]),

and is presented in Figure 6.14.

In the middle of the bottom of the domain, the contaminants ﬂow as

from a permanent source. With the stationary velocity ﬁeld, the contami-

nants are computed with the software package R

T (Radioactive-Reaction-

Retardation-Transport software toolbox, see [73]). The ﬂow ﬁeld transports

the radioactive contaminants up to the top of the domain. The decay chain

162Decomposition Methods for Diﬀerential Equations Theory and Applications

FIGURE 6.14: Flow ﬁeld for a two-dimensional calculation.

is presented with 26 components as follows,

Pu−244 → Pu−240 → U −236 → Th−232 → Ra−228

Cm−244 → Pu−240

U −232

Pu−241 → Am−241 → Np−237 → U −233 → Th−229

Cm−246 → Pu−242 → U −238 → U −234 → Th−230 →

Ra−226 → Pb−210

Am−242 → Pu−238 → U −234

Am−243 → Pu−239 → U −235 → Pa−231 → Ac−227 .

e present the important concentration in this decay chain. In Figure

6.15 the contaminant uranium isotope U-236 is presented after 100[a]. This

movement of isotope is less retarded by adsorption than other nuclides and has

a very long half-life period. The concentration of U-236 in the model domain

is illustrated in Figure 6.15 at two diﬀerent time points in the modeling.

The diﬀusion process spreads the contaminant throughout the left side of the

domain. The impermeable layer is also contaminated. After the time period

of 10000[a], the contaminant is ﬂown up to the top of the domain.

The calculations are performed on uniform grids. The convergence of these

grids is conﬁrmed with adaptive grid calculations, which conﬁrmed the results

of ﬁner and smaller time-steps, cf. Table 6.37. The calculations begin with

explicit methods until the character of the equation is more diﬀusive, at this

point we switch to the implicit methods and use large time-steps. With this

procedure, we can fulﬁll the mandatory maximum calculation time of 1 day.

6.5.1.2 Three-Dimensional Model of a Waste Disposal

In this example, we consider a three-dimensional model, because of the in-

terest in the three-dimensional eﬀects of contaminants in the ﬂowing ground-

Numerical Experiments 163

FIGURE 6.15: Concentration of U-236 at the time point t = 100[a]and

t = 10000[a].

water. We simulate approximately 10000[a] and focus on the important

contaminants ﬂowing farthest with a high concentration. We assume an

anisotropy domain of 6000[m] × 2000[m] × 1500[m] with diﬀerent permeable

layers. We have calculated 26 components as presented in the two-dimensional

case. The parameters for the diﬀusion and dispersion tensor are given as

D =1· 10

−9

/s], α

=4.0[m], α

=0.4[m], |v|

max

=6· 10

−6

[m/s],

ρ =2· 10

, D

= α

|v| and D

= α

|v|, where the longitudinal disper-

sion length is ten times bigger than the transversal dispersion. The source is

suited at the point (4250.0, 2000.0, 1040.0), and the contaminants ﬂow with

a constant rate. The underlying velocity ﬁelds are calculated with D

F,and

we added two sinks at the surface with the coordinates (2000, 2100, 2073) and

(2500, 2000, 2073).

We simulate the transport and the reaction of the contaminants with our soft-

ware package R

T. As a result, we can simulate the needed pumping rate of

the sinks to pump out all of the contaminated groundwater. We present the

velocity ﬁeld in Figure 6.16. The groundwater ﬂows from the right boundary

to the middle of the domain. Due to the impermeable layers, the groundwa-

164Decomposition Methods for Diﬀerential Equations Theory and Applications

TABLE 6.37: Computing the two-dimensional case.

Processors Reﬁnement Number of Number of Time for total

elements time-steps one time-step time

30 uniform 75000 3800 5 sec. 5.5 h.

64 adaptive 350000 3800 14 sec. 14.5 h.

ter ﬂows downward and pools in the middle part of the domain. Because of

the inﬂuence of the salt dome, salt wells up with the groundwater and be-

comes incorporated into the lower middle part of the domain curls. These

parts are interesting for three-dimensional calculations, and because of the

curls, the groundwater pools. We concentrate on the important component

FIGURE 6.16: Flow ﬁeld for a three-dimensional calculation.

U −236. This component is less retarded and is drawn into the sinks. In the

top images of Figures 6.17, we show the initial concentration at time point

t = 100 [a]. We present the data in vertical cut planes. In the next picture,

we present a cut plane through the source term. In the bottom of Figures

6.17, the concentration is presented at the end time point t = 10000 [a]. The

concentration has welled from the bottom up over the impermeable layer and

into the sinks at the top of the domain.

At the beginning of the calculation, we use explicit discretization methods

with respect to the convection-dominant case. After the initializing process,

the contaminants are spread out with the diﬀusion process. We use the im-

plicit methods with larger time-steps and could also calculate the mandatory

time period with a higher-order discretization method.

The computations are performed on a Linux PC-Cluster with 1.6GHz

Athlon processors.

Numerical Experiments 165

FIGURE 6.17: Concentration of U-236 at the time points t = 100[a]and

t = 10000[a].

166Decomposition Methods for Diﬀerential Equations Theory and Applications

In Table 6.38, we report the results of the computations. We begin with

convergence results on uniform reﬁned meshes. We conﬁrm these results with

adaptive reﬁned meshes and obtain the same results with smaller time-steps.

We accomplished these calculations within the 1 day mandatory computa-

tional time limit.

TABLE 6.38: Three-dimensional calculations of a realistic potential

damage event.

Processors Reﬁnement Number of Number of Time for Total

elements time-steps one time-step time

16 uniform 531264 3600 13.0 sec. 13.0 h

72 adaptive 580000 3600 18.5 sec. 18.5 h

REMARK 6.15 Here the performance is obvious, based on the de-

coupling of the equation into convection-reaction and diﬀusion parts. We

could accelerate the solver process with the most adapted discretization and

solver methods. The resources are saved with the direct ODE solvers for the

reaction part and the possibility of large time-steps for the implicit discretized

diﬀusion part. For such problems, a decoupling makes sense with respect to

parallelizing the underlying solver methods.

6.5.2 Sublimation of Crystal Growth

In this next section, we focus on decoupling the complicated sublimation

process in simpler processes. We discretize and solve the simpler equations

by more accurate methods with embedded analytical solutions, see [83] and

[84].

With the test examples, we verify the decomposition methods for the convection-

reaction and heat equation. The models are based on the conduction of heat

in solids, see [38] and a simulation of gas mixture, see [158]. Based on these

results, we present a real-life application. We discuss the heat equation in dif-

ferent material layers and sketch out a treatment to increase computational

eﬃciency based on decomposition methods.

6.5.2.1 Simpliﬁed Models for the Crystal Growth Apparatus:

Convection-Reaction Equation for the Simulation of the

Crystal Growth

We apply the operator-splitting methods for the convection-reaction equa-

tion; see further numerical results in [88] and [89].