Mathematical Modeling and Numerical Simulation of Coupled Mass Transfer and Mechanical Processes in Ceramic Ware at Industrial Drying

An algorithm for numerical simulation of transient moisture content fields and mechanical processes in ceramic ware at drying in industrial aggregates is developed. It is based on mathematical models of the mass transfer and mechanical behavior in the ceramic bodies, data for the drying regime and physical properties of the material as function of water content. The models allow variations of the drying conditions in order to choice the most efficient regime at existing or design dryers. The algorithm is applied for a direct coupled finite element analysis of wet bricks behavior in continuous working drying installation. The shrinkage mode, modulus of elasticity, Poisson ratio, modulus of rupture, effective mass transfer coefficient and critical moisture content are determined by experimental tests of the material. They are used to simulate numerically three-dimensional moisture, stress and strain fields in ceramic bodies at the existing drying regime. Ways for improvement of the models and their application for estimation of the potential for energy savings in industrial dryers are discussed.


Introduction
Industrial drying is a key process in the ceramic manufacturing. The right organization of the drying regime is important for the quality of the products and the associated energy inputs [1]. The advanced ceramic industry is continually evolving in а direction of optimization of geometry and materials of the products in order to increase their functionality. These processes provoke research activities oriented to improving of the energy and technological efficiency of the industrial dryers.
Solutions of such actual tasks are possible by experimental tests in laboratory drying channels or mathematical modelling and numerical analyses of the coupled thermo-structural-diffusion processes in the dried ceramic bodies. The experiments require special equipment and usually are energy and time consuming. They are successfully displaced by modelling investigations in the term of rapid development of the computer technics and software for numerical simulations [2]. Many research efforts are focused to understand, model and predict the complex heat, mass transfer, plastic deformations and shrinkage of the ceramic material in order to prevent the possible failure of the articles [3,4,5]. Such couple field analysis can also be applied for establishment of the initial geometry of the articles that is changed at drying and firing to the final design [6].
The aim of the present study is to develop an algorithm for numerical simulation of the coupled structural-diffusion behavior of ceramic ware at industrial convective drying, suitable for analysis and improvement of the efficiency of the process. Such complex investigations, taking into account the real conditions in industrial dryers have not been reported to now. respectively the drying rate drop/rise (A-B). The change of the moisture content is insignificant and that relatively short period can be examined as a part of the next one [7].
The CDRP begins at temperature of the material, equal to the wet bulb temperature (point B and L in Figure 1). The surface temperature remains constant during the period (line L-M). The CDRP is fully governed by the rates of external heat and mass transfer since a film of free water is always available at the evaporating surface. The moisture migrates from the inside of the porous material to the surface by diffusion [4]. The drying rate remains constant as long as the moisture transport rate from the interior of the material to the exchanging surface is equal to the evaporation rate (B-C). The moisture content at the end of the CDRP is called critical moisture content (W c , point C). The evaporation rate starts to decrease at that point because of the disconnection of the liquid meniscus from the surface. This decreasing trend continues at FDRP until all liquid meniscuses are disconnected. The internal moisture transport rate specifies the drying rate of the FDRP. The temperature of material surpasses the wet bulb temperature and grows up to the drying media temperature (M-N). The drying process stops when the equilibrium moisture content W eq is reached (D).
Several types of strains are expected to exist during the drying the wet body, formed by porous plastic material. They are provoked by the mass transfer (shrinkage), temperature gradients (thermal stresses) and the gravity [4]. Shrinkage is one of the most important factor effecting on drying behavior of highly shrinkable porous media such as clay-like materials. It occurs during the 0 and I period and stops at the critical moisture content.
During the drying process, non-uniform shrinkage results in stresses that may causes cracking, influencing on the quality of dried products. If drying is performed sufficiently slowly, the shrinkage would be closer to uniform and drying induced stresses are not dangerous for the material. By other hand, the increasing of the drying time results in higher embodied energy and lower dryer productivity. For economic reasons the industrial drying processes are usually performed faster. That results in non-uniform shrinkage and subsequent high tensile stresses.
It is important to choose suitable conditions of drying for the different products to prevent failures of the dried articles. The time durations of the CDRP and FDRP (τ 1 and τ 2 on Figure 1) are between the most important factors determining the efficiency of drying. They depend on the potential of drying, expressed by the difference between the dry and wet bulb temperatures of the drying media. As smaller is it as longer is the drying time. To prevent failure of the production due to non-uniform shrinkage, the CDRP is maintain at a relatively small drying potential via keeping of high relative humidity and velocity of the drying gas during that period. The boundary conditions of the models below allow taking into account the influence of these fluid flow parameters on the coupled mass transfer and mechanical processes in the dried solid domain.

Modeling of Structural-Diffusion Processes in the
Wet Ceramic Ware at Drying

