Simulación multi-física del secado de madera en COMSOL Multiphysics 3.4.

 

Multi-Physics Simulation of Wood Drying in COMSOL Multiphysics 3.4.

 

Sandoval-Torres Sadoth¹, Hernández-Bautista Emilio² y Rodríguez-Ramírez Juan³

 

¹ Departamento de Ingeniería, Instituto Politécnico Nacional-CIIDIR Unidad Oaxaca. Correo: sadothsan@yahoo.com

² Instituto Politécnico Nacional-CIIDIR Unidad Oaxaca. Correo: Bautistahe@gmail.com

³ Instituto Politécnico Nacional-CIIDIR Unidad Oaxaca, Departamento de Ingeniería. Correo: jrodrigr@hotmail.com

 

Información del artículo: recibido: febrero de 2012.
Reevaluado: junio de 2012.
Aceptado: agosto de 2012

 

Resumen

En este trabajo se simulan las cinéticas de secado utilizando un modelo matemático de tipo fenomenológico que describe el secado por convección de aire caliente de madera Pinus pseudostrobus. El desarrollo del modelo parte del balance de masa y calor en el elemento de volumen representativo (EVR), que contiene las fases sólida, líquida, y gaseosa. Se obtiene un sistema de ecuaciones diferenciales parciales que es resuelto por factorización numérica utilizando COMSOL Multiphysics 3.4. Se resuelven tres variables primarias: contenido de humedad; temperatura y densidad del aire seco. Los resultados del modelo fueron comparados con los datos experimentales y contra un modelo semi-empírico. Los resultados son satisfactorios; se describen las cinéticas de secado simuladas, los perfiles espaciales de contenido de humedad al interior del material, la evolución de la masa del aire seco y los perfiles de temperatura.

Descriptores: secado, Pinus pseudostrobus, modelo fenomenológico, Comsol multiphysics, separación de fases.

 

Abstract

In this work, the numerical simulation of drying of Pinus pesudostrobus wood is presented by using a phenomenological model. The model is developed by considering the heat and mass balance in the representative elementary volume (REV), which involves the solid, liquid and gas phases. We obtained a system of partial differential equations which was solved by numerical factorization by using COMSOL multiphysics 3.4. Three primary variables were solved: the moisture content, the temperature, and the dry-air mass. The numerical results were compared against experimental data and the semi-empirical model. Our results are satisfactory; we describe the drying kinetics, the moisture distribution within the wood, the dry-air mass evolution and the temperature profiles.

Keywords: Drying, Pinus pseudostrobus, phenomenological model, Comsol Multiphysics, phase separation.

 

Introducción

La eliminación de agua en materiales higroscópicos como la madera es un fenómeno muy complejo, el cual involucra una serie de mecanismos acoplados que ocurren simultáneamente en el proceso de secado. En el proceso de transferencia de masa, el agua libre se elimina primero de la matriz sólida por mecanismos de capilaridad, hasta llegar al punto de saturación de la fibra PSF (aproximadamente en el contenido de humedad de 0.30 kg de agua/kg de sólido seco para coníferas). Al mismo tiempo que el agua es eliminada de los capilares en la madera es reemplazada por el aire que se introduce al material debido a un gradiente de concentraciones además de la diferencia de presiones en el sistema. Finalmente el agua ligada es removida de la matriz sólida por vaporización, y posteriormente por los mecanismos de difusión y presión de vapor de agua hasta llegar al contenido de humedad de equilibrio.

Existen diferentes tipos de modelos de secado de madera; sin embargo, los modelos que involucran en sus ecuaciones los fenómenos de transporte descritos aquí, son los modelos de tipo fenomenológico que parten de las ecuaciones de Whitaker (1977), y que son modificadas por Perré y Turner (1999). Estas ecuaciones involucran balances acoplados de masa, calor y momentum, utilizando propiedades en un volumen promedio EVR (elemento de volumen representativo).

Por lo tanto, un modelo de tipo fenomenológico para el secado de madera se compone de un conjunto de ecuaciones diferenciales parciales EDP, que sólo puede ser resuelto por métodos numéricos. Comsol Multiphysics 3.4 ofrece una forma de incluir la geometría a modelar y el sistema de EDP, evitando la necesidad de programación; además de que el programa incluye los códigos más avanzados para la solución numérica, lo que disminuye tiempos de cómputo, y asegura la convergencia de las ecuaciones.

