Cálculo de las Constantes de Torsión para Vigas de Secciones Transversales Arbitrarias Utilizando Elementos Finitos

Cálculo de las Constantes de Torsión para Vigas de Secciones Transversales Arbitrarias Utilizando Elementos Finitos

Jesús Gerardo Valdés–Vázquez. Departamento de Ingeniería Civil y Ambiental, Universidad de Guanajuato. Av. Juárez 77, Zona Centro. Guanajuato, Gto., México. valdes@ugto.mx

Jesús Fernando Valdés–Vázquez. Departamento de Ingeniería Civil y Ambiental, Universidad de Guanajuato. Av. Juárez 77, Zona Centro. Guanajuato, Gto., México. valdesjfv@ugto.mx

Resumen

En el análisis estructural, es esencial que un ingeniero cuente con un profundo conocimiento de las propiedades de una sección transversal, tales como el área, los momentos de inercia y las constantes de torsión. En particular, el cálculo de las constantes de torsión presenta un desafío significativo, ya que, para secciones conocidas, existen fórmulas exactas o aproximadas que facilitan el proceso. Sin embargo, con las actuales tendencias en construcción que incorporan secciones compuestas o de geometría arbitraria, el cálculo de estas constantes se convierte en un reto considerable. Este trabajo se centra, principalmente, en explorar la obtención de las constantes de torsión mediante métodos numéricos, mejorando su precisión en comparación con las propuestas de otros investigadores.

Palabras clave

Constante de torsión. Torsión de Saint-Venant. Torsión por alabeo.

Introducción

El análisis estructural tridimensional de edificaciones, compuesta por vigas y columnas, demanda un profundo entendimiento de las propiedades geométricas inherentes a las secciones transversales de estos elementos. Entre estas propiedades se encuentran el área, los momentos de inercia y las constantes de torsión. Desde la magnificencia de antiguos templos hasta la imponente altura de los contemporáneos rascacielos, el cálculo de la torsión ha representado un desafío constante para los ingenieros a lo largo de la historia.

En los últimos siglos, el estudio de la torsión ha evolucionado, siendo iniciado por distinguidos pensadores. Johann Bernoulli (1667-1748) destacó como precursor, aportando a la comprensión de las fuerzas internas en estructuras y estableciendo las primeras bases para la teoría de vigas. Posteriormente, Daniel Bernoulli (1700-1782), junto con su hermano Johann, realizó contribuciones fundamentales a la teoría de la elasticidad y la mecánica, cimentando así los pilares para el estudio detallado de la torsión. Casi en paralelo, Leonhard Euler (1707-1783) desplegó una contribución esencial al establecer los cimientos matemáticos para la comprensión de fenómenos torsionales. En el siglo XIX, Adhémar Jean Claude Barre de Saint-Venant (1797-1886) formuló las ecuaciones fundamentales para la torsión, enriqueciendo la teoría de la elasticidad y la mecánica de sólidos deformables. En épocas más recientes, Stephen Timoshenko (1878-1972) emergió como una figura clave, realizando contribuciones significativas, especialmente en la teoría de vigas y torsión. En el ámbito práctico de la ingeniería estructural, (Euler, 1744), (Saint-Venant, 1855) y (Timoshenko, 1934) han dejado un legado influyente, proporcionando herramientas y teorías esenciales que han sido fundamentales para el diseño y análisis de estructuras en la ingeniería moderna.

En el análisis espacial de marcos conformados por vigas y columnas, se emplean predominantemente secciones circulares, rectangulares, tubos circulares, tubos cuadrados y secciones en forma de I. Todas estas tipologías cuentan con fórmulas exactas o aproximadas para el cálculo de la constante de torsión de Saint-Venant ( ${\textstyle J}$), cuyos detalles se pueden encontrar disponibles en (Boresi, Schmidt, & Sidebottom, 1993). Sin embargo, en lo que respecta a la constante de torsión por alabeo o warping ( ${\textstyle {C}_{w}}$), calcularla no resulta tan sencillo, llegando al punto en el cual su influencia en el análisis estructural contemporáneo tiende a ser pasada por alto en la mayoría de los softwares comerciales.

Excluyendo la torsión por alabeo, el análisis de estructuras con secciones arbitrarias o compuestas, como en el caso del análisis de tableros de puentes, plantea un desafío adicional. La determinación de la constante de torsión de Saint-Venant se vuelve todo un reto y, en ocasiones, incluso resulta imposible debido a la complejidad geométrica involucrada. En situaciones de esta índole, se recurre a métodos numéricos para abordar el cálculo de la torsión. Estos métodos pueden incluir la aplicación de técnicas como el método de elementos finitos, tal como se detalla en (Pilkey, 2002).

Propuestas más contemporáneas que emplean el método de los elementos finitos para la determinación de la constante de torsión de Saint-Venant ( ${\textstyle J}$) y la constante de torsión por alabeo ( ${\textstyle {C}_{w}}$) se describen detalladamente en la investigación llevada a cabo por (Haukaas, 2020). Su estudio se centra específicamente en el cálculo mediante elementos finitos de ambas constantes de torsión utilizando triángulos lineales. Nuestro artículo se basa en este último y presenta una mejora en la determinación de las inercias, las cuales son necesarias como paso previo al cálculo de dichas constantes.

Conceptos teóricos y formulación

Comencemos considerando una viga con sección transversal arbitraria que está empotrada en un extremo y sometida a una fuerza de torsión ${\textstyle T}$ en el otro extremo. Según el principio de Saint-Venant, la distribución de esfuerzos en una sección transversal a una distancia ${\textstyle x}$ considerable de los extremos depende principalmente de la fuerza ${\textstyle T}$ y no de la distribución de esfuerzos en dichos extremos. Saint-Venant observó que la deformación por torsión resulta de los desplazamientos de la sección. Siguiendo la formulación de (Boresi, Schmidt, & Sidebottom, 1993), si la fuerza de torsión ${\textstyle T}$ induce una rotación ${\textstyle \theta }$ a una distancia ${\textstyle x}$ de uno de los extremos, las componentes del vector de desplazamiento pueden expresarse mediante:

${\textstyle u=}$$\theta \cdot \Omega (y,z)$, ${\textstyle v=-\theta xz}$, ${\textstyle w=}$$\theta xy$

(1)

Aquí, ${\textstyle \Omega (y,z)}$ representa la función de alabeo. La determinación de esta función se lleva a cabo de manera que satisfaga las ecuaciones de elasticidad. Dado que estamos considerando un tensor de deformaciones infinitesimales, las ecuaciones de compatibilidad se cumplen de manera automática. Las deformaciones en un punto ${\textstyle x}$ de la viga se expresarán entonces como:

${\epsilon }_{x}={\epsilon }_{y}={\epsilon }_{z}={\gamma }_{yz}=0$
${\gamma }_{xy}=\theta \left({\frac {\partial \Omega }{\partial y}}-z\right)$	(2)
${\gamma }_{xz}=\theta \left({\frac {\partial \Omega }{\partial z}}+y\right)$

Al emplear la Ley de Hooke Generalizada, la cual establece la relación entre esfuerzos y deformaciones, derivamos de la ecuación previa las deformaciones distintas de cero de la siguiente manera:

${\tau }_{xy}=G\cdot {\gamma }_{xy}$	(3)
${\tau }_{xz}=G\cdot {\gamma }_{xz}$

Si sustituimos la ecuación (3) en las ecuaciones de equilibrio interno encontradas en la mecánica de medios continuos y despreciamos las fuerzas másicas encontramos:

${\frac {\partial {\tau }_{xy}}{\partial x}}=0$
${\frac {\partial {\tau }_{xz}}{\partial x}}=0$	(4)
${\frac {\partial {\tau }_{yx}}{\partial y}}+{\frac {\partial {\tau }_{zx}}{\partial z}}=$$0$

En virtud de la conservación del momento angular, que conduce a la simetría del tensor de esfuerzos, se demuestra que los esfuerzos cortantes presentes en la ecuación (4) son independientes de la posición ${\textstyle x}$. Como consecuencia, es necesario que exista una función de esfuerzos ${\textstyle \phi (y,z)}$, reconocida como la función de esfuerzos de Prandtl, que satisface la relación:

${\tau }_{xy}={\frac {\partial \phi }{\partial z}}$	(5)
${\tau }_{xz}=-{\frac {\partial \phi }{\partial y}}$

Mediante condiciones de contorno adecuadas y siguiendo los principios de la teoría descrita por (Boresi, Schmidt, & Sidebottom, 1993), es posible establecer la relación entre la fuerza de torsión ${\textstyle T}$, los esfuerzos cortantes y la función de Prandtl ${\textstyle \phi }$ , lo que se traduce en:

$T=\int _{}^{}\int _{}^{}y\cdot {\tau }_{xz}-z\cdot {\tau }_{xy}\,dy\,dz=2\,\int _{}^{}\int _{}^{}\phi dy\,dz$

(6)

En términos sencillos, esta ecuación establece que el torque es equivalente a dos veces el volumen bajo la superficie formada por la función de esfuerzos y el plano de la sección transversal ${\textstyle yz}$. Es crucial destacar que esta expresión fue derivada considerando deformaciones infinitesimales y aplicando relaciones de elasticidad lineal.

Aunque la solución específica para la función de Prandtl nos brinda la capacidad de determinar la constante de torsión en diversas secciones transversales, nuestro propósito principal radica en aplicar la teoría previamente expuesta para alcanzar soluciones numéricas en contextos donde las secciones transversales son de naturaleza arbitraria. En lugar de confinarnos a soluciones analíticas particulares, buscamos la versatilidad y precisión que ofrecen los métodos numéricos al abordar la complejidad inherente a secciones transversales no convencionales. Este enfoque nos permite extender la aplicabilidad de la teoría a una gama más amplia de situaciones prácticas en ingeniería estructural.

Con la finalidad de encontrar las constantes de torsión utilizando el método de los elementos finitos, comencemos con escribir las ecuaciones anteriores en forma matricial. También consideremos que la función de alabeo es debida a una rotación unitaria ${\textstyle \theta =}$$1$, de manera que el desplazamiento longitudinal debido a esta rotación es:

$u=1\cdot \Omega (y,z)$

(7)

La expresión matricial de la ecuación (2) es, por tanto:

${\mathit {\boldsymbol {\epsilon }}}=\left[{\begin{matrix}{\gamma }_{xy}\\{\gamma }_{xz}\end{matrix}}\right]=$$\left[{\begin{matrix}{\frac {\partial \Omega }{\partial y}}\\{\frac {\partial \Omega }{\partial z}}\end{matrix}}\right]+$$\left[{\begin{matrix}-z\\y\end{matrix}}\right]={\mathit {\boldsymbol {\epsilon }}}_{w}+$${\mathit {\boldsymbol {\epsilon }}}_{0}$

(8)

En esta formulación, el vector de deformaciones ${\textstyle {\mathit {\boldsymbol {\epsilon }}}}$ se divide en dos componentes distintas. La primera, ${\textstyle {\mathit {\boldsymbol {\epsilon }}}_{w}}$​, corresponde a las deformaciones generadas por el alabeo, mientras que la segunda, ${\textstyle {\mathit {\boldsymbol {\epsilon }}}_{0}}$​, representa las deformaciones inducidas por este fenómeno, pudiendo ser interpretada como una "carga externa". Por otro lado, la expresión (3) adquiere una forma matricial cuando se expresa de la siguiente manera:

${\mathit {\boldsymbol {\sigma }}}=\,\left[{\begin{matrix}{\tau }_{xy}\\{\tau }_{xz}\end{matrix}}\right]=$$\left[{\begin{matrix}G&0\\0&G\end{matrix}}\right]\cdot \left[{\begin{matrix}{\gamma }_{xy}\\{\gamma }_{xz}\end{matrix}}\right]=$${\boldsymbol {C\cdot }}{\mathit {\boldsymbol {\epsilon }}}$

