Maderas.Ciencia y tecnologia 6(1): 3-18, 2004

ARTICULO

SIMULACIÓN DEL SECADO CONVENCIONAL DE LA MADERA
WOOD DRYING SIMULATION

Carlos SALINAS¹; Rubén A. ANANIAS² ; Mara ALVEAR^3
1Profesor Asistente, Depto. Ing. Mecánica, Universidad del Bío-Bío., Concepción-CHILE. casali@ubiobio.cl
²Profesor Asociado, Depto. Ingeniería en Maderas, Universidad del Bío-Bío, Concepción-CHILE. ananias@ubiobio.cl
³Magister en Ciencia y Tecnología de la Madera. Depto. Ingeniería en Maderas, Universidad del Bío-Bío, Concepción, CHILE. malvear@ubiobio.cl

Autor para correspondencia: casali@ubiobio.cl

Recibido : 26.08.2003. Aceptado: 06.05.2004.

---

RESUMEN

El presente trabajo dice relación con la modelación numérica bidimensional de la cinética del secado de madera basada en el concepto de potencial hídrico. El modelo matemático contempla ecuaciones de transporte del tipo diferencial parcial no lineal de segunda orden resuelto numéricamente a través del Método de Volumen Finitos en coordenadas generalizadas. Resultados de curvas de secado e isoconcentraciones son mostrados para madera de Populus sp. comparados con datos experimentales y/o publicados en la literatura especializada. Se concluye que el método numérico propuesto permite simular satisfactoriamente el secado de la madera obteniendo información transitoria detallada de las distribuciones de humedad al interior de ésta

Palabras Claves: Secado, Potencial hídrico, Simulación, Volúmenes Finitos

---

SUMMARY

The present work says relationship with the two-dimensional numerical model of the wood drying kinetics based on the water potential concept. The mathematical model is based on a set of non-linear second partial differential equations, which is solved by a numerical procedure using the Finite Volume Method in a generalized coordinates system. Results of drying curves and isoconcentrations are shown for Populus sp. wood and compared with experimental data. It has been concluded that the present numerical method allows to simulate the wood drying satisfactorily obtaining detailed information about the wood moisture content distributions.

Keywords: Drying, water potential, simulation, finite volumes method

---

Introducción

La madera es un material heterogéneo, anisotrópico, poroso y no saturado, sobre el cual pueden intervenir esfuerzos en diversas escalas, entre otros: Difusión a nivel molecular, efectos capilares al nivel de escalas intermedias o deformaciones macro escalas, como discutido por Turner y Mujumdar (1997).

Desde el punto de vista del secado, la humedad es su parámetro característico encontrándose el agua que la define en estado liquida y/o gaseosa (Perré et al.,1993), la cual puede estar libre en las cavidades celulares o ligada a las paredes celulares (Siau,1995). Además, en el proceso de secado, se alteran las relaciones y/o equilibrios de esfuerzo mecánicos, lo cual propicia la deformación del material (Lewis et al., 1979, Morgan et al., 1982)

Pese a las diversas complejidades, que devienen de las características de la madera anteriormente expuestas, el transporte de su humedad puede aceptar un tratamiento o modelación a nivel de macro escala asumiendo una homogeneidad (implica seleccionar una madera libre de defectos), despreciando las deformaciones, considerando la anisotropía pero tomando como relevantes la variación de propiedades en sus direcciones principales (ortotropía) y ponderando los fenómenos físicos de diversas escalas en cuanto a los efectos que se producen a grandes escalas (Turner, 1996). Todo lo anterior lleva a considerar, para efectos de modelación, a la madera como un material homogéneo, lo cual permite aplicar las ecuaciones desarrolladas para fenómenos de transporte bajo la hipótesis del continuo.

Por lo anteriormente expuesto, en el presente trabajo se estudia un modelo bidimensional para la simulación del proceso de secado de la madera considerando a ésta como un material poroso, ortotrópico y no deformable, lo cual implica el conocimiento de sus propiedades físicas por lo menos en dos de sus tres direcciones principales (radial, tangencial y longitudinal). La metodología de estudio contempla la caracterización física, la modelación matemática y la solución numérica del proceso de secado (pérdida de humedad) que sufre un trozo de madera sólida cuando es expuesto a una corriente de aire seco. Se trabajará con una topología simple (paralelepipedo de madera sólida), para la cual se disponen datos experimentales y numéricos obtenidos por otros autores. Además, se realiza un estudio sistemático en problemas de difusión unidimensional y bidimensional para efectos de validar los algoritmos al obtener resultados numéricos convergentes y concordantes con la evidencia analítica y/o experimental.

