Introducción

Los fenómenos físicos pueden ser descritos por medio de ecuaciones diferenciales parciales (^{Johnson, 2009}; ^{Boyce, 2000}) para las cuales es imposible encontrar soluciones analíticas en la mayoría de los casos. Cuando esto sucede, es necesario apoyarse en métodos numéricos como el método de elementos finitos (FEM, por sus siglas en inglés) (^{Yuste, 2006}; ^{Moaveni, 1999}). Este método genera sistemas de ecuaciones enormes que se pueden resolver haciendo uso del poder de cálculo computacional. El presente trabajo de investigación se centra en la creación de una primera versión de una herramienta FEM desarrollada en la plataforma RAD Studio XE (cuyo lenguaje base es C++) (^{Deitel & Deitel, 1995}; ^{Joyanes, 2008}). El software desarrollado es llamado SEDEF © (Soluciones de Ecuaciones Diferenciales por Elemento Finito) y se pretende que su primera versión sea capaz de resolver cualquier sistema diferencial de una dimensión (1D) temporal con derivadas espaciales de orden dos. El programa está dotado de una interfaz gráfica modular que permite plantear y resolver problemas de fenómenos de transferencia de calor, modelos de vigas y sistemas generales de ecuaciones diferenciales. El fin que persigue SEDEF es principalmente subsanar las áreas académicas que requieren del uso de la simulación computacional y que poco progresan en México por los altos costos de las herramientas comerciales existentes en el mercado. Los detalles para su descarga y las políticas de uso se pueden consultar en el Sitio Web del proyecto en www.sedef.mx.

En este artículo, se presenta primeramente la interfaz gráfica de SEDEF. Después, se muestran las ecuaciones genéricas que resuelve el software. Luego, se aplica el método de elementos finitos a la resolución de dichas ecuaciones. Posteriormente, se muestran los tipos de estudio y aplicaciones a fenómenos físicos que se pueden hacer con esta primera versión de SEDEF. Finalmente, se presenta la validación de SEDEF mediante la comparación de sus resultados con las del software comercial de elementos finitos COMSOL Multiphysics y con soluciones analíticas.

Interfaz gráfica

El diseño de la interfaz gráfica se basa principalmente en las etapas de configuración que requiere el método del elemento finito. Estas etapas consisten en la definición de parámetros de estudio, definición de la geometría o dominio, la captura de propiedades de materiales y constantes, la aplicación de condiciones de frontera y condiciones iniciales, la construcción del mallado, el cálculo y finalmente el post-procesamiento. Para cada una de estas etapas, el software cuenta con una serie de herramientas que permiten al usuario realizar de manera rápida y precisa la configuración de su proyecto de simulación.

El software SEDEF incluye una interfaz específica para la configuración y simulación de problemas de transferencia de calor (Heat Transfer), de vigas (Beam) y de sistemas de ecuaciones diferenciales parciales generales o PDE. Para cualquiera que sea el tipo de simulación, se puede elegir entre un estudio de tipo estático o de tipo transitorio.

En cuanto a la construcción del dominio, el software está dotado de herramientas tipo CAD (diseño asistido por computadora), con lo cual un usuario puede crear y manipular dominios con facilidad y de la misma manera puede iniciar el condicionamiento de los mismos (Figura 1). Además, SEDEF proporciona un analizador léxico y sintáctico que permite realizar el condicionamiento del problema a través de expresiones matemáticas (funciones de variables de espacio y tiempo) y no solo de valores numéricos constantes.

Figura 1 Ejemplo de CAD para vigas 

Para la etapa del mallado del dominio, las opciones pueden variar entre la imposición de un número de elementos en el dominio, la exigencia de tamaños máximos y mínimos de elementos, hasta poder introducir una tasa de crecimiento del mallado respecto a una extremidad mediante series geométricas (Figura 2).

Figura 2 Mallado de una geometría 

Finalmente, la etapa de post-procesamiento está programada de acuerdo con el tipo de estudio, por lo cual la interfaz puede variar entre simples gráficas XY, hasta un mapa de colores sobre los dominios (Figura 3).

Figura 3 Módulo de post-procesamiento 

Ecuaciones generales

