Whereas the results obtained in [34] were in general quite encouraging, the long term goal of simulating the complete treatment process from when the rods are inserted and until all boron has been washed out, is still not within immediate reach. Apart from the uncertainties regarding the transport of water, which as demonstrated in the previous chapters are significant, the most pressing issue concerns boron diffusion coefficients.
6.5 Conclusions and future work Towards modeling preservation processes
In an experimental study undertaking the determination of these it would be very de-sirable to first determine the diffusion coefficient of boron compound in pure water and then express the boron coefficients in wood as a fraction of this. Not only would that give a direct measure of the resistance offered by wood, but would also give an idea of the continuity of the water phase at low degrees of saturation. Ultimately, it could be hoped that a correlation between boron diffusion coefficients and free water transfer coefficients could be established.
Also, it would be of interest to examine the influence of convective transfer, which would not involve any additional practical difficulties – in the first type of experiment the sam-ples are equilibrated to some moisture content and then sealed, whereas the second type of experiment involves also a transfer of water.
Finally, it is necessary to formulate appropriate boundary conditions and then determine the associated material parameters. As for the first part of this task, the general con-vection type boundary condition which is well–established in the transfer of both heat and moisture, seems promising. Secondly, different physical conditions at the boundary, ranging from pure water to moist soil, should be considered. At the time of writing, a first step in this regard has been taken in an experimental study at the University of Toronto, where the leaching in samples submerged in fresh water is being sought quantified [11].
Department of Civil Engineering - Technical University of Denmark 53
Towards modeling preservation processes 6.5 Conclusions and future work
Chapter 7
Numerical solution procedures
In wood science and in fact, materials science in general, efficient numerical methods are important for at least three reasons.
Firstly, in the development of new models such as those presented in this thesis one is always faced with the problem of supplying material parameters which make the model fit the experimental results as well as possible. In the present this has been done by simple trial and error, and of course with a good feeling for the neighbourhood in which the optimal parameters are to be found. When this is done it is desirable to have an efficient and robust numerical procedure.
The second reason is along the same lines. Because, instead of determining the parameters as sketched in the above, an alternative would be to use some kind of inverse modeling.
Exactly how this should be carried out is uncertain, but will inevitably involve solving the same equations thousands or tens of thousands of times. Thus, a factor of 50 may mean the difference between a week and a year.
Lastly, it must be expected that in the not so far future, inclusion of the statistical varia-tion of the model parameters will become an integral part of any computavaria-tion involving wood. Such computations will again require thousands of runs, each with slightly different parameters, and thus, a factor of 10–100 may mean the difference between making such computations practically feasible or not.
The models presented in previous chapters are all formulated in terms of a set of partial differential equations. In only a few cases can these be solved analytically, or otherwise, by hand calculation methods. Thus, in general some kind of numerical procedure must be applied. For the types of equations dealt with in this thesis, i.e. diffusion type equa-tions, finite difference methods or the more modern equivalents such as the control volume method, are often used, see e.g. [48, 44, 43, 49, 63].
In this work, however, we have exclusively applied the finite element method. Thus, it would not be appropriate to go into a deeper discussion of the relative merits of each method. We will, however, remark that for the type of problems dealt with here, there are no inherent problems with the finite element method, provided that certain trivial provision are taken. These include the use of low order elements, regular grids, and lumped mass matrices. Indeed, if these criteria are met, the two methods are in many cases equivalent. Detailed discussions of the two methods are given in [9, 26, 10].
55
Numerical solution procedures 7.1 Classification of diffusion–type equations
In recent years, coupled heat and mass transfer problems have become quite common, and the primary focus has often been on the numerical solution method rather than on the physics of the problem. Applications include wood drying [49, 48, 63], non–isothermal consolidation of soils [20, 22] and concrete exposed to high temperatures as during fire [21]. However, in our opinion, the numerical methods applied are often far from optimal.
In [20, 22] for example, both eight and nine node plane elements are used, although these are known to cause possibly very severe non–physical oscillations in the computed results.
Moreover, it seems that consistent mass matrices are used, and no particular mention is made of the convergence of the Newton iterations. Finally, the choice of primary state variables is inappropriate and is known to result in non–conservative solutions.
Concurrent with the development in the solution of coupled heat and mass transfer prob-lems, there has been a significant effort, primarily within the soil sciences, to develop efficient numerical methods for the solution of the simpler problem of isothermal single–
phase free water flow, see e.g. [19, 18, 38, 61, 15, 31, 32, 30, 50, 71, 33, 17, 56].
In most coupled heat and mass transfer problems it is the equation governing the flow of water which is by far the most nonlinear and, thus causes the most problems. In the fol-lowing, therefore, rather than solving equations whose validity is questionable by methods whose efficiency is even more questionable, we will concentrate on developing an efficient method for the problem of isothermal flow of water.
The method presented has already been applied successfully to the related problem of heat conduction with phase change [35]. In the following we will demonstrate its applica-bility, first to the problem of free water flow in soils governed by the so–called Richards equation, and then to the problem of water transport in wood where, in addition to free water, also bound water and water vapour is included.
Secondly, some of the more practical problems are discussed. These relate particularly to the appropriate weighting of diffusive and convective terms. Thus, it is well–known that in highly nonlinear diffusion–type problems there is the possibility of encountering solutions with spurious non–physical oscillations. These can, however, be countered by appropriate weighting of the material parameters over each element.
Finally, the issue of computing the tangent matrices required in Newton’s method is treated. Most often these matrices are determined numerically. However, we present a general set of formulas enabling exact analytical determination of the matrices. The use of these will generally be more cost efficient and further, it will reduce or eliminate the possibility of round–off and other numerical errors often encountered when computing tangent matrices numerically.
7.1 Classification of diffusion–type equations
In many cases diffusion–type equations describing the transport of some conservative substance, e.g. heat or mass, are naturally formulated as
∂θ
∂t =∇ ·(λ∇φ) (7.1)
where the left–hand side accounts for the accumulation of the conserved substance per unit volume per unit time. Thus, the variable θcould be enthalpy per unit volume or water
7.1 Classification of diffusion–type equations Numerical solution procedures
-10 -5 0 5 10
0 100 200 300 400
Temperature (oC)
Enthalpy (kJ/kg)
Soild Fluid
(a)
-300 -20 -10 0 10 20
0.2 0.4 0.6 0.8 1
Pressure head (m)
Saturation
Unsaturated Saturated
Clay
Sand
(b)
Figure 7.1Enthalpy–temperature curve for water (a) and pressure–saturation curves for sand and clay (b).
per unit volume. The right hand side describes the diffusive transfer of the substance by proportionality to a gradient in some driving force φ and with the proportionality factor(s) being described throughλ. Thus, for the heat conduction problemφwould be temperature, and for the problem of unsaturated flow in porous media, capillary pressure.
Depending on the particular problem in question additional terms may also be included, e.g. those describing sources and sinks, gravity, etc. The most important point is, however, the appearance of two variables in the governing equation. Thus, in order to solve the equation, a relation betweenθandφmust be established. And this is where the primary source of numerical difficulties originates. Because in many cases, the θ–φ relation is extremely nonlinear. Two examples of this are shown in Figure 7.1. Figure 7.1 (a) shows the enthalpy (θ) of water as function of temperature (φ), and in 7.1 (b) the saturation (θ) of two porous materials are shown as function of the pressure head (φ). In the porous media flow problem also the conductivityλis a strong function of the saturation, but as will be shown, the most critical issue concerns theθ–φrelation.
The governing equation (7.1) supplemented with the appropriateθ–φ relation gives rise to the possibility of three different formulations, which are equivalent in the continuous versions, but which posses very different properties when it comes to their numerical implementation.
7.1.1 φ–form
Theφform of the governing equation can be written as C∂φ
∂t =∇ ·(λ∇φ) (7.2)
where the ‘capacitance’ is given by
C= ∂θ
∂φ (7.3)
Department of Civil Engineering - Technical University of Denmark 57
Numerical solution procedures 7.1 Classification of diffusion–type equations
The discrete finite element formulation with the backward Euler scheme applied to the time derivative ofφ then reads
Cn+1
φ
φφn+1−φφφn
∆t +Kn+1φφφn+1+fn+1=0 (7.4) where the capacitance and conductivity matrices are given by
C= Z
Ω
NTCNdΩ, K= Z
Ω
BTλBdΩ (7.5)
andf defines the boundary conditions which will be dealt with in detail later on.
Of the three possible formulations theφ–form is perhaps the most inefficient. This is due to the representation of the capacitance. Thus, if theθ–φ relation is highly nonlinear as shown e.g. in 7.1 (a) the most important part of it, i.e. the part with the highest slope, may be missed altogether. Thus, the discreteφ–form is generally non–conservative although, of course, in the limit of infinitely small time steps, the method is conservative. In the literature concerned with phase–change problems, i.e. with an enthalpy–temperature relation as shown in Figure 7.1 (a), theφ–form has traditionally been rather popular. In addition to using very small time steps, some approximation to the capacitance is usually applied. Thus, Del Guidice et al. [27] suggest replacing C by a spatially smoothed equivalent given by
C ≈ ∇θ· ∇φ
∇φ· ∇φ (7.6)
whereas Morgan et al. [41] use a temporal average given by Cn+1≈ θn−θn−1
φn−φn−1
(7.7) where n+ 1 refers to the new state,n to the current and n−1 to the previous state.
Other averaging schemes are described in [39]. Although these schemes do improve the enthalpy balance, the time step is still severely restricted and furthermore, non–physical temperature oscillations may occur.
In the case of porous media flow, the φ–form has also been rather widely applied, see e.g. [38] for a survey of commonly applied methods. Here, similar approximations to C as those used in the phase change problem are often applied. This again improves mass balance, but again the magnitude of the time step is rather restricted. However, for porous media flow problems, theφ–form has the distinct advantage that φ, i.e. the water pressure or pressure head, is a unique variable for both partially and fully saturated conditions, such that in this case there seems to be some justification to using theφ–form, at least over theθ–form, which will be discussed subsequently.
7.1.2 θ–form
Theθ–form of the governing equation is given by
∂θ
∂t =∇ ·(D∇θ) (7.8)