4.2 Heat and fluid flow at the 3-phase contact line
4.2.1 Solution and parameterization of the microscale model
pl−pv=−ρlul ul−uint+ρvuv uv−uint
(4.22) This equation corresponds exactly to Eq. (2.16) which has been analytically derived for the normal momentum balance at a moving interface with phase change. Only the pressure difference due to the surface tension forces is not reproduced as the interface of the example is not curved. Hence, it could be shown that there is no additional source term needed for the momentum equation to account for phase change. The recoil pressure is exactly reproduced by balancing the momentum for a control volume containing the liquid-vapor interface and the regions with the source terms. It should be noted that the above integration of the mass and momentum conservation equations is performed in a region around the liquid-vapor interface and not directly at the interface itself. Hence, the recoil pressure is correctly predicted in some distance (several cells) to the interface. Certainly, it is not correctly predicted directly at the liquid-vapor interface.
its thickness can be obtained by setting Eq. (4.27) to zero and by taking into account that the curvature of the adsorbed film is zero. The initial conditions for the remaining variables can then be obtained, too.
δad =
A
ρl∆hV Twall
Tsat −1
1/3
(4.28)
δ0ad = 0 (4.29)
∆pad = A
δ3ad (4.30)
Qad = 0 (4.31)
This system of equations is solved using the softwareMatlab2. A 4th orderRunge-Kuttamethod3with error estimation and step size adaption is used for the integration. As the conditions within the adsorbed film are the trivial solution of the system of ODEs (4.24) to (4.27), they cannot be directly used for the numerical integration. Instead, the film thickness δ and the integrated heat flux Q are slightly perturbed4. The perturbation on the integrated heat flux Q is used as parameter within the shooting method in order to achieve the target value for a macroscopic curvature at the end of the micro region.
A detailed description of the solution procedure is given by Stephan [104].
The system of ODEs (4.24) to (4.27) is integrated from the adsorbed film at ξ=0to a chosen value ξ=ξmic. In principal, the choice ofξmicis arbitrary. However, two aspects must be taken into account.
First, the film thickness at ξmic must be large enough such that the microscale effects (i.e. disjoining pressure and change of thermodynamic equilibrium by extremely high curvature) have decayed. This criterion provides a lower limit for the choice of ξmic. Second, the film thickness and its slope at ξmic
must be small enough such that the assumptions which are made during the derivation of the equations for the contact line heat transfer in section 2.3 are still valid. In particular, the film must still be thin and flat enough to justify the assumptions of 1D heat conduction and of a flow which is nearly parallel to the wall. In the calculations performed for the present thesis, a value ofξmic=0.5µmwas chosen.
A typical result of the contact line evaporation model is shown in Figure 4.4. The film thickness and slope are shown in Figure 4.4(a). Starting at the very small value ofδad which is not visible in the plot, the film thickness starts to grow rapidly at aroundξ=0.1µm. This goes along with a step-like increase of the slopeδ0and thereby of the angleθ =arctanδ0which approaches the apparent contact angle. The sharp rise in contact angle is a result of the rapidly decaying intermolecular forces which are taken into account by the disjoining pressure concept. The corresponding pressures are shown in Figure 4.4(b). The total pressure difference between liquid and vapor is the sum of disjoining pressureA/δ3 and capillary pressureσκ. Although recoil pressure is not taken into account in the modeling equations, it is calculated from the local heat flux (q2/∆h2V 1/ρv−1/ρl
, this equation is valid for a non-moving interface) and also plotted in Figure 4.4(b). It can clearly be seen that the disjoining pressure is highly dominant within the adsorbed film. When the film thickness starts to grow, the disjoining pressure decreases rapidely and is replaced by the capillary pressure. In the presented case, the recoil pressure is very low compared to disjoining and capillary pressure. It was checked that this is also valid for higher wall superheats and different material properties. Thus, under the conditions which are subject of the
2 Technical computing software distributed byThe MathWorks.
3 Theode45function of theMatlabsoftware package is used.
4 The magnitude of the perturbation of the film thickness is chosen as∆δ=δad/1000.
4.2 Heat and fluid flow at the 3-phase contact line 41
0 0.1 0.2 0.3 0.4 0.5 0
0.1 0.2
coordinate,ξ [µm]
filmthickness,δ[µm]
0 0.1 0.2 0.3 0.4 0.50
10 20
slope,arctanδ’[°]
film thickness slope
(a) Film thicknessδand slopearctanδ0
0 0.1 0.2 0.3 0.4 0.5
100 102 104 106 108
coordinate,ξ[µm]
pressure[Pa]
total pressure difference disjoining pressure capillary pressure recoil pressure
(b)Pressure comparison (logarithmic scale)
0 0.1 0.2 0.3 0.4 0.5
0 10 20 30 40 50 60 70 80
coordinate,ξ[µm]
heatflux,q[MW/m2 ]
0 0.1 0.2 0.3 0.4 0.50
2 4 6 8 10 12 14 16
integratedheatflux,Q[W/m]
heat flux
integrated heat flux
(c)Heat fluxqand integrated heat fluxQ
0 0.1 0.2 0.3 0.4 0.5
373 374 375 376 377 378 379
coordinate,ξ[µm]
temperature[K]
saturation temperature of bulk liquid
local saturation temperature
(d)Saturation temperature
Figure 4.4:Example of the results obtained from the microscale model for the evaporation at the 3-phase contact line. The fluid is saturated water (p=1013 mbar) at a wall superheat of∆Twall=5 K.
present thesis, the effect of the recoil pressure can be neglected in the microscopic region at the 3-phase contact line. However, this statement is not necessarily valid in cases with much higher wall superheat, e.g. in a situation close to the critical heat flux. The heat transfer at the contact line is shown in Figure 4.4(c). There is no evaporation from the adsorbed film and therefore the heat flux is zero in this region. At the point where the film thickness starts to grow, the heat flux increases and quickly reaches a maximum value. The position of the maximum heat flux represents the point where the film thickness is still very small (which enhances heat transfer) while the intermolecular forces (which reduce heat transfer) have already become very small. After the maximum, the heat flux continuously decreases due to the growing thickness of the film. The effect of the intermolecular forces can be better understood when looking at Figure 4.4(d) in which the bulk saturation temperature and the local saturation temperature are compared. The bulk saturation temperature of water atp=1013 mbar is Tsat =373.15 K (100◦C). However, the intermolecular forces and the strong capillary pressure lead to a change of the local thermodynamic equilibrium in the vicinity of the wall. The local saturation temperature is the minimum temperature at the interface which is required for the evaporation of liquid.
Its value is equal to the wall temperature within the adsorbed film. Thus, the intermolecular forces preemt evaporation at the adsorbed film. As the film thickness grows, the local saturation temperature drops down to the bulk saturation temperature due to the decay of the microscale effects.
4.2 Heat and fluid flow at the 3-phase contact line 42
0 5 10 15 20 0
50 100
integratedheatflux,Q mic[W/m]
0 5 10 15 200
0.2 0.4
wall superheat, ΔT
wall [K]
filmthickness,δ mic[µm]
integrated heat flux (data points) integrated heat flux (correlation) film thickness (data points) film thickness (correlation)
Figure 4.5:Integrated heat flux and film thickness at the end of the integration domain of the contact line evaporation model for different wall superheats. Results obtained from contact line evap-oration model (data points) and correlations (according to Eq. (4.32) and Eq. (4.33)). The fluid is saturated water (p=1013 mbar).
From the results shown in Figure 4.4(a) to Figure 4.4(d) one can see that the choice ofξmic=0.5µm is appropriate in this case. At the end of the integration of the system of ODEs (4.24) to (4.27) the microscale effects have decayed completely while the thickness and slope of the liquid film are still small enough to justify the assumptions of the model.
The solution algorithm for the integration of the system of ODEs is too time consuming to be used directly in the CFD simulation. The calculation of one solution for a particular fluid and a particular wall superheat only takes several seconds5. However, it would have to be performed at each segment of the contact line during one time step of the CFD simulation (see section 4.2.2 for details) leading to a large number of solutions which must be calculated for the contact line evaporation model. Therefore, the results of the contact line evaporation model are parameterized prior to the CFD simulation. For a given set of material properties, the results of the contact line evaporation model depend mainly on the wall superheat. There is also some influence of the target curvature that is achieved at the end of the integration domain of the shooting method. However, the influence is far inferior to the influence of the wall superheat. Hence, the contact line evaporation model is solved for different wall superheats and a correlation is determined. The results which are required for the coupling to the CFD simulation are the integrated heat flux, the film thickness and prospectively the contact angle at the end of the integration domain (see section 4.2.2 for details). These results are correlated to the wall superheat. A series of root functions was found to fit the data best.
5 These computations were performed on a single core of an Intel Core 2 Duo with 3 GHz usingMatlab7.6.0.324 (R2008a).
4.2 Heat and fluid flow at the 3-phase contact line 43
wall
liquid vapor
reconstructed interface
θ
ξmic
ξmesh
transition region
film thickness heat flux θ
Figure 4.6:Illustration of the coupling between the subgrid scale model for the contact line evaporation and the CFD simulation. A transition region is defined between the different length scales of the subgrid scale model and the CFD simulation.
Qmic=Q ξmic
=
Ni
X
i=1
aQ,i∆Twall1/i (4.32)
δmic=δ ξmic
=
Ni
X
i=1
aδ,i∆Twall1/i (4.33)
θmic=arctanδ0 ξmic
=
Ni
X
i=1
aθ,i∆Twall1/i (4.34)
Typical results of the parameterization procedure are shown in Figure 4.5. It can easily be seen that the correlation fits the data points very accurately. The coefficientsaQ,i, aδ,i and aθ,i are calculated for different fluids. The values of the coefficients are given in appendix C.