Para poder modelar gran cantidad de fenómenos físicos es necesario proponer expresiones que puedan representar cualquier sistema de ecuaciones diferenciales parciales (EDPs) (^{Johnson, 2009}; ^{Boyce, 2000}), con coeficientes tanto dependientes del espacio como del tiempo. En SEDEF, se proponen dos ecuaciones genéricas, ambas consideran un vector incógnita ũ de n funciones escalares de las variables de espacio x y tiempo t (x ϵ Ω_s y t ϵ [t_0'τ], Ω_s es el dominio número s). La primera EDP genérica es:

Para algunos fenómenos físicos que involucran términos de convección y difusión, es más cómodo manejar la ecuación genérica siguiente:

Para las ecuaciones 1 y 2, K~~, L~~, M~~, O~~, N~~ y P~~ son matrices de coeficientes del sistema de dimensión n × n, f~ y r~ son vectores conocidos de dimensión n, ũ es el vector de incógnitas de dimensión n y el punto sobre el vector ũ indica derivada parcial de este respecto al tiempo. Las componentes del vector ũ son:

A continuación, se muestra una interpretación física a los términos de la ecuación (2):

• K~~ = coeficientes de masa

• L~~ = coeficientes de amortiguamiento

• -M~~∂u~∂x-O~~u~+r~ = flujo conservativo

• M~~∂u~∂x,O~~u~u~~'r~ = difusión, la convección y la fuente del flujo, respectivamente

• N~~∂u~∂x = la convección

• P~~u~ = la absorción

• f~ = vector fuente

Aplicación del método de elementos finitos

Las ecuaciones diferenciales parciales propuestas en la sección anterior son válidas en el dominio número s: Ω_s=[a_s, b_s]. Para esta porción de línea y su respectiva ecuación, se aplica el método de elementos finitos (FEM) en una forma débil (^{Felippa, 2003}). Por ejemplo, para la ecuación (2), se obtiene la siguiente ecuación integral escalar:

donde φ~ es un vector de funciones de forma de la variable de espacio y el superíndice t indica la transposición del vector. Se define el vector de flujo:

Para el término de la ecuación integral que involucra al vector de flujo se realiza una integración por partes, la cual hace aparecer las condiciones de frontera sobre el vector flujo. Sin embargo, es posible que en lugar del flujo se conozca directamente el valor de la incógnita en la frontera; por lo que los valores conocidos en las fronteras x = a _s y x = b _s son:

Donde:

vfi=la componente i del vector v~

ui=la componente i del vector uf~

Continuando con el proceso, se divide el dominio en m elementos delimitados por nodos. Se interpola a ũ en cada elemento con polinomios de x de grado p ∈ {1, 2, 3, 4} y se consideran funciones de forma del mismo grado polinomial. Al verificar la ecuación integral para cualquier función de forma se obtiene el sistema de ecuaciones diferenciales ordinarias:

Donde: U~= vector que contiene el valor de cada una de las incógnitas en cada nodo de la malla. σ~~,λ~~, y T~~ = matrices cuadradas de las masas, amortiguamiento y rigideces respectivamente e~= vector fuente.

Es posible que el problema cuente con más de un dominio y además los dominios tengan zonas de unión, para esto se definen fronteras internas y externas como se muestra en la Figura 4. Para asegurar la continuidad en las fronteras internas, se escriben las expresiones 8 y 9.

Figura 4 Ejemplo de un problema con tres dominios conectados 

Donde:

g~ = 0~ si el flujo es continuo, g~≠0≠ si se desea definir una discontinuidad del flujo, u~cses el vector de incógnitas en el nodo de unión en el dominio s, v~fcs es el vector de flujo en el nodo de unión del dominio s, n _Ω es el número de dominios.

Tipos de estudio

Adaptaciones a fenómenos físicos

Vigas en el plano

Para vigas de Euler en el plano (^{Liu, 2003}; ^{Slivker, 2007}), se tienen dos desplazamientos u y v y una rotación θ dentro del plano. Las 3 ecuaciones de equilibrio son:

ρA∂2u∂t2-∂N∂x=fxρA∂2v∂t2-∂Q∂x=fyIzz∂2θ∂t2-∂Mz∂x-Q=m (18)

mientras que el comportamiento mecánico está definido por las 3 ecuaciones siguientes:

∂v∂x-θ=0Mz-EIzz∂θ∂x=Mz0-EIzz∂θ0∂xN-EA∂v∂x=N0-EAε0n (19)

En las ecuaciones (18-19):

• ρ, E, A, I _zz , N, Q y M _z son la densidad del material, el módulo de Young del mismo, el área de la sección, el momento de inercia de la sección, la fuerza normal, la fuerza cortante y el momento, respectivamente.

• F _x´, f _y y m son de las densidades lineales de fuerzas aplicadas en x y y, y la densidad lineal de momento aplicado.

• N _0´, Mz0 _, θ0 y ε0nson valores previos de fuerza normal, momento, ángulo de rotación y deformación normal (pre-esfuerzos, pre-deformaciones).

El vector de incógnitas escogido para la ecuación (2) es:

u~=uvθQMzN (20)

Las matrices de coeficientes y el vector fuente tienen las siguientes expresiones:

K~~=ρA000000ρA000000Izz000000000000000000000,O~~=000001000100000010000000000000000000 (21)

N~~=00000000000000000000000000-EIzz000-EA00000,P~~=000000000000000-10000-1000000010000001 (22)

f~=fxfym0Mz0-EIzz∂θ0∂xN0-EAε0n (23)

Las demás matrices en la ecuación (2) son cero. Estas ecuaciones son válidas en el sistema de coordenadas locales. Para definir un sistema en común para todos los dominios, llamado sistema global (Figura 5), se implementa la matriz de rotación:

R~~=cos⁡θesin⁡θe0000-sin⁡θecos⁡θe0000001000000100000010000001 (24)

Al definir la matriz de rotación R~~, el vector de funciones incógnitas u~g en coordenadas globales (X, Y) es:

u~g=UVθQMzN=R~~uvθQMzN (25)

y el vector de flujo:

v~f=Qsin⁡θe-Ncos⁡θe-Nsin⁡θe-Qcos⁡θe-Mz000 (26)

Por lo tanto, según la ecuación 9, la ecuación de equilibrio de los nodos de conexión (los cuales son fronteras internas) entre dominios distintos se inscribe en:

∑F~i+F~e=0~ (27)

Figura 5 Elemento viga en el plano 

donde F~i es el vector de fuerzas internas aplicadas en el nodo de conexión F~e y es el vector de fuerzas externas aplicadas en dicho nodo.

Validación

Validación 1: Estudio estacionario con dos variables y dos dominios

El problema consiste de los dominios Ω₁ = [0,5] y Ω₂ = [5,7.3] con dos funciones incógnita u ₁ y u ₂. El sistema de ecuaciones, en forma de la ecuación genérica (2), es:

k=1 y 2, ∀xϵΩk:ddx-M~~kdu1dxdu2dx+N~~kdu1dxdu2dx+P~~ku1u2=f~k (28)

Donde:

M~~1=-1-2-x-2,N~~1=00-10,P~~1=1101,f~~1=0x(90x-15x2-90) (29)

M~~2=-x000,N~~2=-1111,P~~2=0000,f~~2=100100-45x2

Para este problema, no se plantea la continuidad del flujo en el nodo x = 5. Las condiciones de frontera son:

• u ₁ (0) = 0, u ₁ (5) = 15 (5)³, u ₁(7.3) = 15 (7.3)³

• u ₂ (0) = u ₂ (5) = 1

La solución analítica del problema es:

∀xϵΩ1,u~=15x390x-15x3+cosx2+sin⁡x21426csc52-cot52 (30)

y

∀xϵΩ2,u~=15x3100x-30x3+3251 (31)

Para el estudio se utilizaron 100 elementos lineales por cada dominio. En la Figura 6 se muestra la comparación entre la solución analítica y la solución numérica de SEDEF. Se aprecia que se obtienen prácticamente los mismos resultados. El error relativo máximo cometido por SEDEF es de 1.38 %. Se puede disminuir este error utilizando un mayor número de elementos o empleando elementos de mayor orden (por ejemplo: cúbicos en lugar de lineales, es decir, funciones de forma de tipo polinomios de grados 3 en lugar de grado 1).

Figura 6 Solución analítica vs solución de SEDEF para un sistema de ecuaciones diferenciales 

Validación 2: Estudio de transferencia de calor dependiente del tiempo

Para esta prueba se desea predecir la distribución de temperatura ∀tϵ[0,1] (en segundos), de una barra de longitud con las siguientes propiedades: κ=1W/(m∙k) = 1 W/(m·k), ρ = 1 kg/m³ y c_p = 1 J/(kg·K)

Las condiciones iniciales están dadas por una temperatura uniforme e igual a 1 °C. Las condiciones de frontera están dadas por la condición que la temperatura en ambos bordes (x = 1m y x = -1m) decrece linealmente con el tiempo:

Tx=1,t=Tx=-1,t=1-γt (32)

donde γ = 1 °C/s. Para validar los resultados obtenidos, se compararán estos contra los resultados del software COMSOL. En el mallado se emplearon 27 elementos de segundo grado.

Las temperaturas obtenidas por ambos software son graficadas en función de la variable de espacio en la Figura 7, y esto en diferentes instantes. El resultado obtenido por ambos software es prácticamente el mismo (el error relativo máximo es 0.113 %); esto significa que el módulo de transferencia de calor de SEDEF fue adaptado correctamente a partir de la ecuación (2) y que los resultados del software son precisos.

Figura 7.  Temperaturas obtenidas por COMSOL y SEDEF 

Validación 3: Conjunto de vigas en el plano

Se desea predecir los desplazamientos en el puente que se muestra en la Figura 8, hecho de acero estructural. En esta figura, los nodos están enumerados del 0 al 15 en color negro mientras que los dominios van del 0 al 22 en color rojo. Los nodos 0 y 4 están empotrados; en los nodos 10 y 15 existen apoyos de articulación y de rodillo, respectivamente.

Figura 8.  Geometría considerada del puente; números negros hacen referencia al nodo, números rojos al dominio 

El puente está sometido a las cargas mostradas en la Tabla 1 y al peso propio de la estructura que está dado por:

W1=gAiρ (33)

Donde:

W_i= peso propio del dominio

i,g = 9.81 m/s²

A_i = área transversal del dominio i y ρ = 7850 kg/m3 = densidad

Nodo	Fx(N)	Fxy(N)
11	3000	3000
12	3000	3000
13	-3000	3000
14	-3000	3000

Se tienen dos tipos secciones transversales diferentes: tipo 1 y tipo 2. La nomenclatura de sus dimensiones se muestra en la Figura 9 y sus valores en la Tabla 2 para cada sección.

Figura 9 Secciones transversales consideradas 

Tabla 2 Dimensiones de la sección transversal de cada dominio 

Dominio(s)	Tipo de sección	h(m)	b(m)	t w (m)	t p (m)
0.1	1		0.4	0.025	0.025
2.3	1		0.4	0.025	0.025
4-12	2	0.25	0.25	0.025	0.025
13-22	1	0.4	0.2	0.0125	0.0125

Los dominios 13 a 22, además de estar sometidos a su peso propio, soportan una densidad lineal de carga en dirección y de grados de -24.7212 kN/m y sufren una variación de temperatura δT = 30 °C. El coeficiente de dilatación térmica del material es α = 12.03 10^-6 C^-1 y el módulo de Young es E = 200 GPa.

La altura de la sección transversal del arco varía linealmente conforme a las siguientes fórmulas:

∀xϵ-8,0:h1X=0.6-0.3L(X)L(0)∀xϵ0,8:h2X=0.31+LX-L(0)L8-L(0) (34)

Donde:

LX=1812X64X2+X2+32lnx8+X264+1+42+ln⁡(1+2) (35)

Se utilizaron 40 elementos por dominio y elementos con funciones de forma de grado 3. En la Figura 10, se muestran las deformadas obtenidas por COMSOL y por SEDEF utilizando un factor de amplificación de 450. No se aprecia diferencia alguna entre las dos deformadas. Para cada nodo, se determinó la magnitud de desplazamiento con SEDEF y se calculó la diferencia relativa respecto al resultado obtenido por COMSOL. En la Figura 11, se muestra un histograma con las diferencias relativas mencionadas para cada nodo. Estas diferencias son muy pequeñas: la máxima diferencia no alcanza 0.5 %. Esta diferencia tiende hacia cero a medida que se refina el mallado. Cabe señalar que la diferencia entre los resultados de SEDEF y COMSOL radica en que COMSOL no maneja el mismo orden de elementos (mismo grado polinomial de funciones de forma) para los desplazamientos y rotaciones mientras que en SEDEF sí, y esto no solo para desplazamientos y rotaciones, sino también para momentos de flexión y fuerzas normales y cortantes. Por todas estas razones, se valida el módulo de vigas de SEDEF y se comprueba su utilidad y precisión para problemas de vigas sometidas a cargas termo-mecánicas.

Figura 10.  Comparación de la geometría deformada entre, a) COMSOL y b) SEDEF 