System of Equations
The transient fields of the moisture content and subsequent mechanical behavior of the materials are obtained by coupled numerical solution of mass transfer equation (1) and stress-strain relationship (2) for 3D finite element mesh, approximating the geometry of the ceramic body.
In a coupled structural -diffusion analysis the total strain is formed of elastic { } and diffusion { } parts: where { } = vector of coefficients of diffusion expansion, m 3 kg -1 ; ∆C = concentration change according reference value C ref : The modulus of elasticity E, Poisson ratio, coefficient of diffusion expansion β and D ef are determined experimentally and used in the models as function of moisture content in the dried material [10,11].

Boundary Conditions
The mass transfer from the boundaries S of the ceramic body to the drying gas can be computed by mass flux ̇, depending on the thermodynamically and fluid flow conditions in the dryer [9]: The intensity of drying in CDRP is expressed by: where ̇ = mass flux of evaporated water by unit surface, kgm -2 s -1 ; h mI = mass transfer coefficient at CDRP, ms -1 ; Т = temperature of the drying media, K; R w = 462 Jkg -1 K -1specific gas constant for water vapor; р w and p s are respectively partial pressures of unsaturated and saturated water vapor, Ра.
The coefficient h mI depends on the fluid velocity, relative humidity and temperature field in the drier and can be computed by Nusselt number at mass transfer: where l = length of the ceramic body on way of the drying gas; D d.m = diffusion coefficient of the drying media, m 2 s -1 .
At isobaric process, it is computed by the temperature of the drying media t, °C: The Nusselt number is given in the common form in the literature [7,8]: where A 1 , m 1 and n 1 depend on Reynolds number (Re). They are known for large limits of Re and can be calibrated on the base of operating conditions in the dryer. Re, Prandtl number at mass transfer and Guhman number, used as non-dimensional drying potential, correspondently are: where w = average velocity, ms -1 ; ν = kinematic viscosity of the drying media, m 2 s -1 ; T w =wet bulb temperature of the drying gas, K. The mass flux, evaporated from ceramic ware at FDRP can be obtained according hypothesis for linear decreasing of the drying rate [8]: where C c = water concentration at the critical moisture content, kgm -3 ; C eq =equilibrium water concentration, kgm -3 ; h mII = mass transfer coefficient in FDRP, ms -1 . The boundary conditions for the structural analysis include zero normal displacements on the contact surface between the transport and the dried body and on the symmetry surface (if there exist any).
The gravity is taken into account.

Analogy between the Heat and Mass Transfer
The advanced software for numerical simulation [2] allow structural -diffusion analysis at limited opportunities for definitions of boundary conditions and material properties for solution of equations (1) and (2). The modeling of the mass flux, mechanical properties and coefficient of mass diffusion as function of the concentration is not directly possible. An alternative approach is to implement the structural-diffusion simulation using the analogies between the mass and heat transfer and between the diffusion and thermal stresses.
The effective diffusion coefficient in (1) can be expressed by: where K art = artificial thermal conductivity, Wm -1 K -1 ; ρ d = density of the dry mass, kgm -3 ; c = specific heat capacity, Jkg -1 K -1 . The concentration can be referred by the dry mass of the material: where u = moisture content of the ceramic mass, kg The boundary conditions for the solution of (14) is obtained taking into account (5), (12) and (13): where: If T replaces W in (14) and (15), the equations are transformed to Fourier equation for heat transfer by conduction and the second kind of boundary conditions for its solution, respectively.
The diffusion strains can be expressed as thermal strains in (2): where α=instantaneous coefficient of thermal expansion, K -1 .

At an analogy between thermal and diffusion stresses α =β
where W in is the initial moisture content, %; l in = initial linear size, ∆l =reduction of the linear size at the shrinkage. The last are obtained experimentally at laboratory drying of test samples.

External Mass Exchange in the Dryer
The proposed models are tested to investigate the structural-diffusion behavior of building bricks in industrial continuously working dryer. Partial results are presented in [12]. The drying media is an air-flue gases mixture with predominant air fraction. The drying gas and the wet ceramic articles are moved in counter flow through the drying tunnel. The change of the fluid flow temperature, relative humidity and velocity are obtained by in situ measurements, material and thermal balances. Detail information about the operating conditions in the dryer is not possible due to confidential rules. Figures 2 and 3 illustrate the change of the non-dimensions numbers Re and Gu with the non-dimensional time and tunnel length. It is obvious that the time duration of the CDRP, corresponding to the smaller Gu number, is a about 55% of the full time duration.
The coefficients in Nu mI number are calibrated to obtain the mass transfer coefficient and mass flux according (6) in order to achieve the evaporated amount of water per article at CDRP according the experimental data: A 1 =0.49; n 1 =0.45 and m 1 =0.135. So obtained average mass flux at the CDRP ̇ is used in the boundary conditions (15). ̇ * is applied on all boundaries, excepting the symmetry ones and the contact surfaces between the ceramic ware to be dried and the floor of the transport in the dryer.

