Carlos Salinas¹, Cristian Chávez², Yerko Gatica³, Rubén Ananias⁴

¹ Ph. D., en Ciencias. Profesor Asociado, Departamento de Ingeniería Mecánica, Universidad del Bio-Bio, Chile. casali@ubiobio.cl

² Candidato a M.Sc. en Ciencia y Tecnología de la Madera, Universidad del Bio-Bio, Chile. crchavez@alumnos.ubiobio.cl

³ Candidato a Ph.D. en Ciencia & Industria de la Madera, Universidad del Bio-Bio, Chile. ygatica@ubiobio.cl

⁴ Ph. D., en Ciencias de la Madera. Professor Asociado, Departamento de Ingeniería en Madera, Universidad del Bio-Bio, Chile. ananias@ubiobio.cl

RESUMEN

El objetivo de este trabajo es el de simular esfuerzos bidimensionales de secado usando el método de volúmenes de control conformado por elementos finitos (CVFEM). Se modelan esfuerzos de secado en una sección transversal de madera sólida, producidos por gradientes de humedad asociados a los fenómenos de contracción y sorción mecánica. En particular, se implementa CVFEM con elementos finitos triangulares, interpolación lineal de la variable independiente al interior de él y programado en lenguaje Fortran 90. El modelo se valida contrastando resultados con trabajos similares disponibles en la literatura especializada. Finalmente, se muestran resultados originales de la modelación aplicada al secado isotérmico (20 ºC) de madera sólida de álamo (Populus tremuloides): distribución bidimensional de esfuerzo/deformación y contenidos de humedad para tiempos de secado igual a 40, 80, 130, 190 y 260 (h) y evolución transitoria de esfuerzos normales ( ^{-2,5<σ_xx<1.2, MPa} ), desde el centro a la superficie de la madera.

Palabras clave: simulación, secado, madera, esfuerzos, creep, CVFEM.

ABSTRACT

The work was aimed at simulating two-dimensional wood drying stress using the control-volume finite element method (CVFEM). Stress/strain was modeled by moisture content gradients regarding shrinkage and mechanical sorption in a cross-section of wood. CVFEM was implemented with triangular finite elements and lineal interpolation of the independent variable which were programmed in Fortran 90 language. The model was validated by contrasting results with similar ones available in the specialised literature. The present model’s results came from isothermal (20ºC) drying of quaking aspen (Populus tremuloides): two-dimensional distribution of stress/strain and water content, 40, 80, 130, 190 and 260 hour drying time and evolution of normal stress ( ^{-2,5<σ_xx<1.2, MPa} ), from the interior to the exterior of wood.

Keywords: simulation, drying, wood, stress, creep, CVFEM.

Recibido: agosto 21 de 2009. Aceptado: febrero 10 de 2011

Introducción

El trabajo se enmarca dentro del área de la modelación y simulación del secado de madera sólida centrado en el fenómeno de transferencia de masa (sorción) y los esfuerzos que dicho proceso ocasiona (creep). Una colección de trabajos relevante en esta área es dada en Turner y Mujumdar (1997). En particular, se abordan dos grandes aspectos: 1) el transporte de humedad, modelado sobre la base del concepto de potencial hídrico, ampliamente discutido en trabajo previo (Salinas et al., 2004); 2) las tensiones de secado asociadas a los gradientes de humedad, objetivo principal del trabajo.

La teoría señala que la contracción que experimenta la madera durante el secado induce el desarrollo de esfuerzos mecánicos a través de su estructura (McMillen, 1963). Estos esfuerzos son afectados por los cambios de humedad, las restricciones mecánicas, la anisotropía y el comportamiento viscoelástico de la madera, entre otros factores (Young, 1957). Se observa que los cambios de humedad favorecen el desarrollo de esfuerzos hidromecánicos, tanto o más importantes que la deformación plástica (Keey et al., 2000).