Figura 11.  Diferencia relativa de desplazamientos nodales de SEDEF respecto a los de COMSOL  

Conclusiones

Se desarrolló la primera versión de un software gratuito llamado SEDEF para la resolución de ecuaciones diferenciales parciales por el método de elementos finitos. En esta versión, SEDEF es capaz de resolver problemas transientes de una dimensión espacial en sus tres módulos: ecuaciones matemáticas genéricas, transferencia de calor en sólidos y vigas en el plano. Se realizaron varias pruebas confrontando los resultados de SEDEF contra soluciones analíticas y numéricas obtenidas por el software comercial COMSOL Multiphysics y se comprobó la precisión de SEDEF para su validación. También, SEDEF resuelve los problemas de manera rápida y utilizando pocos recursos de memoria de la computadora. Con lo anterior, se puede decir que el programa es eficiente y eficaz. La eficiencia conseguida se debe principalmente a la implementación de la librería CS-parse.

En un futuro, el programa se aplicará a más fenómenos físicos, entre ellos se encuentra el de vigas en el espacio. También se desarrollará un algoritmo capaz de obtener los valores y vectores propios de las matrices dispersas. Asimismo, se hará un desarrollo para resolver problemas 2D. Todo esto será implementado en la segunda versión del software gratuito SEDEF.

A pesar de las limitaciones de la primera versión de SEDEF, el software puede ser ya una herramienta de gran valor pedagógico para clases de matemáticas, métodos numéricos, física y mecánica, entre otras. También, el software puede ser de utilidad para científicos e ingenieros, particularmente para ingenieros civiles (Yuste, 2006) por el alcance del módulo de vigas.