Se han realizado diferentes estudios similares para especies de pino y encino europeo (Fernandez y Howell, 1997; Raji et al., 2009; Turner et al., 1998), sin embargo, para especies de pinos mexicanos existen pocos estudios con el enfoque de este trabajo, por lo que al modelar el proceso de secado se obtendrá información que podría utilizarse en la elaboración de modelos de secado más complejos que se apliquen en la optimización del diseño y control automático de secadores para especies endémicas, además de contribuir con el aprovechamiento forestal mexicano.

El objetivo de este trabajo es desarrollar un modelo matemático unidimensional que describa los mecanismos de transporte de calor y masa durante el secado por convección de aire caliente de madera Pinus pseudostrobus y resolverlo numéricamente utilizando COMSOL Multiphysics 3.4

 

Desarrollo

Con el objetivo de comparar las simulaciones del modelo con resultados experimentales se realizaron experimentos de secado basados en un diseño uni-factorial con la temperatura como factor y cuatro niveles: 50, 60, 70 y 80°C, manteniendo la velocidad del aire constante (2.5 m/s), y sin controlar la humedad relativa del aire. Se utilizó un túnel de secado que está diseñado para hacer pasar por la muestra un flujo de aire con velocidad y temperatura uniforme y controlada. El aire se calienta mediante dos resistencias eléctricas de 20Ω, (2.4 kW), controlada por medio de un control PID. La medición de la temperatura se da a través de termopares tipo J, conectados simultáneamente y aislados eléctricamente para reducir el error del sistema. El peso de la muestra durante el proceso es registrado por el sistema de adquisición de datos conectado a una celda de carga (rango de 0-11.34 kg, precisión de 0.1%) que mide la tensión que se produce por el peso de la muestra.

El túnel cuenta con un controlador automático programmable National Instruments de la serie FP-1000, que realiza el control del proceso y la adquisición de datos en tiempo real, conectado a un programa de computadora que registra las cantidades por medio del software LabVIEW.

Las muestras de madera fresca se obtuvieron al azar del almacén del aserradero ubicado en la comunidad de Macuilxóchitl de Artigas Carranza, Oaxaca. Las dimensiones de las muestras fueron de 25 cm de largo, 15 cm de ancho con un espesor de 2.54 cm. Las tablas frescas de madera fueron rociadas con ayuda de un dispersor y cubiertas con plástico impermeable para su traslado al laboratorio, en donde fueron almacenadas a 0°C en el refrigerador, para evitar su degradación.

Por cada experimento se utilizaron dos tablas de madera fresca. Los bordes fueron sellados con silicón con el objetivo de que el transporte de humedad ocurra sólo en la dirección del espesor de la madera. Una tabla fue perforada en los extremos del espesor para la introducción de termopares en el centro y a un milímetro de la superficie, como se muestra en la figura 1. Los experimentos tuvieron una duración promedio de 70 horas, cada uno con su réplica.

 

Formulación matemática

Se desarrolló una formulación matemática, la cual representa el proceso de secado por convección de aire caliente basado en el modelo de (Whitaker, 1977), y posteriormente de (Perré y Turner, 2001). Debido a los arreglos experimentales, configuración de los poros y punteaduras en la madera, además de la temperatura del aire del túnel debajo del punto de ebullición del agua, se puede suponer que el transporte de humedad se da principalmente en la dirección del espesor del material (Turner 1996). La geometría a modelar será el espesor de la madera, una línea recta (figura 1).

Los balances de materia y calor son incorporados en COMSOL Multiphysics 3.4, en el navegador de modelos, en la forma EDP general (ecuación diferencial parcial) para los balances de humedad y aire seco, y en la forma de EDP de coeficientes para el balance de calor, donde se declaran las variables primarias a resolver: contenido de humedad (W); densidad intrínseca del aire (p_a); y temperatura (T). Como se muestran en la tabla 1.

