Model Predictive Control Of Timber Drying

Published at IFAC World Congress, June 30-July 4, San Francisco 1996, Vol. N, 355-360

Model Predictive Control of Timber Drying H. E. Musch, G. Barton, T. Langrish, A. Brooke Department of Chemical Engineering The University of Sydney NSW 2006, Australia

Abstract There is considerable economic incentive for improved timber drying in terms of reduced total drying time, lower wastage levels and improved product quality. This paper presents simulation results in two areas pertinent to improved timber drying productivity. Firstly, it is shown that timber drying can be posed and solved as a nonlinear model predictive control problem. Recognition of two distinct phases (drying to a set average moisture content, followed by flattening of the moisture content profile) reduces the amount of computation required. It is also shown that using only the initial moisture content and an on-line measurement of a timber stack’s weight is sufficient to design a closed-loop control system whose performance is essentially indistinguishable from that achieved using an open-loop system based on perfectly known drying and internal stress models. In both open and closed-loop cases, significant improvements in both total drying time and final product quality were obtained.

1 Introduction The mechanical properties of timber are strongly dependent on its moisture content. When its moisture content is reduced, timber’s strength and durability are increased but shrinkage will start to occur as soon as the moisture content falls below the fibre saturation point. Therefore, the moisture content of timber must be brought into equilibrium with the environment before it can be used as a construction material. This is frequently achieved by cutting the timber into boards which are then stacked and dried using ambient air as the drying medium. However, drying under such conditions is very slow and drying times of up to two years may be necessary. This approach of timber drying necessitates the maintenance of a large (and expensive) inventory. Further, the final moisture can vary from 6% (0.06 kg water/kg dry timber) in hot/ dry ambient conditions to 40% in cool/wet drying conditions. These disadvantages (very slow drying and a highly variable final moisture content) can be overcome by kiln drying where the timber is dried under controlled air temperature and 1

humidity conditions. However, since the drying conditions here are more severe due to increased moisture transport from the timber surface, moisture gradients in the timber can become quite substantial. These gradients cause considerable stress within the timber boards, and degradation by phenomena such as checking (both surface and internal), and collapse can occur, as shown in Figure 1.

Surface checking

Internal checking

Collapse

Figure 1: Timber degradation phenomena

Degradation of timber results in substantial wastage and loss of economic value and, therefore, must be reduced as much as possible. As a result, standard drying schedules have been developed for many species of timber, as each type responds differently under kiln drying. Such a drying schedule is a pair of setpoint profiles (as a function of drying time, or timber moisture content) for the dry and wet bulb temperatures in the kiln. An example of an Australian standard schedule for XXXX hardwood is given in Table 1 [1]. These schedules obviously try to maximize the drying rate at any time. However, they must be conservative to Table 1: A standard drying schedule TYPE??? Moisture change point (% dry basis)

Dry bulb temperature (°C)

Wet bulb depression (°C)

Green

50

4

60

50

5

40

55

8

35

55

8

30

60

10

25

65

15

20

70

20

15 to end

70

20

2

some extent as they also aim to minimize timber losses by degradation. Even between batches of a given timber type the moisture transport within the timber, the thickness of the boards, as well as the initial moisture content can vary widely. Consequently, despite the inherent conservativeness of standard drying schedules, degradation rates of up to 50% of the load are not unknown in kiln drying using these schedules. To take the full economical advantage of kiln drying, a drying schedule for each timber load should really be calculated on-line. This paper presents a solution for the scheduling problem based on nonlinear model predictive control techniques. These methods require models describing both the timber drying and the amount of internal stress, models which are introduced in the following section. Further sections discuss the calculation of drying schedules assuming perfect models (which permit the use of pure feedforward control) as well as the calculation of drying conditions incorporating on-line measurements and closed-loop control.

2 Drying and stress modelling 2.1

The drying model

In recent years, the drying of timber has been studied extensively and many drying models have been proposed (REFERENCES???). In this paper, only an introductory description of the basic modelling principles is given. The drying of timber can be described as coupled heat and mass transport, both to and from the surface and within the timber itself. Since timber boards are usually much longer and wider than they are thick, this three-dimensional transport problem can be reduced to a single dimension with minimal loss of accuracy [reference???]. With this assumption made, is has been shown that moisture transport within hardwood may be adequately described by Fick’s law [2] ∂X ∂ ⎛ ∂X⎞ ------ = D -----∂t ∂ x ⎝ ∂x ⎠

(1)

where D represents the temperature dependent diffusivity coefficient, X is the moisture content and x is the position within the timber board:

D = Dr e

D – ------ET