(9)

en la que ${\textstyle {\boldsymbol {C}}}$ representa la ecuación constitutiva. La ecuación anterior se puede representar por

${\mathit {\boldsymbol {\sigma }}}={\boldsymbol {C\cdot }}{\mathit {\boldsymbol {\epsilon =\,}}}{\boldsymbol {C\cdot }}{\mathit {\boldsymbol {\epsilon }}}_{w}+$$\,{\boldsymbol {C\cdot }}{\mathit {\boldsymbol {\epsilon }}}_{0}$

(10)

En la que hemos sustituido la ecuación (8) en la ecuación (9). En esta última ecuación se aprecia de manera más evidente que la segunda expresión del lado derecho de la igualdad es equivalente a una “carga externa”.

Discretización con elementos finitos

La forma débil de las ecuaciones anteriores se puede expresar mediante el principio del trabajo virtual, en la que:

$\delta {W}_{int}=\delta {W}_{ext}$

(11)

Podemos escribir que el trabajo virtual interno es igual a:

$\delta {W}_{int}=\int _{A}^{}\delta {\mathit {\boldsymbol {\epsilon }}}\cdot {\mathit {\boldsymbol {\sigma }}}\,dA$

(12)

donde la variación de la deformación es igual a:

$\delta {\mathit {\boldsymbol {\epsilon }}}=\delta \,\left({\mathit {\boldsymbol {\epsilon }}}_{w}+\right.$$\left.{\mathit {\boldsymbol {\epsilon }}}_{0}\right)=\delta {\mathit {\boldsymbol {\epsilon }}}_{w}$

(13)

ya que ${\textstyle {\mathit {\boldsymbol {\epsilon }}}_{0}}$ es como una constante. Entonces, el trabajo virtual interno puede expresarse de la siguiente manera:

$\delta {W}_{int}=\int _{A}^{}\delta {\mathit {\boldsymbol {\epsilon }}}_{w}\cdot {\mathit {\boldsymbol {\sigma }}}dA=$$\int _{A}^{}\delta {\mathit {\boldsymbol {\epsilon }}}_{w}\cdot \left({\boldsymbol {C\cdot }}{\mathit {\boldsymbol {\epsilon }}}\right)dA=$$\int _{A}^{}\delta {\mathit {\boldsymbol {\epsilon }}}_{w}\cdot \left({\boldsymbol {C\cdot }}{\mathit {\boldsymbol {\epsilon }}}_{w}\right)dA+$$\,\int _{A}^{}\delta {\mathit {\boldsymbol {\epsilon }}}_{w}\cdot \left({\boldsymbol {C\cdot }}{\mathit {\boldsymbol {\epsilon }}}_{0}\right)dA$

(14)

En la ecuación (14) se puede ver que el último término que involucra a ${\textstyle {\mathit {\boldsymbol {\epsilon }}}_{0}}$ es constante, por lo que no representa un trabajo virtual interno y lo podemos mandar al lado derecho de la ecuación (11) convirtiéndose en trabajo virtual externo. De esta forma, el principio del trabajo virtual se expresa de manera siguiente:

$\left[\delta {W}_{int}=\int _{A}^{}\delta {\mathit {\boldsymbol {\epsilon }}}_{w}\cdot {\boldsymbol {C\cdot }}{\mathit {\boldsymbol {\epsilon }}}_{w}dA\right]=$$\left[\delta {W}_{ext}=-\int _{A}^{}\delta {\mathit {\boldsymbol {\epsilon }}}_{w}\cdot {\boldsymbol {C\cdot }}{\mathit {\boldsymbol {\epsilon }}}_{0}dA\right]$

(15)

La ecuación (8) depende de las variables independientes ${\textstyle \Omega }$ , ${\textstyle y}$ y ${\textstyle z}$, las cuales procederemos a discretizar de la siguiente manera:

$y\approx {y}^{h}=\sum _{J=1}^{{n}_{nodos}}{N}_{J}\left(y,z\right)\cdot {y}_{J}$
$z\approx {z}^{h}=\sum _{J=1}^{{n}_{nodos}}{N}_{J}\left(y,z\right)\cdot {z}_{J}$	(16)
$\Omega \approx {\Omega }^{h}=\sum _{J=1}^{{n}_{nodos}}{N}_{J}\left(y,z\right)\cdot {u}_{J}$

siendo ${\textstyle {N}_{J}(y,z)}$ las funciones de forma del nodo ${\textstyle J}$ del elemento finito bidimensional utilizado en la discretización de la sección transversal, ${\textstyle {y}_{J}}$ la coordenada en dirección del eje ${\textstyle y}$ de cada nodo ${\textstyle J}$ del elemento finito bidimensional y ${\textstyle {z}_{J}}$ su correspondiente coordenada en dirección del eje ${\textstyle z}$. Por otro lado, ya que se esta asumiendo una rotación unitaria ${\textstyle \theta =}$$1$, entonces, de la ecuación (7), los desplazamientos son iguales a la función de alabeo, es decir ${\textstyle u=}$$\Omega (y,z)$, siendo ${\textstyle u}$ el desplazamiento axial. Como se describe en la ecuación (16), ${\textstyle {u}_{J}}$ es el desplazamiento axial del nodo ${\textstyle J}$ en el elemento finito. A continuación, expresaremos matricialmente las deformaciones de la ecuación (8) utilizando los valores ya discretizados en la ecuación (16), obteniendo las siguientes expresiones:

${\mathit {\boldsymbol {\epsilon }}}_{0}\approx {\mathit {\boldsymbol {\epsilon }}}_{0}^{h}=$$\left[{\begin{matrix}-{z}^{h}\\{y}^{h}\end{matrix}}\right]=\sum _{J=1}^{{n}_{nodos}}\left[{\begin{matrix}-{N}_{J}\left(y,z\right)\cdot {z}_{J}\\{N}_{J}\left(y,z\right)\cdot {y}_{J}\end{matrix}}\right]=$$\sum _{J=1}^{{n}_{nodos}}\left[{\begin{matrix}0&-{N}_{J}\left(y,z\right)\\{N}_{J}\left(y,z\right)&0\end{matrix}}\right]\cdot \left[{\begin{matrix}{y}_{J}\\{z}_{J}\end{matrix}}\right]$

(17)

mientras que

${\mathit {\boldsymbol {\epsilon }}}_{w}\approx {\mathit {\boldsymbol {\epsilon }}}_{w}^{h}=$$\left[{\begin{matrix}{\frac {\partial {\Omega }^{h}}{\partial y}}\\{\frac {\partial {\Omega }^{h}}{\partial z}}\end{matrix}}\right]=$$\sum _{J=1}^{{n}_{nodos}}\left[{\begin{matrix}{\frac {\partial {N}_{J}\left(y,z\right)}{\partial y}}\cdot {u}_{J}\\{\frac {\partial {N}_{J}\left(y,z\right)}{\partial z}}\cdot {u}_{J}\end{matrix}}\right]=$$\sum _{J=1}^{{n}_{nodos}}\left[{\begin{matrix}{\frac {\partial {N}_{J}\left(y,z\right)}{\partial y}}\\{\frac {\partial {N}_{J}\left(y,z\right)}{\partial z}}\end{matrix}}\right]\cdot \left[{u}_{J}\right]$

(18)

Sin embargo, la ecuación (18) se puede apreciar que se trata de las derivadas de las funciones de forma, cuya notación habitual en elementos finitos es mediante la matriz ${\textstyle {\mathit {\boldsymbol {B}}}}$. Entonces, la ecuación (18) se puede simplificar y escribirla de la forma:

${\mathit {\boldsymbol {\epsilon }}}_{w}\approx {\mathit {\boldsymbol {\epsilon }}}_{w}^{h}=$$\left[{\begin{matrix}{\frac {\partial {\Omega }^{h}}{\partial y}}\\{\frac {\partial {\Omega }^{h}}{\partial z}}\end{matrix}}\right]=$$\sum _{J=1}^{{n}_{nodos}}\left[{B}_{J}\right]\cdot \left[{u}_{J}\right]$

(19)

La discretización de la ecuación (13) da como resultado la parte virtual, en la que se cambia de notación para el nodo ${\textstyle I}$ en lugar del nodo ${\textstyle J}$ y así evitar confusiones, por lo que:

$\delta {\mathit {\boldsymbol {\epsilon }}}_{w}\approx {\delta {\mathit {\boldsymbol {\epsilon }}}}_{w}^{h}=$$\sum _{I=1}^{{n}_{nodos}}\left[{B}_{I}\right]\cdot \left[\delta {u}_{I}\right]$

(20)

Sustituyendo las expresiones anteriores en el trabajo virtual interno, encontramos:

$\delta {W}_{int}={\left[\delta {u}_{I}\right]}^{T}\cdot \int _{A}^{}\sum _{I=1}^{{n}_{nodos}}\sum _{J=1}^{{n}_{nodos}}{\left[{B}_{I}\right]}^{T}\cdot {\boldsymbol {C\cdot }}\left[{B}_{J}\right]\cdot \left[{u}_{J}\right]\,dA$

(21)

mientras que el trabajo virtual externo es igual a:

$\delta {W}_{ext}={\left[\delta {u}_{I}\right]}^{T}\cdot \int _{A}^{}\sum _{I=1}^{{n}_{nodos}}\sum _{J=1}^{{n}_{nodos}}{\left[{B}_{I}\right]}^{T}\cdot {\boldsymbol {C\cdot }}\left[{\begin{matrix}0&-{N}_{J}\left(y,z\right)\\{N}_{J}\left(y,z\right)&0\end{matrix}}\right]\cdot \left[{\begin{matrix}{y}_{J}\\{z}_{J}\end{matrix}}\right]\,dA$	(22)
$\delta {W}_{ext}={\left[\delta {u}_{I}\right]}^{T}\cdot \int _{A}^{}\sum _{I=1}^{{n}_{nodos}}\sum _{J=1}^{{n}_{nodos}}{\left[{B}_{I}\right]}^{T}\cdot {\boldsymbol {C\cdot }}{\mathit {\boldsymbol {\epsilon }}}_{0}^{h}\,dA$

Ya que los desplazamientos virtuales son arbitrarios, encontramos que las fuerzas internas son iguales a las fuerzas externas, cuya solución nos permite encontrar los desplazamientos ${\textstyle {\mathit {\boldsymbol {u}}}}$ que representa la incógnita del problema denominada función de alabeo. La expresión para resolver resulta ser:

$\int _{A}^{}{\mathit {\boldsymbol {B}}}^{T}\cdot {\mathit {\boldsymbol {C}}}\cdot {\mathit {\boldsymbol {B}}}\,dA\cdot {\mathit {\boldsymbol {u=-}}}\int _{A}^{}{\mathit {\boldsymbol {B}}}^{T}\cdot {\mathit {\boldsymbol {C}}}\cdot {\mathit {\boldsymbol {\epsilon }}}_{o}^{h}\,dA$

(23)

Hay que destacar que, en planteamiento anterior, ${\textstyle {\mathit {\boldsymbol {B}}}=}$$\sum _{}^{}\left[{\mathit {\boldsymbol {B}}}_{I}\right]$ . La ecuación (23) resulta muy común para aquellos familiarizados con el análisis matricial de estructuras donde, la ecuación anterior resulta ser:

${\mathit {\boldsymbol {K}}}\cdot {\mathit {\boldsymbol {u}}}={\mathit {\boldsymbol {f}}}^{ext}$

(24)