El modelo matemático implementado sigue la línea de investigación fundamentada por Luikov (1966) y desarrollada en el tiempo por diversos autores, entre otros: Liu y Cheng (1989), Gui et al. (1994) y Cloutier y Fortin (1994), Tremblay et al. (2000) y hasta llegar a la formulación simplificada de Defo et al. (2000). Procediendo de esta forma, se obtiene un sistema de ecuaciones diferenciales parciales no lineales de segundo orden para los procesos de transporte por conducción y convección de energía y concentración descritas para las variables de humedad, temperatura y presión. En particular, para el caso del secado convencional, que se estudia aquí, las variables de temperatura y presión pueden ser omitidas, restando entonces sólo la ecuación de transporte de concentración de humedad (Cloutier y Fortin, 1994).

Siendo así, esta ecuación de transporte del tipo diferencial parcial no lineal de segunda orden, es descrita numéricamente en términos del Método de Volúmenes Finitos (Patankar, 1980), en coordenadas generalizadas (Thompson et al., 1985 y Hirsch, 1990) con representación de segundo orden para términos difusivos, con formulación implícita para el avance en el tiempo y considerando “Update” linearización en los términos no lineales (Lapidus y Pinder, 1982). El sistema lineal de ecuaciones algebraicas resultante de la aplicación del método numérico es resuelto a través del esquema iterativo Gauss Seidel incluido factores de relajación. Las mallas son generadas en base a una ecuación de Poison como descrita en Thompson et al., 1985, con el algoritmo desarrollado en Salinas (1996) y en Moraga y Salinas (2000). Aplicaciones del modelo para la simulación de difusión isotrópica y ortotrópica unidimensional y bidimensional son realizados, siendo los resultados discutidos y comparados con datos experimentalmente o numéricos disponibles en la literatura.

Modelo Matemático

Considerando la variación local de la concentración de humedad equivalente a la divergencia del flujo se puede escribir, de acuerdo con el modelo de Luikov, la siguiente ecuación diferencial de transporte de concentración:

(1)

Donde:
C Concentración de humedad (kg_agua/m³_{madera-húmeda})
q_m flujo de humedad (kg_agua/m²_{madera-húmeda})

Suponiendo pequeñas variaciones de temperatura y equilibrio entre las fases del agua en la madera, el flujo de humedad es descrito en función de la conductividad y el potencial hídrico como:

(2)

De esta forma, reemplazando (2) en (1) y asumiendo la madera como un medio poroso indeformable se tiene:

(3)

Donde:
M Contenido de humedad (100 kg_agua/kg_madera-seca)
P Presión total (Pa)
T Temperatura (°K)
K Tensor de conductividad efectiva, función de M, T y P.(kg²_aguam s J)
Ψ Potencial hídrico (J/kg)
G_m Gravedad especifica (kg_madera-seca m³_agua/m³_{madera-húmeda} kg_agua)
ρ_w Densidad del agua (kg_agua/m³_agua)

Modelo Numérico

La ecuación de transporte de masa (1), puede ser representada por conveniencia de la siguiente forma:

(4)

Donde:
Φ Variable dependiente (Humedad)
c_Φ Capacidades(concentración de humedad)
q_Φ Flujo de Φ

operador de divergencia

S_Φ Fuente generadora de Φ

En particular, para el transporte de humedad Φ =M (porcentaje de contenido de humedad (kg_agua/kg_madera-seca)) se define el vector de flujo q_Φ y la capacidad c_Φ como:

=

(5)

(6)

Siendo k_ii(i=1,3) difusibilidades en las direcciones principales (radial, tangencial y longitudinal) y ∂M/∂ψ variación de humedad en relación al potencial, parámetros físicos de transporte a ser determinados experimentalmente (Cloutier y Fortin, 1993).

Por otra parte, integrando la ecuación (4) de acuerdo al Método de Volúmenes Finitos (MVF) resulta:

(7)

Asumiendo que el volumen finito esta conformado por seis caras, en las cuales el flujo a través de cada una de ellas es constante y considerando el valor medio de las integrales de volumen, se puede escribir la ecuación anterior, en coordenadas de cuerpo como (Salinas 1996), esto es:

(8)

Donde:

con

	base contravariante de vectores;
	jacobiano de la transformación.

De las relaciones de transformación (Thompson et al., 1985) se obtiene que la divergencia de un escalar (gradiente) es:

Por lo tanto, el flujo en cada cara del VF queda dado por:

donde

En particular, para el caso de coordenadas mutuamente ortogonales la expresión anterior para el flujo en cada cara puede ser substantivamente simplificada, esto es:

Evaluando Φ_ε en términos de diferencia central se puede escribir:

(9)

y definiendo

y

se obtiene:

(10)

luego representando el término temporal a través de un esquema explícito de Euler:

(11)

Donde:

Incorporando las definiciones de la ecuación (11) en (10) y omitiendo el superíndice n para mayor simplicidad:

(12)

Los coeficientes a_t,, a_ia y a_ir son dependientes de Φ lo cual incorpora la no linearidad en la ecuación anterior. Siendo así, se requiere de un esquema iterativo para su solución. En particular, si s indica un nivel iterativo, se puede plantear la siguiente ecuación para el ciclo iterativo s+1:

(13)

Esta es la ecuación de transporte Φ en términos numéricos del tipo implícita en el espacio y explícita en el tiempo a ser resuelta para cada VF considerado en el dominio discretizado lo que implica resolver un sistema de N_VF x N_VF ecuaciones, siendo N_VF el numero de volúmenes finitos. La expresión entrega el valor de la variable en el centroide P del VF (Φ_P) como función explícita de sus adyacentes inmediatos en cada dirección principal (Φ_iA;Φ_iR), lo cual conlleva una representación discreta a través de un esquema padrón de 7 puntos para el caso tridimensional, como es mostrado en la figura (1)

Figura 1. Disposición de la variable Φ

 

Resultados y Discusión

Tendientes a validar los algoritmos computacionales, se simulan tres problemas de transporte de diversa complejidad: Primero, difusión unidimensional isotrópica, la que sirve como prueba básica para el algoritmo, en cuanto refleja o no el comportamiento físico de difusión, corroborado con la solución analítica del problema planteado. En segundo lugar, se considera la solución de un problema bidimensional isotrópico, cuyos datos son comparados similares obtenidos por otros autores mediante diversos métodos. En ambos casos se incluyen efectos de convección externa y se realiza un análisis de convergencia de las variables utilizando diferentes mallas (10x10, 12x12, 20x20, 40x40, 60x60) a fin de determinar la sensibilidad del método. Y por último, se analiza un caso de secado convencional bidimensional ortotrópico de madera, aquí los resultados son descritos en relación con datos experimentales recogidos en la literatura (Cloutier et al 1992).

Difusión Unidimensional Isotrópica

El problema plantea la transferencia de humedad transitoria por difusión unidimensional en una cavidad de 0.5 x 1 (m²) (figura 2), considerando 0% de humedad en X₁=0, aislada en X₂ = 0 y X₂=H=0.5, y transferencia de humedad por convección en X₁=1. Se considera una capacidad c_Φ =1, un coeficiente de difusión isotrópico k=10 (m²/s), coeficiente de convección de materia h_Φ =100, esquematizado en la figura 2. Cuyas propiedades, condiciones iniciales y de contorno de acuerdo con la ecuación (4) para las diversas aplicaciones son resumidas en la tabla 1.

Figura 2. Difusión unidimensional isotrópica

Tabla 1. Propiedades1, condiciones iniciales y de contorno

1Valores de

son determinados experimental para cada especie. En particular son usados valores obtenidos por Cloutier et al. (1992).

2Unidades de las propiedades para el caso isotrópico 1D y ortotrópico 2D son dadas en la ecuación (3).

La figura 3 presenta los resultados en el centroide de la cavidad tendiente a analizar la sensibilidad del algoritmo con relación a la forma de malla. Se observa la similitud entre las curvas en relación con la cantidad de elementos, siendo convergentes las soluciones para mallas superiores 30x30. Esta aplicación sencilla, sirve al propósito de probar la funcionalidad y convergencia de los algoritmos de cálculo.

Figura 3. Análisis de convergencia para el caso de difusión unidimensional

Difusión Bidimensional Isotrópica

El problema dice relación con la transferencia de calor por difusión en una cavidad bidimensional de 0.01x0.02 m (figura 4)considerando una temperatura impuesta en los costados de la cavidad T_h=300 (°C), aislada en el fondo y expuesta convección forzada en su cara superior con coeficiente de convección h=200(W/m²). En la tabla 2 se comparan las temperaturas en los puntos indicados en la figura 4, obtenidos utilizando el presente algoritmo, aquellos entregados por el software ALGOR y los publicados por Holman (1989). Además se considera una temperatura ambiente T_amb=50, conductividad térmica k=3.0 (W/mºC), densidad ρ =1600 (kg/m³) y calor específico c_p=0.8(J/kgºC) (ver detalles en tabla 1)

Figura 4. Difusión bidimensional isotrópica (L=0.02 m; H=0.01 m)