La ecuación de conservación del aire seco (2) permite una correcta descripción de los gradientes de presión, los cuales se hacen presentes durante el secado. La ecuación de transferencia de humedad (1) es una versión más elaborada del modelo difusivo. Sin embargo, la principal diferencia entre las dos aproximaciones es la separación de fases, así como las relaciones constitutivas que comprenden la contribución de cada una de las fases (Krabbenhøft, 2003). La ecuación de conservación de energía (3) involucra la conducción de calor, y la transferencia de calor convectivo, además de los cambios de entalpía debidos a los cambios de fase (m_bv y m_iv)

Los flujos de humedad, calor y aire seco, son descritos en COMSOL en ajustes del subdominio; el flujo de agua libre se supone que sigue un comportamiento generalizado que se puede describir con la ley de Darcy, para la cual la velocidad promedio de masa está dada por la ecuación 4. Además se ha demostrado que las fuerzas que ejerce el cuerpo debido a la gravedad son muy pequeñas y por lo tanto despreciables (Plumb y Prat, 1992; Turner et al., 2010).

La transferencia de vapor de agua y aire seco puede ser descrita por la combinación de la ley de Fick y la ley de Darcy. A través de esto podemos representar el transporte de vapor de agua con la siguiente formulación matemática de la ecuación 6. El flujo de aire seco (ecuación 7) tiene un comportamiento parecido al del vapor, pero en sentido opuesto (tabla 2).

Por último, el mecanismo de transporte de agua ligada en el dominio higroscópico se describe en la ecuación 5, donde el flujo de difusión-desorción del agua ligada toma en cuenta dos mecanismos, el primero es difusión debido a los gradientes de contenido de humedad, y la segunda parte a los gradientes de temperatura (efecto Soret).

En la figura 1 se puede observar al lado derecho la geometría a modelar. Existen dos condiciones de frontera CF1 y CF2, las cuales están a la misma temperatura y por las que pasa el mismo flujo de masa (mezcla aire-vapor), a la misma velocidad de flujo de aire. Por lo tanto, los ajustes de contorno de la geometría propuesta para las superficies expuestas al secado son las mostradas en la tabla 3.

La ecuación 8 representa el flujo de agua que se está perdiendo en la superficie de la madera por evaporación; debido a las presiones atmosféricas la cantidad de aire en los alrededores no cambiará, por lo tanto, la condición de frontera para la ecuación de conservación total de aire seco (9) se mantendrá constante en todo el proceso de secado (condición de contorno tipo Dirichlet). El calor en la superficie se remueve por convección forzada que se representa con la ley del enfriamiento de Newton en la ecuación 10. Asimismo, es necesario considerar el calor que se pierde en la superficie, debido a la cantidad de agua que se está evaporando. h y k_m son los coeficientes de transferencia de calor y masa, respectivamente, cuyos valores son determinados experimentalmente a las condiciones de secado.

Para la simulación de las cinéticas de secado es necesario contar con las propiedades termo-físicas de la madera de la especie Pinus pseudostrobus, la cual ha sido poco estudiada. Por lo tanto, existen parámetros que aún se desconocen, sin embargo, se han utilizado propiedades de maderas similares en densidad. Algunas de las propiedades utilizadas se muestran en la tabla 4.

El sistema de ecuaciones fue resuelto utilizando una computadora con procesador AMD Athlon(tm)X2 Dual-Core, 2100 Mhz, con un tiempo de cómputo de 20 s.

 

Discusión y análisis de resultados

A continuación se presentan los resultados de la simulación, en la que se considera una temperatura del aire en el túnel de 60°C, con temperatura inicial de 25°C y contenido de humedad inicial de la madera de 96%.

Las cinéticas de contenido de humedad promedio son simuladas con el modelo fenomenológico y comparadas con datos experimentales, así como con el modelo de la curva característica de secado previamente publicados en Hernández et al. (2010). La figura 2 muestra los ajustes de los modelos comparados con datos experimentales, además se incluye un eje del lado derecho para mostrar los perfiles de temperatura y su evolución durante el proceso de secado simulados y experimentales en el centro de la madera y en la superficie.

