Conventional wood drying simulation using CVFEM 

Simulación del secado convencional de madera usando CVFEM 

Carlos Salinas, Cristian Chávez, Yerko Gatica, Rubén Ananias 

Universidad del Bío-Bío. Avda. Collao 1202, Concepción, Chile. Telf. (56) 041 2-731509.
Fax (56) 041 2-731022. casali@ubiobio.cl, crchavez@alumnos.ubiobio.cl, ygatica@ubiobio.cl, ananias@ubiobio.cl 

Abstract 

The main objective of the present work is the development of an original two-dimensional numerical simulation of wood drying, based on the water potential concept using the Control Volume Finite Element Method (CVFEM). The mathematical model was based on a set of non-linear second partial differential equations, which were solved by a numerical procedure using CVFEM. Results of drying curves and isoconcentrations are shown for aspen (Populus tremuloides) sapwood and compared with numerical and experimental data. Furthermore, the impact of the variation of various numerical parameters was studied to evaluate the algorithms consistence. It had been concluded that the current numerical method was consistent and allows to simulate wood drying satisfactorily resulting in detailed information about wood moisture content distribution. 

Key words: drying, wood, simulation, water potential. 

Resumen 

El principal objetivo del presente trabajo es realizar una simulación bidimensional original del transporte de humedad en la madera, modelado en base al potencial hídrico usando el Método de Volúmenes Finitos conformado por Elementos Finitos (CVFEM). El modelo matemático contempló ecuaciones de transporte del tipo diferencial parcial no lineal de segundo orden integrado numéricamente a través CVFEM. Resultados de simulaciones de curvas de secado e isoconcentraciones fueron obtenidos para madera de álamo (Populus tremuloides), comparados con datos experimentales y numéricos publicados en la literatura especializada. Además, se presentan resultados de parámetros numéricos que permitieron evaluar la consistencia de los algoritmos. Se concluye que el método numérico propuesto permite simular consistentemente el secado de la madera, obteniendo información detallada de las distribuciones de humedad al interior de ésta. 

Palabras clave: secado, madera, simulación, potencial hídrico. 

Recibido el 5 de Agosto de 2009 

En forma revisada el 6 de diciembre de 2010 

Introducción 

La madera es un material no homogéneo, anisotrópico, poroso y no saturado sobre el cual intervienen mecanismos de transporte de humedad, en diversas escalas, por ejemplo: Difusión a nivel molecular, efectos capilares al nivel de escalas intermedias y deformaciones a macro escalas [1]. Las diversas complejidades, inherentes a la madera, exigen fuertes simplificaciones para modelar y simular el secado de la madera. En términos generales implica aceptar la hipótesis del medio continuo, eliminar la heterogeneidad, considerar la anisotropía sólo relevante en las direcciones principales, asumirla como indeformable para efectos de fenómenos de transporte y en equilibrio termodinámico para balances de calor y masa. La humedad es uno de los aspectos críticos del secado, encontrándose el contenido de agua que la define en estado líquido, gaseoso y ligada en las paredes celulares [8]. Además, en el proceso de secado se alteran los equilibrios de esfuerzo mecánicos ocasionando la deformación del material [3], lo cual también ocurre durante su acondicionamiento [4]. 

En base a lo anterior se han desarrollado diversos modelos y estrategias para la simulación del transporte de calor y masa bajo los preceptos clásicos difusivos [5] y aquellos de transportes desarrollados por Luikov [6] y Whitaker [7]. Una revisión de estos modelos es dada por Hernández y Quinto [8]. 

El modelo matemático implementado en este trabajo sigue la línea de investigación fundamentada por Luikov [6] y desarrollada en el tiempo por autores de la Universidad Laval, entre otros, Cloutier et al. [9], Cloutier y Fortin [10] y [11], Tremblay et al. [12] y [13] y Defo et al. [14]. Se establece que el potencial hídrico es la fuerza motriz que induce el flujo de humedad en la madera bajo condiciones isotérmicas. Se obtiene así un sistema de ecuaciones diferenciales parciales no lineales de segundo orden para los procesos de transporte energía y concentración en función de la humedad, la temperatura y la presión. En particular, para el caso del secado convencional que se estudia en este trabajo, las variables de temperatura y presión pueden ser omitidas, restando sólo la ecuación de transporte de concentración [11]. 

