6 Numerical Applications
6.1 The Heat-Pipe Problem
qenergy pl= 1 bar
T = 70◦C
s
Figure 6.1: The spatial domain of the heat-pipe problem: The experiment consists of a tube of2 mlength filled by a fully water-saturated porous medium. On the left side, the tube is open and exposed to atmospheric pressure; on the right, it is closed but heated. Initially the temperature within the domain isT = 70◦C.
Figure 6.2: Uniform extrusion of a one-dimensional grid to three spatial dimensions.
also only be transported by means of heat conduction. The most interesting part of this setup lies in between these two extremes: There, liquid water as well as steam is present. The liquid is drawn from left to right by the capillary pressure where the liquid is evaporated. Then, the material gets transported back to the left by a pressure gradient in the gas phase and it condenses to liquid again.
6.1.1 Uniform Domain Extrusion
Since the mass and energy conservation equations (2.10, 2.23) are only valid for three spatial dimensions, we first have to become aware of their interpretation in the case of one space-dimension. The answer turns out to be quite simple: In order to apply Equations 2.10 and 2.46, the spatial domain can be extruded uniformly to three dimensions as depicted in Figure 6.2 and the usual conservation equations thus apply.
6.1.2 Semi-Analytical Solution
Besides this small glitch in the interpretation of the results, the quasi one-dimensional nature of this setup allows us to derive a semi-analytical solution. To obtain it, we first suppose the system to be in steady-state. This means that there cannot be any net flux of mass along the axissof the tube,i.e.,
ρlkr,l µl
K∂ pl
∂s =−ρgkr,g µg
K∂ pg
∂s (6.1)
The Heat-Pipe Problem 79
holds.
Since the intrinsic permeabilityK is equal on both sides, it can be eliminated. If we also include the definition of capillary pressurepg =pl+pc,gl, and if we use kinematic viscosi-tiesνα=µα/ραinstead of dynamic ones, we can transform Equation 6.1 to
kr,l νl
∂ pl
∂s =−kr,g νg
∂(pl+pc,gl)
∂s .
Reordering, we get
∂ pl
∂s =−
kr,g
νg
∂ pc,gl
∂s kr,l
νl +kνr,g
g
. (6.2)
We can now multiply the numerator as well as the denominator of Equation 6.2 withνg/kr,g, and apply the chain rule. We get
∂pl
∂Sl
∂Sl
∂s =−
∂ pc,gl
∂Sl
∂Sl
∂s
1 +kkr,g
r,l
νl
νg
. (6.3)
Finally, multiplying both sides of Equation 6.3 with∂s/∂Sl, we get
∂pl
∂Sl =−
∂pc,gl
∂Sl
1 +kkr,g
r,l
νl νg
. (6.4)
For the next part of the semi-analytical solution we take advantage of the fact that the energy flux is a given constant for the whole length of the tube. This means that
hlρlvl+hgρgvg+λpm∂T
∂s =qenergy (6.5)
holds, whereλpmis the overall heat conductivity of the medium. In this context, we also know that
ρlvl=−ρgvg
holds, since the net mass flux must be zero. This allows us to transform Equation 6.5 to ρgvg(hg−hl) +λpm∂T
∂s =qenergy .
Further, using DARCY’s relation (2.39) for the velocity of the liquid phase, and using kinematic instead of dynamic viscosities, we get
hvap
kr,g νg
K∂pg
∂s +λpm
∂T
∂s =qenergy ,
wherehvap=hg−hlis thespecific enthalpy of vaporizationof liquid water. Next, we use the
chain rule and the definition of capillary pressure, so we get hvapkr,g
νg K pl
∂s +∂pc,gl
∂Sl
∂Sl
∂Sl
∂Sl
∂s
+λpm∂T
∂s =qenergy . Multiplying both sides with∂s/∂Slleads us to
hvap
kr,g
νg
K ∂pl
∂Sl +∂pc,gl
∂Sl
+λpm
∂T
∂Sl =qenergy
∂s
∂Sl . (6.6)
In Equation 6.6 we now observe that, since we expressedpg in terms of pl and pc,gl, the term∂pl/∂Slis identical to zero, so that we get
hvapkr,g
νg K∂pc,gl
∂Sl +λpm∂T
∂Sl =qenergy ∂s
∂Sl . Now, we apply the chain rule to the temperature term. We get
hvapkr,g νg
K∂pc,gl
∂Sl
+λpm∂T
∂pg
∂(pl+pc,gl)
∂Sl
=qenergy ∂s
∂Sl
.
After some reordering, we get
∂s
∂Sl = ∂pc,gl
∂Sl
hvapKkνr,g
g +λpm∂p∂T
g
qenergy
.
Finally, we take advantage of
pg =pvap(T),
i.e., the fact that the pressure of the gas phase is identical to the vapor pressure of water [88]
at a given temperature in the two-phase, one-component case. This leads us to the second equation of the semi-analytical solution
∂s
∂Sl = ∂pc,gl
∂Sl
hvapKkνr,g
g +λpm
∂p
vap
∂T
−1
qenergy
. (6.7)
Together, Equations 6.4 and 6.7 form a coupled system of two ordinary differential equations which can be solved using the methods presented in Section 4.2.
6.1.3 Results
For the results presented in this section, the capillary pressure relation proposed by LEV
-ERETT[51]
pc,gl=pc,gl,0γ(1.263Sg3−2.120Sg2+ 1.417Sg),
The Heat-Pipe Problem 81
s/cm]
Sl
100 75
50 25
0 1
0.75
0.5
0.25
0
Analytic NCPPVS 50 elements 100 elements 200 elements 400 elements 800 elements 1600 elements
(a) Liquid Saturation (b) Temperature
s/cm T/K
100 75
50 25
0 377
376
375
374
373
Figure 6.3: Comparison of the semi-analytical solution calculated using Equations 6.4 and 6.7 (dotted) with results obtained by the numerical models using primary variable switching (PVS, dashed) and linear complementary functions (NCP, non-dashed). The numerical models were run untilt= 109swith a uniform spatial resolution along the tube axis of50,100 200,400,800, and1600elements.
was applied, whereγ = 0.05878N/mrepresents the surface tension, andpc,gl,0=p
φ/Kis the scaling pressure. Further, the porosityφis0.4, and the intrinsic permeabilityK is10−12m2. The overall heat conductivityλpmis approximated by the SOMERTON[76] relation
λpm=λpm,dry+p
Sl(λpm,wet−λpm,dry) ,
where λpm,dry is the heat conductivity of the porous medium fully saturated by steam and λpm,wet is the heat conductivity of the porous medium when it is fully saturated by liquid water. For these parameters, we assume
λpm,dry = 1.85W/K m2 and λpm,wet = 2.67W/K m2 .
Finally, the energy influx at the right side of the domain isqenergy = 200W/m2.
Figure 6.3 shows the obtained results for the steady-state saturation of liquid water and the temperature for the semi-analytical solution calculated using Equations 6.4 and 6.7 in comparison to the results of the finite-volume based numerical models. In this figure, the origin of the spatial domain is set to the leftmost point where the saturation of the gas phase deviates from zero.
Generally, it can be observed that both numerical solutions are identical up to rounding errors, and that both converge to the semi-analytical solution when the grid gets refined.
From these results, we can conclude that both numerical models and their implementation are correct in the sense that the discretization is convergent and that both models lead to
NCP
(200 elements)
NCP
(1600 elements) PVS
(200 elements)
PVS
(1600 elements)
n[−] 48 232 344 1629
tCPU[s] 1.056 79.99 12.212 469.603
nNEWTON[−] 234 2161 3858 17750
nNEWTON/tCPU[1/s] 221.59 27.01 315.91 37.79
Table 6.1: Comparison of the performance of the numerical models using primary variable switching (PVS) and non-linear complementarity functions (NCP) for the heat-pipe problem. Here,ncorresponds to the number of time steps required to reach steady state, tCPUis the CPU time required to reach steady state andnNEWTONis the total number of iterations of the NEWTON-RAPHSONmethod required for the whole simulation.
Figure 6.4: The spatial domain of the radially symmetric CO2injection problem.
physically meaningful results.
Table 6.1 compares the number of time steps, NEWTON-RAPHSONiterations, and the CPU time which both numerical models required to reach the simulation time of109son a com-puter equipped with an IntelR CoreTMi7-3770K processor at3.5 GHz. These results clearly demonstrate that the model which directly incorporates the non-linear complementarity functions exhibits much better stability,i.e., the number of time steps required is about an order of magnitude lower compared to the model which uses primary variable switching.
This effect more than amortizes the slightly higher computational cost required for each iteration of the NEWTON-RAPHSONscheme.