En la ecuación anterior resulta que la matriz ${\textstyle {\mathit {\boldsymbol {K}}}}$ es singular, ya que el sistema no cuenta con condiciones de contorno. La solución propuesta por (Haukaas, 2020) consiste en igualar a uno algún grado de libertad de los desplazamientos de la ecuación y después resolver para el resto de los desplazamientos. Las discretizaciones anteriores nos permiten encontrar el momento de torsión ${\textstyle T}$ definido en la ecuación (6), que escrita en forma matricial resulta:

$T=\,\int _{A}^{}\left[{\begin{matrix}-z&y\end{matrix}}\right]\cdot \left[{\begin{matrix}{\tau }_{xy}\\{\tau }_{xz}\end{matrix}}\right]\,dA=\,\int _{A}^{}\left[{\begin{matrix}-z&y\end{matrix}}\right]\cdot \left[{\begin{matrix}G&0\\0&G\end{matrix}}\right]\cdot \left[{\begin{matrix}{\gamma }_{xy}\\{\gamma }_{xz}\end{matrix}}\right]\,dA$	(25)
$T=\,\int _{A}^{}\left[{\begin{matrix}-{z}^{h}&{y}^{h}\end{matrix}}\right]\cdot \left[{\begin{matrix}G&0\\0&G\end{matrix}}\right]\cdot \left\{{\mathit {\boldsymbol {B}}}\cdot {\mathit {\boldsymbol {u}}}+\left[{\begin{matrix}-{z}^{h}\\{y}^{h}\end{matrix}}\right]\right\}\,dA$

Aspectos computacionales

Una vez encontrado el vector de los desplazamientos ${\textstyle {\mathit {\boldsymbol {u}}}}$, se tiene que normalizar el sistema mediante la ecuación:

$\int _{A}^{}{\Omega }^{h}\left(y,z\right)\,dA=0$

(26)

De esta manera, se calcula la constante de corrección por alabeo para corregir el campo de los desplazamientos. En nuestro caso por haber considerado una rotación unitaria, los desplazamientos coinciden con el resultado de nuestra función de alabeo ${\textstyle \Omega \left(y,z\right)}$ . Otra corrección al campo de los desplazamientos viene dada por el centro de cortante. De acuerdo con la torsión de St.-Venant y lo referido considerando la formulación por elementos finitos de (Haukaas, 2020) las coordenadas del centro de cortante se pueden calcular mediante:

${y}_{sc}={z}_{sc}\cdot {\frac {{I}_{yz}}{{I}_{y}}}-\,{\frac {{\left({z}_{K}\right)}^{T}\int _{A}^{}{N}_{I}^{T}\cdot {N}_{J}\,dA\,{u}_{K}}{{I}_{y}}}$	(27)
${z}_{sc}={y}_{sc}\cdot {\frac {{I}_{yz}}{{I}_{z}}}+\,{\frac {{\left({y}_{K}\right)}^{T}\int _{A}^{}{N}_{I}^{T}\cdot {N}_{J}\,dA\,{u}_{K}}{{I}_{z}}}$

Recordando que el producto de inercia vale cero sobre los ejes principales, se pueden calcular los centros de cortante desacoplando las ecuaciones anteriores de la siguiente manera. Si la sección presenta algún eje de simetría paralelo al eje ${\textstyle y}$ o al eje ${\textstyle z}$, entonces ${\textstyle {I}_{yz}=}$$0$ y los centros de cortante valen:

${y}_{sc}=-\,{\frac {{\left({z}_{K}\right)}^{T}\int _{A}^{}{N}_{I}^{T}\cdot {N}_{J}\,dA\,{u}_{K}}{{I}_{y}}}$	(28)
${z}_{sc}=\,{\frac {{\left({y}_{K}\right)}^{T}\int _{A}^{}{N}_{I}^{T}\cdot {N}_{J}\,dA\,{u}_{K}}{{I}_{z}}}$

En cualquier otro caso, los centros de cortante se pueden calcular si trabajamos con las inercias principales ${\textstyle {I}_{y}^{P}}$, ${\textstyle {I}_{Z}^{P}}$ donde el producto principal de inercia ${\textstyle {I}_{yz}^{P}=}$$0$. De esta manera, los centros de cortante se calculan de acuerdo con:

${y}_{sc}=-\,{\frac {{\left({z}_{K}\right)}^{T}\int _{A}^{}{N}_{I}^{T}\cdot {N}_{J}\,dA\,{u}_{K}}{{I}_{y}^{P}}}$	(29)
${z}_{sc}=\,{\frac {{\left({y}_{K}\right)}^{T}\int _{A}^{}{N}_{I}^{T}\cdot {N}_{J}\,dA\,{u}_{K}}{{I}_{z}^{P}}}$

Las ecuaciones anteriores revelan la necesidad de calcular los momentos de inercia para determinar los centros de cortante. Una de las primeras propuestas para este cálculo se atribuye a (Owen & Hinton, 1980), quienes sugirieron que la inercia podría obtenerse multiplicando el área de cada elemento finito por su distancia al cuadrado con respecto al eje de referencia. Aunque esta aproximación proporciona resultados cercanos al valor analítico, una metodología más precisa fue desarrollada por (Haukaas, 2020). Este último investigador consideró que las coordenadas de la sección transversal no estaban referidas al centroide de la sección, logrando así resultados prácticamente idénticos al valor analítico. No obstante, una desventaja de este enfoque surge cuando la malla no es simétrica, generando valores distintos de cero del producto de inercia para secciones transversales con un eje de simetría paralelo a uno de los ejes coordenados. Es importante señalar que la magnitud de esta discrepancia puede variar según el sistema de unidades elegido, lo que puede conducir a resultados precisos en algunas ocasiones (dentro de una tolerancia aceptable) y, en otras, introducir errores significativos no deseados.

Con el objetivo de superar esta limitación y lograr que el cálculo de las inercias no esté sujeto ni al sistema de unidades ni a mallas estructuradas, se sugiere calcularlas siguiendo el enfoque presentado por (Hally, 1987). Este autor propuso expresiones universales aplicables a cualquier polígono, proporcionando así una solución más generalizada y libre de dependencias respecto a unidades o estructuras de malla específicas, obteniéndose:

${I}_{z}={\frac {1}{12}}\sum _{}^{}\left({y}_{i}{z}_{i+1}-{y}_{i+1}{z}_{i}\right)\left({y}_{i}^{2}+\right.$$\left.{y}_{i}{y}_{i+1}+{y}_{i+1}^{2}\right)$	(30)
${I}_{y}={\frac {1}{12}}\sum _{}^{}\left({y}_{i}{z}_{i+1}-{y}_{i+1}{z}_{i}\right)\left({z}_{i}^{2}+\right.$$\left.{z}_{i}{z}_{i+1}+{z}_{i+1}^{2}\right)$
${I}_{yz}={\frac {1}{24}}\sum _{}^{}\left({y}_{i}{z}_{i+1}-{y}_{i+1}{z}_{i}\right)\left({y}_{i}{z}_{i+1}+\right.$$\left.2{y}_{i}{z}_{i}+2{y}_{i+1}{z}_{i+1}+{y}_{i+1}{z}_{i}\right)$

Dado que el polígono que describe el contorno de la sección transversal en cuestión puede presentar huecos o intersecciones entre líneas, la aplicación de las ecuaciones (30) en estos escenarios queda excluida. En nuestro caso, el uso apropiado de las ecuaciones (30) implica aplicarlas individualmente a cada elemento finito de la sección transversal y luego sumar las contribuciones de cada uno para obtener el momento de inercia correcto. Es crucial destacar que las coordenadas utilizadas deben hacer referencia al centroide geométrico de la sección para garantizar la precisión buscada.

Con los cálculos anteriores se encuentra la función final de alabeo en cada nodo sumando a los desplazamientos la constante de corrección por alabeo más el centro de cortante ${\textstyle {y}_{sc}}$ multiplicado por las coordenadas en ${\textstyle z}$ del nodo referidas al centroide menos el centro de cortante ${\textstyle {z}_{sc}}$ multiplicado por las coordenadas en ${\textstyle y}$ del nodo también referidas al centroide.

Cálculo de las constantes de torsión

Una vez que se ha obtenido la función final por alabeo ${\textstyle \Omega (y,z)}$, el cálculo de la constante de torsión por alabeo ${\textstyle {C}_{w}}$ se obtiene fácilmente mediante:

${C}_{w}=\int _{A}^{}{\Omega }^{2}dA={u}_{K}\,\int _{A}^{}{N}_{I}^{T}\cdot {N}_{J}dA\,\,{u}_{K}$

(31)

Por otro lado, para encontrar la constante de torsión de Saint-Venant ${\textstyle J}$, recordemos que la relación existente entre el momento de torsión ${\textstyle T}$ y el giro de la sección ${\textstyle \theta }$ es mediante:

$T=G\cdot J\cdot \theta$

(32)

Ya que hemos supuesto que el giro es unitario, entonces podemos igualar la ecuación (25) con la ecuación (32), de donde despejamos ${\textstyle J}$ obteniedo:

$J=\,\int _{A}^{}\left[{\begin{matrix}-{z}^{h}&{y}^{h}\end{matrix}}\right]\cdot \left\{{\mathit {\boldsymbol {B}}}\cdot {\mathit {\boldsymbol {u}}}+\left[{\begin{matrix}-{z}^{h}\\{y}^{h}\end{matrix}}\right]\right\}\,dA$

(33)

La determinación precisa de las constantes de torsión, representadas por ${\textstyle {C}_{w}}$ y ${\textstyle J}$, para cualquier sección transversal abre las puertas a la innovación en el diseño estructural. Estas constantes proporcionan un entendimiento profundo de cómo una estructura resiste y transmite fuerzas torsionales. Este conocimiento no solo es esencial para garantizar la integridad estructural, sino que también brinda la oportunidad de explorar y crear formas geométricas únicas y estéticas que pueden embellecer significativamente las edificaciones modernas. Al comprender las características torsionales de diferentes secciones transversales, los ingenieros civiles y estructurales pueden experimentar con nuevas configuraciones, contribuyendo así a la creación de espacios habitacionales innovadores y visualmente atractivos que cumplen con los estándares contemporáneos de funcionalidad y estética. La investigación y aplicación cuidadosa de estas constantes no solo mejoran la eficiencia estructural, sino que también permiten la incorporación de elementos distintivos que enriquecen la identidad visual de las edificaciones, contribuyendo a un paisaje urbano más diverso y emocionante.

Estás constantes también permiten mejorar el comportamiento mecánico de los elementos estructurales, ya que particularmente hablando de ${\textstyle {C}_{w}}$, se agregan nuevos miembros a la matriz de rigidez convencional de un elemento, que de acuerdo con (Haukaas, 2020) es:

$K=\left[{\begin{matrix}{k}_{aa}&{k}_{ab}\\{k}_{ba}&{k}_{bb}\end{matrix}}\right]$

(34)

En cada una de estas cuatro submatrices ${\textstyle ({k}_{aa},\,{k}_{ab},\,{k}_{ba},\,{k}_{bb})}$ aparece un nuevo término que involucra la constante por alabeo ${\textstyle {C}_{w}}$ de la sección transversal y que da lugar a un nuevo grado de libertad conocido como bimomento en la sección. Cada una de estas cuatro matrices es igual a:

${k}_{aa}=\left[{\begin{matrix}{\frac {EA}{L}}&0&0&0&0&0&0\\0&{\frac {12E{I}_{z}}{{L}^{3}}}&0&0&0&{\frac {6E{I}_{z}}{{L}^{2}}}&0\\0&0&{\frac {12E{I}_{y}}{{L}^{3}}}&0&-{\frac {6E{I}_{y}}{{L}^{2}}}&0&0\\0&0&0&{\frac {12E{C}_{w}}{{L}^{3}}}+{\frac {6GJ}{5L}}&0&0&{\frac {6E{C}_{w}}{{L}^{2}}}+{\frac {GJ}{10}}\\0&0&-{\frac {6E{I}_{y}}{{L}^{2}}}&0&{\frac {4E{I}_{y}}{L}}&0&0\\0&{\frac {6E{I}_{z}}{{L}^{2}}}&0&0&0&{\frac {4E{I}_{z}}{L}}&0\\0&0&0&{\frac {6E{C}_{w}}{{L}^{2}}}+{\frac {GJ}{10}}&0&0&{\frac {4E{C}_{w}}{{L}^{3}}}+{\frac {2GJL}{15}}\end{matrix}}\right]$

(35)

${k}_{ab}=\left[{\begin{matrix}-{\frac {EA}{L}}&0&0&0&0&0&0\\0&-{\frac {12E{I}_{z}}{{L}^{3}}}&0&0&0&-{\frac {6E{I}_{z}}{{L}^{2}}}&0\\0&0&-{\frac {12E{I}_{y}}{{L}^{3}}}&0&-{\frac {6E{I}_{y}}{{L}^{2}}}&0&0\\0&0&0&-{\frac {12E{C}_{w}}{{L}^{3}}}-{\frac {6GJ}{5L}}&0&0&{\frac {6E{C}_{w}}{{L}^{2}}}+{\frac {GJ}{10}}\\0&0&{\frac {6E{I}_{y}}{{L}^{2}}}&0&{\frac {2E{I}_{y}}{L}}&0&0\\0&-{\frac {6E{I}_{z}}{{L}^{2}}}&0&0&0&{\frac {2E{I}_{z}}{L}}&0\\0&0&0&-{\frac {6E{C}_{w}}{{L}^{2}}}-{\frac {GJ}{10}}&0&0&{\frac {2E{C}_{w}}{L}}-{\frac {GJL}{30}}\end{matrix}}\right]$

(36)

${k}_{ba}=\left[{\begin{matrix}-{\frac {EA}{L}}&0&0&0&0&0&0\\0&-{\frac {12E{I}_{z}}{{L}^{3}}}&0&0&0&{\frac {6E{I}_{z}}{{L}^{2}}}&0\\0&0&-{\frac {12E{I}_{y}}{{L}^{3}}}&0&{\frac {6E{I}_{y}}{{L}^{2}}}&0&0\\0&0&0&-{\frac {12E{C}_{w}}{{L}^{3}}}-{\frac {6GJ}{5L}}&0&0&-{\frac {6E{C}_{w}}{{L}^{2}}}-{\frac {GJ}{10}}\\0&0&-{\frac {6E{I}_{y}}{{L}^{2}}}&0&{\frac {2E{I}_{y}}{L}}&0&0\\0&{\frac {6E{I}_{z}}{{L}^{2}}}&0&0&0&{\frac {2E{I}_{z}}{L}}&0\\0&0&0&{\frac {6E{C}_{w}}{{L}^{2}}}+{\frac {GJ}{10}}&0&0&{\frac {2E{C}_{w}}{L}}-{\frac {GJL}{30}}\end{matrix}}\right]$

(37)

${k}_{bb}=\left[{\begin{matrix}{\frac {EA}{L}}&0&0&0&0&0&0\\0&{\frac {12E{I}_{z}}{{L}^{3}}}&0&0&0&-{\frac {6E{I}_{z}}{{L}^{2}}}&0\\0&0&{\frac {12E{I}_{y}}{{L}^{3}}}&0&{\frac {6E{I}_{y}}{{L}^{2}}}&0&0\\0&0&0&{\frac {12E{C}_{w}}{{L}^{3}}}+{\frac {6GJ}{5L}}&0&0&-{\frac {6E{C}_{w}}{{L}^{2}}}-{\frac {GJ}{10}}\\0&0&{\frac {6E{I}_{y}}{{L}^{2}}}&0&{\frac {4E{I}_{y}}{L}}&0&0\\0&-{\frac {6E{I}_{z}}{{L}^{2}}}&0&0&0&{\frac {4E{I}_{z}}{L}}&0\\0&0&0&-{\frac {6E{C}_{w}}{{L}^{2}}}-{\frac {GJ}{10}}&0&0&{\frac {4E{C}_{w}}{{L}^{3}}}+{\frac {2GJL}{15}}\end{matrix}}\right]$

(38)

El lector podrá advertir que todo lo concerniente a ${\textstyle {C}_{w}}$ representa un concepto novedoso, y se observará la presencia tanto de un renglón adicional como de una columna extra en cada submatriz, ampliando así la matriz local de la viga de 12x12 a una matriz local de 14x14. Este tipo específico de matriz rara vez se encuentra incluido en los softwares comerciales convencionales, lo cual genera un espacio propicio para el avance del conocimiento en el campo de la ingeniería estructural. La inclusión de estas dimensiones adicionales en la matriz local no solo expande las posibilidades de modelado, sino que también plantea desafíos que, al ser abordados y comprendidos, contribuirán al desarrollo y perfeccionamiento de herramientas y enfoques analíticos en el ámbito estructural.

Ejemplos de validación

Sección circular

En este ejemplo se determinarán las constantes de torsión de una sección circular que tiene un diámetro de ${\textstyle 1\,m}$. La respuesta analítica para la constante de torsión de Saint-Venant es ${\textstyle J=}$$\pi {D}^{4}/32$, que en nuestro caso valdrá ${\textstyle J=}$$0.09817\,{m}^{4}$ mientras que la constante de torsión por alabeo es ${\textstyle {C}_{w}=}$$0$. Se añade la comparación del producto de inercia ${\textstyle {I}_{yz}}$ ya que es el valor que tiene más variación e importancia porque cuando no es cero, genera una rotación de los ejes. El resto de los valores obtenidos (ej.: área ${\textstyle A}$, momentos de inercia ${\textstyle {I}_{y}}$, ${\textstyle {I}_{z}}$, centros de cortante ${\textstyle {y}_{sc}}$, ${\textstyle {z}_{sc}}$) aunque se calculan no se reportan por ser bastante parecidos. Los resultados numéricos con tres diferentes mallados (Fig. 1) y comparando con los resultados de (Haukaas, 2020) se muestran a continuación.

	Malla Gruesa			Malla Media			Malla Fina
	${I}_{yz}$	$J$	${C}_{w}$	${I}_{yz}$	$J$	${C}_{w}$	${I}_{yz}$	$J$	${C}_{w}$