En este sentido, se han propuesto modelos para la madera enfocados principalmente a la deformación provocada por el transporte de masa (humedad): Perre et al. (1993), Chen et al.(1997), Stevensson and Martensson (2002), Kang et al. (2004) y Pang (2007), han propuestos modelos unidimensionales para determinar los esfuerzos desarrollados durante el secado debido principalmente a los fenómenos de contracción y sorción mecánica. De la misma forma, modelos bidimensionales para la deformación han sido propuestos por autores como: Turner y Ferguson (1995), Martensson y Stevensson (1996), Ferguson (1998), Perre y Turner (2001), Thuvander et al. (2002), Cheng et al. (2007); como también tridimensionales, realizados por Ormarsson et al. (2003), Remond et al. (2006), Keuneke y Niemz(2008), entre otros.

En este trabajo se modelan las deformaciones recurrentes de los gradientes de humedad durante el proceso de secado sobre la base del trabajo hecho por Turner y Ferguson (1995). Desde el punto de vista numérico, se postula integrar las ecuaciones de transporte y de tensión/deformación usando el método de volumen de control conformado por elementos finitos (Control Volume Finite Element Method, CVFEM) (Baliga y Patankar, 1983). En particular, se efectúan variaciones de parámetros geométricos y físicos para el análisis de estabilidad y consistencia de los algoritmos. Dichos algoritmos son validados comparando resultados con trabajos similares disponibles en la literatura y se aplican al secado isotérmico de madera sólida de álamo (Populus tremuloides), generando datos numéricos originales de distribución bidimensional transitoria de esfuerzo/deformación.

Modelo físico

Se presenta un modelo para el fenómeno físico esfuerzo/ deformaciones ocasionado en una pieza de madera sólida, la cual es sometida a un proceso de secado. Se modelan los efectos transitorios no uniforme inducidos por la variación del contenido de humedad (M), esto es: esfuerzos (s_ij), deformaciones (e_ij) y desplazamientos (u = (u, v) o u_i con (i = 1,2)).

En particular, se estudia una sección transversal bidimensional de madera en el plano radial tangencial, esquematizada en la figura 1, cuyas propiedades son dadas en la tabla 1. Se asume un proceso isotérmico para el transporte de ortotrópico de humedad durante el proceso de secado y una deformación plana para el problema esfuerzo-deformación, producto de una deformación libre (contracción) y del fenómeno de sorción mecánica, ocasionados por gradientes de humedad bajo el punto de saturación de la fibra (PSF).

Las condiciones iníciales y de contorno son: a) para el problema de transporte de humedad se considera una humedad inicial M= M_ini, con condiciones de contorno del tipo Neumann de 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; b) para el problema tensión/deformación se asume un estado inicial no deformado (σ_ij = ϒ_ij = u_i = 0), con condiciones de contorno del tipo Dirichlet en los ejes de simetría (u = 0 en x = L/2 y v = 0 en y = 0) y libre en las superficies x = 0 e y = H/2.

El modelo reológico está inspirado en los trabajos desarrollados por Ferguson (1997 y 1998), donde el creep es descrito como:

(1.1)

Lo cual implica que está compuesto por deformaciones del tipo: elástica ( ), contracción libre ( ) y sorción mecánica .

Modelo matemático

El modelo matemático de transporte de humedad usado en este trabajo es basado en el concepto de potencial hídrico ( ψ), según modelo propuesto en Cloutier y Fortin (1994) y en la forma numérica implementada en Salinas et al. (2004). Lo anterior, tomando como variable independiente el contenido de humedad (M), resumido en la siguiente expresión:

(1.2)

donde conductividades de M capacidad .

donde C_M ,K_xx y K_yy dependen del contenido de humedad M y la temperatura T.

El modelo matemático para la deformación bidimensional resultado del equilibrio mecánico aplicado a cada punto, en ausencia de fuerza de cuerpo, puede ser descrito como:

(1.3)

donde .