La ecuación de transporte de concentración del tipo diferencial parcial no lineal de segundo orden, es descrita numéricamente en términos del Método de Volúmenes Control conformado por Elementos Finitos [15], reconocido por las siglas CVFEM. En términos generales el método consiste en definir Volúmenes de Control o Volúmenes Finitos en base a contribuciones parciales de Elementos Finitos. En este sentido, CVFEM es un método de integración numérica que puede ser explicado por la combinación del Método de Volúmenes Finitos (FVM) y del Método de Elementos Finitos (FEM). En consecuencia, CVFEM se caracteriza, por preservar la simplicidad y conservatividad de los esquemas FVM, y al mismo tiempo, puede escribirse bajo una forma variacional similar a los esquemas FEM. Esto último, permite un mejor análisis matemático de estimación de errores, mediante el método de la energía. Por lo tanto, la principal característica de este método CVFEM, es conjugar las bondades de estos dos métodos clásicos (FVM) y (FEM): La conservación de los flujos, aportado por FVM, y la formulación variacional, característica del FEM. 

En particular, se seleccionan elementos finitos triangulares con distribución lineal de la variable dependiente al interior de este, debido a su versatilidad topológica y a la simplicidad de los algoritmos desarrollados. 

El objetivo del presente trabajo es realizar una simulación original del transporte de humedad en la madera, modelado en base al potencial hídrico y usando CVFEM. En consecuencia, el modelo es programado en lenguaje Fortran y aplicado a la simulación bidimensional isotérmica de difusión de humedad en madera de álamo considerando ésta como un continuo de propiedades ortotrópicas variables, en función del contenido de humedad. Se detallan aspectos transitorios del transporte de humedad y se comparan resultados con datos experimentales y numéricos dados en la literatura. 

Modelo físico 

Se considera una pieza bidimensional de madera sólida de álamo sometida a un proceso de secado convectivo (Figura 1), donde se asume un medio homogéneo ortotrópico de propiedades variables con el contenido de humedad M. Esto es: Conductividades k_ii y potencial hídrico y en función de contenidos de humedad M (k_ii(M) y y(M) respectivamente. 

La simetría respecto al flujo y forma de la madera permite reducir a un cuarto el dominio de cálculo del transporte bidimensional de humedad en el plano transversal (radial/tangencial), caracterizado por un largo L y un ancho H iguales a 0.045 (m). 

Las condiciones iniciales y de contorno son de humedad inicial M = M_ini y condición de contorno del tipo Neumann de flujo nulo (q_m = 0) en los ejes de simetría (x = L/2 e y = 0.) y 3) y convección en las superficies x = 0 e y = H/2. 

Modelo matemático 

El modelo matemático considera que la variación local de la concentración de humedad, es equivalente a la divergencia del flujo y se puede escribir de acuerdo con el modelo propuesto por Cloutier et al., en la forma implementada por Salinas et al. [16]. Esto es: 

        (1) 

donde C es la concentración (kg_agua/m³_{maderahúmeda})  y el flujo (kg_agua/m²_{maderahúmeda}s). 

Suponiendo pequeñas variaciones de temperatura y equilibrio entre las fases del agua en la madera, el flujo de humedad q_m puede ser descrito en función del tensor de conductividad y el potencial hídrico y (J/kg) como: 

.         (2) 

En particular, el transporte de humedad es determinado en forma indirecta a través del potencial hídrico, considerando la madera como un medio ortotrópico bidimensional en el plano x,y. Este puede ser expresado como¹: 

    (3) 

donde y conductividades ortogonales (kg²_agua/m²_{maderahúmeda}s J), 

Capacidad (kg²_agua/J m³_{maderahúmeda}). 

