Format And Type Fonts CCHHEEMMIICCAALL EENNGGIINNEEEERRIINNGG TTRRAANNSSAACCTTIIOONNSS VOL. 45, 2015 A publication of The Italian Association of Chemical Engineering www.aidic.it/cet Guest Editors: Petar Sabev Varbanov, Jiří Jaromír Klemeš, Sharifah Rafidah Wan Alwi, Jun Yow Yong, Xia Liu Copyright © 2015, AIDIC Servizi S.r.l., ISBN 978-88-95608-36-5; ISSN 2283-9216 DOI:10.3303/CET1545235 Please cite this article as: Harun Z., Ong T.C., 2015, A complete mathematical model of convective drying in membrane fabrication, Chemical Engineering Transactions, 45, 1405-1410 DOI:10.3303/CET1545235 1405 A Complete Mathematical Model of Convective Drying in Membrane Fabrication Zawati Harun, Tze Ching Ong* Faculty of Mechanical and Manufacturing Engineering, University Tun Hussien Onn,Parit Raja, 86400, Batu Pahat, Johor, Malaysia alex_ongtc@yahoo.com This article presents a fully coupled heat and mass transfer model for convective drying of ceramic porous material. A two dimensional mathematical model that coupled mass, heat and gas transfer was developed and integrated with finite element method and the backward difference method to solve the spatial discretization and temporal discretization of the governing equations respectively. The numerical was computed using Skyline solver to capture highly nonlinear and transient process. Two different meshes (coarse and fine) were applied for the rectangular domain as a representative of the drying object. Results show two typical periods during drying namely the constant rate period (CRP) and falling rate period (FRP) which demonstrated a very small difference between both meshes. However, computational time was increased significantly for the finer mesh as compared to coarser mesh due to their high computational complexity. Thus, we can deduce that this model is sufficient to describe the mechanisms controlling the drying process in membrane fabrication. 1. Introduction Membrane drying is a process that involves removing moisture or solvent during precursor formation before the sintering process. Sometimes it is also known as the most energy intensive processes in any industrial as it involves a lot of variables that change concurrently especially in the mixed pore phase condition (Mujumdar, 2006). The losses due to convective drying defects alone can be up to the order of 10 % in industrial production (Tarantino et al., 2010). Hence, optimization of the drying process in terms of energy usage and production time, without compromising the defects of the product quality such as cracking or warping required a comprehensive study (Gariev et al., 2014) and better understanding of controlling transport mechanisms that evolve during the drying process. Generally, drying means initially wet body evaporates from the drying surface at a constant rate corresponding to liquid flux (Perré et al., 2007). This stage of drying is denoted a constant rate period (CRP) (Scherer, 1990). As drying progress, the fraction of the wet area at the surface decreases with decreasing free water concentration. When the water content reaches critical free water content, the surface evaporation will form discontinue wet patches. This period of drying is denoted as the falling rate period (FRP) (Harun and Ong, 2014a). As the drying surface recedes further into the interior of solid matrix, intense drying emerges the formation of a wet region and a dry region with vapour flux substitute the previous liquid flux as a dominant driving force for moisture transfer. Finally, drying ceases as water content decreases to zero and subsequently vapour flux decays toward atmosphere condition as internal vapour pressure and external vapour pressure is almost equilibrium. Further complications arise if drying involves materials in a multilayer structure such as membrane. Ceramic membrane layers have a top layer with hygroscopic material behaviour which bound water is strongly attached to the capillary wall (Harun et al., 2014a). Meanwhile, the bottom layers consist of non- hygroscopic material characteristic where moisture removal is easier due to lesser water retention (Harun et al., 2014c). Typically, drying of the membrane layering structure not just involve the compatible issue of drying parameter but the different materials that acts with different properties also strongly influence the consistency of shrinkage geometry which is often associated to failure of membranes due to the improper 1406 of drying and sintering processes. The capability of the unique modeling approach based on fundamental heat and mass transfer relationships and thermodynamic equilibrium on the desired drying technique and specified material of interest will provides a tool in predicting and exploring the drying behaviour as well as understanding the impact of important variables towards the drying process (Harun and Ong, 2014b). In this way, the impact of the many variables on the drying behaviour of a specific material of interest and applied drying technique can be examined and interpreted without having an intensive experimental task. In fact the kinetic and dynamic during drying involve a very complex mechanism and investigation via measurement is difficult to apply. Therefore, modeling approach seems to be the only solution to eliminate all those unnecessary troubles. Hence, understanding of various stages and transport parameters that evolve dynamically with time throughout the drying is extremely important to ensure a good quality of the final products. Normally, mathematical modeling results a set of governing equations and numerically solving them using a specialized solver and an adaptive mesh scheme. Thus a domain meshing sensitivity towards the model in term of mesh dependency to give accurate results also needs to be checked first. Thus this work is devoted to investigate this problem and saturation variable during drying using a mathematical model with coupled mass, heat and gas transfer in two dimensional porous materials for convective drying process. This case study is an extended work from Harun et al. (2014b). The significant differences between two mesh used, namely a coarse and fine mesh are the details of elements and nodes layout in a domain as a representative of dried object. Fine mesh contains more elements and nodes respectively when compared to coarse mesh. The process of derivation of the transport governing equations and numerical process would only be discussed briefly in the next section. 2. Mathematical modeling 2.1 Theoretical formulation for phase transport Unsaturated membrane is a three-phase system namely liquid, solid membrane particles and gas. In this formulation, the liquid phase is considered to be water with dissolved vapour and air, whilst, binary mixture of vapour and dry air was considered a gas phase. Capillarity for liquid phases, diffusion by vapour and air, gas transport and also heat transport are the transport phenomena integrated in the model with the governing equations were expressed into three primary variables; water pressure, Pl, temperature, T and gas pressure, Pg. The liquid mass conservation equation is formulated as Eq(1). )(.)(. )(.)(. )()( b V lg VvvVvl V l t gSv t l S l         (1) The water velocity, V1 and gas velocity, Vg can be easily derived from Darcy’s law. Meanwhile, the vapour velocity by diffusion is defined in Kanno et al. (1996). For a hygroscopic material, when the water content reaches the maximum irreducible level the bound water flux is taken into account. The bound water velocity was derived from Zhang et al. (1999). Full derivation of those velocity accordingly to the measured variables can be found in subsequent work by Harun et al. (2014a). By applying conservation of mass for the flow of dry air within the pores of the material body dictates that the time derivative of the dry air content is equal to the spatial derivative of the dry air flux and is written as Eq(2). )( )( vVvρgVaρ t gSaρ     (2) The heat transfer formulation employed considers the effect of conduction, latent heat and convection and therefore given as Eq(3). )iVa,l,vi pici)ρr(T-T()L-gVvρvVv(ρT)λ t a,l,vi iciρis s ρp)c1-      ( )((  (3) 2.2 Thermodynamic relationship The existence of a local equilibrium at any point within the porous matrix is assumed. Kelvins law is applied to give the Eq(4). 1407           T v R l ρ g P w P exph (4) The vapour partial pressure can be defined as a function of local temperature and relative humidity as Eq(5). h o ρ v ρ  (5) where the saturation vapour pressure, ρo is estimated using the equation in Mayhew and Rogers (1976) as a function of temperature. The degree of saturation, S is an experimentally determined function of capillary pressure and temperature (Harun et al., 2014c) as Eq(6) . m n (T))(α1 1 r θ s θ r θθ l S               (6) where the parameters α, n and m are dependent on the porous material properties and these influence the shape of the water retention curve. The permeability of water and gas are based on Muelem’s model (Mualem, 1976). 2.3 Solution of governing equations and numerical method The coupled heat and mass transfer equations described above, in 2-dimensions can be written into a matrix form as Eq(7). )(}){)]([)](([}{ t )[C( ZR y i cy K x i cx K     (7) where {Φ}={Pw, T, Pg} is the column of unknowns;[C],[Kcx]and [Kcy] are 3x3 matrices. Each element of the matrix is a coefficient for the unknown {Φ}; ix and iy are the unit direction vectors. In order to discretize this simplified second order non-linear coupled partial differential equation, finite element method was used. Afterwards, Galerkin method was used to minimize the residual error before the application of Greens theorem, to the dispersive term involving second order derivatives. The transient matrix and nonlinear second order differential equations above were then solved by using a fully implicit backward time stepping scheme along with a Picard iterative method take into account for non-linearity. 3. Material data In this work, a two layers system mapped with a coarse mesh (10 elements) and a fine mesh (20 elements) was observed and investigated. The domains are assumed to be a porous medium that is homogeneous, isotropic and composed of solid phase, water and vapour phase, gas phase and dry air phase. Both domains have the same initial saturation of S0=0.6 and the initial temperature T0=298 K. The two domains that have been created are shown in Figure 1. Table 1: Physical properties of ceramic body and transport parameter Layer Symbol Unit Bottom Top Density of porous matrix ρs kg/m 3 2,000 3,970 Porosity Ø 1 0.4 0.1 Intrinsic permeability K m 2 2.9 x 10 -15 5.13 x 10 -17 Thermal conductivity of the porous matrix λ W/m K 1.8 39 Specific heat capacity of porous matrix Cp J/kg K 925 775 Critical saturation Scri 1 0.3 0.3 Irreducible saturation Sirr 1 0.09 0.09 Equilibrium saturation Seq 1 - 0.0499 The temperature was set for slow convective drying with the temperature at 300 K and heat and mass coefficient were 0.7 W/m 2 K and 0.00075 ms -1 respectively. Both domains are exposed to the convective boundary condition at one side (coarse-node 1,12,18,29 and 35; fine-node 1, 22, 33, 54 and 65) while the 1408 other sides are insulated and impermeable. The results for multilayer system are presented in the form of saturation contours and all figures are presented based on the same scale of contour colour in order to compare the mesh sensitivity scale. The material properties and transport parameter used as input data for simulation are listed in Table 1 with top layer consists of hygroscopic material whilst bottom layer is a non-hygroscopic material. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 X-coord (m) Y -c o o rd (m ) 0 0.00025 0.0005 0.00075 0.001 0 0.00025 0.0005 0.00075 0.001 0.00125 0.0015 0.00175 0.002 (a) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 X-coord (m) Y -c o o rd (m ) 0 0.00025 0.0005 0.00075 0.001 0 0.00025 0.0005 0.00075 0.001 0.00125 0.0015 0.00175 0.002 (b) Figure 1: Domain meshing (a) coarse; (b) fine 4. Results and discussion (a) (b) Figure 2: Saturation curve for validation (a); between fine mesh and coarse mesh at node 22, 54(fine) and 12, 29(coarse) over 8 hour drying time for two layer structure Validation of the model was done by comparing the saturation profiles gained from experiments and simulation. Figure 2(a) shows a comparison of the saturation profile obt ained by experiment and computed results by Stanish et al. (1986) model together with our proposed model. The saturation curve obtained is in closed agreement at all times. The typical drying stages can be clearly distinguished from the figure. The CRP which is indicated by a straight line from point A to B and FRP as a gradual curve from point B to C. Detailed explanations of all these periods of drying will be discussed later. After instilled confidence in the proposed model efficiency to generate the right evolution of the drying curve, case study to explore the meshing sensitivity in membrane structure is performed. Figure 3 depicts the saturation evolution with time in contour presentation for both meshes. Generally, the drying saturation curve is characterized with the linear line in the CRP (refer Figure 2(a), point A to B) and gradually decreasing curve in the FRP (refer Figure 2(a), point B to C). In natural convective drying, capillary mechanism is the dominant force for moisture transport. Initial moisture content and hydraulic diffusivity of porous materials are sufficiently high enough and the entire surface area is covered by the 1409 film of free water (Shahidzadeh-Bonn et al., 2007). This lead to liquid continuously evaporates through the surface of the heating boundary without being replaced by air and area of contact between the air and medium is constant. Liquid migrates from regions of high concentration of moisture content in large pores towards the low concentration of moisture content in the small pores. Drying rate is greater at CRP and depends only on the external conditions such as temperature, relative humidity, velocity and others, rather than the rate of transport within the material itself (Perré et al., 2007). As drying proceeds, moisture content reaches a critical value of 0.3 for most porous materials (Zhang et al., 1999). This is stage denotes as the FRP. FRP happens once insufficient free water passes through randomly the flow paths in a medium leading to a formation of percolation threshold. Consequently, water phase will continue when the free water content is above than the critical value and the formation of discontinuous wet patches on surface occurs once the free water level is below the critical regardless of the rate of internal moisture transfer. Those dry patches form at surface still contains bound water in their pores and it offers a resistance to heat transport, resulting in a reduction of evaporation rate as well as drying rate. At this moment, capillary mechanism which is dominated before FRP is slowly diminishing as the water is strongly bonded to the porous matrix. Vapour movement by diffusion will control the much slower drying rate at this moment associated with the decreasing in moisture content. Finally, drying process eventually ceases when moisture content reaches equilibrium moisture content denotes the final value at the end of drying (Kowalski and Mierzwa, 2013). The end of this drying process at this stage is generally only referring to non-hygroscopic material (denotes as nhg in Figure 2(b)). As for hygroscopic material (denotes as hg in Figure 2(b)) which exhibit bound water mechanism, drying stage will proceed to the next stage. Bound water transport mechanism only effective when saturation irreducible is reached for hygroscopic materials. The movement of bound water in hygroscopic materials is also known as liquid moisture transfer near dryness or sorption diffusion with driving force of vapour transport and without liquid transport (Zhang et al., 1999). Also noted in Figure 2(b) and Figure 3, top layer (hygroscopic) has a higher saturation all times due to higher water retention during CRP and FRP before the bound water mechanism. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 X-coord (m) Y -c o o rd (m ) 0 0.00025 0.0005 0.00075 0.001 0 0.00025 0.0005 0.00075 0.001 0.00125 0.0015 0.00175 0.002 Level Sl 15 38.7584 14 38.5727 13 38.387 12 38.2013 11 38.0156 10 37.8299 9 37.6442 8 37.4585 7 37.2728 6 37.0871 5 36.9013 4 36.7156 3 36.5299 2 36.3442 1 36.1585 (a) Level : 15 = 38.9441 1= 35.9728 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 X-coord (m) Y -c o o rd (m ) 0 0.00025 0.0005 0.00075 0.001 0 0.00025 0.0005 0.00075 0.001 0.00125 0.0015 0.00175 0.002 Level Sl 15 38.7584 14 38.5727 13 38.387 12 38.2013 11 38.0156 10 37.8299 9 37.6442 8 37.4585 7 37.2728 6 37.0871 5 36.9014 4 36.7157 3 36.53 2 36.3442 1 36.1585 (b) Level : 15 = 38.9441 1= 35.9728 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 X-coord (m) Y -c o o rd (m ) 0 0.00025 0.0005 0.00075 0.001 0 0.00025 0.0005 0.00075 0.001 0.00125 0.0015 0.00175 0.002 Level Sl 15 6.9 14 6.85 13 6.83 12 6.8 11 6.74672 10 6.7 9 6.65 8 6.6 7 6.55 6 6.5 5 6.49636 4 6.45 3 6.4 2 6.3 1 6.246 (c) Level : 15 = 9 1= 4.99418 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 X-coord (m) Y -c o o rd (m ) 0 0.00025 0.0005 0.00075 0.001 0 0.00025 0.0005 0.00075 0.001 0.00125 0.0015 0.00175 0.002 Level Sl 15 6.99627 14 6.97 13 6.85 12 6.8 11 6.78 10 6.7458 9 6.7 8 6.6 7 6.5 6 6.49533 5 6.45 4 6.4 3 6.35 2 6.3 1 6.24486 (d) Level : 15 = 9 1= 4.99253 Figure 3: Comparison of saturation level at 1 hour((a),(b)); 8 hour((c),(d)) for both meshes 1410 Meanwhile, in our case study of meshing sensitivity, both meshes show very little different (less than 0.02 %)(denotes as diff in Figure 2(b)) over the period of drying (refer at Figure 2(b) and Figure 3). However, the significant difference in computing time is observed as fine mesh need almost 15 times computational times than coarse mesh. This might due to the complexity of in the numerical as more nodes and elements to take into consideration when obtain solution. Thus, it can be concluded that our model is mesh independence and a coarse mesh is always adequate to be used in future study. 5. Conclusions In this present work, a novel model with fully coupled heat, mass transfer and gas transport of the dynamic drying process together with thermodynamic of different phases of the porous multilayer structure was proposed. Validation work on this model and meshing dependency work were conducted using the complex and dynamic of the drying curve. One of the great advantages of our fully coupled model is able to offer a less mesh dependency in simulating highly dynamic changes of drying variables. Two simulated results of the mesh dependence test show only a slight difference between a coarse and fine mesh. Model is also able to simulate precisely the drying curve for both meshes using a coarse mesh. Such a model also enables construction processes to be optimised with respect to drying time, material selection and energy saving as well as reducing the environmental effect via less waste energy losses. Acknowledgements The authors would like to thank MOHE and UTHM for the financial support. References Gariev, A.O., Klemeš J.J., Kusakov S.K., Tovazhnyansky L.L., Anokhin P., Kapustenko P.O., Arsenyeva O.P., Čuček L., 2014, The Development of Heat Substation for Drying Waste Heat Utilisation, Chemical Engineering Transactions, 39, 1405-1410. Harun Z., Ong T.C., 2014a, Comparative Studies of Heat and Mass Transfer by Convective and Microwave-Convective Drying for Nonhygroscopic Ceramic, Aust. J. Basic Appl. Sci., 8, 218–225. Harun Z., Ong T.C., 2014b, Material parameters sensitivity in modeling drying of porous materials, Adv. Mater. Inf. Technol. Process., 87, 91–98. Harun Z., Ong T.C., Fauzi A.F., 2014a, Simulation of Drying for Multilayer Membranes, J. Teknol., 69, 107- 115. Harun Z., Ong T.C., Rosli A., 2014b, Drying Comparison of Nonhygroscopic and Hygroscopic Materials, Adv. Mater. Res., 465, 637–641. Harun Z., Ong T.C., Rosli A., 2014c, Estimation of Water Desorption in Drying of Membrane Structure System, Appl. Mech. Mater., 1004-1005, 405–408. Kanno T., Kato K., Yamagata J., 1996, Moisture movement under a temperature gradient in highly compacted bentonite, Eng. Geol., 41, 287–300. Kowalski S.J., Mierzwa D., 2013, Numerical analysis of drying kinetics for shrinkable products such as fruits and vegetables, J. Food Eng., 114, 522–529. Mayhew Y.R., Rogers G.G.C., 1976, Thermodynamic and Transport Properties of Fluids. Blackwell, Oxford, England. Mualem Y.A., 1976, A new model for predicting the hydraulic conductivity of unsaturated porous media, Water Resour. Res., 12, 513–522. Mujumdar A.S., Eds., 2006, Handbook of Industrial Drying. Taylor & Francis Group, Boca Raton, USA. Perré P., Remond R., Turner I.W., 2007, Comprehensive drying models based on Volume Averaging: Background, Application and Perspective, vol. 1: Modern Drying Technology, Eds. Tsotsas E., Mujumdar A.S., Wiley VCH Verlag GmbH & Co. KGaA, Weinheim, Germany. Scherer G.W., 1990, Theory of Drying, J. Am. Ceram. Soc., 73, 3-14. Tarantino A., Sacchet A., Dal Maschio R., Francescon F., 2010, A Hydromechanical Approach to Model Shrinkage of Air-Dried Green Bodies, J. Am. Ceram. Soc., 93, 662-670. Shahidzadeh-Bonn N., Azouni A., Coussot P., 2007, Effect of wetting properties on the kinetics of drying of porous media, J. Phys. Condens. Matter, 19, DOI: 10.1088/0953-8984/19/11/112101. Stanish M.A., Schajer G.S., Kayihan F., 1986, A mathematical model of drying for hygroscopic porous media, AIChE J., 32, 1301–1311. Zhang Z., Yang S., Liu D., 1999, Mechanism and mathematical model of heat and mass transfer during convective drying of porous materials, Heat Transf. Asian Res., 28, 337–351.