Producto de la distribución no uniforme de humedad, se producen esfuerzos inducidos por la deformación libre y por sorción mecánica. La modelación de dichos esfuerzos es realizada por una función implícita de cinco parámetros: deformación inicial , contracción ∝ &sorción mecánica m, esfuerzo σ, variación de concentración ΔC y razón de Poisson v:

(1.4)

Siendo así, los esfuerzos son determinados a partir de siguiente expresión:

(1.5)

donde σ^o representa los valores de los esfuerzos iníciales y D es la matriz constitutiva.

(1.6)

El problema es integrado en el tiempo de acuerdo al método de deformación inicial (Zienkiewicz y Taylor, 2000). Para esto, se define ΔC como la variación de humedad entre el tiempo y , donde dt es el intervalo de tiempo considerado. Esto es: y de forma similar y . Siendo así:

(1.7)

Modelo numérico

En este ítem se mostrará en detalle el modelo numérico aplicado a la ecuación de desplazamientos directamente relacionado con el problema esfuerzo/deformación. Una discusión similar fue realizada para el problema de transporte de humedad, ecuación, en Salinas et al. (2011).

Considerando el dominio ∀conformado por nl volúmenes finitos (VF) ( ) de contorno Ω^l y de forma similar, cada VF es conformado por nk contribuciones parciales de elementos finitos (EF) de contorno (figura 2), se puede escribir:

(1.7)

En consecuencia, la integral de (1.3) en los subdominios , de acuerdo con el teorema de Green, es igual a la divergencia del flujo en el contorno Ω¹. Esto es:

(1.9)

donde:

σ_i es el esfuerzo en la dirección i = 1,2 (x e y, respectivamente) y el vector unitario normal externo a Ω¹.

cuya solución conjunta (sistema de ecuaciones) permite obtener valores discretos de la variable independientes involucradas (desplazamientos u).

Cada VF está formado por contribuciones parciales de ^nk EF, cual es mostrado en la figura 2. Siendo así, la contribución parcial del EF k al VF centrado en el nodo local 1 (global l) es igual a:

(1.10)

donde los subíndices a y b representan los segmentos de contorno, respectivamente.

Si consideramos el valor de σ_i en el punto medio de cada segmento de contorno igual a como predominante en el mencionado segmento, e incorporando las definiciones para σ_i dadas en (1.3), se puede escribir:

(1.11)

donde Por otro lado, redefiniendo las tensiones σ_i en funciones de las deformaciones lineales ^ε_ii y angulares ϒ_ij como:

(1.12)

donde son esfuerzos por deformación inicial.

siendo las deformaciones ^ε_ii y ϒ_ij descritas en función del desplazamiento como:

(1.13)

considerando, para la integración anterior, una variación lineal ^u al interior del EF, es decir:

(1.14)

donde ^A,B,C son constantes definidas en función de los desplazamientos nodales ^{u1,u2 y u3} , iguales a:

se tiene que la variación de los desplazamientos son iguales a:

(1.15)

Introduciendo las definiciones dadas por las ecuaciones (1.12) a (1.15), la contribución puede ser expresada en función de los valores nodales de desplazamiento como:

(1.16)

Una expresión similar se puede obtener para .

Expresiones análogas a las contribuciones del segmento iguales a , dada por (1.16), de los segmentos y iguales a y pueden ser obtenidas simplemente intercambiando el vector normal ⁿ_a por ⁿ_b o , ⁿ_c según corresponda.

Ahora, recuperando las definiciones dadas en (1.10) e incorporando los desarrollos de (1.16), se tiene:

(1.17)

(1.18)

Siendo así, se puede escribir el siguiente sistema de ecuaciones:

De esta forma, para cada valor l se obtendrán dos ecuaciones algebraicas, una para cada componente del desplazamiento medio , en cada subdominio ∀^l, lo cual configura un sistema lineal de 2nl x 2nl ecuaciones de la forma

Resultados

A continuación se presentan resultados de dos simulaciones: 1) difusión de calor y tensiones térmicas en una placa de acero, realizada para efectos de validación; y 2) transporte de humedad y tensiones de secado en madera de álamo (Populus tremuloides).