La tabla 2 muestra resultados para el estado transitorio en las posiciones espaciales T_1, T₂ y T₃ obtenidos por Holman a través del método de Diferencias Finitas, por el software ALGOR que utiliza el método de elementos finitos y el presente algoritmo que es basado en el método de volúmenes finitos. Los datos son obtenidos de una malla uniforme 40x40 a tiempo de 12 s, tomando como referencia de comparación los resultados de Holman (1989), se aprecia que las diferencias de valores para los diversos métodos son inferiores al 5%, lo cual indica una buena concordancia.

La figura 5 muestra el análisis de sensibilidad para diversas mallas en el dominio bidimensional isotrópico, por el método Volúmenes Finitos. Se observa que entre 20x20 y 40x40 existe una mínima diferencia, llegando a interceptarse en determinado momento, lo que indica que existe convergencia para mallas uniformes iguales o superiores a 20x20 elementos.

La figura 6 (a,b y c), muestra resultados numéricos de distribución de temperatura en la mitad del dominio bidimensional isotrópico para diversos tiempos de evolución, los cuales son cualitativamente concordantes cuando comparados con los aportados por el software ALGOR (Algor Tutorial,versión 2002).

Tabla 2. Comparación de Temperaturas en los puntos T1, T2 y T3 para t=12 (s)

Posición	Holman (1981)	ALGOR		Programa
	T(°C)	T(°C)	Dif. %	T(°C)	Dif. %
T3	243.32	240	1.36	245.4	0.85
T2	279.87	270	3.52	276.0	1.4
T1	289.71	280	3.35	286.1	1.26

Figura 5. Análisis de convergencia de temperaturas en el centroide

a) 1 segundo

b) 5 segundos

c)10 segundos

Figura 6. Distribución de temperatura para diferentes tiempos

 

Difusión Bidimensional Ortotrópica

El planteamiento es dirigido a la transferencia de masa bidimensional, modelada en base al potencial hídrico en una pieza de madera sólida de 9x9 (cm²), de la especie Populus sp. expuesta a un flujo de humedad (q_m) por convección en sus caras exteriores. Naturalmente, el problema tiene una doble simetría (Ejes x₁ y x₂, ver figura 7), lo cual permite trabajar con un cuarto del dominio. Para el problema se considera un coeficiente convectivo de transferencia de masa de k_m= 9.36*10^-10 (kg².m²/s.J) (velocidad del aire ν=2.5 (m/s)), gravedad especνfica G_m=0.419 (1), densidad del agua =1000 (kg/m³) y potencial de equilibrio ψ_w=-119400 (J/kg) (implica una humedad de equilibrio ambiental CHE=9 (%) a una temperatura T_amb=20 (°C)) (Ver detalles en Tabla 1).

Figura 7. Esquema del problema (difusión bidimensional ortotrópica)

Los valores experimentales del coeficiente de difusión y potencial hídrico, recogidos de Cloutier et al. (1992), son presentados en figura (A1 – A3) del apéndice A.

La figura 8 muestra las curvas de secado experimental y numéricas obtenidas por volúmenes finitos para temperaturas de 20 (°C), con mallas de 10x10, 12x12, 20x20 y 40x40, para el secado convencional de Populus sp. se aprecia la influencia de la ampliación del numero de elementos en la curva de secado observándose la similitud de resultados entre mallas 20x20 y 40x40, lo que permite afirmar que para mallas iguales o superiores 20x20 elementos se obtiene una solución convergente.

Las figuras 9 (a, b y c) muestran la distribución espacial transitoria de humedad modelada. Se observa el carácter ortotrópico del material en la asimetría de la difusión de masa conforme las direcciones x₁ y x₂, existiendo zonas de marcados gradientes de concentración (región inferior de la Figuras 9 a y b), lo cual puede dar indicios para apoyar estudios de defectos en la madera.

Figura 8. Curvas de secado experimental y numéricas (20 (°C) dt=36 (s), malla 40x40)

Comparando estos resultados con los obtenidos experimentalmente y numéricamente (Método de Elementos Finitos) por Cloutier et al (1992), 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 de Volúmenes Finitos, con la ventaja este último de representar con mayor consistencia fenómenos de transporte debido a sus características intrínsecas de conservatividad.

a) 1 (h)

b) 10 (h)

c) 40 (h)

Figura 9. Curvas numéricas de isoconcentraciones (20 (°C) Ñt=36 (s) malla 40x40).

Conclusiones

Se puede concluir que la aplicación del Método Volúmenes Finitos permite modelar satisfactoriamente el secado convencional de madera considerando el enfoque de Luikov y apoyado en la existencia de un potencial hídrico total. En particular, validado para casos de difusión bidimensional ortotrópica de humedad en madera de Populus sp. Además, la predicción de distribuciones transitorias de humedad considerando la anisotropía del material, permite observar regiones de marcados gradientes de concentraciones que favorecen un mejor análisis de las cualidades y consecuencias de un determinado secado.