.

(2)

Similarly, heat transport within the timber is described by a second partial differential equation 3

∂T ∂ ⎛ k ∂T⎞ ------ = -------- -----∂t ∂ x ⎝ ρc p ∂x ⎠

(3)

Since the thermal conductivity k , the density ρ , and the heat capacity c p are functions of moisture content X , and the moisture diffusivity D is a function of temperature T , Equations (1) and (3) must be solved simultaneously. At the timber surface, the mass flux of water from the board ( J ) is proportional to the difference between the gas humidities just above the timber surface Y S and in the bulk air stream YG : J = ρ Air R ( Y S – Y G )

(4)

While the humidity of the air Y G can be calculated from the dry and wet bulb temperatures ( T G , T W ) in the kiln, Y S is a function of the moisture content at the timber surface X S and thus Equation (4) forms a boundary condition to Equation (1). At the board surface, both convective heat transport, as well as heat transfer associated with the evaporation of moisture must be considered: Q = h F ( T G – T S ) – ( λ W + λ S )J

(5)

In the above equation h F is the heat transfer coefficient, λ W is the latent heat of vaporization, λ S the heat of sorption, and T S the surface temperature. This equation forms a boundary condition for Equation (3). If we assume that the same amount of heat and mass transport is occurring at the upper and lower surfaces, then the temperature and moisture profiles will be symmetric about the centre of the timber board. Consequently, moisture and heat transport cannot occur across the centre, a fact which defines another two boundary conditions: dX ------dx dT -----dx

= 0

(6)

= 0

(7)

C

C

If the moisture and temperature gradients are approximated by finite differences (Reference for Crank Nicholson?), then the drying model reduces to a stiff system of ordinary differential equations: Surface 4

dX S 1 D ( X1 – XS ) --------- = -------- --------------------------- – J dt Δx S Δx S

(8)

k ( T1 – TS ) dT S 1 --------- = ----------------- ------------------------- + Q dt Δx S ρc p Δx S

(9)

dX i D Xi + 1 – Xi Xi – Xi – 1 -------- = ------ ----------------------- – ----------------------dt Δx Δx Δx

(10)

dT i Ti + 1 – Ti Ti – Ti – 1 k -------- = --------------- ---------------------- – ---------------------dt Δxρc p Δx Δx

(11)

dX C –D XC – XC – 1 ---------- = --------- -------------------------dt Δx C Δx C

(12)

dT C TC – TC – 1 –k --------- = ------------------ ------------------------Δx C dt Δx C ρc p

(13)

Timber

Centre

In the present case, the board was divided into 16 finite elements which led to a system of 34 differential equations (30 describing conditions within the board, and two each for the centre and the surface) which can be easily solved using an integrator suitable for stiff systems. 2.2

The stress model

As timber dries out, it will begin to shrink as soon as the moisture content falls below the fibre saturation point X fsp . Since we are using a one-dimensional drying model, it is reasonable to assume that the moisture content in every finite layer is uniform (see Figure 2). If any one layer i could move independently of all other layers, its shrinkage ε could be well approximated by the linear relation [3], ε i = βM i

(14)

where the moisture change is given by M i = min ( X i, X fsp ) and β is the moisture swelling coefficient. We could then calculate the stress σ which was required to keep the layer at original length using Hooke’s law together with a temperature and moisture dependent Young’s module E σ i = E i βM i

(15)

However, none of the layers can move independently as any board will keep its general shape as shrinkage occurs. Therefore, the internal forces must be calculated in relation to the average 5

shrinkage of the whole board. The sum of these internal forces must be zero, which can be expressed in terms of the finite elements as follows:

∑ Ei ( ε – βMi )Δxi

= 0

(16)

From (16) the average shrinkage of timber ε is easily calculated according to

∑ Ei βMi Δxi ε = ------------------------------∑ Ei Δxi

(17)

Consequently, for each thin layer the strain can be approximated as ε i = ε – βM i .

(18)

Figure 2 illustrates these relations. It should be noted that this simple stress model takes only Tension Layer Compression

Initial length

Shrinkage ε

Figure 2: Stress development in timber during drying

elastic components of the stress into account. However, particularly wet timber shows plastic behaviour which causes a slow relaxation of the local stress. In such cases, the stress levels calculated with this stress model will be higher than the true values. The complexity of the stress model used, however, dos not affect our approach to calculating optimal drying schedules.

6