Difusión de calor y tensiones térmicas en una placa de acero.

Se considera la simulación transitoria del esfuerzo plano en una placa bidimensional que experimenta un calentamiento transitorio no uniforme, esquematizada en la figura 3, cuyas propiedades son dadas en la tabla 2. Sus magnitudes son: largo L = 0,6 (m) y altura H = 0,2 (m). Se considera para efectos de condiciones de contorno: desplazamiento restringido en el sentido Y en el tramo AD (v = 0 en y = 0), restringido en el sentido X en el tramo BC (u = 0 en x = L) y libre en los tramos CD y DA. Naturalmente, el punto B está restringido en X e Y ((u,v) = (0,0) en B). Al mencionado problema se le conoce solución analítica dada por Boley y Weiner (1960).

La distribución transitoria no uniforme de temperatura T (ºC), de forma exponencial en función de la coordenada y (m) y el tiempo t(h), es dada por:

(1.19)

Producto de la distribución no uniforme de temperaturas se producen esfuerzos inducidos por la dilatación y también debidos a esfuerzos sostenidos en el tiempo (creep). La modelación de dichos esfuerzos puede ser realizada de forma equivalente al modelo de tensiones presentado para la variación de humedad. En particular, la deformación inicial ε⁰ es función de los siguientes parámetros: dilatación térmica ∝, coeficiente de creep m, esfuerzo σ, variación de temperatura ΔT y, para el caso de deformación plana, la razón de Poisson v.

La figura 4 muestra mapeamientos del estado de deformación y tensiones para t = 1 (h) (estado permanente). En dichos mapeamientos se puede observar la consistencia en cuanto a la imposición de condiciones de contorno (desplazamientos dados en la figura 4a). Los efectos de las restricciones a la deformación libre en el desarrollo de tensiones normales pueden ser apreciados en la figura 4b.

El análisis de convergencia y consistencia, con relación al tipo de malla, contrastado con la solución analítica, es mostrado en la figura 5. Se aprecia que resultados transitorios convergentes pueden ser logrados con mallas uniformes de 52 x 20.

Por último, la figura 6 muestra el análisis de consistencia en cuanto al paso del tiempo versus solución convergente en el tramo BC. Se aprecia una buena convergencia para pasos de tiempo de 0,001 (h).

Transporte de humedad y tensiones de secado en madera de álamo.

Se considera la simulación transitoria no uniforme del transporte de humedad en la madera, que induce esfuerzos y deformaciones. Para esto se determina en cada tiempo de integración la distribución de humedades, según lo descrito en Salinas et al. (2004), lo cual permite posteriormente calcular las deformaciones libres (ε⁰ ) que motivan desplazamiento y esfuerzos, calculados estos últimos en forma desacoplada, de acuerdo al método de deformación inicial (Zienkiewicz y Taylor, 2000). Análisis de consistencia revelaron que mallas con refinamiento exponencial hacia superficies de convección permiten captar gradientes de humedad con un menor número de elementos discretos.

La figura 7 muestra distribuciones espaciales para diversos tiempos de secado de los parámetros de cálculo del problema esfuerzo/ deformación. En la figura 7a se pueden apreciar las distribuciones transitorias y espaciales del contenido de humedad, donde se destacan sus mayores gradientes (concentración de isoconcentraciones) en alrededor del 60% de contenido de humedad (CH) y el comportamiento ortotrópico del secado, reflejado en la falta de simetría de los frentes de secado. La figura 7b muestra desplazamientos en la dirección radial (u), apreciándose la correcta implementación de las condiciones de contorno y concentrándose los mayores desplazamientos en los bordes restringidos en un grado de libertad.

La figura 8a presenta en detalle los esfuerzos normales residuales de secado en la dirección radial. Se aprecia cómo se concentran los esfuerzos principales en A, B y C (figura 1), estando en tensión en A y C (superficie) y compresión en B (centro). Por otro lado, los esfuerzos de corte mostrados en la figura 8b se concentran en la diagonal levemente desplazado hacia el punto D (superficie).

