Two-dimensional wood drying stress simulation using control-volume mixed finite element methods (CVFEM)

Simulación bidimensional de esfuerzos de secado en la madera usando CVFEM

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.

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.

• HTML (es)
• HTML (en)

Simulación bidimensional de esfuerzos de secado en la madera fuerzos usando CVFEM

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

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

ABSTRACT

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.

---

Two Two-dimensional wood drying stress simulation using control control-volume mixed finite element methods ( CVFEM)

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

¹Ph. D., in Science. Associate Professor, Department of Mechanical Engineering, Universidad del Bio-Bio, Chile. casali@ubiobio.cl

² M.Sc. Candidate in Science and Technology of Wood, Universidad del Bio-Bio, Chile. crchavez@alumnos.ubiobio.cl

³Ph.D. Candidate Science & Industry Wood, Universidad del Bio-Bio, Chile. ygatica@ubiobio.cl

⁴Ph. D. in Wood Science. Associate Professor, Department of Wood Engineering, Universidad del BioBio, Chile. ananias@ubiobio.cl

---

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 isothermp><b>al (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.

Received: August 21th 2009. Accepted: Feuary 10th 2011

---

Introduction

The work describes modelling and simulating solid wood drying focused on mass transfer (sorption) and stresses caused by drying. A collection of relevant work was found in Turner and Mujumdar (1997). The paper deals with two major aspects: moisture transport, modelled based on the water potential concept which has been widely discussed in previous work (Salinas et al., 2004 and 2011) and drying stress distribution associated with moisture content (this work’s main objective).

Theory says that wood shrinkage during drying induces mechanical stress developing throughout the structure (McMillen, 1963); such stress is affected by changes in moisture content, mechanical restraint, anisotropy and the viscoelastic behaviour of wood (Young, 1957). It has been noted that moisture gradients have favour the development of hydro-mechanicl stress, this bing as much or more important than plastic deformation (Keey et al., 2000).

The proposed model has thus been focused primarily on strain caused by mass transport (moisture content) in wood. Perre et al., (1993), Chen et al., (1997), Stevensson and Martensson (2002), Kang et al., (2004) and Pang (2007) have all proposed one-dimensional models for determining the drying stress caused by shrinkage and mechanical sorption. Likewise, twodimensional models for the strains have been proposed by authors such as Turner and Ferguson (1995), Stevensson and Martensson (1996), Ferguson (1998), Perre and Turner (2001), Thuvander et al., (2002), Cheng et al., (2007), as well as threedimensional models by Ormarsson et al., (2003), Remond et al., (2006), and Niemz Keuneke (2008).

This paper studied the strain and stress caused by moisture gradients during drying based on work by Turner Ferguson (1995). It was thus proposed to integrate the partial differential equation model using the control-volume finite element method (CVFEM) developed by Baliga Patankar (1983). In particular, changes were made to the geometrical and physical parameters for analysing stability and consistency. These algorithms were validated by comparing this study’s results with similar studies in the pertinent literature applied to the isothermal drying of solid quaking aspen wood (Populus tremuloides). Numerical data for transient twodimensional distribution of stress/strain is given.

Physical model

A model for the strain/stress produced when drying solid wood is presented. Non-uniform transient effects induced by variation of moisture content (M) was modelled: stress (s_ij), strain (e_ij) and displacement (u = (u, v) or u_i con (i = 1,2)).

A two-dimensional cross-section of wood was studied in transverse plane (Figure 1); its properties are given in Table 1. An isothermal process was assumed for the orthotropic transport of moisture during drying and plane strain for strain/stress, the product of free deformation (shrinkage) and mechanical sorption caused by moisture gradients below fie saturation point (FSP).

Initial and boundary conditions were as follows. An initial moisture M= M_ini, with Neumann-type boundary conditions with non-flow on the symmetry axis (x = L/2 and y = 0) and convection on surfaces x = 0 and y = H/2 were considered for moisture transport. A non-deformed initial state (s_i =g_ij=u_i =0) was assumed with Dirichlet boundary conditions on the symmetry axis (u=0 at x=L/2 and v=0 at y=0) and free on surfaces x=0 and y=H/2 for tstrain/stress.

The rheological model was based on work performed by Ferguson (1997 1998) where creep was described as:

(1.1)

leading to the next kind of deformation: elastic ( ), free shrinkage ( ) and mechanical sorption ( ).

Mathematical model

The mathematical moisture transport model used in this work was based on the water potential concept (ψ), like the model proposed by Cloutier Fortin (1994) and implemented in numerical form in Salinas et al., (2004). Taking moisture content (M) as independent variable in the above can be summarised in the following expression:

(1.2)

where conductivities of M, capacity

C_M ,K_xx y K_yy were functions of moisture content M and temperature T.

The mathematical model for two-dimensional deformation resulting from applied mechanical equiliium at every point in the absence of body force can be described as:

(1.3)

where

The products of non-uniform distribution of moisture were induced stress due free deformation and mechanical sorption. Such stresses was modelled by five implicit-function parameters: initial strain , shrinkage a, mechanical sorption m, stress s, change in concentration DC and Poisson ratio n:

(1.4)

Thus, stresses were determined from the following expression:

(1.5)

where σ^o represented initial stress data and D the constitutive matrix.

(1.6)

The problem was integrated using the Initial Strain Method, Zienkiewicz Taylor (2000). DC was thus defined as moisture variation between time and , , where dt was the time interval considered, i.e. and similarly and . Thus:

(1.7)

Numerical model

The numerical model applied to the displacement equation for the problem directly related to stress/strain is shown below; a similar procedure has been performed for the moisture transport equation (1.2) given in Salinas et al., (2011).

Considering that domain ∀ formed nl finite volume (FV)( ) contour Ω^l and similarly, each VF was composed by nk partial contributions of finite element (FE) with boundary (see Figure 2), then

(1.7)

Consequently, the integral of (1.3) in subdomain , according to Green's theorem, was equal to flux divergence on boundary Ω¹, i.e.

(1.9)

Where

were stresses with i=1,2: stresses in x and y direction, respectively and the external normal unit vector for Ω¹.

Joint solution (system of equations) led to discrete values for the independent variable involved (displacement u).

Each VF was formed by partial contributions of ^nk EF (Figure 2). Thus, the partial contribution of EF k to VF focused on local node 1 (global l) was:

(1.10)

Where subscript a and b represented contour segments , respectively.

Considering the value of σ_i at the midpoint for each segment was equal to as dominant in that segment and incorporating the definitions given in (1.3), then:

(1.11)

where On the other hand, redefining stresses σ_i as functions of linear strains ^ε_ii and angular stress ϒ_ij as:

(1.12)

where were initial stresses.

Deformation ^ε_ii and ϒ_ij described in terms of displacement were:

(1.13)

Where as, for the previous integration, a linear variation of ^u within the EF, of the form:

(1.14)

where ^A,B,C were constants defined in terms of nodal displacements ^{u1,u2 y u3}, equal to:

Displacement variation would thus be:

(1.15)

Introducing the definitions given by equations (1.12) to (1.15), contribution may be expressed in terms of nodal values of displacement as:

(1.16)

A similar expression may be obtained for .

Similar expressions for segment contributions equal to , given by (1.16), for segments and equal to and can be obtained simply by exchanging normal vector ⁿ_a by ⁿ_b or ⁿ_c as appropriate.

Taking the definitions given in (1.10) and incorporating developments (1.16):

(1.17)

(1.18)

The following system of equations was written:

Two algeaic equations were thus obtained for each l value; one for each component of the average displacement in each subdomain ∀^l, all of which constitutes a 2nL x 2nL type linear system of equations.

Results

Diffusion of heat and thermal stress in a steel plate

The transient simulation of plane stress was considered in a twodimensional plate undergoing non-uniform transient heating (Figure 3); properties are given in Table 2. Magnitudes were L = 0.6 (m) length and H = 0.2 (m) height. The following were considered for boundary condition purposes: constrained displacement in direction Y in segment AD (v = 0 at y = 0), constrained in direction X on section BC (u = 0 at x = L) and free CD and LD segments. Naturally, point B was constrained in X and Y ((u,v) = (0,0) in B). The above problem has been dealt with by Boley and Weiner (1960).

Figure 3. Two-dimensional plate with transient non-uniform temperature variation

Transient non-uniform temperature distribution T (ºC) as an exponential function of y (m) coordinate and time t (h) given by:

(1.19)

Product of non-uniform distribution of temperature occurring with stress induced by sustained expansion (creep). Such effect can be modelled with the same model and stress presented for moisture. Initial deformation ε⁰ was based on the following parameters: thermal expansion a, creep coefficient m, stress s, temperature variation DT and for plane strain Poisson ratio n.

Figure 4 shows the mapping of the strain/stress state for t = 1 (h) (steady state). Consistency in imposing boundary conditions can be seen in such mapping (displacement given in Figure 4a). The effects of restriction on free deformation in the development of normal stresses can be seen in Figure 4b

Analysis of convergence and consistency regarding mesh type contrasted with analytical solution is shown in Figure 5, showning that transient convergent results can be achieved with 52x20 uniform mesh.

Figure 6 shows analysis of consistency in terms of time compared to convergent solution in segment BC. Good convergence was observed for 0.001 (h) time steps.

Moisture transport and drying stress in quaking aspen wood

Non-uniform transient simulation of moisture transport in wood inducing stress and strain was considered. Each integration time for moisture distribution was thus determined (as described in Salinas et al., 2004), giving free deformation, and then (ε⁰) encouraging displacement and stress was calculated (the latter was then calculated as being uncoupled, according to Zienkiewicz and Taylor’s initial deformation method, 2000). Consistency analysis revealed that meshes having refinement exponential to surface convection could capture moisture gradients having a smaller number of discrete elements.

Figure 7 shows spatial distribution for drying times regarding stress/strain calculation. Figure 7a shows moisture content transient and spatial distribution having the greatest gradients (concentration of iso-lines) around 60% M. Orthotropic drying pattern was reflected in the lack of drying front symmetry. Figure 7b shows radial displacement (u), showing the correct implementation of boundary conditions, largest displacements being concentrated at the edges, constrained to one degree of freedom.

Figure 8a shows in detail the residual normal stress of radial drying out. It can been observed that the main efforts were concentrated on A, B and C (see Figure 1), being in tension in A and C (surface) and compression in B (centre). The shear stresses shown in Figure 8b were concentrated on the diagonal, slightly shifted towards point D (surface).

Figure 9 highlights in detail the transient evolution of normal radial stress on section BC where its complex evolution according to simulated drying stress stages can been appreciated. Figure 9a shows the first stage of drying characterised by increasing development of tension on the wood surface by contrast with compression in the centre of it. As a complement, Figure 9b shows the second stage characterised by increasing reversal stress

Figure 10 shows the transient changes in kind of stress in the centre and at the surface of the domain of calculation (point B and C, Figure 1); it shows the dynamics of change in normal radial stress intensity correlated with moisture content. There was marked tension and compression at the beginning of drying, on the surface and in the centre, respectively. This contrasted with final stress at the end of drying in which normal stresses were reversed. The latter is an important qualitative feature of the physical problem being simulated. It highlighted the singular strong transient stress developed during drying at t = 195 (h) with CH = 10%.

Conclusions

The results presented in Figures 4-7 led to concluding that the algorithm developed here generated effective simulations of stress/strain caused by free deformation and sustained stress. Consistency analysis of mesh size and time steps inferred that simulation can be achieved with 20x52 mesh and 0001 time steps (h).

The results presented in Figures 8-10 led to concluding that the qualitative features present in wood drying had been simulated (i.e. regions subjected to tension and compression during transitional drying, resulting in residual drying stress characterised by tension on the surface and compression in the centre of the wood). Moreover, accepting the validation of numerical analysis on the thermal plate, the quantitative values of simulated wood drying should be very similar to physical properties. Naturally, it would be desirable to corroborate experimental data, which has been left for future work.

---

References

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, ail 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.

---