3 Calculation of optimal drying schedules From an economic point of view, a drying schedule is optimal if some weighted combination of the overall drying time and the amount of timber degradations is kept as small as possible. Short drying times require severe drying conditions (high temperature, large wet bulb depression) which frequently causes high degradation rates. Therefore, the optimal schedule will involve some compromise between drying time and degradation rate. Since timber is a natural material with somewhat different characteristics from load to load, it is not possible to derive a very reliable correlation between degradation rate and stress (or other phenomena within the timber). Nevertheless it is well known that degradation becomes highly probable if a certain stress level within the timber is exceeded [reference?]. Timber also becomes very stiff at low moisture contents, a condition which increases the risk of degradation. Thus, excessive drying of the timber surface during kiln drying should also be avoided. Consequently, a reasonable description of an optimal drying schedule is one that involves the minimum possible drying time subject to the constraints of the local strain being below some upper limit ε max and the surface moisture content being above some minimum level X Smin during the entire drying period. It is worth noting that usually the main objective is to dry the timber to a specified average moisture content X req . Since drying requires a moisture gradient, once this objective has been achieved the moisture profile within the timber will not be flat with the moisture content at the centre exceeding that at the surface. Therefore, a second objective can be defined as flattening the moisture profile so that the difference between the maximum un minimum local moisture contents is below some nominated upper limit ΔX req . 3.1

The optimization problem

Initially we will assume that the drying and stress models perfectly describe the real drying process. This assumption is obviously not justified as the timber properties (even for a given species) will vary from load to load. At best, these models can only provide a reasonable approximation to reality. In Section 4, however, this problem will be considered further. The optimal drying schedule determination can be posed as an optimization problem (in terms of minimizing the drying time t f ) using the controller inputs u ( t ) , representing here the time dependent dry and wet-bulb temperatures, as the free (or optimization) variables. In Equation 7

(19), the parameters p describe the timber properties while X ( t ) and T ( t ) describe the moisture and temperature profiles within the timber:

min u(t)

subject to

tf dX ⎧ ------- – f [ X ( t ), T ( t ), u ( t ), p ] = 0 ⎪ dt ⎪ dT ⎪ ------- – g [ X ( t ), T ( t ), u ( t ), p ] = 0 ⎪ dt ⎪ ε i ( t ) < ε max ⎪ ⎪ X S ( t ) ≥ X Smin ⎪ ⎨ X ( t f ) = X req ⎪ ⎪ ⎪ max [ X ( t f ) ] – min [ X ( t f ) ] ≤ ΔX req ⎪ ⎪ T G ≤ T Gmax ⎪ ⎪ X ( t=0 ) = X 0 ⎪ ⎩ T ( t=0 ) = T 0

(19)

The solution of such an optimal control problem would require considerable computation time. Thus, some computationally less expensive approach is required. If the total drying time is split up into a “drying” phase and a subsequent “moisture profile flattening” phase, then the overall optimization problem can be solved as a series of smaller optimization problems. 3.2

The drying phase

In the first phase, the timber should be dried to the specified average moisture content without ever exceeding the maximum strain level and minimum surface moisture content constraints. If the total drying time is split up into a series of time increments Δt , the it is obvious that the first drying phase is minimized if the drying rate is maximized (subject to the above constraints) during each time interval. An alternative (but equivalent) objective here is the minimization of the average moisture content within the entire board at the final time t i + 1 for each time step i . If the time increments are kept small compared to the overall drying time, then constant dry and wet bulb temperatures can be used during each interval without significantly compromising the accuracy of the solution. If such time discretisation is used, then it is in fact sufficient to formulated the constraints as functions of the stress level and surface moisture content at the end point t i + 1 of each step. The resulting nonlinear programming problem to be solved at each step can now be stated as

8

min ui

subject to

X ( ti + 1 ) ⎧ dX ------- – f [ X ( t ), T ( t ), u i, p ] = 0 ⎪ dt ⎪ ⎪ dT ------- – g [ X ( t ), T ( t ), u i, p ] = 0 ⎪ dt . ⎪ ε j ( t i + 1 ) < ε max ⎪ ⎨ X S ( t i + 1 ) ≥ X Smin ⎪ ⎪ T G ≤ T Gmax ⎪ ⎪ X ( t=t i ) = X i ⎪ ⎪ T ( t=t i ) = T i ⎩

(20)

This approach is equivalent to nonlinear model predictive (or receding horizon) control using a one step prediction of the process behaviour [5]. Within the context of this study, Equation (20) was solved sequentially using Gear’s method for the solution of the differential equations, and a sequential quadratic programming method [Reference?] for the solution of the optimization problem. Figure 3 illustrates the results obtained using a one-step prediction horizon of four hours, a minimum surface moisture of 6% (dry basis), and a maximum absolute strain level of 0.022 . The initial conditions were a uniform timber temperature of 15°C and a uniform moisture content of 60%. The board thickness was 30 mm. The optimized schedule is considerable more aggressive than the standard drying schedule. However, despite the optimized schedule taking a considerable shorter time to achieve 12% moisture content, both the maximum and average stress levels are less than those encountered in the standard drying schedule. 3.3