La figura 9 destaca en detalle la evolución transitoria de las tensiones normales radiales en el tramo BC, donde se puede apreciar su compleja evolución conforme las etapas de las esfuerzos de secado simuladas. La figura 9a muestra la primera etapa de secado, caracterizada por el desarrollo creciente de tensiones en la superficie de la madera en contrapartida a la compresión del centro de ésta. Como complemento, la figura 9b expone la se gunda etapa de esfuerzos, caracterizada por una reversión creciente de los esfuerzos de tensión por compresión en la superficie y viceversa en el centro de la madera.

La figura 0 despliega las evoluciones transitorias de las tensiones en el centro y en la superficie del dominio de cálculo, puntos B y C de la figura 1. En ella se puede apreciar la dinámica de la variación de las intensidades del esfuerzo normal radial correlacionado con el contenido de humedad. Básicamente, al inicio del secado, en el superpie y en el centro, se nota una marcada tensión y leve compresión, respectivamente. Esto contrasta con el estado final de esfuerzos, al final del proceso de secado, en el cual se invierten los esfuerzos normales. Esto último es una de las características cualitativas relevantes del problema físico simulado. Se destaca la marcada singularidad transitoria en el desarrollo de esfuerzos para t = 195 (h), donde CH = 10%.

Conclusiones

Los resultados presentados en las figuras 4-7 permiten concluir que el algoritmo desarrollado genera simulaciones efectivas del fenómeno esfuerzo/deformación ocasionado por deformaciones libres y esfuerzos sostenidos en el tiempo. Los análisis de consistencia realizados, en cuanto al tamaño de la malla y el tiempo, permiten inferir que simulaciones convergentes pueden ser logradas con malla del orden de 20 x 52 y pasos de tiempo de 0,001 (h).

Los resultados presentados en las figuras 8-10 permiten concluir que se han simulado efectivamente los característicos presentes en el secado de madera, esto es: regiones sometidas indistintamente a tensión y compresión durante el proceso transitorio de secado, resultando esfuerzos residuales de tensión en la superficie y compresión en el centro de la madera. Además, aceptada la validación del modelo en base al análisis realizado a la placa térmica, los valores cuantitativos simulados deberían ser muy similares a los físicos. Naturalmente que sería deseable contar con datos experimentales para corroborar resultados numéricos, lo cual es dejado para futuros trabajos.

Referencias

Baliga B. R., Patankar S. V., A new finite element formulation for convection diffusion problems., Numerical Heat Transfer, Vol. 3, 1980, pp. 393-409.

Boley B. A., Weiner J. H, Theory of Thermal Stresses., John Wiley and Sons, Inc., New York, 1960.

Chen G, Keey, R.B, Walker J.F.C., The drying stress and check development on high-temperature kiln seasoning of sapwood Pinus radiate boards., Holz als Roh-und Werkstoff, Vol. 55, Nº 2, 1997, pp. 59-64.

Cheng, W., Morooka, T., Wu, Q., Liu, Y., Characterization of tangential shrinkage stresses of wood during drying under heated steam above 100 ºC., Forest Prod. J, Vol. 57, Nº 11, 2007, pp. 39-43.

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, pp. 1151-1181.

Cloutier A., Fortin Y., Wood drying modeling based on the water potential concept: Effect of the hysteresis in the M-y relationship, Drying Tech., Vol. 12, Nº 8, 1994, pp. 1793-1814.

Ferguson, W.J., A control volume finite element numerical solution of creep problems., Int. J. Num. Meth. Eng, Vol. 40, No. 18, 1997, pp. 3463-3475

Ferguson, W.J., A numerical prediction of the effect of airflow and wet bulb temperature on the stress development during convective wood drying., In mathematical modeling and numerical techniques in drying technology, Ed. I. Turner, A.S. Mujumdar, 1997, pp. 259-277.