La figura 3 muestra el contenido de humedad en el espesor de la madera y su evolución durante el tiempo en intervalos de dos horas. Durante el secado, a medida que se elimina el agua, es reemplazada con aire. En las superficies la densidad del aire se incrementará primero (figura 4). Estos perfiles tienen un comportamiento contrario a la densidad del vapor, ya que el desplazamiento de vapor de agua es reemplazado por aire seco. Recordemos que la mezcla gaseosa es una combinación de vapor de agua y aire seco, la cual es descrita por la ley de Dalton, que dice que la presión de una mezcla de gases, que no reaccionan químicamente, es igual a la suma de las presiones parciales de cada una de ellas.

La temperatura del material es la tercera variable que fue resuelta con el modelo. La figura 5 muestra los perfiles de temperatura en el espesor de una tabla de madera cada media hora durante el proceso.

 

Discusiones

El modelo matemático del tipo fenomenológico implementado muestra una asertiva simulación de las curvas de secado, lo cual es expuesto en la figura 2. Se puede observar un mejor ajuste utilizando el modelo de la curva característica (CCS); esto es de esperarse debido a que el modelo fenomenológico no parte de datos experimentales, ni es un ajuste de datos. Muchas de las propiedades de transporte requeridas por el modelo corresponden a especies de pino similares en densidad, por lo cual las variaciones en las cinéticas simuladas pueden disminuir si se tiene acceso a esta información (Perré y Turner, 2001).

En la figura 2, se muestra una parte del calentamiento al principio de la cinética de secado, que es mucho más notable a altas temperaturas debido a la solución numérica, incluso mostrada en otros modelos, por ejemplo en Perré (1999). Al principio se muestra siempre el fenómeno de incremento en el contenido de humedad en los primeros segundos. Por otra parte, el contenido de humedad de equilibrio se simula adecuadamente en la mayoría de las cinéticas, dado que los datos utilizados sobre las isotermas de desorción son específicos para esta especie de madera (Pinus pseudostrobus). Por lo tanto, para contenidos de humedad debajo del PSF se simula adecuadamente en comparación con la fase capilar.

En cuanto a los perfiles de temperatura que se muestran en la figura 2, éstos exponen una diferencia notable entre la temperatura superficial del material, la temperatura del centro y las simuladas por el modelo. Estas discrepancias pueden ser atribuidas al desconocimiento de las funciones exactas de la conductividad térmica del material, de la permeabilidad y de la presión capilar, entre otras propiedades utilizadas en el modelo; en este trabajo se tomaron las funciones de otros estudios de pino previamente realizados en Hernández y Puiggali (1994).

La cantidad de calor que se suministra a la madera también depende de la velocidad y de la humedad relativa del aire; estos parámetros son considerados en el coeficiente de transferencia de calor y de masa, pues entre más grande sea este coeficiente mayor será la cantidad de calor transferida en la interfase. Si el coeficiente es muy grande, la temperatura superficial aumenta rápidamente.

En la figura 3 se muestra la variación espacial del contenido de humedad en el espesor de una muestra de madera y su evolución durante el tiempo en intervalos de dos horas. A tiempo cero se tiene un contenido de humedad homogéneo (0.9 kg de agua/kg de sólido seco) en todo el espesor, el cual comienza a decrecer rápidamente en las orillas, producto del transporte de masa en la interfase. En dos horas de secado los extremos casi llegan al contenido de humedad de equilibrio. Por otra parte, se observa que al contenido de humedad cercano al PSF, existe una fluctuación en la parábola que forma el perfil de contenido de humedad en cada determinado tiempo. Esto se debe a la transición de la fase capilar a la fase higroscópica. Esta fluctuación depende principalmente de la permeabilidad relativa del líquido que decrece durante el secado. Krabbenhøft (2003) describe la influencia de la permeabilidad relativa en la simulación del secado y explica cómo la ecuación en función de la saturación suaviza las parábolas con funciones de saturación al cuadrado. Los perfiles presentados en este trabajo son similares a los presentados por Krabbenhøft (2003).

Al incrementar la densidad del aire en el interior de la madera también incrementa la presión en el gas. En el secado por convección de aire caliente a bajas temperaturas, la ecuación del aire seco describe con precisión el desarrollo de los gradientes de presión de gas en el material. El lugar en donde primero se elimina el agua libre es en la superficie (figura 3). Es ahí donde la densidad del aire incrementa primero, al tener también un perfil de densidad del aire dentro del material, que varia con el tiempo de secado (figura 4). Nótese cómo la masa de aire en el medio poroso es también alterada por la transición de la zona capilar a la zona higroscópica, producto de una fuerte interacción entre las fases.