The flattening phase

Once the required average moisture content in the timber has been achieved, the moisture profile within the timber should be flattened to improve the final quality of the product. In contrast to the drying phase, however, where the optimal schedule is obtained by simply combining the optima from each time interval, optimization of the flattening phase required a more “global” approach. If we assume N stepwise constant dry and wet-bulb temperatures during the flattening phase, then the optimization problem can be stated as

9

0.025 0.02

40

Max. strain

Av. moisture (% dry basis)

60

20

0.015 0.01

0 0

100

200 Time (h)

300

0.005 0

400

100

200 Time (h)

300

400

100

200 Time (h)

300

400

Wet bulb depression (°C)

Dry bulb temperature (°C)

0 80

60

40

20 0

100

200 Time (h)

300

-10

-20

-30 0

400

Figure 3: Comparison of optimal and standard drying schedules. Solid lines: Optimized schedule Dotted lines: Standard schedule

min u

t flat

⎧ ti + 1 – ti = tf ⁄ N ⎪ ⎪ dX ⎪ ------- – f [ X ( t ), T ( t ), u i, p ] = 0 ⎪ dt ⎪ dT ⎪ ------- – g [ X ( t ), T ( t ), u i, p ] = 0 ⎪ dt ⎪ ε j ( t i ) < ε max ⎪ (21) ⎪ subject to X S ( t i ) ≥ X Smin ⎨ ⎪ X ( t f ) = X req ⎪ ⎪ ⎪ max [ X ( t f ) ] – min [ X ( t f ) ] ≤ ΔX req ⎪ ⎪ T G ≤ T Gmax ⎪ ⎪ X ( t=t 0 ) = X 0 ⎪ ⎪ T ( t=t 0 ) = T 0 ⎩ This optimization problem can be solved using the same numerical methods as were employed for the drying phase. The results presented here are for a maximum moisture difference of

10

Wet bulb depression (°C)

Moisture (% dry basis)

0 15 12.5 Final profiles

10

Initial profile 7.5 5 0

10

20

-10

-20

-30

30

0

10

x (mm)

20 Time (h)

30

Figure 4 Optimum results for the flattening phase. Dotted lines: N = 1

Dash-dotted lines: N = 2

Dashed lines: N = 3

Solid lines: N = 4

ΔX req = 3% together with the same constraints for the stress and surface moisture content as imposed during the drying phase. As shown in Table 2, the time required for the flattening of Table 2: Time required for flattening of the moisture profile Number of steps (N)

Flattening time (h)

1

32.5

2

19.1

3

19.3

4

19.1

the moisture content profile depends on the number of steps used N . The time required is significantly shortened if the number of step is increased from one to two. However, further increases in the number of steps show negligible improvement. Since internal moisture transport is faster at the higher temperatures, it was not surprising that all calculated optimal results for this flattening phase employ the maximum dry bulb temperature (80°C). The initial and final moisture profiles as well as the optimum wet-bulb temperatures are shown in Figure 4.

4 Closed-loop control of the drying phase The solution to the scheduling problem presented in the previous section was based on the assumption of perfect models being available. In practice, the best that can be expected is that the drying and stress models are reasonable approximation to reality, in which case some form of feedback must be introduced to permit the moisture profile within the timber to be esti11

Objectives Controller

T G, T W

Timber in kiln

) )

X m, T G, T W

D r, X , T

Gross error detection

Observer

Figure 5 Structure of the closed-loop nonlinear model predictive controller.

mated as closely as possible. Unfortunately, on-line measurement of local moisture contents during kiln drying is not possible. The only feasible on-line measurement is the weight of a sample of the timber within the kiln. Since the initial average moisture content of the timber can be measured quite accurately prior to drying, measurement of the weight of a timber stack can be used to calculated the actual average moisture content of the timber X m . In timber drying, heat transport (in general) as well as moisture transport away from the timber surface are much faster processes than the transport of moisture within the timber. Consequently, internal moisture diffusion limits the drying rate and the parameters D r and D E describing the temperature dependency of the diffusivity coefficient are the most critical parameters in the overall drying model. Using a history horizon approach (Reference?), these two parameters could theoretically be estimated on-line by minimizing the error between the measured average moisture content and that predicted by the drying model. However, numerical experiments showed that the temperature range in the data used was not large enough (since the time horizon has to be kept relatively short) to obtain reliable estimates of both parameters. Fortunately, as the temperature setpoints during the drying phase are changed in relatively small steps, a precise knowledge of the activation energy D E is not necessary for a sufficiently accurate prediction of the internal moisture profiles to be calculated. It would appear from our results to be adequate to use an estimate of the activation energy that had been calculated (for any given species of timber) by some form of off-line regression using a complete set of drying data. The structure of the resulting closed-loop control system is illustrated by Figure 5. The controller contains the solution of the open-loop drying phase problem (see Equation 20) but 12

) )

