Dual-Scale Modelling of the Vacuum Drying Process for Transformer Cellulose-Based Insulation

1. Introduction

It is a well-established fact that the moisture inside the cellulose-based insulation increases the aging rate and shortens the service life of the apparatus [1]. To attain acceptable level of residual moisture content, before being put into service, the transformers are subjected to a thorough drying process which typically involves cyclic heating and vacuum application [2,3].

When it comes to establishing the control strategy for cellulose insulation drying process, the manufacturers frequently rely on experience [3] which is the reason why some recent works have been focused on mathematical modelling of moisture dynamics inside the non-impregnated cellulose insulation. Namely, by having at disposal a reliable and flexible model of vacuum drying process, one can improve it in terms of drying time or energy consumption.

A common approach in describing the moisture transport inside the cellulose insulation is by using Fick’s law in which diffusion coefficient is determined empirically. Notable work in this area was done by Du et al. [4] who determined diffusion coefficient for non-impregnated pressboard. Afterwards, this work was extended by García et al. [5] and Villarroel et al. [6] who experimentally determined the moisture diffusion coefficient in the transformer paper as a function of temperature, moisture content and thickness of the insulation. The similar methods were more recently used by Wang et al. [7] to additionally allow for the effects of insulation ageing on the value of moisture diffusion coefficient.

While results of those models are in good agreement with experimental values, as an intrinsic property of the material, it is not physically based for the diffusion coefficient to be dependent on the dimensions of the sample such as thickness. As stated in [8,9,10] this is a clear indicator that the drying process cannot solely be described by the Fick’s law which was also noticed by García et al. [11] who demonstrated that the thickness dependence of the diffusion factor is not a consequence of interlayer spaces or voids but rather a compensation for inaccuracy of the used diffusion model.

This shortcoming of the classical Fickian description of moisture transport was also recognized by Kang et al. [12] who used a Langmuir model to study the moisture dynamics. This model is based on assumption that the moisture contained within the insulation can be separated into free and bound moisture which are mutually convertible. To extract the bound moisture, it first needs to be transformed into a free moisture which is then transported according to the Fick’s law. Additional parameters of this model are adsorption and desorption coefficients which were determined by fitting the analytical solution of the Langmuir model to the experimentally obtained drying curves. The comparison between the models showed that the Fick model significantly underestimates the necessary drying time of the pressboard.

On the other hand, Brahami et al. [3] adopted an empirical model to describe the drying of the transformer insulation. This model additionally considers degradation of the cellulose insulation due to its exposure to oxygen and elevated temperatures thus providing a mathematical tool for manufacturers to determine the desired compromise between energy consumption and drying time of the process.

All the above-mentioned works adopt assumption of the uniform temperature distribution inside the cellulose insulation. Although heat diffuses through insulation faster than moisture does, there are indications that considerable temperature gradients through cellulose insulation might exist even in later stages of the drying process. They might be caused by local cooling due to the evaporation of water vapor from the cellulose fibres [13] or uneven heating [14]. Furthermore, as drying process consists of series of cycles of heating and vacuum application, transient periods will occur in-between the phases, making it unclear on how boundary conditions on the interface between cellulose insulation and chamber atmosphere should be imposed in those periods. The vacuum chamber atmosphere also interacts with the cellulose insulation through the exchange of heat and mass thus making the whole problem strongly coupled. Accurately capturing this behaviour requires a modelling framework that accounts for both internal dynamics of the cellulose insulation and the surrounding chamber atmosphere. This complexity has been previously addressed in other fields such as wood vacuum drying [15,16,17,18] and pharmaceutical freeze drying [19] where dual-scale modelling approach has been used. The basic idea is to introduce a computational fluid dynamics (CFD) [20] or simplified zero-dimensional (0D) models at the scale of the dryer thus extending the existing (macroscale) model of heat and mass transfer through the porous media. Such models enable more accurate predictions of energy consumption and drying time thus supporting the process analysis and optimization.

In this study, a comprehensive mathematical model that describes the vacuum drying process at both the scale of the cellulose insulation and at the scale of vacuum chamber is presented. The model is implemented in Python and employed to simulate various process configurations. Numerical experiments conducted for a simplified drying system highlight the influence of key model parameters and demonstrate the model’s potential as an effective tool for analysing and optimizing the vacuum drying process of transformer cellulose-based insulation.

2. Mathematical Model

The mathematical model that describes the vacuum drying process of transformer cellulose-based insulation is structured to capture the relevant physical phenomena occurring within the drying system. The following subsections first provide a description of the observed drying system, followed by the formulation of the governing equations, interface conditions between the cellulose and the chamber atmosphere, initial conditions and numerical solution procedure.

2.1. Drying System Description

In this study, a simplified yet sufficiently sophisticated drying system is considered, designed to enable control of the two primary process variables: drying temperature and pressure. As illustrated in Figure 1, the system consists of a thermally insulated vacuum chamber with two ports. The first port is connected to a vacuum pump through the suction line, with a throttle valve installed to reduce or completely shut off the pumping speed when necessary. Additionally, to ensure that no significant over-pressure occurs and to enable the venting process at the end of each vacuum phase, a vacuum chamber is equipped with a valve that has both breathing and venting function. Lastly, the heating of the cellulose insulation is achieved with infra-red (IR) heaters mounted on the inner walls of the vacuum chamber shell. In addition to radiation, the heat is also being transferred by natural convection since no fan or forced air circulation is employed.

Inside the chamber, transformer samples with simple geometry are placed. Each sample consists of cellulose insulation paper wrapped around a cylindrical metal core, serving as a reasonable approximation for actual transformer. In that case, one can assume that the outer side of the sample is completely exposed to the vacuum chamber conditions and inner side is thermally insulated and impermeable to mass transfer. The number, size and spatial arrangement of the samples can vary depending on the available space inside the vacuum chamber.

2.2. Governing Equations at the Scale of the Cellulose Insulation

Unlike previous studies, different approach is considered in modelling the heat and mass transport through the cellulose insulation. Knowing that the transformer insulation is composed of the cellulose fibres containing moisture and pores filled with humid air (Figure 2), it is assumed that the moisture inside the fibres is practically immobile and that its transport is achieved in the form of water vapor through the pore space by means of diffusion and convection. Furthermore, a local non-equilibrium is assumed between cellulose fibres and humid air inside the pores. The moisture content of the fibres will tend to reach the equilibrium state which will then result in local moisture transport between the fibres and humid air inside the pores. That is, moisture inside cellulose fibres will either be absorbed or desorbed, depending on the difference between current moisture content and equilibrium moisture content.

The cellulose insulation usually has cylindrical shape with radial dimensions much smaller than height. Additionally, the conditions on the outer and inner boundaries are assumed to be uniform along the perimeter and height therefore making it reasonable to assume that variations of physical fields in axial and azimuthal directions can be neglected. In that case, the transport of heat and mass is considered only in the radial direction, as shown on Figure 1.

The pore space of the insulation can be filled with dry air and water vapor so those species must be conserved. The conservation law for dry air has the form of convection-diffusion equation:

$${\epsilon }_{\mathrm{p}}{\displaystyle \frac{\partial C_{\mathrm{a}}}{\partial t}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(rv{C}_{\mathrm{a}}\right)+{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(r{D}_{\mathrm{e}\mathrm{q}}{C}_{\mathrm{g}}{\displaystyle \frac{\partial y_{\mathrm{v}}}{\partial r}}\right)=0,$$

(1)

whereas, in conservation law for water vapor additional source term emerges due to the absorption or desorption of moisture from the cellulose fibres:

$${\epsilon }_{\mathrm{p}}{\displaystyle \frac{\partial C_{\mathrm{v}}}{\partial t}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(rv{C}_{\mathrm{v}}\right)-{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(r{D}_{\mathrm{e}\mathrm{q}}{C}_{\mathrm{g}}{\displaystyle \frac{\partial y_{\mathrm{v}}}{\partial r}}\right)+{\displaystyle \frac{{\rho }_{\mathrm{c}}}{M_{\mathrm{v}}}}{\displaystyle \frac{\partial X}{\partial t}}=0.$$

(2)

In equations (1) and (2), $C_{\mathrm{a}}$ (kmol/m³) and $C_{\mathrm{v}}$ (kmol/m³) are molar densities of air and water vapor respectively. Both water vapor and dry air are assumed to behave like ideal gases, so following relations hold:

$$p=p_{\mathrm{v}}+p_{\mathrm{a}}=C_{\mathrm{g}}{R}_{\mathrm{m}}T,$$

(3)

$$C_g=C_{\mathrm{v}}+C_{\mathrm{a}},$$

(4)

$$y_{\mathrm{v}}={\displaystyle \frac{C_{\mathrm{v}}}{C_{\mathrm{g}}}}={\displaystyle \frac{p_{\mathrm{v}}}{p}},$$

(5)

where $p$ is the total pressure of the gaseous phase with $p_{\mathrm{v}}$ and $p_{\mathrm{a}}$ being partial pressures of water vapor and dry air, respectively.

The equivalent binary diffusivity $D_{\mathrm{e}\mathrm{q}}$ in equations (1) and (2) includes Knudsen effect that could arise during vacuum phases and is given by the following relation [18]:

$${\displaystyle \frac{1}{D_{\mathrm{e}\mathrm{q}}}}={\displaystyle \frac{1}{D_{\mathrm{e}\mathrm{f}}}}+{\displaystyle \frac{1}{D_{\mathrm{K}}}},$$

(6)

where, $D_{\mathrm{e}\mathrm{f}}$ is an effective binary diffusivity of water vapor and dry air mixture, expressed as a function of absolute temperature and pressure [21]:

$$D_{\mathrm{e}\mathrm{f}}={\displaystyle \frac{{\epsilon }_{\mathrm{p}}}{{\tau }_{\mathrm{p}}}}1.8947775·{10}^{-5}{\displaystyle \frac{T^{2.072}}{p}},$$

(7)

and $D_{\mathrm{K}}$ (m²/s) is a Knudsen diffusivity that should be empirically determined during the model calibration [18].

The velocity of the gaseous phase is calculated according to Darcy’s law:

$$v=-{\displaystyle \frac{k_{\mathrm{e}\mathrm{f}}}{{\mu }_{\mathrm{g}}}}{\displaystyle \frac{\partial p}{\partial r}},$$

(8)

where, $K_{\mathrm{e}\mathrm{f}}$ (m²) is an effective permeability coefficient of the cellulose insulation obtained with Klinkenberg correction to account for the slip effects at low pressures [22]:

$$k_{\mathrm{e}\mathrm{f}}=k_0\left(1+{\displaystyle \frac{b}{p}}\right).$$

(9)

In the equation (9) $k_0$ (m²) is an absolute permeability that is measured at large pressures and $b$ (Pa) is a Klinkenberg correction factor that increases with decrease of absolute permeability $k_0$.

Moisture inside the cellulose fibres can be absorbed or desorbed depending on the equilibrium moisture content which it tends to reach. The mass conservation law for moisture contained in the cellulose fibres has the following form:

$${\displaystyle \frac{\partial X}{\partial t}}=K\left(X_{\mathrm{e}\mathrm{q}}-X\right),$$

(10)

where $K$ (1/s) is an internal overall mass transfer coefficient whose value depends on the resistance of the moisture transport from the inner parts of the cellulose fibres to the pore-fibre interface. It also strongly depends on the total pore-fibre interface surface area per unit of volume of the insulation. It is not an easy task to accurately calculate $K$ so a simplified approach is undertaken in this work by assuming that it has a constant value that should be empirically determined.

Dry-basis equilibrium moisture content $X_{\mathrm{e}\mathrm{q}}$ in equation (10) is a function of temperature and water activity of the surrounding media. An analytical form with temperature and water vapor partial pressure as independent variables was early given by Fessler et al. [24] and later corrected by Du et al. [25]. This expression, although accurate for low moisture contents, might give non-physical results for water activities closer to unity. Przybylek et al. [23] calculated the parameters of the GAB (Guggenheim, Anderson and de Boer) sorption isotherm model based on the experimental results for Kraft paper samples with different ageing degrees. This model gives accurate results in a wider range of relative humidities (Figure 3) and was used in calculating equilibrium moisture content in this model.

Since water activity was not chosen as a dependent variable of the model, it is expressed in terms of molar density and temperature. This follows from the definition of water activity and the assumption that dry air and water vapor are ideal gases.

$$a_{\mathrm{w}}(C_{\mathrm{v}},T)={\displaystyle \frac{p_{\mathrm{v}}}{p_{\mathrm{s}\mathrm{a}\mathrm{t}}\left(T\right)}}={\displaystyle \frac{C_{\mathrm{v}}{R}_{\mathrm{m}}T}{p_{\mathrm{s}\mathrm{a}\mathrm{t}}\left(T\right)}}.$$

(11)

In the above equation, $p_{\mathrm{s}\mathrm{a}\mathrm{t}}\left(T\right)$ is water vapor saturation pressure that can easily be calculated for a known temperature [26].

Equilibrium moisture content according to the GAB model is calculated from the following expressions:

$$X_{\mathrm{e}\mathrm{q}}\left(C_{\mathrm{v}},T\right)={\displaystyle \frac{X_m·C\left(T\right)·k·a_{\mathrm{w}}\left(C_{\mathrm{v}},T\right)}{\left[1-k·a_{\mathrm{w}}\left(C_{\mathrm{v}},T\right)\right]·\left[1+k·a_{\mathrm{w}}\left(C_{\mathrm{v}},T\right)·\left(C\left(T\right)-1\right)\right]}},$$

(12)

$$C\left(T\right)=C_0·\mathrm{e}\mathrm{x}\mathrm{p}\left[{\displaystyle \frac{Q}{R_{\mathrm{m}}}}\left({\displaystyle \frac{1}{T}}-{\displaystyle \frac{1}{T_{\mathrm{r}\mathrm{e}\mathrm{f}}}}\right)\right].$$

(13)

For more detailed explanation of GAB model and physical meaning behind parameters $X_m$, $k$, $C_0$, $Q$ and $T_{\mathrm{r}\mathrm{e}\mathrm{f}}$, one can refer to the work of Staudt et al. [27]. The values of the above parameters used for this model are taken from [23] and are given in Table A1 in the Appendix.

To close the system of the governing equations at the scale of the cellulose insulation, conservation of energy needs to be applied. Unlike moisture concentration, an assumption of local thermal equilibrium between cellulose fibres and gaseous phase inside the pores is assumed. The specific heat of the insulation is assumed to be equal to the specific heat of the solid cellulose fibres, so the conservation of energy reads:

$${\epsilon }_{\mathrm{f}}{\rho }_{\mathrm{f}}{c}_{\mathrm{f}}{\displaystyle \frac{\partial T}{\partial t}}-{\lambda }_{\mathrm{e}\mathrm{f}}{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(r{\displaystyle \frac{\partial T}{\partial r}}\right)-{\mathsf{\Delta }h}_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{p}}{\epsilon }_{\mathrm{f}}{\rho }_{\mathrm{f}}{\displaystyle \frac{\partial X}{\partial t}}=0,$$

(14)

where ${\lambda }_{\mathrm{e}\mathrm{f}}$ (W/(m∙K)) is an effective thermal conductivity of the cellulose insulation and ${\mathsf{\Delta }h}_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{p}}$ (J/kg) is specific heat of water evaporation.

2.3. Governing Equations at the Scale of the Vacuum Chamber

The governing equations at the scale of the vacuum chamber are derived under the assumption that temperature, pressure and molar densities of water vapor and dry air are uniformly distributed fields within the available volume $V$ (m³) inside the vacuum chamber. Moreover, the temperature fields of the chamber walls and IR heaters are also considered uniform and dependent only on time variable.

The gaseous phase (water vapor and dry air) can enter the chamber through the venting valve and leave the chamber through the pressure relief valve or during the pump down process. At the interface between the transformer cellulose insulation and vacuum chamber atmosphere, the molar flow rate of gaseous phase can be directed either inward or outward. The conservation laws are thus expressed accordingly:

$$V{\displaystyle \frac{d{C}_{\mathrm{v},\mathrm{c}\mathrm{h}}}{dt}}={\dot{n}}_{\mathrm{v}\mathrm{e}\mathrm{n}\mathrm{t},\mathrm{v}}-{\dot{n}}_{\mathrm{r}\mathrm{e}\mathrm{l}\mathrm{i}\mathrm{e}\mathrm{f},\mathrm{v}}-{\dot{n}}_{\mathrm{p}\mathrm{u}\mathrm{m}\mathrm{p},\mathrm{v}}+{\dot{n}}_{\mathrm{t},\mathrm{v}},$$

(15)

$$V{\displaystyle \frac{d{C}_{\mathrm{a},\mathrm{c}\mathrm{h}}}{dt}}={\dot{n}}_{\mathrm{v}\mathrm{e}\mathrm{n}\mathrm{t},\mathrm{a}}-{\dot{n}}_{\mathrm{r}\mathrm{e}\mathrm{l}\mathrm{i}\mathrm{e}\mathrm{f},\mathrm{a}}-{\dot{n}}_{\mathrm{p}\mathrm{u}\mathrm{m}\mathrm{p},\mathrm{a}}+{\dot{n}}_{\mathrm{t},\mathrm{a}},$$

(16)

After expanding and rearranging the terms, the following equations are obtained:

$${\displaystyle \frac{d{C}_{\mathrm{v},\mathrm{c}\mathrm{h}}}{dt}}={\displaystyle \frac{C_{\mathrm{v},\mathrm{a}\mathrm{m}\mathrm{b}}{Q}_{\mathrm{v}\mathrm{e}\mathrm{n}\mathrm{t}}-C_{\mathrm{v},\mathrm{c}\mathrm{h}}{Q}_{\mathrm{r}\mathrm{e}\mathrm{l}\mathrm{i}\mathrm{e}\mathrm{f}}-C_{\mathrm{v},\mathrm{c}\mathrm{h}}{Q}_{\mathrm{p}\mathrm{u}\mathrm{m}\mathrm{p}}+A_{\mathrm{t}}\sum _{i=1}^{N_{\mathrm{t}}}{J}_{\mathrm{v},t, i}}{V}},$$

(17)

$${\displaystyle \frac{d{C}_{\mathrm{a},\mathrm{c}\mathrm{h}}}{dt}}={\displaystyle \frac{C_{\mathrm{a},\mathrm{a}\mathrm{m}\mathrm{b}}{Q}_{\mathrm{v}\mathrm{e}\mathrm{n}\mathrm{t}}-C_{\mathrm{a},\mathrm{c}\mathrm{h}}{Q}_{\mathrm{r}\mathrm{e}\mathrm{l}\mathrm{i}\mathrm{e}\mathrm{f}}-C_{\mathrm{a},\mathrm{c}\mathrm{h}}{Q}_{\mathrm{p}\mathrm{u}\mathrm{m}\mathrm{p}}+A_{\mathrm{t}}\sum _{i=1}^{N_{\mathrm{t}}}{J}_{\mathrm{a},t, i}}{V}}.$$

(18)

In the above equations $Q_{\mathrm{v}\mathrm{e}\mathrm{n}\mathrm{t}}$ and $Q_{\mathrm{r}\mathrm{e}\mathrm{l}\mathrm{i}\mathrm{e}\mathrm{f}}$ are volume flow rates through venting valve and pressure relief valve, respectively. Since a single valve is used for both venting and pressure relief, volume flow rates $Q_{\mathrm{v}\mathrm{e}\mathrm{n}\mathrm{t}}$ and $Q_{\mathrm{r}\mathrm{e}\mathrm{l}\mathrm{i}\mathrm{e}\mathrm{f}}$ are mutually excluding (only one at the time can have non-zero value). In view of that, $Q_{\mathrm{v}\mathrm{e}\mathrm{n}\mathrm{t}}$ and $Q_{\mathrm{r}\mathrm{e}\mathrm{l}\mathrm{i}\mathrm{e}\mathrm{f}}$ are calculated as follows:

$$Q_{vent}=\left\{\begin{array}{l}{K}_{\mathrm{v}}\sqrt{p_{\mathrm{a}\mathrm{m}\mathrm{b}}-p_{\mathrm{c}\mathrm{h}}}, p_{\mathrm{a}\mathrm{m}\mathrm{b}}>p_{\mathrm{c}\mathrm{h}}\\ 0, p_{\mathrm{a}\mathrm{m}\mathrm{b}}\le p_{ch}\end{array}\right.,$$

(19)

$$Q_{\mathrm{r}\mathrm{e}\mathrm{l}\mathrm{i}\mathrm{e}\mathrm{f}}=\left\{\begin{array}{l}{K}_{\mathrm{v}}\sqrt{p_{\mathrm{c}\mathrm{h}}-p_{\mathrm{a}\mathrm{m}\mathrm{b}}}, p_{\mathrm{c}\mathrm{h}}>p_{\mathrm{a}\mathrm{m}\mathrm{b}}\\ 0, p_{\mathrm{c}\mathrm{h}}\le p_{ch}\end{array}\right.,$$

(20)

where $p_{\mathrm{a}\mathrm{m}\mathrm{b}}$ (Pa) and $p_{\mathrm{c}\mathrm{h}}$ (Pa) are ambient and chamber pressures respectively. During the venting process, a choked flow might occur, however, since venting is relatively rapid part of the process in comparison to other phases, the chocked flow phenomena is not considered in this work.

The pumping speed at the outlet port of the vacuum chamber (effective pumping speed), $Q_{\mathrm{p}\mathrm{u}\mathrm{m}\mathrm{p}}$ (m³/s), is generally different from the pumping speed at the inlet of the vacuum pump. The inlet pumping speed $Q_{\mathrm{p}\mathrm{u}\mathrm{m}\mathrm{p},\mathrm{i}\mathrm{n}}$ can be found in the manufacturer’s datasheet as a function of the pressure at the vacuum pump inlet $p_{\mathrm{i}\mathrm{n}}$. To calculate effective pumping speed, throughput continuity equation under the assumption of isothermal flow in the suction line is used:

$$Q_{\mathrm{p}\mathrm{u}\mathrm{m}\mathrm{p}}={\displaystyle \frac{p_{\mathrm{i}\mathrm{n}}·Q_{\mathrm{p}\mathrm{u}\mathrm{m}\mathrm{p},\mathrm{i}\mathrm{n}}\left(p_{\mathrm{i}\mathrm{n}}\right)}{p_{\mathrm{c}\mathrm{h}}}}.$$

(21)

To obtain the pressure at the inlet of the vacuum pump, effective pumping speed from the equation (21) is substituted into the equation for the conductance of the suction line:

$$C_{\mathrm{s}}={\displaystyle \frac{p_{\mathrm{i}\mathrm{n}}·Q_{\mathrm{p}\mathrm{u}\mathrm{m}\mathrm{p},\mathrm{i}\mathrm{n}}\left(p_{\mathrm{i}\mathrm{n}}\right)}{p_{\mathrm{c}\mathrm{h}}-p_{\mathrm{i}\mathrm{n}}}}.$$

(22)

For a given conductance $C_{\mathrm{s}}$ (m³/s) and pressure inside the vacuum chamber $p_{\mathrm{c}\mathrm{h}}$, equation (22) is solved for pressure at the inlet $p_{\mathrm{i}\mathrm{n}}$. This pressure is then substituted into equation (21) to obtain effective pumping speed $Q_{\mathrm{p}\mathrm{u}\mathrm{m}\mathrm{p}}$ which is then used in the governing equations. More details regarding the calculation of the effective pumping speed can be found in [28].

The energy conservation laws for dry air and water vapor can be incorporated into a single equation since both species are at the same temperature $T_{\mathrm{c}\mathrm{h}}$. Accumulation of the internal energy of the open system is due to the convective heat flow rate $\dot{Q}$ and enthalpy flows through the boundary:

$$V{\displaystyle \frac{d\left[\left(C_{\mathrm{m},v,\mathrm{v}}{C}_{\mathrm{v},\mathrm{c}\mathrm{h}}+{C_{\mathrm{m},v,\mathrm{a}}C}_{\mathrm{a},\mathrm{c}\mathrm{h}}\right)T_{\mathrm{c}\mathrm{h}}\right]}{dt}}=\dot{Q}+{\dot{H}}_{\mathrm{i}\mathrm{n}}-{\dot{H}}_{\mathrm{o}\mathrm{u}\mathrm{t}}.$$

(23)

The enthalpy difference can be written as:

$${\dot{H}}_{\mathrm{i}\mathrm{n}}-{\dot{H}}_{\mathrm{o}\mathrm{u}\mathrm{t}}=Q_{\mathrm{v}\mathrm{e}\mathrm{n}\mathrm{t}}{C}_{\mathrm{g},\mathrm{a}\mathrm{m}\mathrm{b}}{C}_{\mathrm{m},p,\mathrm{g}}{T}_{\mathrm{a}\mathrm{m}\mathrm{b}}-\left(Q_{\mathrm{r}\mathrm{e}\mathrm{l}\mathrm{i}\mathrm{e}\mathrm{f}}+Q_{\mathrm{p}\mathrm{u}\mathrm{m}\mathrm{p}}\right)C_{\mathrm{g},\mathrm{c}\mathrm{h}}{C}_{\mathrm{m},p,\mathrm{g}}{T}_{\mathrm{c}\mathrm{h}}+A_{\mathrm{t}}\sum _{i=1}^{N_{\mathrm{t}}}{J}_{\mathrm{v},t, i}{C}_{\mathrm{m},p,\mathrm{v}}{T}_{\mathrm{s},i}+A_{\mathrm{t}}\sum _{i=1}^{N_{\mathrm{t}}}{J}_{\mathrm{a},t, i}{C}_{\mathrm{m},p,\mathrm{a}}{T}_{\mathrm{s},i},$$

(24)

where $T_{\mathrm{s},i}$ is cellulose insulation surface temperature of the i-th transformer.

Atmosphere in the chamber can convectively exchange heat with chamber walls, IR heaters and cellulose insulation of the transformers. The term $\dot{Q}$ in equation (23) can therefore be written as:

$$\dot{Q}=A_{\mathrm{t}}\sum _{i=1}^{N_{\mathrm{t}}}{h}_{\mathrm{t},i}\left(T_{\mathrm{s},i}-T_{\mathrm{c}\mathrm{h}}\right)+A_{\mathrm{h}}{h}_{\mathrm{h}}\left(T_{\mathrm{h}}-T_{\mathrm{c}\mathrm{h}}\right)+A_{\mathrm{w}}{h}_{\mathrm{w}}\left(T_{\mathrm{w}}-T_{\mathrm{c}\mathrm{h}}\right),$$

(25)

where $h_{\mathrm{t},i}$, $h_{\mathrm{h}}$ and $h_{\mathrm{w}}$ are convective heat transfer coefficients at the interface between chamber atmosphere and other components: transformer cellulose insulation, IR heaters and chamber walls.

After combining equations (23), (24), (25) and rearranging, the final form of the energy equation is obtained:

$$\begin{array}{l}{\displaystyle \frac{d{T}_{\mathrm{c}\mathrm{h}}}{dt}}={\displaystyle \frac{A_{\mathrm{t}}\sum _{i=1}^{N_{\mathrm{t}}}{h}_{\mathrm{t},i}\left(T_{\mathrm{s},i}-T_{\mathrm{c}\mathrm{h}}\right)+A_{\mathrm{h}}{h}_{\mathrm{h}}\left(T_{\mathrm{h}}-T_{\mathrm{c}\mathrm{h}}\right)+A_{\mathrm{w}}{h}_{\mathrm{w}}\left(T_{\mathrm{w}}-T_{\mathrm{c}\mathrm{h}}\right)}{V\left(C_{\mathrm{m},v,\mathrm{v}}{C}_{\mathrm{v},\mathrm{c}\mathrm{h}}+C_{\mathrm{m},v,\mathrm{a}}{C}_{\mathrm{a},\mathrm{c}\mathrm{h}}\right)}}\\ +{\displaystyle \frac{Q_{\mathrm{v}\mathrm{e}\mathrm{n}\mathrm{t}}{C}_{\mathrm{g},\mathrm{a}\mathrm{m}\mathrm{b}}}{V\left(C_{\mathrm{v},\mathrm{c}\mathrm{h}}+C_{\mathrm{a},\mathrm{c}\mathrm{h}}\right)}}\left(T_{\mathrm{a}\mathrm{m}\mathrm{b}}{\gamma }_{\mathrm{g}}-T_{\mathrm{c}\mathrm{h}}\right)+{\displaystyle \frac{Q_{\mathrm{p}\mathrm{u}\mathrm{m}\mathrm{p}}+Q_{\mathrm{r}\mathrm{e}\mathrm{l}\mathrm{i}\mathrm{e}\mathrm{f}}}{V}}{T}_{\mathrm{c}\mathrm{h}}\left(1-{\gamma }_{\mathrm{g}}\right)\\ +{\displaystyle \frac{C_{\mathrm{m},v,\mathrm{v}}{A}_{\mathrm{t}}\sum _{i=1}^{N_{\mathrm{t}}}{J}_{\mathrm{v},t, i}\left(T_{\mathrm{s},i}{\gamma }_{\mathrm{v}}-T_{\mathrm{c}\mathrm{h}}\right)}{V\left(C_{\mathrm{m},v,\mathrm{v}}{C}_{\mathrm{v},\mathrm{c}\mathrm{h}}+C_{\mathrm{m},v,\mathrm{a}}{C}_{\mathrm{a},\mathrm{c}\mathrm{h}}\right)}}+{\displaystyle \frac{C_{\mathrm{m},v,\mathrm{a}}{A}_{\mathrm{t}}\sum _{i=1}^{N_{\mathrm{t}}}{J}_{\mathrm{a},t, i}\left(T_{\mathrm{s},i}{\gamma }_{\mathrm{a}}-T_{\mathrm{c}\mathrm{h}}\right)}{\left(C_{\mathrm{m},v,\mathrm{v}}{C}_{\mathrm{v},\mathrm{c}\mathrm{h}}+C_{\mathrm{m},v,\mathrm{a}}{C}_{\mathrm{a},\mathrm{c}\mathrm{h}}\right)}}.\end{array}$$

(26)

where $C_{\mathrm{g},\mathrm{a}\mathrm{m}\mathrm{b}}$ (kmol/m³) is molar density of ambient air and ${\gamma }_{\mathrm{g}}$ is specific heat ratio of the gaseous phase inside the vacuum chamber.

Energy equation for the IR heaters, derived in a similar fashion, reads:

$${\displaystyle \frac{d{T}_{\mathrm{h}}}{dt}}={\displaystyle \frac{1}{m_{\mathrm{h}}{c}_{\mathrm{h}}}}\left(A_{\mathrm{h}}{h}_{\mathrm{h}}\left(T_{\mathrm{c}\mathrm{h}}-T_{\mathrm{h}}\right)+{\dot{Q}}_{\mathrm{r}\mathrm{a}\mathrm{d},\mathrm{h}}+{\dot{Q}}_{\mathrm{e}\mathrm{l}}\right),$$

(27)

where ${\dot{Q}}_{\mathrm{r}\mathrm{a}\mathrm{d},\mathrm{h}}$ (W) and ${\dot{Q}}_{\mathrm{e}\mathrm{l}}$ (W) are radiative heat flow rate from heaters surface and input electrical power respectively.

Vacuum chamber walls energy equation reads:

$${\displaystyle \frac{d{T}_{\mathrm{w}}}{dt}}={\displaystyle \frac{1}{m_{\mathrm{w}}{c}_{\mathrm{w}}}}\left(A_{\mathrm{w}}{h}_{\mathrm{w}}\left(T_{\mathrm{c}\mathrm{h}}-T_{\mathrm{w}}\right)+A_{\mathrm{w}}{\displaystyle \frac{{\lambda }_{\mathrm{i}\mathrm{n}\mathrm{s}}}{{\delta }_{\mathrm{i}\mathrm{n}\mathrm{s}}}}\left(T_{\mathrm{a}\mathrm{m}\mathrm{b}}-T_{\mathrm{w}}\right)+{\dot{Q}}_{\mathrm{r}\mathrm{a}\mathrm{d},\mathrm{w}}\right),$$

(28)

where ${\dot{Q}}_{\mathrm{r}\mathrm{a}\mathrm{d},\mathrm{w}}$ is radiative heat flow rate at wall surface, ${\lambda }_{\mathrm{i}\mathrm{n}\mathrm{s}}$ is thermal conductivity of the vacuum chamber thermal insulation and ${\delta }_{\mathrm{i}\mathrm{n}\mathrm{s}}$ is thermal insulation thickness.

2.3.1. Convective Heat Transfer Coefficients Calculation

It would be wrong to assume that convective heat transfer coefficients don’t change during the vacuum drying process. Rarified atmosphere at low pressures causes convective heat transfer to diminish. To calculate the heat transfer coefficients at low pressures, a model from Gonzales et al. [29] which gives relation between Nusselt and Rayleigh number for natural convection along a vertical plate is used:

$$Nu=0.737\left[{Ra}^{1/4}+5.725·{Ra}^{0.019}\right],$$

(29)

where $Nu$ and $Ra$ are dimensionless quantities defined as:

$$Nu={\displaystyle \frac{h_j{H}_j}{{\lambda }_{\mathrm{g}}}},$$

(30)

$$Ra=Pr·Gr={\displaystyle \frac{{\mu }_{\mathrm{g}}{c}_{p,\mathrm{g}}}{{\lambda }_{\mathrm{g}}}}\left({\rho }_{\mathrm{g}}-{\rho }_{\mathrm{g},j}\right){\rho }_{\mathrm{g},j}{\displaystyle \frac{g{H}_j^3}{{\mu }_{\mathrm{g},j}^2}}.$$

(31)

where $H_j$ is a height of the j-th surface. Every surface is considered to behave like vertical plate, even the transformers (vertical cylinders) since their height is much larger than outer diameter. The properties with j in index are evaluated for the temperature of the j-th surface.

Equation (29) is valid for Rayleigh numbers in the range from 10^-2 to 10⁵, a condition that is satisfied throughout this type of the drying process.

2.3.2. Radiative Heat Flow Rates Calculation

Under the assumption that all surfaces in the system are grey and diffuse, balance of the radiative het flow rates can be written [30]:

$$\sum _{j=1}^{N_{\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{f}}}\left[\left({\displaystyle \frac{{\delta }_{kj}}{{\epsilon }_j}}-F_{k-j}{\displaystyle \frac{1-{\epsilon }_j}{{\epsilon }_j}}\right){\displaystyle \frac{{\dot{Q}}_{\mathrm{r}\mathrm{a}\mathrm{d},j}}{A_j}}\right]=\sum _{j=1}^{N_{\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{f}}}\left[F_{k-j}\sigma \left(T_k^4-T_j^4\right)\right], k=\mathrm{1,2},\dots ,N_{\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{f}}$$

(32)

where $N_{\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{f}}$ is the number of surfaces that exchange heat by radiation, ${\epsilon }_j$ is total emissivity of the j-th surface, $F_{k-j}$ is the view factor from surface k to surface j, $A_j$ is the surface area of the j-th surface and ${\delta }_{kj}$ is Kronecker delta operator that is equal to unity when indices k and j are the same and zero otherwise.

Equation (32) written down for all the surfaces results in linear system of equations which is then solved for radiative heat flow vector ${\dot{Q}}_{\mathrm{r}\mathrm{a}\mathrm{d},j}$.

2.4. The Interface Between Transformer Cellulose Insulation and Vacuum Chamber Atmosphere

To close the system of the governing equations, the interface between transformer cellulose insulation and chamber atmosphere should be properly modelled.

It is assumed that the partial pressures of water vapor and dry air at the interface are the same as in the chamber atmosphere. Therefore, for all transformers, the following equations hold:

$${\left.p_{\mathrm{v},i}\right|}_{r=r_{\mathrm{o}\mathrm{u}\mathrm{t}}}=p_{\mathrm{v},\mathrm{c}\mathrm{h}},$$

(33)

$${\left.p_{\mathrm{a},i}\right|}_{r=r_{\mathrm{o}\mathrm{u}\mathrm{t}}}=p_{\mathrm{a},\mathrm{c}\mathrm{h}},$$

(34)

however, since every transformer can have different temperature field, the following Dirichlet boundary conditions are used for i-th transformer:

$${\left.C_{\mathrm{v},i}\right|}_{r=r_{\mathrm{o}\mathrm{u}\mathrm{t}}}={\displaystyle \frac{C_{\mathrm{v},\mathrm{c}\mathrm{h}}{R}_{\mathrm{m}}{T}_{\mathrm{c}\mathrm{h}}}{R_{\mathrm{m}}{\left.T_i\right|}_{r=r_{\mathrm{o}\mathrm{u}\mathrm{t}}}}},$$

(35)

$${\left.C_{\mathrm{a},i}\right|}_{r=r_{\mathrm{o}\mathrm{u}\mathrm{t}}}={\displaystyle \frac{C_{\mathrm{a},\mathrm{c}\mathrm{h}}{R}_{\mathrm{m}}{T}_{\mathrm{c}\mathrm{h}}}{R_{\mathrm{m}}{\left.T_i\right|}_{r=r_{\mathrm{o}\mathrm{u}\mathrm{t}}}}}.$$

(36)

From equations (33) and (34) it consequently follows that the total pressure at the interface is equal to the pressure in the vacuum chamber:

$${\left.p_i\right|}_{r=r_{\mathrm{o}\mathrm{u}\mathrm{t}}}=p_{\mathrm{c}\mathrm{h}},$$

(37)

On the inner radius, zero molar flux boundary conditions are imposed:

$${\left.\left(v_i{C}_{\mathrm{v},i}-C_{\mathrm{g},i}{D}_{\mathrm{e}\mathrm{f},i}{\displaystyle \frac{d{y}_{\mathrm{v},i}}{dr}}\right)\right|}_{r=r_{\mathrm{i}\mathrm{n}}}=0,$$

(38)

$${\left.\left(v_i{C}_{\mathrm{a},i}+C_{\mathrm{g},i}{D}_{\mathrm{e}\mathrm{f},i}{\displaystyle \frac{d{y}_{\mathrm{v},i}}{dr}}\right)\right|}_{r=r_{\mathrm{i}\mathrm{n}}}=0,$$

(39)

The boundary condition for energy equation (14) is a Robin-type and includes both convective and radiative heat fluxes:

$$-{\lambda }_{\mathrm{e}\mathrm{f}}{\left.{\displaystyle \frac{\partial T_i}{\partial r}}\right|}_{r=r_{\mathrm{o}\mathrm{u}\mathrm{t}}}=h_{\mathrm{t},i}\left({\left.T_i\right|}_{r=r_{\mathrm{o}\mathrm{u}\mathrm{t}}}-T_{\mathrm{c}\mathrm{h}}\right)-{\displaystyle \frac{{\dot{Q}}_{\mathrm{r}\mathrm{a}\mathrm{d},\mathrm{t},i}}{A_{\mathrm{t}}}},$$

(40)

where ${\dot{Q}}_{\mathrm{r}\mathrm{a}\mathrm{d},\mathrm{t},i}$ is radiative heat flow rate at the i-th transformer surface which can either be transfer from the surface (negative) or onto the surface (positive).

On the inner radius, zero heat flux is imposed:

$$-{\lambda }_{\mathrm{e}\mathrm{f}}{\left.{\displaystyle \frac{\partial T_i}{\partial r}}\right|}_{r=r_{\mathrm{i}\mathrm{n}}}=0,$$

(41)

For known fields of molar densities and temperature, molar fluxes of water vapor and air at the interface between the transformer cellulose insulation and vacuum chamber atmosphere can be calculated:

$$J_{\mathrm{v},t, i}={\left.\left(v_i{C}_{\mathrm{v},i}-C_{\mathrm{g},i}{D}_{\mathrm{e}\mathrm{f},i}{\displaystyle \frac{d{y}_{\mathrm{v},i}}{dr}}\right)\right|}_{r=r_{\mathrm{o}\mathrm{u}\mathrm{t}}},$$

(42)

$$J_{\mathrm{a},t, i}={\left.\left(v_i{C}_{\mathrm{a},i}+C_{\mathrm{g},i}{D}_{\mathrm{e}\mathrm{f},i}{\displaystyle \frac{d{y}_{\mathrm{v},i}}{dr}}\right)\right|}_{r=r_{\mathrm{o}\mathrm{u}\mathrm{t}}}.$$

(43)

These fluxes are then substituted into the governing equations (17), (18) and (26).

2.5 Initial Conditions

All the components of the vacuum drying system as well as cellulose insulation are assumed to be at ambient temperature. The gaseous phase inside the pores of the cellulose insulation and atmosphere of the vacuum chamber are under ambient pressure and ambient relative humidity. In that case, the initial moisture content of the cellulose insulation is equilibrium moisture content for ambient air. The initial conditions are as follows:

$$T_0=293.15 \mathrm{K},$$

(44)

$$p_0=100000 \mathrm{P}\mathrm{a},$$

(45)

$$a_{\mathrm{w},0}=0.6,$$

(46)

$$C_{\mathrm{v},0}={\displaystyle \frac{a_{\mathrm{w},0}{p}_{\mathrm{s}\mathrm{a}\mathrm{t}}\left(T_0\right)}{R_{\mathrm{m}}{T}_0}}\approx 5.7566·{10}^{-4} {\displaystyle \frac{\mathrm{k}\mathrm{m}\mathrm{o}\mathrm{l}}{{\mathrm{m}}^3}}$$

(47)

$$C_{\mathrm{a},0}={\displaystyle \frac{p_0}{R_{\mathrm{m}}{T}_0}}-C_{\mathrm{v},0}\approx 4.045223·{10}^{-2} {\displaystyle \frac{\mathrm{k}\mathrm{m}\mathrm{o}\mathrm{l}}{{\mathrm{m}}^3}},$$

(48)

$$X_0=X_{\mathrm{e}\mathrm{q}}\left(C_{\mathrm{v},0},T_0\right)\approx 0.081465=8.1465 \%$$

(49)

2.5 Numerical Solution Procedure

The governing equations (1), (2), (10) and (14) at the scale of cellulose insulation and governing equations (17), (18), (26), (27) and (28) at the scale of vacuum chamber, along with all the boundary conditions and constitutive relations form a nonlinear coupled PDE-ODE system.

One way to decouple the system and find the solution would be to solve the ODE system with guessed values of cellulose insulation fields. The ODE solutions can then be used to define necessary boundary conditions for which PDE system is solved and the better guess for cellulose insulation fields is obtained. This is the iterative procedure that is repeated until convergence is achieved. In this study, however, solving a coupled PDE-ODE system was addressed using an alternative approach.

Firstly, all the equations were scaled so that variables of the model have order of magnitude around unity or less. Then spatial derivatives in governing equations ((1), (2), (14)) and in interface equations ((38), (39), (40), (41), (42), (43)) were discretized using second order finite difference scheme. This leads to a large system of ODEs which is then solved using the LSODA method implemented in in the SciPy package for Python programming language [31]. This method automatically switches between Adams (non-stiff part of solution) and implicit BDF (stiff part of the solution) methods. Relative and absolute tolerances were set to 10^-5 and 10^-7 respectively. A similar methodology for solving PDE-ODE systems was investigated by Filipov et al. [32] so readers interested in further details are referred to their work.

3. Results and Discussion

Numerical simulations for all the cases studied below were run with the parameter values and relations listed in the Table A1, Table A2, Table A3 and Table A4 given in the Appendix section. The exceptions are two parameters: absolute permeability $k_0$ and internal overall mass transfer coefficient $K$ for which series of values were used to study their influence on the moisture dynamics.

The vacuum dryer dimensions correspond to those of a laboratory-scale vacuum chamber that is used at our institution. Likewise, the choice of vacuum pump and its characteristic pumping speed correspond to a commercially available unit (VD 65, Leybold) currently in use.

Apart from the parameter values and relations, to solve the system of governing equations, the values of the view factors need to be calculated beforehand. This was done firstly by generating rudimentary geometry of the surfaces that exchange heat by radiation, namely, the IR heaters, transformer samples and vacuum chamber walls, as shown on Figure 4. The generated mesh is a pre-requisite for computation of view factors; however, it is not necessary to calculate them for all cells since reciprocity and summation rules can be applied [30].

Furthermore, due to the spatial layout of the transformer samples, three representative ones, corner, side and center transformer sample, as pointed out on the Figure 4 can be identified out of the total of nine. Since these samples are exposed to a distinct radiative conditions, it is sufficient to compute view factors for only these three and thus significantly save computational resources.

The calculation of view factors was done with pyviewfactor [34], a Python package that utilizes PyVista [33] functionalities, such as ray tracing that was used for checking the shadowing of the cells. The calculated view factors subsequently used in equation (32) are listed in Table 1.

The drying process described by this model can be controlled by on/off switching of three main components:

• IR heaters (when off, ${\dot{Q}}_{\mathrm{e}\mathrm{l}}=0$ W)
• Vacuum pump (when off, $C_{\mathrm{s}}=0$ m³/s)
• Venting/pressure relief valve (when off, $K_{\mathrm{v}}=0$)

Depending on the component type and timing of the switching, the process can discretely by controlled in numerous ways. In this work, test case and drying case were studied.

3.1. Test Case Results

Firstly, to analyze the results in terms of spatio-temporal physical fields, a short-duration test case was simulated. The main component states and their corresponding durations are listed in Table 2 and as such provide enough information to define the test case process. The states of the components are defined in such way that all the commonly encountered conditions, such as IR heating, vacuum pump-down and venting are encompassed in this numerical simulation. The test case therefore consists of 2.5-hour IR heating period with the pressure relief valve open to the atmosphere, followed by 2.5 hours of vacuum pump-down with heaters turned off, ending with a 2.5-hour venting period during which the heaters remain off.

As heaters are activated, the temperature of cellulose insulation increases. This reduces the equilibrium moisture content of the cellulose fibers according to equation (12). Since cellulose fibers tend to reach a lower moisture content state, desorption of moisture occurs consequently making the pores of the cellulose insulation more saturated with water vapor, as seen on Figure 5a. The developed gradient of the water vapor molar fraction acts as a driving force for equimolar diffusion of dry air and water vapor within the pores. The diffusive flux during the heat-up period, however, is not unidirectional; it changes direction at certain radius which is the consequence of water vapor molar ratio maxima at that point. This is not surprising; since heat is transferred to the outer surface only, the molar density of water vapor is highest near the surface and lower at the very interface. A similar thing happens with the pressure field, as seen on Figure 6b. On the same plots, an under-pressure region can also be identified near the internal surface. This is a consequence of the combined effect of dry air diffusing outwards and absorption of water vapor that diffuses inwards, toward the parts that are still unheated.

The temporal variations of physical fields for corner and side transformer samples are nearly identical, which is expected given the relatively similar values of the view factors from IR heaters to these samples. In contrast, the center transformer sample is completely obstructed by the surrounding samples, so the radiative flux reaching its surface originates from neighboring samples and not directly from the heaters. This manifests as a slower increase in temperature (Figure 5c) and slower decrease of moisture content (Figure 5d). Moreover, the moisture content of the center transformer even shows a slight increase during the heating period near the outer surface which has to do with the variation of molar density in the vacuum chamber atmosphere.

Water vapor is convectively and diffusively transferred at the surface of corner and side transformer samples to the vacuum chamber atmosphere according to the interface equation (42). Furthermore, both dry air and water vapor exit the chamber through the venting/over pressure valve due to the rising pressure. This results in an increase in water vapor and a decrease in dry air molar density as shown in Figure 7a. The elevated water vapor molar fraction and lower surface temperature of the center transformer sample correspond to higher equilibrium moisture content causing a slight absorption of water vapor at the surface of the center transformer sample.

When the vacuum pump down period starts, dry air molar density within the pores rapidly drops (Figure 5b) and the similar decrease is observed in the vacuum chamber (Figure 7a). In contrast, water vapor molar density remains nonzero due to continued desorption from the cellulose fibers (Figure 5a). As a result, pressure gradient develops along thickness of the cellulose insulation (Figure 6b) which gives rise to a more intensive convective transport and consequently to an increase in the drying rate. The diffusive flux in this period is zero since water vapor molar fraction is uniform and equal to unity (Figure 6a).

During the vacuum chamber venting period, dry air diffuses into the cellulose insulation leading to an increase in its molar density. The elevated dry air concentration combined with a still relatively high temperature causes desorption of moisture in the inner layers and in consequence an increase in pressure shortly after opening the venting valve.

Figure 7b shows temporal variations of surface temperatures inside the vacuum chamber. Among the transformer sample surfaces, the side transformer sample surface exhibits the most rapid increase in temperature that is attributed to the highest view factor from the heaters (Table 1). A similar trend is observed for the corner sample surface, though its temperature is slightly lower due to a lower value of view factor. On the other hand, a center transformer sample surface shows significantly slower increase in temperature and additionally maintains lower temperatures throughout the entire heating period. This behavior is clearly a result of the surrounding samples obstructing the center one from direct exposure to the heaters.

These results suggest that the surface temperature of the sample most exposed to the IR heaters (in this case, side transformer sample) is a suitable candidate for use as a control variable in terms of temperature control. The choice of highest temperature is due to the upper temperature limit imposed to prevent rapid degradation of the cellulose insulation.

The temperature of the chamber atmosphere initially increases mostly due to convective heat transfer from the surfaces within the chamber. A sudden drop is observed when the pump down period starts, which can be attributed to the expansion of the gaseous phase. The opposite effect occurs when the venting valve opens. As ambient air rushes inside the chamber it compresses the gaseous phase, resulting in sharp increase in temperature. Lastly, vacuum chamber walls exhibit the slowest temperature increase due to their relatively high heat capacity.

3.2. Impact of the $k_0$ and $K$ Parameters

The parameters whose impact is investigated are absolute permeability $k_0$ and internal overall mass transfer coefficient $K$. All other parameters are either held constant or can be calculated from the relations, both of which are provided in the Table A1 in the Appendix. Moreover, the states of the main components for all simulations are the same as in the test case (Table 2). Lastly, given the small difference in physical fields between the corner and side transformer samples, only the center and side samples are observed in this analysis.

To evaluate the influence of the above-mentioned parameters, an average moisture content of the cellulose insulation was computed using the following equation:

$$X_{\mathrm{a}\mathrm{v}\mathrm{g}}\left(t\right)={\displaystyle \frac{1}{r_{\mathrm{o}\mathrm{u}\mathrm{t}}-r_{\mathrm{i}\mathrm{n}}}}\underset{r_{\mathrm{i}\mathrm{n}}}{\overset{r_{\mathrm{o}\mathrm{u}\mathrm{t}}}{\int }}X\left(r,t\right)dr.$$

(50)

Figure 8a illustrates the influence of the internal overall mass transfer coefficient $K$ on the evolution of average moisture content during the test case process. As $K$ increases, the characteristic time for equilibration decreases, causing local absorption and desorption processes to become significantly faster than diffusive and convective transport. In that case, local equilibrium can be assumed to occur instantaneously, allowing the moisture content field to be computed directly using equation (12) for given fields of water vapor molar density and temperature. In contrast, lower values of $K$ lead to more inert moisture dynamics characterized by reduced drying rates.

The influence of absolute permeability $k_0$ is illustrated in Figure 8b. As shown, higher permeability results in an increased drying rate during the vacuum period, which is a direct consequence of enhanced convective transport, consistent with Darcy’s law (equation (8)). Conversely, lower value of $k_0$ result in lower drying rate in vacuum period. It is also observed that value of $k_0$ has a negligible effect on the drying rate during the heating period, at least while pressure gradients remain low. However, as pressure gradients develop, the convective transport becomes more prominent, thereby increasing the sensitivity of the results to the $k_0$ value.

Both $k_0$ and $K$ parameters significantly affect the drying rate, yet their values are not easily determined through theoretical means as both are highly dependent on the microstructure of the insulation. The most time-efficient approach would therefore be to estimate these parameters through model calibration, by comparing the simulation results with experimental data.

The Knudsen diffusivity parameter $D_{\mathrm{K}}$ was also examined as a part of this analysis and it was found to have no significant influence on the drying rate. This observation is consistent with the model formulation. Specifically, equation (6) which indicates that the contribution of $D_{\mathrm{K}}$ becomes relevant at low pressures. However, under vacuum conditions, the molar fraction of water vapor approaches unity and becomes uniform across the domain, eliminating concentration gradients and, consequently, diffusive transport.

3.3. Drying Case Results

The drying case was defined to resemble a process with conditions typically met in transformer manufacturing procedure, consisting of a sequence of alternating heating and vacuum phases [35,36]. To maintain the specified temperature and pressure levels throughout the process, a temperature and pressure two-point (on/off) control was implemented in the simulation. The operational states of the main components and their respective durations are summarized in Table 3.

To ensure that the cellulose insulation does not exceed the upper temperature limit of 130 °C, the controlled temperature was chosen to be the surface temperature of the side transformer sample, as explained in the previous section. Figure 8 presents the temporal variations of the surface temperatures and the innermost layer temperatures of the corresponding transformer samples. The same figure also displays the variation in pressure of the chamber atmosphere over time. Significant temperature gradients across the thickness of the cellulose insulation are evident. Their negative signature is characteristic of one-sided heating and is unfavorable in terms of drying rate [10].

Furthermore, the obstructed center sample remains at somewhat lower temperature throughout the entire drying process, a condition that is expected to extend the overall drying time. The observed decrease in temperature of the transformer samples during the vacuum phases results from radiative cooling due to the turned off IR heaters as well as cooling due to the water vapor evaporation. This is evident from Figure 9, where each successive vacuum phase is characterized by a smaller drop in surface temperatures, indicating a decline in evaporative cooling as the drying rate diminishes over time.

Convective and radiative heat fluxes at the surface of side and center transformer samples are illustrated in Figure 10. The results demonstrate that convective heat flux diminishes during the vacuum phases which is consistent with the reduction of the convective heat transfer coefficient at low Rayleigh numbers, as described by equation (29). For all transformer samples, convective heat flux remains generally negative due to the lower temperature of the vacuum chamber atmosphere relative to the sample surfaces. Peaks in convective heat flux appearing at the onset of each vacuum phase result from the rapid cooling of the vacuum chamber atmosphere caused by its expansion during pump-down. On the other hand, radiative heat flux is characterized by periodic oscillations which arise from the on/off temperature control achieved by changing the operating states of the heaters and thus the emitted radiative flux.

Figure 11a illustrates the molar flux at the surface of side and center transformer throughout the drying process. Initially, a diffusion dominates, with diffusive flux reaching a peak before gradually decreasing. This behavior results from the rising water vapor molar fraction in the chamber, which eventually causes the diffusive flux to vanish entirely. Not long after the start of the process, as insulation continues to heat up and water vapor gets desorbed from the cellulose fibers into the pore space, pressure gradients begin to develop across the insulation thickness, giving rise to convective molar transport. The convective flux reaches its peak after the diffusive flux and then declines as the pressure field gets more uniform.

To extract the water vapor from the pore space of the insulation, at some point, chamber pressure must be reduced to regain an increase in convective molar flux. As these gradients again start to diminish, the chamber is vented and temperature control is reactivated. At the beginning of the next heating phase, the molar fraction of the water vapor is no longer uniformly distributed since ambient air is introduced which then results in subsequent increase in diffusive flux.

Periods of increased molar flux correspond to a higher drying rate intervals as shown in the drying rate time evolution plot for side and center transformer sample in Figure 11b. For both transformer samples, each successive phase exhibits a lower drying rate values. Interestingly, though, the center transformer sample showing a lower drying rate at the beginning of the process eventually surpasses the side sample. This behavior can be attributed to its higher moisture content throughout the process which makes it more responsive to subsequent heating and vacuum phases.

The drying curves of all three transformer samples are shown on Figure 11c. The side and corner samples, both exposed to similar radiative conditions, exhibit nearly identical drying behavior. In contrast, the center transformer sample initially experiences a slight increase in moisture content. This is attributed to a rapid increase in water vapor molar fraction of the vacuum chamber atmosphere (primarily due to moisture desorption from the surrounding samples) and comparatively slow heating of the center sample. As the center sample heats up, its drying rate changes the signature, however, its moisture content remains elevated throughout the entire process. This observation clearly suggests that the overall drying duration is dictated by the moisture dynamics of the radiatively obstructed regions and cold spots within the cellulose insulation. These challenges are among the reasons why alternative drying technologies, such as low-frequency heating and vapor phase drying [35], are often employed in industrial practice.

The drying case results presented should be interpreted qualitatively, as the underlying model has not been calibrated with respect to key parameter values. Nevertheless, the simulation results suggest that by modifying the controlled values, as well as the timing and duration of heating and vacuum phases, it is possible to influence the drying rate trend and, consequently, the total drying time. As the model accounts for the process at the vacuum chamber scale, it also enables estimation of the total energy consumption of the IR heaters. This allows for process optimization with respect to drying time, energy consumption or a trade-off between the two.

4. Conclusions

This study has presented a detailed mathematical model of the vacuum drying process for transformer cellulose-based insulation accounting for coupled heat and mass transfer phenomena at both the cellulose insulation scale and the vacuum chamber scale. Although the simulations were carried out for a laboratory-scale dryer with a simplified geometry, implementing real dimensions and geometry is largely a matter of parameter adjustment. The presented model could also be adapted to alternative heating methods, such as hot air circulation, by modifying the constitutive relations used to compute the convective heat transfer coefficient or drying with internal heating, by modifying the boundary conditions at the cellulose insulation interface.

By incorporating radiative heat transfer along with view factor computation, the model results have highlighted the influence of the spatial layout of transformer samples, an important practical consideration when setting up the drying process.

Unlike previous studies discussed in the introduction section, the results of this model demonstrate the important role of combining heating and vacuum phases in the drying procedure. While this strategy is already common in transformer manufacturing industry, the present study provides clear physical explanation for its effectiveness.

The model, with its dual-scale formulation, offers a valuable framework for optimizing the drying process with respect to both drying time and energy consumption. However, a major limitation of this model lies in the uncertainty of several cellulose insulation properties that have substantial effect on moisture dynamics. Such are absolute permeability and internal overall mass transfer coefficient as well as sorption isotherm data. These limitations can be addressed through experimental calibration of the model; an essential step planned for our future work.