La figura 5 muestra los perfiles de temperatura en el espesor de una tabla de madera cada media hora, es importante ver cómo la temperatura del material aumenta rápidamente. Aproximadamente en 4 horas se alcanza una temperatura homogénea de 60°C en el material, desde el centro hasta la superficie, contrariamente a lo que pasa en la realidad. La estabilización de temperaturas ocurre en 10 horas y el gradiente de temperatura entre la superficie y el centro es aproximadamente de 5°C. La diferencia entre la temperatura del centro y la superficie disminuye con el tiempo de secado (figura 2). En la simulación la diferencia es aproximadamente de un grado centígrado o menor entre el centro y la superficie de la madera. Como se explicó anteriormente esto se debe al flujo de calor en la madera, expresado en la conductividad térmica. En este caso, dicha conductividad térmica varía con el contenido de humedad solamente (tabla 4), y ésta va disminuyendo a lo largo del secado linealmente, ya que el agua tiene un coeficiente de conductividad térmico más grande que el del aire o el vapor, por lo que la conductividad térmica de la madera es alterada por la presencia de las fases líquida y gaseosa. Al final del secado, sólo se tiene la matriz sólida y aire seco, los cuales son menos conductivos. Es necesario mencionar que en la realidad el coeficiente de conductividad térmica no sólo varía con respecto al contenido de humedad: también lo hace con la temperatura, pero dicho parámetro no siempre está disponible para todos los materiales.

 

Conclusiones

Se desarrolló un modelo matemático unidimensional que describe los mecanismos de transporte de calor y masa durante el secado convencional de madera de Pinus pseudostrobus. El modelo fue comparado con el modelo de la curva característica de secado y con datos experimentales, describiendo correctamente las cinéticas de secado.

En las simulaciones realizadas se encontraron fluctuaciones en los perfiles espaciales de masa, en la transición de la fase capilar a la higroscópica, y diferencias en los perfiles térmicos; por lo tanto, es necesario un análisis más detallado de esta transición, y el cálculo de las propiedades termofísicas de esta madera.

Sin embargo el modelo fenomenológico desarrollado toma en cuenta automáticamente las modificaciones en las condiciones del proceso y los fenómenos de transporte de calor y masa que existen dentro del espesor de la madera, describiendo con mayor detalle el transporte de humedad y de energía al interior del material.

Se espera que la información generada por las simulaciones pueda ofrecer un preámbulo para propósitos de análisis de procesos de secado o diseño de nuevos modelos que se aproximen mejor a la realidad. Estos resultados pueden ser utilizados para la optimización y diseño de procesos de secado.

 

 

Agradecimientos

Al Instituto Politécnico Nacional.

 

Referencias

Baronasa R., Ivanauskasa F., Juodeikienéc I., Kajalaviciusc A. Modelling of Moisture Movement in Wood during Outdoor Storage. Nonlinear Analysis: Modelling and Control, volumen 6 (número 2), 2001.         [ Links ]

Colakoglium M.H. Determination of Change in Moisture Ratios of Some Woods During Air-Drying by Finite Element Analysis. Journal of Applied Sciences, volumen 9 (número 22), 2009: 4091-4094.         [ Links ]

Fernandez M., Howell J. Convective Drying Model of Southern Pine. Drying Technology, volumen 15 (número10), 1997: 2343-2375.         [ Links ]

Fuentes-Salinas M. Estimación del punto de saturación de la fibra (PSF) de las maderas. Revista Chapingo Serie Ciencias Forestales y del Ambiente, volumen 6 (número 1), 2000: 79-81.         [ Links ]

Hernández-Bautista E., Sandoval-Torres S., Rodríguez-Ramírez J., Velasco-Cruz O. A. Modelización del secado de madera de Pinus pseudostrobus utilizando el método de la curva característica. Bois et forêts des tropiques, volumen 306 (número 4), 2010.         [ Links ]

Hernández J.M., Puiggali J.R. Simulation of Drying of Coniferous Wood Using Various Processes. International Chemical Engineering, volumen 34 (número 3), 1994: 339-350.         [ Links ]