Ferguson, W. J., The control volume finite element numerical solution technique applied to creep in softwoods., Int. J. Solid Structures, Vol. 35, No. 13, 1998, pp.1325-1338.

Kang, W., Lee, N., Jung, H., Simple analytical methods to predict one-an two-dimensional drying stresses and deformations in lumber., Wood Sci. Technol., Vol. 38, Nº 6, 2004, pp. 417- 428.

Keey, R., Langrish, T., Walker, J., Kiln-drying of lumber., Springer-verlag, N.Y. 2000.

Keunecke, D.S.H, Niemz, P., Three-dimensional elastic behaviour of common yew and Norway spruce., Wood Science and Technology, Vol. 42, Nº 8, 2008, pp. 633-647.

Mcmillen, J., Stresses in wood during drying., Res. Rap. 1652 USDA Forest Service. Forest Products Lab. Madison, WI., 1963.

Martensson, A., Stevensson, S., Application of a material model describing drying stresses in Word., 5th IUFRO International wood drying conference, 1996, pp. 93-102.

Ormarsson, S., Cown, D., Dahlblom, O., Finite element simulations de moisture related distortion in laminated timber products of norway spruce and radiata pine., 8th IUFRO International wood drying conference, 2003, pp. 27-33.

Pang, S., Mathematical modeling of kiln drying of softwood timber: Model development, validation and practical application., Drying Technology, Vol. 25, 2007, pp. 421-431.

Perre, P., Moser M., Matin, M., Advances in transport phenomena during convective drying with superheated steam and moist air., Int. J. Heat and Mass Transfer, Vol 36, Nº. 11, 1993, pp. 2725-2746.

Perre P., Turner, I., Determination of the material property variations across the growth ring of softwood for use in a heterogeneous drying model., Holzforschung, Vol. 55, Nº. 4, 2001, pp. 417-425.

Remond, R., Passard, J., Perre, P., The effect of temperature and moisture content on the mechanical behaviours of wood: a compressive model applied to drying and bending., European Journal of Mechanical Solid, Vol. 26, 2006, pp 558- 575.

Salinas, C., Ananias, R. A., Alvear, M., Simulación del secado convencional de la madera., Maderas Ciencia y Tecnología, Vol. 6, Nº. 1, 2004, pp. 3-18.

Salinas, C., Chávez, C., Gatica, Y., Ananias, R., Simulación del secado convencional de madera usando CVFEM., Revista Técnica de Ingeniería, Universidad de Zulia, vol. 34. Nº1, abril 2011.

Stevensson, S., Martensson, A., Simulation of drying stresses in wood Part III Convective air drying of sawn timber., Holz Roh Werkst, Vol. 60, 2002, pp. 72-80.

Thuvander, F., Kifetew, G., Berglund, L.A., Modeling of cell wall drying stresses in wood., Wood Sci. Technol., Vol. 36, 2002, pp.241-254.

Turner, I.W., Ferguson, W., An unstructured mesh cell-centered control volume method for simulating heat and mass transfer in porous media: Application to softwood drying, part I: The isotropic model., Appl. Math. Modeling, Vol. 19, 1995a, pp. 654-667.

Turner, I.W., Ferguson W., An unstructured mesh cell-centered control volume method for simulating heat and mass transfer in porous media: Application to softwood drying. Part II. The anisotropic model., Appl. Math. Modeling, Vol. 20, 1995b, pp. 669-674.

Turner, I., Mujumdar, A.S., Mathematical modeling and numerical techniques in drying technology., Marcel Dekker Inc., New York, ISBN 0-8247-9818-X, 1997, pp. 1-20.

USDA., Wood Handbook: Wood as un engineering material., report: FPL-GTR-190, 2010.

Young, R. The perpendicular to grain mechanical properties of red oak as related to temperature, moisture content and time., USDA, FPL-2079, Madison, USA,1957.

Zienkiewicz, O.C., Taylor R.L., The Finite Element Method., Fifth edition, published by Butterworth-Heinemann, Vol. 2. Solid Mechanics, 2000, pp. 365.