Las conductividades y y la variación de humedad en relación al potencial, , son parámetros físicos de transporte determinados experimentalmente (ver Figura 1A). 

Alternativamente, (3) puede ser expresada en función del contenido de humedad M (100 kg_agua/kg_maderaseca): 

    (4) 

donde , y . 

Modelo numérico 

El modelo matemático dado por la ecuación (1, 3 ó 4) para el caso de transporte de humedad; puede ser presentado en forma genérica por una Ecuación de Difusión Transiente no lineal de segunda orden. Esto es: 

, con     (5) 

cuyas variables, propiedades y condiciones de contorno específicas son dadas en la Tabla 1. 

Al integrar (5), de acuerdo con el Teorema de Green, se tiene que: 

    (6) 

donde V es dominio de contorno S, cuya normal externa es n. 

El dominio V es conformado por nl Volúmenes Finitos (VF) ( con ) de contorno . De forma similar, cada VF es constituido por nk contribuciones parciales de Elementos Finitos (EF) de volumen con (Figura 2). 

con y
con .        (7) 

El término temporal es evaluando en forma implícita de primera orden (Euler Implícito): 

con .        (8) 

El término difusivo es integrado espacialmente sobre el contorno del EF k igual a  con para el tiempo  (omitido por simplicidad). Siendo así, las contribuciones parciales difusivas del EF k a los VF centrados en los nodos locales l_i (i = 1,2 y 3), se pueden escribir como: 

    (9) 

Expresiones similares pueden ser escritas para las contribuciones centradas en l₂ y l₃. 

Los subíndices a y b representan segmentos de contorno y respectivamente. 

Para la integración anterior se define una variación lineal de f al interior del EF, es decir: 

        (10) 

donde A, B, C son constantes definidas en función de valores nodales , y iguales a: 

    (11) 

    (12) 

              (13) 

      .        (14) 

De acuerdo con la notación de (9), la contribución integral difusiva del contorno definido por el segmento , es dada por: 

.    (15) 

Y en consideración a (10), el gradiente de f es dado por: 

.        (16) 

Los flujos difusivos y en direcciones ortotrópicas x e y, respectivamente, son dados por: 

.    (17) 

Asumiendo el gradiente y la normal como constantes en el segmento de contorno , la integral resultante es igual a: 

             (18) 

siendo igual a la longitud del segmento . 

Similarmente, son determinadas las contribuciones del complemento de segmentos que conforman el contorno del VF V^l. 

Las ecuaciones (8) y (18) tienen implícito el requerimiento de evaluar numéricamente los coeficientes , y , dependientes de la variable transportada, según modelo matemático dado por (4). En este sentido se modela considerando como predominante el valor de estos coeficientes en los centroides de cada VF. Para hacer esto viable, los valores experimentales discretos de los mencionados coeficientes (Figura A1-A3) son representados en forma continua por medio de una interpolación, de tipo exponencial, para rangos entre valores discretos disponibles. 

De esta forma, al considerar las contribuciones transiente, dada por la ecuación (8), sumado a los términos difusivos del tipo (18), se obtiene la ecuación de difusión transiente bidimensional discreta, planteada en forma implícita para f en en un volumen genérico l. Esto es: 

        (19) 

válida para ecuaciones, donde y . 

Siendo así, para cada valor l se obtendrá una ecuación algebraica formulada para el valor medio de f en cada subdominio , configurando un sistema de nl × nl ecuaciones de la forma 

.         (20) 

Según lo observado anteriormente, la dependencia de los coeficientes , y del contenido de humedad, induce la no linealidad del sistema de ecuaciones resultante. Es decir, la matriz A y el vector B de (20) son función de f. Siendo así, al mencionado sistema se le aplica un esquema de linealización update (update method: valores de los coeficientes evaluados ). Para hacer esto viable, se define una función continua de dichos coeficiente en base a una interpolación logarítmica de sus valores discretos experimentales, dados en el Apéndice A. Estos, varían linealmente al interior de cada EF. 

Resultados y discusión 