AGRADECIMIENTOS
Agradecemos el financiamiento de la Dirección de Investigación de la Universidad del Bío-Bío a través del proyecto 031110 3/R.

REFERENCIAS

CLOUTIER, A.; FORTIN, Y.; DHATT, G. 1992. A wood drying finite element model based on the water potential concept. Drying Technology 10(5): 1151-1181.

CLOUTIER, A.; FORTIN, Y. 1993. A model of moisture movement in wood based on water potential and the determination of the effective water conductivity. Wood Sci. Technol. 27: 95-114.

CLOUTIER, A.; FORTIN, Y. 1994. Wood drying modeling based on the water potential concept: Effect of the hysteresis in the M-y relationship. Drying Technology 12(8):1793-1814.

DEFO, M.; CLOUTIER, A.; FORTIN, Y. 2000. Modeling vacuum-contact drying of wood: The water potential approach. Drying Technology 18(8): 1737-1778.

GUI, Y.Q.; JONES, E.W.; TAYLOR, F.W.; ISSA, C.A. 1994. An applications of finite element analysis to wood drying. Wood Fiber Sci, 26(2): 281-293.

HIRSCH, C. 1990. Numerical computation of internal and external flows. Vol. 2, chap. 2, John Wiley & Sonc, Inc., New York.

HOLMAN, J.P. 1989. Heat transfer. McGraw-Hill, N.Y, USA, SI Metric Edition, ISBN 0-07-100487-4, pag. 173.

LAPIDUS, L.; PINDER, G.F. 1982. Numerical solution of partial differential equations in science and engineering. John Wiley & Sons, Inc[ STANDARDIZEDENDPARAG]

LEWIS, R.W.; MORGAN, K.; THOMAS, H.R. 1979. Drying induced stresses in porous bodies – An elastoviscoplastic model. Computer Methods in Applied Mechanics an Engineering, 20: 291-301.

LIU, J. Y.; CHENG, S. 1989. Solution of Luikov equations of heat and mass transfer in capillary porous bodies. Int. J. Heat and Mass Transfer 34(7):1745-1754.

LUIKOV, A.V. 1966. Heat and mass transfer in capillary porous bodies. Pergamon Press, Oxford.

MORAGA, N.O.; SALINAS, C. H. 2000. Numerical study of unsteady 2d natural convection and solidification of a food inside a freezing chamber. Numerical Heat Transfer, Part A, 37: 755-777.

MORGAN, K.; THOMAS, H.R.; LEWIS, R.W. 1982. Numerical modeling of stress reversal in timber drying. Wood Science, 15 (2): 139-149.

PATANKAR, S.V. 1980. Numerical heat transfer and fluid flow, Hemisphere Publishing Corporation, Washington, DC.

PERRE, P.; MOSER, M.; MARTIN, M. 1993. Advances in transport phenomena during convective drying with superheated steam and moist air Int. J. Heat and Mass Transfer 36(11):2725-2746.

SALINAS, C. H. 1996. Modelação de escoamentos tridimensionais em geometrias complexas, D. Sc. Thesis, COPPE/PENO, UFRJ, RJ-Brasil.

SIAU, J.F. 1995. Wood: Influence of moisture on physical properties. Virginia Tech. USA.

THOMPSON, J.F.; WARSI, Z.U.A.; MASTIN, C.W. 1985. Numerical grid generation, Elsevier Science. Publishing.

TREMBLAY, C.; CLOUTIER, A.; GRANDJEAN, B. 1999. Experimental determination of the ratio of vapor diffusion to the total water movement in wood during drying. Wood Fiber Sci. 31(3):235-248.

TREMBLAY, C.; CLOUTIER, A.; FORTIN, Y. 2000. Determination of the effective water conductivity of red pine sapwood. Wood Sci. Technol. 34(2): 109-124.

TURNER, I.; MUJUMDAR, A.S. 1997. Mathematical modeling and numerical techniques in drying technology, Marcel Dekker Inc., New York, ISBN 0-8247-9818-X.

TURNER, I.W. 1996. A two dimensional orthotropic model for simulating wood drying process, Applied Mathematical Modeling, 20: 60-81.

APENDICE A: Datos experimentales de potencial y conductividad a 20 (°C)
(Fuente: Cloutier et al. 1992)

Figura A1: Potencial hídrico

Figura A2: Conductividad radial

Figura A3: Conductividad tangencial.