Material Properties of the Wet Mass
The wet mass to be dried is a mixture of clay and combustible supplements. Standard geometry samples (plates and cylinders) made by a material after the extruder in the plant were dried consecutively in a laboratory dryer to different moisture contents. The relative linear deformation at CDRP, critical moisture content and the variations of the modulus of elasticity, Poisson ratio, compressive and flexure strength (modulus of rupture) of the material with the water content were obtained at series of experiments. They are obvious from  The effective diffusion coefficients are computed based on the information, obtained at drying of plates, insolated at the smaller sides in order to imitate infinite plates [7]: where: u in = initial moisture content, u c = critical moisture content, � = average moisture content, kg water /kg dry mass; � = average drying rate, referred to the dry mass, kgkg -1 s -1 at CDRP. Г=3 for infinite plate with thickness 2R. R v is:

Numerical Simulation of Structural -Diffusion Fields in the Dried Bodies
Numerical simulation of the transient mass transfer and mechanical behavior in the domain of a typical article ( Figure 10) is implemented. The geometrical model includes a brick, divided by a symmetry plane. It is discretized by a finite element swept mesh with higher-order hexahedral element Solid 226 (20 nodes), available at ANSYS Mechanical APDL for thermalstructural-diffusion analysis (Figure 11). The nonlinear relations between the degrees of freedom and coordinates within Solid 226 allow precise analysis at a relatively coarse mesh. The computation procedure takes a relatively long CPU time (19 days at workstation with 128 GB of RAM and processor Intel Xeon E5).
Transient fields of the moisture content, stresses and strains in the body during the drying process are obtained.  The models of the mechanical processes are validated comparing the final sizes after the deformations, obtained at the simulation and the measured real ones. Differences below 4 % are established. That prove the adequacy of the mechanical models. Additional validation on the models of the mass transfer is performed by comparison of the mass of the evaporated water computed by the final average moisture content, established at the real conditions and the same, obtained at the simulation. The difference of 10 % between the computed and measured water mass is established. The ignoring of the surface reduction of the body at the shrinkage in the models of external mass exchange is accepted as one of the reasons for that difference. Second numerical simulation of the processes is implemented in order to check the influence of that ignoring on the results. The mass flow in (16) is divided into the square of the average relative linear reduction of the sizes at that step. The final average moisture content is smaller according the same at the first simulation and the difference between the computed and measured masses of evaporated water is decreased to 7 %.
The change of the volume averaged moisture content with the time is shown on Figure 12  Fields of moisture content, first principal stress and displacements at a moment of the CDRP with higher local stresses in the body are shown on Figures 14 and 15. The maximal values of the first principal stress does not exceed the modulus of rupture and the minimal values are higher than the compressive strength of the material. So fractures of the products due to non-uniform shrinkage are not expected at the simulated drying regime. But it proceeds at lower drying rate in comparison with the maintained one in the dryer due to the imperfection of the models of the mass transfer. So the information, obtained at that stage of finite elements analyses gives approximate idea of the real processes. It is expecting to be more accurate after corrections of the models of the external mass exchange. Possible ways for refining of model of the mass transfer are:  increasing of the mass flux ̇ at the assumption that CDRP is shorter than 55 % of the full time duration of the draying;  approximation of ̇ with nonlinear dependence on the moisture content at FGRP in equation (16).
These corrections of the models have to be made in order to achieve differences in the computed quantity of the evaporated water per article and the maintained in the factory under 1 %. So calibrated models can be used to simulate mass transfer and subsequent mechanical behavior in ceramic ware with different geometries at changes in the drying regime. A possibility for improvement of the efficiency of the investigated dryer is a reduction of the time duration of the process. This also will increase the productiveness of the installation. The time reduction can be realized by a higher speed of the ceramic ware trough the drying channel and additional corrections of the velocity, temperature and relative humidity variations with the time, consistent with the technical possibilities of the drying installation. The mechanical stresses, obtained as results of numerical simulations of the coupled fields, taking into account these changes in the models, have to be smaller than the experimentally obtained strength of the material at the correspondent moisture content (Figures 7 and 8) in order to avoid failures of the production. The energy effects at the changes of the drying regime can be estimate by energy balances of the installation.

Conclusions
An algorithm for numerical investigation of the structural-diffusion processes in ceramic ware at drying in industrial conditions is composed and tested on a real object. Detail information about the mass transfer and mechanical behavior in articles with complex geometry during convective drying is obtained. The comparison of the results with experimental data showed necessity of improvements of the models of the mass transfer. The corrected algorithm can be used for investigation of the possibilities for improvement of the efficiency of the dryer installations by reduction of the time duration of the process. That is possible if the local stresses, obtained at the numerical simulations at different drying rate don't exceed the strength of the material at the correspondent moisture content.