Se simula el transporte bidimensional transiente de humedad en una sección transversal en el plano radial-tangencial de madera de álamo. El problema esquematizado en la Figura 1 y sus propiedades dadas en el apéndice A, Tabla A1 y Figuras A1-A3. En particular, su magnitud característica se define por L = H = 0.045 m. Como condición inicial se considera M_ini = 135 % y como condición de contorno: flujo nulo en los ejes de simetría (x = L/2 e y = 0.) y convección en las superficies x = 0 e y = H/2. Lo anterior según discutido en el ítem de modelo matemático del presente trabajo. 

Todas las simulaciones son realizadas resolviendo los sistemas de lineales resultantes a través del Método de Gauss Seidel con SOR (Successive over relaxation) usando un factor de relajación w = 0.85, con un criterio de convergencia basado en normal residual menor que 10^–9. 

La Figura 3, muestra una malla característica utilizada en las modelaciones, la cual presenta un refinamiento de tipo logarítmico hacia las superficies con convección. En dichas regiones se tienen grandes gradientes de humedad, los cuales pueden ser óptimamente captados por mallas de este tipo. 

La Figura 4, muestra un análisis de consistencia y convergencia en relación al tipo de malla (distribución de elementos: uniforme, cosenoidal y logarítmica) para valores de Contenido de Humedad (CH) en el tramo BC. Se observa que las distribuciones que concentran elementos hacia las regiones de convección se aproximan mejor a la solución convergente, lo cual se consigue con una malla logarítmica 30×30. 

La Figura 5 muestra curvas de secado (contenido de humedad media v/s tiempo) para diversos pasos de tiempo. Se puede apreciar soluciones convergentes para dt = 0.1 h al utilizar una malla logarítmica 30×30. Este resultado sumado a los mostrados en la figura 4 permiten concluir que una malla logarítmica de 30×30 y un paso de tiempo de 0.1 h son adecuados para obtener simulaciones convergentes. Esta malla y paso de tiempo son utilizados en la generación de resultados mostrados en las Figuras 5-9. 

La Figura 6 muestra una curva de secado en donde se compararan simulaciones de acuerdo a tipo de variable dependiente (potencial o contenido de humedad). Se pueden apreciar leves diferencias, con una mejor aproximación de la formulación de contenido de humedad como variable dependiente. 

La Figura 7, muestra detalles de distribución de humedades simuladas para t igual 2 y 50 (h) de secado. Se puede observar el avance del frente de secado: Mayor concentración de iso-líneas de humedad, lo cual refleja mayores gradientes en el entorno de 70%. El comportamiento ortotrópico puede observarse en la falta de simetría polar de las iso-líneas de humedad. 

En la Figura 8 se comparan distribuciones de humedad simuladas por el presente trabajo (CVFEM) y las obtenidas en [9] a través de FEM. Se puede apreciar una buena concordancia cualitativa de resultados, incluidas las regiones de mayores gradientes en el entorno del contenido de humedad crítico (65% aprox.), que caracteriza el frente de evaporación. Este último fenómeno también puede ser apreciado visualmente en la Figura 7. Las diferencias medias entre niveles de CH no superan el 5%. 

Por último, la Figura 9 muestra resultados de curvas de secado obtenidas por la presente modelación, contrastada con resultados experimentales obtenidos por Cloutier et al. [9]. Se aprecia una buena concordancia, lo cual refleja la captación de la esencia del fenómeno físico en estudio. 

Comparando estos resultados con los obtenidos experimentalmente y numéricamente (Método de Elementos Finitos) por Cloutier et al. [9], se observa similitud tanto en la forma de las curvas de distribución de humedad, como en los valores simulados. Este comportamiento permite inferir que la modelación bidimensional del secado basado en el potencial hídrico, como fuerza inductora de la humedad interna, se ha modelado satisfactoriamente a través del presente algoritmo basado en el Método CVFEM, con la ventaja de que este último representa con mayor consistencia los fenómenos de transporte debido a sus características intrínsecas de conservación; sumado a la ventaja de su carácter no estructurado el cual le permite potencialmente abordar topologías complejas. 