Kang W. y Chung W.Y. Liquid Water Diffusivity of Wood from the Capillary Pressure-Moisture Relation. The Japan Wood Research Society, volumen 55 (número 1), 2009: 91-99.         [ Links ]

Krabbenhøft K. Moisture Transport in Wood. A Study of Physical-Mathematical Models and their Numerical Implementation, Technical University of Denmark, Denmark, 2003,169 p.         [ Links ]

Perré P. How to Get a Relevant Material Model for Wood Drying Simulation?, Advances in Drying of Wood, 1999.         [ Links ]

Perré P. y Turner I. Transpore: A Generic Heat and Mass Transfer Computational Model for Understanding and Visualising the Drying of Porous Media. Drying Technology, volumen 17 (números 7 y 8), 1999: 1273-1289.         [ Links ]

Perré P. y Turner I. Determination of the Material Property Variations Across the Growth Ring of Softwood for Use in a Heterogeneous Drying Model, Part 2-Use of Homogenisation to Predict Bound Liquid Diffusivity and Thermal Conductivity. Holzforschung, volumen 55 (número 4), 2001: 417-425.         [ Links ]

Perré P. y Turner I.W. A Dual-Scale Model for Describing Drier and Porous Medium Interactions. AIChE Journal, volumen 52 (número 9), 2006: 3109-3117.         [ Links ]

Plumb O. y Prat M. Microscopic Model for the Study of Drying of Capillary Porous Media. Drying '92, 1992, pp. 397-406.         [ Links ]

Raji S., Jannot Y., Lagière P., Puiggali J. Thermophysical Characterization of a Laminated Solid-Wood Pine Wall. Construction and Building Materials, volumen 23, 2009: 3189-3195.         [ Links ]

Sandoval-Torres S., Rodríguez-Ramírez J., Mendez-Lagunas L., Carrillo-Parra A. Experimental Determination and Modeling of Equilibrium Moisture Content from the Sapwood of Mexican Pine (Pinus pseudostrobus), Forestry Studies in China, 2011.         [ Links ]

Turner I., Puiggali J., Jomaa W. A Numerical Investigation of Combined Microwave and Convective Drying of a Hygroscopic Porous Material: A Study Based on Pine Wood. Trans IChemE, volumen 76 (número A), 1998: 193-209.         [ Links ]

Turner I., Rousset P., Rémond R., Perré P. An Experimental and Theoretical Investigation of the Thermal Treatment of Wood (Fagus Sylvatica L.) in the Range 200-260 C. International Journal of Heat and Mass Transfer, volumen 53, 2010: 715-725.         [ Links ]

Whitaker S. Simultaneous Heat, Mass and Momentum Transfer in Porous Media: a Theory of Drying. Advances in Heat Transfer, volumen 13 (número 1), 1977: 119-203.         [ Links ]

Wook K., Chung W.Y., Eom C.D., Yeo H. Some Considerations in Heterogeneous Nonisothermal Transport Models for Wood: a Numerical Study. The Japan Wood Research Society, volumen 54 (número 1), 2008: 267-277.         [ Links ]

 

Semblanza de los autores

Sadoth Sandoval-Torres. Es ingeniero químico egresado del Instituto Tecnológico de Oaxaca, realizó los estudios de doctorado en Francia en la Universidad de Bordeaux, 2008 y realizó el pos-doctorado en el Instituto de la Filtración y de Técnicas Separativas de Francia. Es miembro de la red de talentos mexicanos en el extranjero y actualmente profesor colegiado en el Instituto Politécnico Nacional.

Emilio Hernández-Bautista. Es ingeniero químico egresado del Instituto Tecnológico de Oaxaca, maestro en ciencias en conservación y aprovechamiento de recursos naturales por el Instituto Politécnico Nacional y actualmente candidato a ingresar al doctorado en ciencias en el Centro Interdisciplinario de Investigación para el Desarrollo Integral Regional Unidad Oaxaca.

Juan Rodríguez-Ramírez. Doctor en ciencias en ingeniería química por el Instituto Tecnológico de Celaya, es miembro del SNI Nivel I y Premio Estatal de Ciencia y Tecnología Oaxaca 2011. Semblanza de los autores.