Haukaas	5.09e-05	0.09595	2.14e-14	4.45e-05	0.09737	1.81e-14	1.01e-06	0.09790	1.54e-14
Valdés	7.58e-11	0.09595	2.14e-14	3.83e-18	0.09737	1.81e-14	5.42e-19	0.09790	1.54e-14

Fig. 1 Mallas utilizadas

Si el mismo ejemplo se repite, pero cambiando las unidades a centímetros, el diámetro será de ${\textstyle 100\,cm}$. La constante de torsión de Saint-Venant valdrá ${\textstyle J=}$$9,817,477\,c{m}^{4}$. Se hará la comparación utilizando solamente una malla fina entre nuestra propuesta y la de Haukaas.

	Malla Fina
	${I}_{yz}$	$J$	${C}_{w}$
Haukaas	101.02	9,789,572	4.41e-06
Valdés	5.09e-11	9,789,572	4.41e-06

Como se puede ver en este último ejemplo, la propuesta de Haukaas genera productos de inercia muy diferentes de cero mientras que nuestra propuesta mejora este cálculo. Ya que se trata de una sección circular, la rotación a que da lugar el cálculo de las constantes de torsión no se ve afectada en este ejemplo debido a que, aunque los ejes se roten, las inercias no cambian. Por último, el error entre la solución analítica y la solución numérica con la malla más fina de este ejemplo es de solo el 0.28%.

Sección rectangular

En el segundo ejemplo se determinarán las constantes de torsión de una sección rectangular con unidades en centímetros y que tiene una base ${\textstyle b=}$$30\,cm.$ y una altura ${\textstyle h=60\,cm.}$ La fórmula para calcular la torsión de Saint-Venant para este ejemplo es: ${\textstyle J=}$$0.229\cdot h\cdot {b}^{3}=370,980$. Esta fórmula no es exacta y considera que puede tener errores del 4%. Los resultados que obtenidos son:

	Malla Fina
	${I}_{yz}$	$J$	${C}_{w}$
Haukaas	7.4860	375,969	14,863,370
Valdés	1.81e-11	375,969	14,863,370

La malla utilizada en este ejemplo aparece en la siguiente figura.

Fig. 2 Malla utilizada

En este ejemplo, el error en el cálculo de la constante de torsión de Saint-Venant con la solución analítica (aproximada) y la solución numérica es del 1.34%. Se aprecia además que hay una importante diferencia en el cálculo del producto de inercia ${\textstyle {I}_{yz}}$, resultando más exacto con nuestra propuesta.

Sección arbitraria de puente

Con la finalidad de mostrar las capacidades de nuestra propuesta, se combinado con el software de procesamiento de elementos finitos (GiD – The personal pre and post processor, 2023). Los resultados de una sección transversal que pudiera servir para un puente se muestran en las figuras 3-7. Por cuestión de programación, los ejes ${\textstyle {I}_{y},\,{I}_{z}}$ se cambian de nombre por ${\textstyle {I}_{x},\,{I}_{y}}$.

Fig. 3 Área de la sección transversal ${\textstyle A}$

Fig. 4 Momento de inercia ${\textstyle {I}_{x}}$

Fig. 5 Momento de inercia ${\textstyle {I}_{y}}$

Fig. 6 Constante de torsión ${\textstyle J}$

Fig. 7 Constante de torsión ${\textstyle {C}_{w}}$

Conclusiones

El análisis de marcos espaciales implica la necesidad de comprender a fondo las propiedades de las secciones transversales bajo estudio, destacando especialmente el área, los momentos de inercia y las constantes de torsión. A medida que la construcción moderna incorpora cada vez más secciones compuestas o con geometrías caprichosas, el cálculo preciso de estas propiedades se convierte en un desafío fundamental para los ingenieros que se embarcan en el análisis de dichas estructuras. Para abordar esta complejidad, diversos investigadores han desarrollado técnicas numéricas avanzadas capaces de calcular con precisión las propiedades de cualquier tipo de sección transversal.

En este contexto, el presente trabajo no solo se suma a esta línea de investigación, sino que también destaca al presentar una mejora significativa en la precisión del cálculo de dichas propiedades. Esta mejora se respalda con ejemplos concretos que demuestran la eficacia del método propuesto. De este modo, se contribuye no solo a la comprensión teórica de las propiedades de las secciones transversales en marcos espaciales, sino también a la aplicación práctica de herramientas más precisas y avanzadas en el análisis estructural.

Referencias

Boresi, A. P., Schmidt, R. J., & Sidebottom, O. M. (1993). Advanced Mechanis of Materials. John Wiley & Sons, Inc.

Euler, L. (1744). Methodus Inveniendi Lineas Curvas Maximi Minimive Proprietate Gaudentes.

GiD – The personal pre and post processor. (2023). Barcelona (Spain); CIMNE. accessed 2023 Dec 13.

Hally, D. (1987). Calculation of the Moments of Polygons. Technical Memorandum 87/209. Canadian National Defense.

Haukaas, T. (2020). Computational Cross-section Analysis.

Owen, D. R., & Hinton, E. (1980). Finite Elements in Plasticity. Theory and Practice. Swansea, U.K.: Pineridge Press Limited.

Pilkey, W. D. (2002). Analysis and Design of Elastic Beams. Computational Methods. John Wiley & Sons.

Saint-Venant, A. C. (1855). Memoire sur la Torsion des Prismes.

Timoshenko, S. (1934). Theory of Elasticity.