Conclusiones 

Se puede concluir que el presente modelo basado en el Método CVFEM permite modelar satisfactoriamente el secado convencional de madera considerando el enfoque de Luikov y apoyado en la existencia de un potencial hídrico total, validado para casos de difusión bidimensional ortotrópica de humedad en madera de álamo. 

La predicción de distribuciones transitorias de humedad considerando la anisotropía del material permite observar regiones de marcados gradientes de concentración que pueden auxiliar un mejor análisis de las cualidades y consecuencias de un determinado secado. 

Por otro lado, los resultados del análisis de consistencia en relación al tamaño, forma de la malla y paso de tiempo, denotan un comportamiento consistente del algoritmo de cálculo desarrollado. 

Por último, no se observaron dificultades ni variaciones significativas en la curva de secado al utilizar como variable dependiente el contenido de humedad o potencial hídrico. 

Referencias 

1. Turner, I.; Mujumdar, A.S. “Mathematical modeling and numerical techniques in drying technology”, Marcel Dekker Inc., New York, ISBN 0-8247-9818-X. 1997.         [ Links ] 

2. Siau, J.F. “Wood: Influence of moisture on physical properties”. Virginia Tech. USA. 1995.         [ Links ] 

3. Lewis, R.W.; Morgan, K.; Thomas, H.R. “Drying induced stresses in porous bodies – An elastoviscoplastic model”. Computer Methods in Applied Mechanics an Engineering Vol. 20 (1979) 291-301. 

4. Morgan, K.; Thomas, H.R.; Lewis, R.W. “Numerical modeling of stress reversal in timber drying”. Wood Science Vol. 15 No 2 (1982)139-149. 

5. Chen, X. D. “Moisture Diffusivity in food and biological materials”. Drying Technology. Vol. 25 (2007) 1203-1213. 

6. Luikov, A.V. “Heat and mass transfer in capillary porous bodies”. Pergamon Press, Oxford. 1996.         [ Links ] 

7. Whitaker, S. “Simultaneous heat, mass, and momentum transfer in porous media: A theory of drying”. Advances in Heat Transfer Vol. 13 (1977)119- 203. 

8. Hernández, R.; Quinto P. “Secado en medios porosos: Una revisión a las teorías actualmente en uso. Cinética” Revista Científica. Vol. 9 No 2 (2005) 63-71. 

9. Cloutier, A.; Fortin, Y.; Dhatt, G. “A wood drying finite element model based on the water potential concept”. Drying Technology Vol. 10 No 5 (1992) 1151-1181. 

10. Cloutier, A.; Fortin, Y. “A model of moisture movement in wood based on water potential and the determination of the effective water conductivity”. Wood Sci. Technol. Vol. 27 (1993) 95-114. 

11. Cloutier, A.; Fortin, Y. “Wood drying modeling based on the water potential concept: Effect of the hysteresis in the M-y relationship”. Drying Technology Vol. 12 No 8 (1994)1793-1814. 

12. Tremblay, C.; Cloutier, A.; Grandjean, B. “Experimental determination of the ratio of vapor diffusion to the total water movement in wood during drying”. Wood Fiber Sci. Vol. 31 No 3 (1999) 235-248. 

13. Tremblay, C.; Cloutier, A.; Fortin, Y. “Determination of the effective water conductivity of red pine sapwood”. Wood Sci. Technol. Vol.34 No 2 (2000)109-124. 

14. Defo, M.; Cloutier, A.; Fortin, Y. “Modeling vacuum-contact drying of wood: The water potential approach”. Drying Technology. Vol.18 No 8 (2000) 1737-1778. 

15. Baliga, B.R.; Patankar, S.V. “A new finite element formulation for convection diffusion problems”. Numerical Heat Transfer Vol.3 (1980) 393-409. 

16. Salinas, C.; Ananias, R.A.; Alvear, M. “Simulación del secado convencional de la madera”. Maderas. Ciencia y tecnología, Vol. 6, No 1 (2004) 3-18.