using the estimated moisture and temperature profiles ( X , T ) as well as the estimated parameter D r . This parameter was estimated “on-line” by solving the following nonlinear programming problem:

)

∑

X ( tj ) – Xm ( tj )

)

min Dr

2

)

= 0 (22)

)

= 0

)

subject to

)

)

)

)

⎧dX - – f [ X ( t ), T ( t ), u i, p ] ⎪ -------dt ⎪ ⎪ T ⎪ d-------- – g [ X ( t ), T ( t ), u i, p ] ⎨ dt ⎪ ⎪ X ( t=t i ) ⎪ ⎪ ⎩ T ( t=t i )

)

)

= Xi = Ti

Here X m ( t j ) represent the “measured” average moisture contents values (that is, the average moisture contents that would be calculated from the weight balance readings), filtered to remove gross errors. For this simulation of closed-loop behaviour, an error of +15% was introduced into the diffusion parameter D E and an error of –40% into the surface mass transfer coefficient. Further, white noise with? was added to the “measured” average moisture content values. Drying phase results were calculated using the same timber initial conditions and constrains as before. A one-step control horizon of four hours and a three-step observation horizon (of 12 hours) were used (Figures 6 and7. Despite the introduction of quite large errors in both model parameters and moisture content measurements, the estimated moisture profile were very accurate indeed. Consequently, the prediction of the maximum local strain was very close to the “true” value with the overall result that the time required for the (closed-loop) drying phase was just one (four hour) step more than that required for the open-loop solution based on perfect drying and stress models. Although is was not carried out as part of this study, the previously described open-loop approach could be used for the (relatively short duration) flattening phase using the latest estimated value for the parameter D r .

5 Conclusions This study has clearly shown that the use of nonlinear model predictive control offers considerable economic benefits for the kiln drying of timber. In comparison with a standard hard13

60

Moisture (% dry basis)

50

40

30

20

t(Δt=8h) 10 0

10

20

x (mm)

30

60

0.025 0.02 Max. strain

40

20

0.015 0.01

0 0

50

100 150 Time (h)

200

0.005 0

250

60

40

20 0

100 Time (h)

200

100 150 Time (h)

200

0

80

Wet bulb depression (°C)

Dry bulb temperature (°C)

Av. moisture (% dry basis)

Figure 6: Estimated and “true” moisture profiles during closed-loop control of the kiln drying. Solid lines: Estimated profiles Dotted lines: “True” profiles

50

100 150 Time (h)

200

-10

-20

-30 0

250

50

Figure 7: Closed-loop simulation results for the drying phase. Upper right: Estimated strain: solid line, true strain: dotted line

14

250

wood drying schedule, the use of such a controller offers significantly reduced drying times. In addition, as constraints for maximum local strain and minimum surface moisture content can be explicitly included in the controller design specification, decreased degradation rates would also be expected. Most importantly though, the on-line measurements necessary to implement a closed-loop version of such a control scheme and obtain these economic benefits for an operational kiln are modest indeed.

6 References [1]

Campbell, G. S.: “Index of Kiln Dying Schedules for Timber Dried in Australia,” CSIRO Building Research Division, (1980)

[2]

Collignan, A., J. P. Nadeau, J. R. Puiggali: “Description and Analysis of Timber Drying Kinetics,” Drying Technology, 11, 489-506 (1993)

[3]

Johnson, J. A.: “Stress Development,” IUFRO 1989 Wood Drying Symposium, Seattle, Washington, Appendix B (1989)

[4]

Little, R. L., R. L. Toennison: “Drying Hardwood Lumber Using Computer Controller Mini-Step Schedules,” IUFRO 1989 Wood Drying Symposium, Seattle, Washington, 203-212 (1989)

[5]

Patwardhan, A. A., J. B. Rawlings, and T. F. Edgar: “Nonlinear Model Predictive Control,” Chem. Eng. Comm., 87, 123-141 (1990)

15