Analisis No Lineal
Analisis No Lineal
Analisis No Lineal
INGENIERIA
Introducción al Análisis No Lineal de Estructuras
Texto del curso Análisis No Lineal de Estructuras
Este documento es publicado bajo una licencia Creative Commons Attribution-ShareAlike 4.0
International License. Ver detalles de la licencia en creativecommons.org/licenses/by-sa/4.0.
Documento producido usando software libre: LATEX, TeXstudio, Geany, Inkscape y GNU-Octave.
Prefacio
En este documento se presentan los conceptos impartidos en la primer edición
del curso Análisis No Lineal de Estructuras, incluido en el Posgrado en Ingeniería
Estructural de la Facultad de Ingeniería de la Universidad de la República. A lo largo
del texto se hace referencia a bibliografía especializada en la cual se pueden encon-
trar desarrollos detallados que complementan lo presentado, no obstante, cualquier
error es responsabilidad de los autores.
Con el objetivo de facilitar la asimilación de conceptos, los autores desarrolla-
ron una herramienta informática abierta (usando GNU-Octave) que permite resolver
problemas de análisis no lineal de estructuras de barras articuladas. La herramienta,
llamada Open Nonlinear Structural Analysis System (ONSAS), fue desarrollada con el
objetivo de mostrar al lector una implementación de algunos de los métodos de aná-
lisis descritos en el texto, así como también poder brindar una herramienta gratuita
que permita resolver una amplia gama de problemas estructurales de interés.
El Capítulo 1 comienza con una descripción de problemas estructurales que re-
quieren modelos no lineales, continuando con la introducción a los procedimientos
numéricos necesarios para resolver las ecuaciones no lineales correspondientes.
En el Capítulo 2 se presentan herramientas analíticas para el análisis de estruc-
turas geométricamente no lineales y su aplicación en ejemplos canónicos. También
se presenta una formulación para la resolución de problemas estructurales con no
linealidad geométrica basada en el Método de los Elementos Finitos.
En el Capítulo 3 se aborda el análisis de estructuras formadas por materiales en
los cuales la relación entre tensión y deformación es no lineal. Se presentan modelos
de comportamiento constitutivo: hiperelástico, lineal a tramos y elasto-plástico.
El Capítulo 4 introduce al lector al análisis dinámico de estructuras con com-
portamiento no lineal. Para ello, se realiza una introducción al análisis dinámico en
modelos estructurales lineales para luego pasar a considerar la dinámica de modelos
estructurales no lineales.
Los autores confían en que este texto puede resultar de utilidad para toda persona
interesada en comprender los conceptos básicos del análisis no lineal de estructuras
siguiendo literatura en español. El material, condensa una parte de la extensa biblio-
grafía que existe sobre el tema en cuestión, por lo que se invita al lector a continuar
la lectura en las referencias citadas a lo largo del texto. Los autores agradecen a los
colegas, amigos y familiares que de diversas formas han apoyado el desarrollo de
este material.
Montevideo, Uruguay, 22 de diciembre de 2017.
A los Profesores Adam Sadowski y Julio Borghi.
Lista de contenidos
Prefacio iii
2 No Linealidad Geométrica 25
2.1 Principio de Mínima Energía Potencial Total . . . . . . . . . . 25
2.1.1 Definición de Energía Potencial Total . . . . . . . . . . 25
2.1.2 Equilibrio y Estabilidad . . . . . . . . . . . . . . 27
2.1.3 Clasificación de Puntos Críticos . . . . . . . . . . . 29
2.1.4 Ejemplo de Análisis de Equilibrio y Estabilidad . . . . . . . 33
2.2 Principio de los Trabajos Virtuales . . . . . . . . . . . . . 34
2.2.1 PTV para Partículas y Sistemas de Partículas . . . . . . . 35
2.2.2 PTV para Cuerpos Rígidos . . . . . . . . . . . . . 36
2.2.3 PTV para Estructuras con Componentes Deformables. . . . . 38
2.2.4 PTV para Elementos de Tipo Biela. . . . . . . . . . . 41
2.3 Medidas No Lineales de Deformación . . . . . . . . . . . . 42
2.3.1 Estructura de Referencia . . . . . . . . . . . . . 42
2.3.2 Deformación Unitaria Ingenieril Rotada . . . . . . . . . 43
2.3.3 Deformación Unitaria de Green . . . . . . . . . . . 44
2.3.4 Deformación Unitaria de Logarítmica Rotada . . . . . . . 45
2.3.5 Comparación de las Medidas de Deformación . . . . . . . 46
2.4 Método de Elementos Finitos Incremental. . . . . . . . . . . 48
2.4.1 Método de los Elementos Finitos en análisis lineal de reticulados . 49
2.4.2 Método de Newton-Raphson aplicado al PTV . . . . . . . 51
2.4.3 Método de Newton-Raphson Incremental . . . . . . . . 63
2.4.4 Otras técnicas avanzadas . . . . . . . . . . . . . 65
2.4.5 Ejemplo: Arco bi-articulado . . . . . . . . . . . . . 66
2.5 Estabilidad estructural . . . . . . . . . . . . . . . . 69
2.5.1 Clasificación de equilibrio . . . . . . . . . . . . . 69
2.5.2 Análisis lineal de pandeo (LBA) . . . . . . . . . . . 71
2.5.3 Análisis no lineal de pandeo . . . . . . . . . . . . 72
2.5.4 Ejemplo: Análisis de estabilidad de cercha de Von Mises . . . . 72
2.5.5 Ejemplo: modo de pandeo asimétrico de arco bi-articulado . . . 75
3 No Linealidad Material 77
3.1 Elasticidad no lineal . . . . . . . . . . . . . . . . . 77
3.1.1 Hiperelasticidad en sólidos . . . . . . . . . . . . . 77
3.1.2 Hiperelasticidad en barras . . . . . . . . . . . . . 82
3.1.3 Análisis de reticulados bi-modulares . . . . . . . . . . 83
3.2 Introducción a la Teoría de Plasticidad . . . . . . . . . . . . 84
3.2.1 Teoría unidimensional de plasticidad . . . . . . . . . . 85
3.2.2 Método numérico de análisis elastoplástico . . . . . . . . 90
3.2.3 Implementación y ejemplos . . . . . . . . . . . . 93
3.2.4 Aplicación a reticulados . . . . . . . . . . . . . . 96
3.3 Solicitaciones últimas en sección de hormigón armado . . . . . . . 96
B ONSAS 137
B.1 Descripción de código . . . . . . . . . . . . . . . . . 137
B.2 Ejemplos de archivos de entrada . . . . . . . . . . . . . 157
B.2.1 Arco bi-articulado . . . . . . . . . . . . . . . 157
B.2.2 Ejemplo de cercha de Von Mises . . . . . . . . . . . 161
Capítulo 1: Motivación y Métodos Numéricos
Carga
1 LA
2 LBA 3 GNA
5 GNIA
4 MNA
6 GMNIA
Desplazamiento
Enfoque
Métodos Incrementales,
Métodos Iterativos,
g(x, λ) = 0. (1.1)
f (x) − λv = 0. (1.2)
Las soluciones de los sistemas (1.1) o (1.2) determinan el conjunto de todas las
parejas de valores (λ, x) que satisfacen dichos sistemas no lineales. De forma equi-
valente, las soluciones se pueden considerar como trayectorias que pasan por (0, 0)
en el espacio de configuraciones de la estructura S ≡ R+ × Rn y que satisfacen la
Ecuación (1.2) en todo punto de dichas trayectorias. Aquí y en los capítulos siguien-
tes se llamará a dichas trayectorias como Curvas de Carga-Desplazamiento.
A continuación, se describen métodos para hallar soluciones numéricas de (1.1)
o en su defecto (1.2).
f (x(λ)) − λv = 0. (1.3)
Fx (x(λ))ẋ(λ) − v = 0, (1.4)
donde Fx es una matriz cuadrada (llamada matriz tangente) cuya entrada ij está
dada por (Fx )i,j = ∂x
∂fi
j
. Asumiendo que Fx es invertible se obtiene la siguiente
E.D.O. de primer orden:
ẋ(λ) = F−1
x (x(λ))v. (1.5)
Dada una E.D.O. de primer orden genérica con variable independiente t, ẋ(t) =
h(x(t), t) y una aproximación numérica aceptable de la solución de la E.D.O. en tk ,
xk ≈ x(tk ), se puede escribir el siguiente desarrollo de Taylor de x(t) en tk (con
ξ ∈ (tk , t)):
1
x(t) = x(tk ) + ẋ(tk )(t − tk ) + ẍ(ξ)(t − tk )2 . (1.7)
2
Evaluando la expresión dada por la Ecuación (1.7) en t = tk + ∆t y usando que
ẋ = h(x, t) se tiene:
1
x(tk + ∆t) = x(tk ) + h(x(tk ), tk )∆t + ẍ(ξ)∆t2 . (1.8)
2
Despreciando el término de segundo orden y llamando xk+1 y xk a las aproxi-
maciones de x(tk + ∆t) y x(tk ) respectivamente, se obtiene:
(
xk+1 = xk + ∆t · h(xk , tk )
F.Euler (1.9)
x0 = x(0)
0.2 0.4
Solución Numérica
Solución Exacta
0.3
0.15
0.2
0.1 0.1
λ
0
0.05
-0.1
Solución Numérica
Solución Exacta
0 -0.2
0 0.1 0.2 0.3 0.4 0.5 -0.5 0 0.5 1 1.5 2 2.5
x x
Figura 1.2: Soluciones obtenidas usando Euler Hacia Adelante. Izquierda: solución
para λ∗ = 0,18. Derecha: solución para λ∗ = 0,26.
Código 1.1: Implementación del Método de Euler para resolución del Ejemplo Nu-
mérico 1.
1 %% Ejemplo Numerico 1 : Solucion Incremental % %
2 clc, clear all, close all
3
4 % Condiciones Iniciales
5 x0 = 0; L0 = 0;
6
7 % Definicion de v
8 v = 1;
9
10 % Definicion de f(x)
11 f = @(x) x−3/2*x.^2+1/2*x.^3;
12
13 % Definicion de F_x(x)
14 Fx = @(x) 1−3*x+3/2*x.^2;
15
16 % Definicion de h(x,lambda)
17 h = @(x) Fx(x)\v;
18
19 %% % %Solucion Incremental con Forward Euler % % % %
20 % Defino hasta que Lambda incrementamos
21 Lmax = .25;
22
23 x(:,1) = x0 ; L(1)=0;
24
25 % Defino Incremento Delta Lambda
26 DL = .02;
27
28 % Comienza secuencia de incrementos
29 i=1;
30 while L(end)<Lmax
31 x(:,i+1) = x(:,i) + DL * h(x(:,i)); % Forward Euler
32 L(i+1) = L(i)+DL;
33 i = i+1;
34 end
35
36 % Genero graficas de las soluciones
37 subplot(1,2,1) % Solucion Por debajo de punto limite
38 plot(x(1,1:end−2),L(1:end−2),'*−k',0:.01:.5,f(0:.01:.5),'−r')
39 xlabel('x'), ylabel('$\lambda$')
40 legend('Solucion Numerica','Solucion Exacta')
41
42 subplot(1,2,2) % Solucion Por encima de punto limite
43 plot(x(1,:),L,'*−k',0:.01:2.25,f(0:.01:2.25),'−r')
44 xlabel('x'), ylabel('$\lambda$')
45 legend('Solucion Numerica','Solucion Exacta')
10 Capítulo 1. Motivación y Métodos Numéricos
Método de Newton-Raphson
λv Recta tangente
λk v
(0, 0) (0)
xk
(1)
xk x
(0)
La ecuación de recta de la tangente en el punto xk está dada por:
(0) (0) (0)
λv = f (xk ) + f 0 (xk )(x − xk ), (1.14)
Es claro a partir de la Ecuación (1.18) que este método no puede iterar si la de-
12 Capítulo 1. Motivación y Métodos Numéricos
rivada primera de f (x) es igual a cero en el punto donde se traza la tangente. Esta
limitación también estará presente en el caso de Newton-Raphson para varias varia-
bles.
Actividad
Verificar que si x ∈ R el método dado por la Ecuación (1.21) se
reduce al presentado en la Ecuación (1.18).
(0)
Asumiendo que el punto de arranque (xk ) está suficientemente cerca de la so-
lución buscada, la sucesión de iterados que genera N-R converge con orden 2. Esto
significa que las distancias entre cada iterado y la solución exacta (xk ) verifican
(i+1) (i)
kxk − xk k ≈ βkxk − xk k2 , siendo β > 0 la velocidad de convergencia. En este
documento no se presentarán formalmente los conceptos de orden y velocidad de
convergencia, ver (Quarteroni et al., 2007) por dichas definiciones.
(i)
El método de N-R presentado en (1.21) requiere que Fx (xk ) sea invertible para
poder calcular el paso iterativo. Más en general, para poder tener pasos iterativos
con precisión confiable se deberá cumplir que el número de condición de la matriz
(i)
Fx (xk ) sea pequeño. Ver la Sección 7.6 del libro (Björk, A. and Dahlquist, 2008) por
Sección 1.2. Métodos Numéricos para Ecuaciones No Lineales 13
El método de N-R modificado presenta una alternativa para no tener que evaluar
y factorizar la matriz Fx en cada iteración. Precisamente esto es lo que define al
(0)
método en cuestión, el mismo deja fija la matriz Fx (xk ) a lo largo de las iteraciones.
Algunas herramientas computacionales, como por ejemplo el software SAP 2000® ,
resuelven problemas con no linealidad geométrica realizando un cierto número de
iteraciones de NR modificado para luego pasar a realizar iteraciones Newton-Raphson
si no se alcanza convergencia.
La iteración de Newton-Raphson modificado puede ser escrita como:
(0) (i) (i)
F x (x k )∆x k = − (f (x k ) − λ k v
N-R Modif. x(i+1) = x(i) + ∆x(i) (1.23)
k k k
x(0)
k
(0)
Esto permite factorizar la matriz Fx (xk ) al comienzo de las iteraciones y alma-
cenar dichos factores. Luego, en cada paso iterativo se deberá resolver una sustitu-
(i)
ción hacia adelante y una sustitución hacia atrás para hallar ∆xk . El costo compu-
tacional de dichas sustituciones es un orden de magnitud menor (en n) que el de la
(0)
factorización de la matriz Fx (xk ).
14 Capítulo 1. Motivación y Métodos Numéricos
Métodos Cuasi-Newton
λv Recta secante
λk v
(1)
(0, 0) xk
(0)
xk x(2)
k x
Es claro que se deben especificar dos puntos de arranque para poder comenzar
a iterar con el Método de la Secante.
Sección 1.2. Métodos Numéricos para Ecuaciones No Lineales 15
Actividad
Verificar que la expresión de la iteración del Método de la secante
está dada por:
(i) (i−1)
(i+1) (i) xk − xk
(i)
xk = xk − (i) (i−1)
f (xk ) − λk v (1.24)
f (xk ) − f (xk )
√
El Método de la Secante tiene orden de convergencia super-lineal (Orden =
2 ), lo cual lo posiciona como más lento que N-R pero más rápido que N-R Mo-
1+ 5
dificado.
La generalización del Método de la Secante a sistemas de ecuaciones no lineales
corresponde al Método de Broyden. En ese método, en lugar de la matriz tangente
Fx utilizada en N-R, se utiliza una matriz secante que verifica una relación en los
incrementos de f (x) y x de la forma: ∆fk = Bk ∆xk . La condición de matriz se-
cante no es suficiente para definir la completamente y es ahí donde surgen distintas
variantes. En la Sección 7.1.4 de (Quarteroni et al., 2007) se presenta una definición
más detallada del método.
Otro método cuasi-Newton que deriva de la condición secante mencionada an-
teriormente es el Método BFGS, cuya formulación puede encontrarse en la Sección
8.4.2 de (Bathe, 2014).
El motivo por el cual estos métodos son atractivos es que en ellos las inversas
de las matrices secantes se actualizan de forma computacionalmente económica a
partir de las anteriores. Se obtiene, por lo tanto, un compromiso entre velocidad y
economía computacional.
mientras que la iteración del Método de Newton-Raphson Modificado está dada por:
0.2 0.2
0.15 0.15
λ
0.1 0.1
0.05 0.05
0 0
0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5
x x
Actividad
Modificar los parámetros de la implementación presentada y resol-
ver el ejemplo para λk = 0,26 y otros dos valores superiores a 0,2.
Analizar los resultados y comparar con los obtenidos por el método
incremental.
Es claro a partir de lo anterior que estos métodos son esenciales para el análisis
computacional del colapso de estructuras. Algunas herramientas computacionales
como: ABAQUS® , LUSAS® , ANSYS® o ADINA® tienen la capacidad de llevar a cabo
análisis estáticos mediante este tipo de procedimientos.
En lo que sigue se presenta una versión básica del método de Longitud de Arco,
la cual no es una implementación eficiente pero tiene como ventaja la claridad en
su formulación. En la Sección 9.3 de (Crisfield, 1996) se presenta una exposición
completa de estos métodos y su historia.
La idea básica de estos métodos es considerar a x ∈ Rn y λ ∈ R+ como incógni-
tas. Se requiere, naturalmente, que esas incógnitas satisfagan la Ecuación (1.3). Sin
embargo, ésta condición por sí sola no es suficiente para tener un problema deter-
minado, ya que se puede comprobar que en ese caso se tienen n + 1 incógnitas y n
ecuaciones. Por lo tanto, los métodos de Longitud de Arco imponen una restricción
adicional de manera de obtener un problema determinado.
Dependiendo de cuál sea dicha restricción adicional y la implementación elegida,
se obtienen los distintos métodos que forman parte de la familia de Métodos de
longitud de arco. Ver en (Crisfield, 1996) Métodos de Longitud de Arco Linealizado,
Cilíndrico y Esférico. La restricción adicional está, tal como lo indica el nombre de
cada método, directamente relacionada al concepto de longitud del arco de la curva
carga-desplazamiento. En la Figura 1.6 se presenta un esquema geométrico de la
restricción de longitud de arco.
El incremento diferencial de la coordenada de longitud de arco s puede escribirse
en función de los incrementos de carga y desplazamiento como:
λv radio: ∆`
f (x)
λk v
sk
(0, 0) xk x
y que el arco de radio ∆l está centrado en ese punto. Con lo cual: ∆xk = xk − xk−1
y ∆λk = λk − λk−1 .
De esta manera el sistema de ecuaciones no lineales que se deberán satisfacer en
el método de Longitud de Arco propuesto son:
(
f (xk ) − λk v = 0
(1.30)
∆xTk ∆xk + ∆λ2k ψ 2 vT v − ∆l2 = 0
Un punto adicional por comentar que surge de analizar la Figura 1.6, es que el
círculo que queda definido por la restricción de Longitud de Arco corta a la curva de
carga-desplazamiento en dos puntos, uno avanza la solución y el otro la retrocede.
Dado lo anterior, es posible que la iteración converja al punto que hace que la so-
lución retroceda en lugar de avanzar. Una manera de controlar esto es estudiando a
medida que avanza la solución numérica si se pasó un punto crítico (Det(Fx ) = 0)
y en ese caso estudiar la necesidad de cambiar el signo del punto de arranque dado
por Euler Hacia Adelante.
Una observación final del método presentado es que la matriz del sistema lineal
dado en la Ecuación (1.34) que debe ser resuelto en cada iteración es no simétrica
y eso presenta una desventaja económica en cuanto a la resolución numérica del
sistema.
2.5
Solución Exacta
Solución Numérica
1.5
1
λ
0.5
-0.5
0 0.5 1 1.5 2 2.5 3
x
Notar que el método de Longitud de Arco utilizado logra pasar por dos puntos
críticos, uno en x ' 0,5 y otro en x ' 1,5.
También se puede observar en la Figura 1.7 que los puntos que componen la
solución mediante el procedimiento de Longitud de Arco se encuentran equiespa-
ciadas en la longitud de la curva correspondiente a la solución exacta. Se debe notar
que para que dichos puntos parezcan equiespaciados, el gráfico debe tener el mis-
mo escalado de ejes que la métrica usada como definición de longitud de arco en la
Ecuación (1.29).
Sección 1.2. Métodos Numéricos para Ecuaciones No Lineales 23
Actividad
Estudie cómo cambiar el código dado para poder resolver un sistema
de ecuaciones no lineales en Rn .
Imperial College London, titulada Sandwiches, snakes and stays: interacting with instabilities, disponible en:
https://youtu.be/NSFqwxW7gRo (último acceso diciembre de 2017).
26 Capítulo 2. No Linealidad Geométrica
donde U (q) es la energía elástica interna del sistema en la configuración dada por q y
WE (q) es el trabajo realizado por las fuerzas externas conservativas aplicadas sobre
el sistema desde una configuración de referencia arbitraria hasta la configuración q.
Un ejemplo sencillo mediante el cual se pueden fijar estos conceptos es el de un
sistema de una masa puntal (m) y un resorte elástico lineal de constante (K) sujeto
a la acción externa de la gravedad (g), como el que se muestra en la Figura 2.1.
q
L0
m
g
K
Kq 2
U (q) = + C1 y WE (q) = mgq + C2 , (2.2)
2
respectivamente. Las constantes C1 y C2 corresponden a las referencias que se to-
men tanto para definir energía elástica interna como trabajo de fuerzas conservati-
vas externas. Como se verá más adelante, estas constantes son intrascendentes para
estudiar el equilibrio y la estabilidad y pueden ser consideradas nulas por conve-
niencia.
La expresión de energía elástica interna del resorte se puede obtener consideran-
do el trabajo que se debe realizar sobre el resorte para llevarlo a una configuración
comprimida (o extendida), en particular evaluando la siguiente integral:
Z q Z q
Kq 2
U (q) = F (u) du = Ku du = + C1 (2.3)
q1 q1 2
Sección 2.1. Principio de Mínima Energía Potencial Total 27
Kq 2
V (q) = − mgq. (2.4)
2
Axioma 1
Axioma 2
1 00 1 000
V (q̂ + δ) − V (q̂) = V (q)|q̂ δ 2 + V (q)|q̂ δ 3 + ... ∀δ < (2.9)
2! 3!
con lo cual, se observa que si V 00 (q)|q̂ > 0 (y δ es pequeño) entonces se satisface la
desigualdad de la Ecuación (2.7) y por lo tanto q̂ es un punto estable.
Notar que podría ocurrir que V 00 (q)|q̂ = 0 y se podría tener estabilidad si
V (q)|q̂ = 0 y V (4) (q)q̂ > 0. Con lo cual, estabilidad no necesariamente im-
000
λ
λ2 Vλ2 (qi )
λc Vλc (qi )
λ1 Vλ1 (qi )
qi
Punto Límite
P punto límite
P Pc salto dinámico
θ α (0, 0) α−θ
K
(a) Esquema de estructura. (b) Curva carga-desplazamiento.
Una bifurcación simétrica estable corresponde a una situación en la cual una es-
tructura alcanza un punto crítico y a partir de ese punto, su curva de carga-desplazamiento
se bifurca en varias ramas: una inestable y otras dos estables.
El pandeo de columnas elásticas pertenece a esta categoría. En la Figura 2.4a se
puede ver una versión discreta de ese tipo de estructura. En estos casos, al alcanzar-
se la bifurcación las estructuras logran mantener la carga aunque sea con grandes
deformaciones (ej. regla elástica en compresión). En la Figura 2.4b se muestra la
curva carga-desplazamiento correspondiente. Notar que estas estructuras tampoco
son sensibles a imperfecciones, siempre y cuando el comportamiento permanezca
elástico.
Sección 2.1. Principio de Mínima Energía Potencial Total 31
Pc
θ con imperf.
(0, 0) θ
(a) Esquema de estructura. (b) Curva carga-desplazamiento.
P punto de bifurcación
θ sin imperf.
Pc
K
con imperf.
θ
(0, 0) θ
(a) Esquema de estructura. (b) Curva carga-desplazamiento.
Bifurcación Asimétrica
P
P punto de bifurcación
Pc
L K θ sin imperf.
con imperf.
L
(0, 0) θ
(a) Esquema de estructura. (b) Curva carga-desplazamiento.
V 00 (θ)|F = 4K − 2P L. (2.15)
Actividad
milar al de (McGuire et al., 1999), sin embargo se puede encontrar otros desarrollos
en libros como (Reddy, 2002) o (Bathe, 2014).
f2
f3
f1
Actividad
Mostrar el recíproco del resultado anterior, es decir, mostrar que el
PTV implica la segunda ley de Newton.
Se propone avanzar, usando el resultado anterior, hacia el planteo del PTV para
un cuerpo rígido.
La hipótesis central considerada es que el cuerpo que estudiado es rígido, por lo
tanto incluso ante desplazamientos virtuales el cuerpo permanece rígido. Se ilustra
la situación mediante el ejemplo de una barra rígida simplemente apoyada en sus
extremos y sometida a una fuerza lateral en la mitad de su largo, ver Figura 2.8a.
Para poder realizar un desplazamiento virtual de la barra rígida, se cambia el apo-
yo deslizante derecho por una fuerza estáticamente equivalente (la reacción vertical),
ver Figura 2.8b. Una vez liberado ese vínculo, se puede considerar un desplazamien-
to virtual δx de la barra rígida. Dicho desplazamiento debe satisfacer la condición de
borde en desplazamiento dada en el apoyo simple izquierdo y también la condición
de cuerpo rígido de la barra. Se concluye que el desplazamiento virtual debe ser un
giro rígido en torno al apoyo izquierdo.
Sección 2.2. Principio de los Trabajos Virtuales 37
L/2 L/2
P
δu
δu/2
L/2 L/2
RB
(b) Movimiento virtual de cuerpo rígido.
Dado que el PTV para una partícula también se aplica a un sistema de partículas,
se puede considerar el cuerpo rígido como un sistema de partículas con la propiedad
que las distancias relativas entre ellas se mantienen constantes. De lo anterior se
concluye que las fuerzas entre partículas no realizan trabajo virtual ya que durante
el desplazamiento virtual no hay cambios en las distancias relativas entre ellas y
por consiguiente no hay trabajo virtual. Por otro lado, la fuerza externa aplicada
sí realiza trabajo sobre el desplazamiento virtual y lo mismo vale para la reacción
derecha. A partir del razonamiento anterior y el PTV para sistemas de partículas, se
tiene:
δu
δWext = −P + RB δu = 0 ∀δu (2.22)
2
Se obtiene la siguiente ecuación:
P
− + RB δu = 0 ∀δu, (2.23)
2
P
RB = . (2.24)
2
La reacción vertical del apoyo izquierdo se puede determinar mediante un pro-
cedimiento idéntico al presentado para la reacción del apoyo derecho.
En conclusión, el PTV para cuerpos rígidos tiene la misma forma que el dado
para sistemas de partículas. Notar que se deben liberar vínculos de manera de poder
realizar los desplazamientos virtuales necesarios para determinar las ecuaciones de
equilibrio.
38 Capítulo 2. No Linealidad Geométrica
Actividad
Determinar todas las reacciones del ejemplo visto en esta sección,
liberando todos los vínculos de apoyo al mismo tiempo y generando
así un sistema de ecuaciones para determinar dichas reacciones.
Dados los conceptos presentados hasta aquí, se puede pasar a estudiar el PTV
aplicado a estructuras con componentes deformables. Puntualmente, se considerará
como ejemplo la cercha de Von Mises, la cual fue introducida en la Sección 2.1.3.
La estructura está compuesta por dos barras de largo L, rígidas, unidas por una
articulación perfecta. El extremo libre izquierdo de las barras se soporta mediante un
apoyo fijo y el extremo libre derecho mediante un apoyo deslizante. Ambos apoyos
se unen con un resorte, es decir, un cuerpo deformable (ver Figura 2.3a). Notar que no
se especifica la relación constitutiva del resorte. La estructura se somete a una carga
externa vertical descendente P y el ángulo que forma una de las barras inclinadas
con la horizontal para P = 0 es igual a α.
El paso clave para determinar el equilibrio del sistema estructural consiste en
sustituir el cuerpo deformable (resorte) por un conjunto de fuerzas estáticamente
equivalente y equilibrado aplicado al resto de la estructura.
Cuando la estructura está en una configuración de equilibrio, el resorte puede
ser sustituido por dos fuerzas opuestas de valor Fint aplicadas en los nodos de los
extremos. De esta forma se vuelve a tener una estructura compuesta de cuerpos
rígidos sometida a fuerzas externas (ver Figura 2.9), con lo cual se puede aplicar el
PTV tal como fue presentado en la sección anterior, es decir:
δWext = 0. (2.25)
P
L
w
Fint
θ Fint α
P δw = Fint δu
P L cos(θ)δθ = Fint 2L sin(θ)δθ
(2Fint sin(θ) − P cos(θ)) δθ = 0 (2.29)
El análisis realizado permite escribir el PTV para estructuras con cuerpos defor-
mables de manera general como se enuncia a continuación.
∂U
Fint = (2.32)
∂u
δWext − Fint δu = 0
∂U
δWext − δu = 0
∂u
δWext − δU = 0
δV = 0
Actividad
Determinar la ecuación de equilibrio de la cercha de Mises bajo la
hipótesis de resorte elástico lineal de constante (K). Comparar la
solución hallada con el resultado mostrado en la Sección 2.1.3.
En base a los casos del PTV presentados en las secciones anteriores, se puede pa-
sar a enunciar el PTV para el caso de elementos sólidos deformables de tipo Biela. Es
decir, se estudia el PTV para barras que solamente pueden resistir esfuerzos axiales.
Se planteará el PTV para una barra aislada y esto permitirá estudiar estructuras
compuestas de ensamblajes de barras de este tipo, simplemente sumando los trabajos
virtuales de todos los elementos considerados.
A partir del PTV para estructuras con componentes deformables visto en la sec-
ción anterior, se puede obtener una expresión para barras tipo biela. La expresión
del PTV es
δWint = δWext ∀δu ∈ V, (2.33)
considerando a la barra como una colección de componentes deformables de largo
Fint dδu
dx y definiendo σ = y δε = se puede llega a expresar el trabajo virtual
A dx
interno como una integral en el largo L de la barra en la configuración de equilibrio.
Se obtiene entonces:
Z
σδε dV = δuT fext ∀δu ∈ V. (2.34)
V
σAn
ln , An w
σAn z
l0 , A0
ln − lo
εE = (2.36)
lo
Notar que la deformación unitaria se mide respecto de una dirección que rota
en conjunto con la barra y no respecto de ejes globales cartesianos fijos. Usando la
geometría dada en la Figura 2.10:
1/2
ln = (z + w)2 + x2 (2.37)
Actividad
Notar que la expresión de la Ecuación (2.44) es bastante más simple que la ob-
tenida con la deformación unitaria ingenieril rotada. Además, se puede ver que si
ln ' lo entonces tanto las deformaciones unitarias como las relaciones de equilibrio
coinciden. En conclusión, en hipótesis de pequeñas deformaciones unitarias los re-
sultados obtenidos con ambas deformaciones unitarias coincidirán de manera muy
aproximada.
Actividad
Mostrar que las deformaciones unitarias vistas están relacionadas
según la identidad:
1
εL = ln(1 + εE ) = ln(1 + 2εG ).
2
Se asume ahora que se trabaja bajo pequeñas deformaciones y por lo tanto: Ao '
An y también lo ' ln , lo cual permite asumir que el volumen de la barra no cambia
al llegar a la configuración deformada: lo Ao ' ln An .
46 Capítulo 2. No Linealidad Geométrica
Se asume una relación constitutiva con módulo elástico constante E (σL = EεL )
y luego se expresa todo en función de la variable w:
p !
E(z + w)Ao lo (z + w)2 + x2
PL = ln (2.49)
(z + w)2 + x2 lo
15
PE : Def. Unit. Ing. Rotada
-5
-10
0 1000 2000 3000 4000 5000 6000
Desplazamiento: (-w)
cos(φe ) − sin(φe )
0 0
sin(φe ) cos(φe ) 0 0
ue = Qe ueloc con Qe = , (2.53)
e e
0 0 cos(φ ) − sin(φ )
0 0 sin(φe ) cos(φe )
donde se puede ver que φe es el ángulo medido en sentido antihorario desde el eje
x del sistema global al eje x local y Qe es la matriz de cambio de base del sistema
50 Capítulo 2. No Linealidad Geométrica
donde la matriz KeL es llamada matriz de rigidez del elemento e y está dada por:
Z
T
e
KL = (beL ) EAbeL dx. (2.58)
`e
El trabajo virtual interno del conjunto de la estructura puede ser escrito como la
suma del trabajo de todas las barras:
ne
T T
X
δWint (u) = (δue ) KeL ue = (δu) KL u, (2.59)
e=1
Tal como se vió en la Sección 2.2.4, al analizar bielas la única componente del
tensor de tensiones que realiza trabajo es la axial. Considerando la Ecuación (2.35) y
las hipótesis de pequeños cambios de área A ≈ A0 y de longitud ` ≈ `0 , se obtiene
la siguiente expresión del PTV para bielas:
Z
A0 σδε dx = δuT fext ∀δu ∈ V (2.63)
`0
donde σ representa la tensión de Green en cada punto y fext las fuerzas externas
nodales, que en el caso de un único elemento estaría dada por:
T
fext = [fext,1,x , fext,1,y , fext,2,x , fext,2,y ] . (2.64)
fext,2
δu2
y
u2
fext,1 2
δu1
u1
`0
x
Figura 2.12: Esquema de desplazamientos y fuerzas de elemento de barra.
en su expresión matricial:
−1 0 1 0 e −1 0 1 0 e
u2 − u1 = u x2 − x1 = x (2.68)
0 −1 0 1 0 −1 0 1
donde
1 0 −1 0
0 1 0 −1
Ge = (2.71)
.
−1 0 1 0
0 −1 0 1
∂g
g(u + δu) ≈ g(u) + (u)δu (2.75)
∂u
54 Capítulo 2. No Linealidad Geométrica
siendo ∂g
(u) unamatriz tal que la entrada de la i-ésima fila y j-ésima columna está
∂u
dada por ∂u∂g
(u) ∂gi
= ∂u j
(u). Dado que el vector δu es un vector columna, el
ij
producto del segundo término del miembro derecho está bien definido.
Operando se obtiene:
∂ε
(u) = be1 + be2 (u), (2.76)
∂u
por lo que la variación de la deformación está dada por:
Actividad
Demostrar la identidad de la Ecuación (2.77) de dos formas:
utilizando la definición y propiedades de la derivada,
realizando una aproximación de primer orden de la diferencia
ε(ue + δue ) − ε(ue ).
El miembro del trabajo virtual interno puede ser escrito en función del vector de
fuerzas internas definido como:
Ensamblado
que la barra 2 une los nodos 3 y 2 (en ambos casos se consideró una orientación
definiendo ejes locales). Ambas barras tienen sección transversal de área A0 y largo
`0 .
Sea un vector columna de desplazamientos virtuales de la estructura, formado
por los vectores columna de desplazamientos de los tres nodos: δuT = δuT1 , δuT2 , δuT3 .
El trabajo virtual de las fuerzas internas está dado por la suma de los trabajos vir-
tuales de cada barra:
1 2
T T
fint,1 T T
fint,1
δWint (u) = δu1 , δu2 1 + δu3 , δu2 2 (2.82)
fint,2 fint,2
donde fint,i
e
representa el vector columna de fuerzas internas correspondiente al nodo
i del elemento e. Desarrollando y factorizando se tiene
1
fint,1
δWint (u) = δu1 , δuT2 , δuT3 fint,2
T 1 2
(2.83)
+ fint,2 .
2
fint,1
siendo `e0 y Ae0 los largos y áreas de las secciones transversales del elemento e-ésimo.
En esta última ecuación se hizo explícita la dependencia de σ respecto a los despla-
zamientos y también se introdujo la notación abreviada
Iteración de Newton-Raphson
∂fint (k)
u ∆u(k) = − fint (u(k) ) − fext . (2.92)
∂u
Se cuenta entonces con las ecuaciones (2.89) y (2.92), que junto con la condición
de desplazamiento nulo para carga nula definen el procedimiento numérico a aplicar,
similar al descrito en la Ecuación (1.21).
Se deben calcular las derivadas del vector de fuerzas internas, lo cual puede ser
realizado para cada elemento dado que el ensamblado es un proceso de suma de
términos. Para calcular la derivada se utiliza la siguiente identidad para derivada de
producto de funciones:
∂fint (k)
(u )∆u(k) = KT (u(k) ) ∆u(k) , (2.97)
∂u
donde la matriz KT es llamada matriz tangente, la cual está dada por:
σ (k) A0
... + G, (2.98)
`
| 0{z }
Kσ
Criterios de parada
k∆u(k) k
< tolu
ku(k) k
kr(u(k+1) )k
< tolf
kfext k
En esta sección se presenta una implementación del método presentado para cal-
cular los desplazamientos de estructuras reticuladas planas. Como ejemplo se consi-
Sección 2.4. Método de Elementos Finitos Incremental 59
700
3000
650
2000 600
|uk (4)|
550
z
1000
500
0 450
Figura 2.13: Resultados de ejecución del Código 2.2: geometrías de referencia en trazo
discontinuo en azul y deformada en rojo continuo (izquierda), gráfico de valores de
desplazamiento vertical en cada iteración (derecha).
Tabla 2.1: Valores obtenidos por el método para el problema de Von Mises.
Actividad
Modificar los parámetros del Código 2.2 para reproducir los resul-
tados presentados en el Ejemplo 1 del artículo (Li and Khandelwal,
2017).
Apoyos elásticos
kx
fext,j
ky
Las fuerzas correspondientes a los resortes pueden ser consideradas como fuer-
zas externas dadas por:
Actividad
Mostrar que el sistema lineal asociado al método de Newton-
Raphson considerando apoyos elásticos está dado por:
Las cargas son también consideradas como aplicadas de forma estática, o con veloci-
dad de variación pequeña en comparación con la inercia y la cantidad de movimiento
de la estructura. Es decir que el tiempo identifica los incrementos pero no se consi-
dera la dinámica.
Método Incremental
u2,t+∆t
y u2,t
2
u1,t+∆t
u1,t
`0
x
Figura 2.15: Esquema de deformadas incrementales de elemento de barra.
para cada valor de fuerza considerado, sin embargo presenta limitaciones para ob-
tener soluciones aceptables para valores de carga de puntos límite.
1e+07
2000
8e+06
1000
Load factors
0
6e+06
y
-1000 4e+06
-2000
2e+06
-3000
0
0 2000 4000 6000 0 1000 2000 3000 4000 5000 6000
x Displacements
Figura 2.16: Gráficos obtenidos por el ONSAS al resolver una cercha de Von Mises.
valo de fuerzas. En la sección siguiente se estudiarán las causas del fenómeno del
salto. Este salto no está provocado por el pandeo local de las barras, fenómeno que
no fue considerado en este análisis.
Lo mencionado anteriormente representa una limitación considerable del méto-
do si se desea realizar un análisis para grandes desplazamientos de una estructura
con comportamiento de punto límite. Para poder obtener la curva completa de des-
plazamientos es necesario utilizar otras técnicas avanzadas que serán brevemente
descritas a continuación.
1e+07
2000
1000 5e+06
Load factors
0
0
y
-1000
-2000 -5e+06
-3000
-1e+07
0 2000 4000 6000 0 1000 2000 3000 4000 5000 6000
x Displacements
Figura 2.17: Ejemplo de cercha de Von Mises usando método de longitud de arco.
Se observa cómo este método permite obtener las fuerzas de equilibrio corres-
pondientes a cada configuración intermedia del movimiento de la estructura.
Es importante destacar que los gráficos carga-desplazamiento generados por la
versión actual del ONSAS no presentan los segmentos de puntos de equilibrio ines-
table punteados como se mostró al inicio del capítulo.
Los arcos han sido utilizados desde la antigüedad en estructuras de gran rele-
vancia como por ejemplo el Pont du Gard, estructura que sostiene un acueducto
construido por los romanos hace casi 2000 años y que sigue en pié al día de hoy
(Timoshenko, 1953). En la actualidad los arcos o estructuras con geometría de ar-
co suelen ser utilizadas en soluciones originales de estructuras de puentes o como
soporte de cubiertas de grandes luces. Frecuentemente se plantean soluciones con
esbelteces considerables volviendo necesario considerar deformaciones de segundo
orden.
En la Figura 2.18 se muestra un ejemplo de un arco reticulado utilizado como
soporte del techo de la estación de Retiro en Buenos Aires, Argentina. Se puede ob-
Sección 2.4. Método de Elementos Finitos Incremental 67
Figura 2.18: Estación de trenes de Retiro, Buenos Aires, Argentina (foto: JPZ).
250000
60
200000
50
Load factors
150000
y
40 100000
50000
30
-20 -10 0 10 20 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
x
Displacements
Figura 2.19: Gráficos obtenidos en ejemplo de Arco para carga de punto límite.
400000
300000
50
200000
Load factors
40 100000
y
0
30
-100000
-200000
-20 -10 0 10 20 0 2 4 6 8 10 12 14
x Displacements
Figura 2.20: Gráficos obtenidos en ejemplo de Arco para carga superior a punto lí-
mite.
estabilidad antes de alcanzar dicho factor de carga (esto será visto más adelante).
También se puede verificar que no se supera el 5 % de deformación axial lo que,
recordando la Sección 2.3, representa una verificación importante para el modelo y
la validez de los resultados usando la medida de deformación de Green y material
elástico lineal.
Es de interés observar los gráficos mostrados en la Figura 2.21, donde a la iz-
quierda se puede ver la configuración deformada para carga nula, es decir corres-
pondiente al punto de intersección de la curva carga-desplazamiento con la recta
λ = 0. La máxima magnitud de deformación axial es 4,71 %.
Esta configuración es claramente inestable ya que cualquier perturbación lleva-
ría a la estructura a adoptar otra configuración. En la siguiente sección se aborda
esta temática.
Se observa también que la deformación de la estructura es totalmente simétri-
ca, lo cual está asociado al hecho de que no existen imperfecciones presentes en la
geometría de la estructura. La inclusión de imperfecciones en el análisis es algo que
brindaría resultados más realistas y será realizado más adelante.
Sección 2.5. Estabilidad estructural 69
250000
60
200000
50
Load factors
150000
y
40
100000
30 50000
0
-20 -10 0 10 20 0 2 4 6 8 10
x Displacements
Sea una estructura cuyos desplazamientos nodales están dados por el vector u,
la función de energía potencial total está dada por:
T 1
δV ≈ (δu) (fint (u) − fext ) + (δu)T KT (u)δu, (2.109)
2
donde fint es el vector de fuerzas internas definido en la sección anterior. Se pue-
de observar que la condición de equilibrio establecida por el principio de trabajos
virtuales consiste en anular el término de primer orden de la variación de energía
potencial.
Para poder distinguir diferentes tipos de soluciones de equilibrio se pasa a ana-
lizar el término de segundo orden. Se establece que una configuración dada por u es
de equilibrio estable si se cumple el PTV y la energía potencial total aumenta para
cualquier perturbación virtual impuesta compatible con los vínculos.
Dado que el signo de la variación de energía potencial está dado por el signo de
la forma cuadrática del segundo término, es posible y conveniente establecer las de-
finiciones o clasificaciones de configuraciones de equilibrio en términos de la matriz
KT :
En los casos de inestabilidad, los vectores propios asociados a los valores propios
nulos de KT corresponden a los modos de pandeo de la estructura. Existen otras
definiciones o estudios más completos que pueden ser realizados utilizando términos
de mayor orden, los cuales no serán abordados en este material.
Es de interés entonces clasificar los equilibrios que se obtienen al resolver un
problema de forma incremental. Para eso es posible calcular los valores propios de
la matriz tangente asociada a los desplazamientos luego de obtenida la convergencia.
La herramienta ONSAS calcula, para cada incremento de carga, los valores propios
de la matriz KT mostrando en pantalla la cantidad de valores propios positivos y
negativos.
El cálculo de valores propios puede resultar costoso (en el caso de que no se
consideren técnicas de factorización particulares en la resolución del sistema lineal
(Bathe, 1982)) y en algunos casos poco eficiente ya que tendrá principal interés pa-
ra valores de carga cercanos a la inestabilidad, y no necesariamente para valores
previos.
Para estimar los valores de carga de inestabilidad, es posible realizar análisis
aproximados de los modos de pandeo. En las siguientes secciones se presentan y
realizan comentarios sobre dos de los métodos aproximados más frecuentemente
utilizados.
El análisis lineal de pandeo, también llamado LBA por su sigla en inglés (Linear
Buckling Analysis), es frecuentemente realizado en la práctica profesional para esti-
mar la carga crítica de pandeo de una estructura y observar sus modos de colapso
por inestabilidad (pandeo).
Se considera una estructura con un estado de fuerzas de referencia fext y un factor
de carga λ ∈ R. Partiendo de la expresión de la matriz tangente, asumiendo que la
estructura no se deforma de manera significativa previo a pandear, se desprecian los
términos que dependen directamente de u (matriz de desplazamiento inicial)
(2.111)
KT (u) ≈ KT 1 + λK̄σ ,
donde la matriz K̄σ puede ser fácilmente calculada obteniendo las tensiones de cada
elemento mediante un análisis lineal considerando la carga de referencia y realizan-
do el respectivo ensamblado de las matrices elementales.
El problema de LBA consiste en encontrar el factor de carga correspondiente a un
equilibrio inestable utilizando la matriz aproximada. Esto es equivalente a encontrar
72 Capítulo 2. No Linealidad Geométrica
los valores λ para los cuales la matriz KT se vuelve singular, es decir que pasa a tener
algún valor propio nulo. Encontrar valor propios nulos de KT consiste en resolver
el siguiente problema
KT (u)u = 0, (2.112)
donde sustituyendo la expresión aproximada se puede ver que esto es equivalente a
resolver el siguiente problema de valores propios generalizados:
KT 1 u = −λK̄σ u. (2.113)
Es necesario tener en cuenta que para realizar un análisis realista de una estruc-
tura se debe poder estimar las cargas vinculadas a los equilibrios inestables de forma
suficientemente precisa de forma de garantizar la seguridad de la estructura. En tra-
bajos como (Wriggers and Simo, 1990) se proponen métodos aproximados para el
cálculo de puntos de equilibrio inestable. Considerando literatura más reciente, se
recomienda seguir el abordaje presentado en (Bathe, 2014) para este tipo de cálculos.
La herramienta ONSAS realiza una estimación de la carga crítica de pandeo para ca-
da valor de factor de carga a través de dicho tipo de análisis aproximado. El método
utilizado también es presentado esquemáticamente en (Li and Khandelwal, 2017).
20000
Punto de bifuración
15000
Fuerza aplicada (kN)
10000
5000
Ramas inestables
0
0 50 100 150 200
Desplazamiento vertical (cm)
1.4e+07
1.5 1.2e+07
1e+07
Load factors
1 8e+06
y
6e+06
0.5 4e+06
2e+06
0 0
0 0.5 1 1.5 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4
x Displacements
Figura 2.23: Gráficos obtenidos para ejemplo de cercha de Von Mises con imperfec-
ción.
La carga máxima soportada por la estructura es 1,27 × 107 lo cual coincide con
el punto de bifurcación obtenido en el artículo citado.
Sección 2.5. Estabilidad estructural 75
Para finalizar la sección se presentan los elementos básicos del análisis del modo
de colapso del arco presentado en la Sección 2.4.5. Tal como se vió, al aumentar la
carga aplicada al arco se produce una deformación simétrica hasta alcanzar el punto
límite provocando el snap-through. Sin embargo en este problema existe un valor de
factor de carga menor al del punto límite para el cual se produce un pandeo según
un modo no simétrico. En (Timoshenko and Gere, 2009) se presenta el cálculo del
valor de carga distribuida crítica
EI π 2
qcr = 3 −1 , (2.114)
R α2
muestra la deformación posterior al colapso por pandeo así como también la curva
carga-desplazamiento.
120000
60
100000
50 80000
Load factors
60000
y
40
40000
30 20000
0
-20 -10 0 10 20 0 1 2 3 4 5 6
x Displacements
60
50
40
y
30
20
-20 0 20
x
Figura 2.25: Geometría de modelo simplificado de puente arco con tablero inferior.
Aspectos básicos
ut (x0 ) = xt − x0 . (3.2)
La dependencia de t está dada de forma implícita a través del subíndice. El útil calcu-
lar el gradiente de dicho campo denotado por D cuya expresión se obtiene derivando
ambos miembros respecto a x0 :
∂xt
F(x0 ) = (x0 ).
∂x0
Dado que se considera la deformada en un único instante, se omitió el tiempo t así
como también se podrá omitir la evaluación en el punto x de ahora en adelante.
Sean dos puntos x0 y x̂0 cercanos de la configuración de referencia y sus corres-
pondientes puntos deformados xt y x̂t . Utilizando la definición de χ y realizando un
desarrollo de primer orden respecto a x0 se obtiene la relación:
donde se consideró que la distancia entre los puntos x0 y x̂0 es pequeña. Esta iden-
tidad permite tener una relación entre los segmentos en las configuraciones de re-
ferencia y deformada. La magnitud escalar J(x0 ) = |F(x0 )| representa la razón de
variación volumétrica local en un entorno de la partícula que ocupa la posición x0
en la configuración de referencia. La función de deformación χ debe cumplir cier-
tas hipótesis de continuidad, como por ejemplo verificar que la razón de variación
volumétrica sea positiva durante todo el movimiento.
El estado local de deformaciones en un entorno de x0 está definido por lo tanto
por el tensor F, sin embargo existen diferentes definiciones de tensores que permiten
representar el estado de deformaciones local. Una herramienta útil es el tensor de
deformaciones de Green E el cual será definido de forma simplificada como una
extensión de la definición considerada para elementos de barra.
Se consideran los segmentos dx0 y dxt dados por:
Por otra parte los segmentos dxt y dx0 están relacionados a través de la Ecua-
ción (3.4) y las expresiones la Ecuación (3.3), por lo que se obtiene:
1 T T
(3.7)
εdx0 = dx0 F F − I dx0 .
2
Por lo que definiendo el tensor de deformaciones de Green como:
1 T
(3.8)
E= F F−I ,
2
se obtiene la siguiente expresión simplificada de la deformación unitaria:
Se puede ver que la matriz asociada de E es simétrica por lo que tendrá valo-
res propios reales. Estos valores propios representan las deformaciones unitarias de
Green correspondientes a los segmentos dados por las direcciones principales.
Sustituyendo la expresión de F en función del gradiente de desplazamiento dada
por la Ecuación (3.3) se obtiene
1
D + DT + DT D . (3.10)
E=
2
Se puede ver que tal como ocurre en el caso de barras, al considerar pequeñas de-
formaciones se obtienen equivalencias entre distintas medidas de deformación. En
el caso de la medida de Green, al considerar las componentes de primer orden se
obtiene:
1
D + DT = ε, (3.11)
E'
2
donde ε es el tensor de deformaciones infinitesimales asociado a la medida de de-
formación presentada en el capítulo anterior como deformación ingenieril.
Material Hiperelástico
λ
Ψ(E) = Tr(E)2 + µ Tr(E2 ), (3.16)
2
Sección 3.1. Elasticidad no lineal 81
donde λ y µ son dos parámetros reales del modelo que caracterizan el comporta-
miento del material. Esta función establece entonces una expresión para las tensio-
nes en todo punto a partir de la deformación, dada por:
Pequeñas deformaciones
∂ψ
σ= (ε). (3.20)
∂ε
No serán abordadas en este documento las condiciones que debe cumplir la función
ψ, sin embargo se considerará que deberán ser funciones que garanticen al menos
continuidad de la función σ(ε). Esta función determina la nolinealidad material para
comportamiento elástico en pequeñas deformaciones.
En el caso del material de Saint-Venant-Kirchhoff para pequeñas deformaciones
se obtiene una expresión de ψ dada por:
λ 2
(Tr(ε)) + µ Tr(ε2 ) (3.21)
ψ(ε) =
2
y la relación del tensor de tensiones está dada por:
∂ψ
σ= (ε) = λTr(ε) + 2µε, (3.22)
∂ε
lo que coincide con la ecuación constitutiva del material elástico-lineal.
∂ψ
σx = (εx ) = σ(εx ), (3.23)
εx
por simplicidad se considerará que la tensión está dada por la función σ(ε).
En el capítulo anterior se presentaron las ecuaciones del PTV, cuya aplicación
lleva al sistema de ecuaciones no lineales:
∂σ ∂ε
(u) = E tan (ε(u)) (u) (3.26)
∂u ∂u
donde E tan es el llamado módulo tangente, dado por:
∂σ
E tan (ε) = (ε). (3.27)
∂ε
ET ε si ε > 0 ET si ε > 0
∂σ
σ(ε) = (ε) = (3.29)
EC ε si ε < 0 ∂ε EC si ε < 0
Actividad
Modificar el Código 2.2 para resolver el problema de la cercha de
Von Mises considerando un comportamiento bi-modular para am-
bas barras donde la rigidez a compresión es la mitad de la rigidez a
tracción. Validar los resultados del código comparando con solucio-
nes analíticas para pequeños desplazamientos.
Formulaciones alternativas
Se considera una barra sometida a una historia de deformaciones tal que parte
de la deformación axial es remanente, es decir que se mantendrá luego de retirar las
cargas externas. Esta deformación remanente es producida cuando la magnitud de
la tensión alcanza un cierto valor de tensión de fluencia σY .
ε = εe + εp . (3.30)
Por otra parte se considera la ley elástica que establece que la relación entre la
tensión y la deformación elástica es lineal y dada por:
Función de fluencia
Para encontrar las condiciones que deben cumplir las tensiones solución del pro-
blema se plantean las condiciones de optimalidad de Karush-Kuhn-Tucker (KKT). A
continuación se presentan de forma esquemática las ecuaciones de KKT (ver (Luen-
berger and Ye, 2008) por detalles del desarrollo).
Sección 3.2. Introducción a la Teoría de Plasticidad 87
todo punto x que satisface las condiciones de optimalidad de primer orden (”punto
crítico”) debe necesariamente cumplir
∂Φ
ε̇p = η (σ, σY ) (3.40)
∂σ
Φ(σ, σY ) 6 0 (3.41)
η > 0 (3.42)
η Φ(σ, σY ) = 0. (3.43)
Estas condiciones permiten derivar las siguientes leyes sobre el proceso de plastifi-
cación.
La condición dada por la Ecuación (3.40) es llamada ley de flujo plástico. Calcu-
lando la derivada e introduciendo la variable γ dada por γ̇ = η se tiene:
γ̇ > 0 (3.45)
Φ(σ, σY ) 6 0 (3.46)
γ̇Φ(σ, σY ) = 0 (3.47)
Función de endurecimiento
σY = σY (ε¯p ), (3.48)
Esta nueva magnitud representa la memoria del material, siendo este el elemento
conceptual más relavante que define al material como no elástico.
Derivando respecto al tiempo ambos miembros de la Ecuación (3.49) y utilizando
la Ecuación (3.44) se obtiene:
Φ=0 y Φ̇ = 0. (3.52)
Sección 3.2. Introducción a la Teoría de Plasticidad 89
∂σY ¯˙p
Φ̇ = signo(σ)σ̇ − ε , (3.53)
∂ ε¯p
por lo tanto, igualando a cero se tiene que durante la plastificación se cumple:
∂σY ¯˙p
signo(σ)σ̇ = ε . (3.54)
∂ ε¯p
Por otra parte, derivando la ley elástica dada por la Ecuación (3.31) se tiene la
relación:
σ̇ = E ε̇ − ε˙p . (3.56)
∂σY p
E (ε¯ )
E tan
(ε¯p ) = ∂ ε¯p . (3.58)
∂σY p
E+ ¯
(ε )
∂ ε¯p
A partir de esta relación se puede también obtener una expresión para la pendien-
te de la función de endurecimiento en función de E tan . Es importante resaltar que
existe una relación implícita entre la deformación total y la deformación plástica
acumulada cuyo análisis no debe ser despreciado.
Relaciones algebraicas
Para obtener un método numérico para resolver las ecuaciones que gobiernan el
proceso de plastificación, es necesario utilizar alguna regla para integrar las ecua-
ciones diferenciales del modelo.
En esta sección se presenta el desarrollo de un método numérico para resolver
las ecuaciones correspondientes a una barra sometida a una deformación axial, asu-
miendo que en cada paso se conoce la deformación total aplicada. Se considera que
el rango temporal del proceso estudiado es discretizado en intervalos uniformes. Por
simplicidad y sin pérdida de generalidad, se considera ∆t = 1. Nuevamente, el tiem-
po no representa que las cargas sean aplicadas dinámicamente, sino simplemente
una referencia a factores de carga o deformación impuesta.
Sea un instante de tiempo n para el cual se conocen los valores de σn , εn , ε¯p n y
γn . Se aplica un incremento en la deformación total por lo que se conoce el valor de
deformación en el instante de tiempo siguiente εn+1 = εn + ∆ε. Se desea determi-
nar las otras magnitudes. Se comienza operando con la derivada de la ley elástica,
dada por la Ecuación (3.56), donde se sustituye la Ley de flujo plástico dada por la
Ecuación (3.44). Se obtiene por lo tanto:
Se puede considerar una regla de integración de tipo Euler hacia atras o Backward
Euler, donde se integra en el tiempo numéricamente, evaluando el integrando en el
instante n + 1, obteniendo:
Actividad
Mostrar que se verifica la siguiente igualdad
Para poder obtener una expresión del valor absoluto de la tensión se continúa
desarrollando la identidad de la Ecuación (3.62),
Dado que ∆γ > 0, los resultados de la función signo de cada miembro son mul-
tiplicados por valores no negativos por lo tanto deben tener el mismo resultado,
obteniendo:
Método numérico
σ
σE,n+1
∆γE E
σn+1
σn
E
εn εn+1 ε
donde σY,0 y K son parámetros del modelo, que determinan la tensión de fluencia
inicial y la pendiente de la función, respectivamente.
Sustituyendo la función dada por este modela en la Ecuación (3.71) se obtiene:
Implementación en GNU-Octave
3e+08 3e+08
2e+08 2e+08
1e+08
1e+08
0
Tension
Tension
0
-1e+08
-1e+08
-2e+08
-2e+08 -3e+08
-3e+08 -4e+08
-0.002 -0.001 0 0.001 0.002 -0.002 -0.001 0 0.001 0.002
Deformacion Deformacion
As rm
As
ε(y) = −κ (y − yN ) (3.78)
si ε ∈ [−3,5‰, −2‰]
−fcd
1000 . fcd . ε . (250 . ε + 1) si ε ∈ [−2‰, 0]
σ C (ε) = (3.79)
si ε > 0
0
mientras que para el acero se tiene que la tensión de límite elástico es εy = 2,17‰
y por lo tanto la ecuación está dada por
Actividad
Comparar el diagrama de interacción obtenido con el correspon-
diente presentado en (Jimenez Montoya et al., 2009).
Sección 3.3. Solicitaciones últimas en sección de hormigón armado 99
σD
−σ C
fyd
−0,85 . fcd
εy ε
−2‰ −3,5‰ −ε
M
µ=
bh2 fcd
N
ν=
bhfcd
Afyd
ω=
bhfcd
Las cargas dinámicas externas están dadas por el vector fext,t y pueden ser perió-
dicas (harmónicas, no harmónicas) o aperiódicas (impulsivas, transitorias). El tipo
de carga influye en la definición del modelo estructural y a su vez en el método de
resolución de las ecuaciones de movimiento resultantes.
El vector de fuerzas internas fint incorpora los efectos de no linealidad geométrica
y eventualmente otros tipos de comportamiento no lineal de los elementos de barra,
como puede ser la no linealidad material.
La matriz M es llamada Matriz de Masa. Si ésta se deduce usando las funciones de
Sección 4.2. Dinámica Lineal 103
Dado lo anterior, el vector de fuerzas internas fint (ut ) está dado por Kut , por
lo que la ecuación de movimiento para el análisis dinámico lineal de la estructura
resulta:
Müt + Cu̇t + Kut = fext,t (4.6)
La Ecuación (4.6) representa un sistema de ecuaciones diferenciales ordinarias
de segundo grado, con coeficientes constantes y no homogéneas. Para determinar
una solución, se deben dar condiciones iniciales en el instante t0 de desplazamiento
(ut0 ) y velocidad (u̇t0 ).
Las definiciones de las matrices M y C, junto con comentarios generales sobre
las mismas, ya fueron presentadas en la sección anterior y mantienen su validez.
Formulación y Pseudo-código
por esto que al aplicar este método se utiliza la matriz de Masa Concentrada y no la
matriz de Masa Consistente. Adicionalmente, la matriz de efectos viscosos puede ser
despreciada en fenómenos rápidos de muy corta duración. En su defecto, se utiliza
una matriz C diagonal de manera de no aumentar el costo computacional.
Estabilidad Numérica
En el Algoritmo 4 se indica que el valor del paso temporal elegido debe satisfacer
la condición: ∆t < Tmin /π, donde Tmin es el mínimo período de vibración natural
del modelo de elementos finitos. Para entender esta restricción, se realiza el análisis
de la estabilidad numérica del Método de Diferencia Centrada.
El análisis de estabilidad numérica se lleva a cabo en una estructura sin amor-
tiguamiento, es decir con C = 0. Sin embargo, es posible analizar la estabilidad
tomando en cuenta el amortiguamiento, en particular bajo la hipótesis de amorti-
guamiento proporcional. Se realizará algún comentario adicional a este punto más
adelante en esta sección.
Previo a realizar el análisis de estabilidad numérica, se debe reconocer que es
suficiente estudiar la estabilidad de la solución numérica de un sistema con un grado
de libertad. La justificación de este hecho se basa en la descomposición modal de la
respuesta de la estructura y se describe brevemente a continuación.
Se definen los modos normales de vibración de la estructura como aquellas so-
luciones de vibración libre no forzada en las que todos los puntos de la estructura
se desplazan de forma sinusoidal respecto del tiempo, en fase y a una misma fre-
cuencia; es decir: ut = sin(ω(t − t0 ))φ, donde φ es un vector de desplazamientos
Sección 4.2. Dinámica Lineal 107
ω 2 Mφ = Kφ. (4.11)
La ecuación anterior define un problema de valores propios para una matriz si-
métrica definida positiva. Esto implica que existen tantos modos normales {φi , ωi2 }
como grados de libertad tenga la estructura (n), es decir tantos como la dimensión
de las matrices M o K (reducidas).
Actividad
Mostrar que los modos normales de vibración satisfacen la relación:
Agrupando cada uno de los n modos normales φi como columnas de una matriz
Φ, y sus respectivas frecuencias al cuadrado (ωi2 ) en la diagonal de una matriz Λ, se
tiene:
MΦΛ = KΦ (4.12)
La Ecuación (4.12) define la relación matricial que verifican los modos normales
de la estructura. Queda definida así una base de vectores dada por Φ en la cual se
pueden escribir los vectores de desplazamiento a lo largo del tiempo.
A partir del desarrollo anterior, se define ut = Φxt . Es decir, xt son las coorde-
nadas que representan a ut en la base de modos normales Φ. La ecuación de movi-
miento resulta:
MΦẍt + KΦxt = fext,t . (4.13)
Si se multiplica la Ecuación (4.13) a la izquierda por ΦT se obtiene,
2 − ωi2 ∆t2
−1
A= . (4.21)
1 0
El vector x̂0 es un vector constante, no nulo, que queda definido a partir de las
condiciones iniciales del Problema Test. Se puede suponer sin perder generalidad
que t = n∆t con n ∈ N. Por lo tanto, la aplicación repetida de la Ecuación (4.20)
permite obtener:
x̂t = An x̂0 (4.22)
El miembro derecho de la Ecuación (4.23) está acotado si, se cumple que el radio
espectral de A es menor a 1, es decir ρ(A) < 1. El radio espectral de A se define
como el máximo de los módulos de los valores propios de la matriz A, es decir:
ρ(A) = máxi |λi |. Si la multiplicidad de los valores propios es uno, la condición se
puede relajar a ρ(A) 6 1.
La justificación formal de este punto se puede buscar en las notas del curso Mé-
todos Numéricos (Facultad de Ingeniería, Universidad de la República). Alternativa-
mente, un idea intuitiva de la validez de este resultado puede obtenerse pensando
en que si A es diagonalizable (A = P−1 DP), entonces An = P−1 Dn P. Si los
módulos de todos los valores propios que hay en la diagonal de D son menores a 1,
se tiene que Dn → 0.
Por lo tanto, para evaluar la condición de estabilidad del método en cuestión
debemos calcular los valores propios de A, verificar que tienen multiplicidad igual
a uno, e imponer que el módulo de éstos sea menor o igual a 1.
Se define z = ω∆t > 0 con lo cual los valores propios resultan:
r
2 − z2 (2 − z 2 )2
λ1,2 = ± − 1. (4.24)
2 4
Se puede ver que los valores propios tienen multiplicidad 1, incluso cuando z = 2.
En la Figura 4.1 se presenta la gráfica de la función ρ(A)(z) = máx{|λ1 (z)|, |λ2 (z)|}
para valores z ∈ [0, 3].
4
ρ
0
0 0.5 1 1.5 2 2.5 3
ωi∆t
En una estructura completa, todos los grados de libertad deben satisfacer esta
condición. Se observa que la condición más estricta para ∆t corresponde al grado
de libertad con el período natural más corto (es decir, el de mayor frecuencia), con
lo cual el criterio de estabilidad numérica para una estructura completa es:
Tmin
∆t 6 , (4.25)
π
siendo Tmin el menor período natural.
u1
0.01
500 u2
u3
400
0.005
desplazamiento [m]
presion [kPa]
300
0
200
100 -0.005
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0 0.1 0.2 0.3 0.4 0.5
t [s] t [s]
Actividad
Obtener los valores de la fuerza cortante en un pilar en función del
tiempo.
11
12 a = 0.7; %m Lado de Columna
13 E = 30e9; %N/m2 Mod. Young Hormigon
14 I = a^4/12; %m4
15 H = 3.5; %m Altura Entre Pisos
16
17 % Definicion de Matriz de Rigidez
18 k = ncols * 12*E*I/H^3; %N/m
19 %u1 u2 u3
20 K = k*[2 −1 0 ; −1 2 −1 ; 0 −1 1]; %N/m
21
22 % Definicion de Matriz de Masa
23 Mpiso = 2500*(e*A^2+hv*bv*A*8); %kg
24 Mpilar = 2500*ncols*a^2*H; %kg
25 Mmuro = muro*H*A*4*0.6; %kg
26 Mint = Mpiso + Mpilar + Mmuro;
27 Msup = Mpiso + Mpilar/2 + Mmuro/2;
28
29 M = diag([Mint, Mint, Msup]); %kg
30
31 % Definicion de Matriz Amortiguamiento
32 % Usamos Rayleigh Damping: C = eta*M + delta*K
33 % Amortiguamiento: 3 % critico para 25rad/s y 106rad/s
34 C = 1.21*M + 4.6e−4*K;
35
36 % Calculo de Modos Normales
37 [PHI, w2] = eig(K,M);
38
39 w = sqrt(diag(w2)); %rad/s
40 f = w/2/pi; %Hz
41 T = 1./f %sec
42
43 % Calculo de Paso Critico para Diferencia Centrada
44 dtcr = min(T)/pi
45
46 % Definicion de Historia de Presiones y Fuerzas Laterales
47 Aint = A*H; %m2
48 Asup = A*H/2; %m2
49
50 % funcion analitica para evaluar presion triangular en el tiempo
51 press = @(t,ta,te,pr1,pr2) (pr1−(pr1−pr2)/te*(t−ta)).*and(t>ta,t<ta+te);
52
53 % historia de presion en el tiempo como suma de tramos lineales (triangulos)
54 pt = @(t) press(t,15.6e−3,6.2e−3,0.55e6,0)+ ...
55 press(t,66.0e−3,11.9e−3,0,−0.035e6)+...
56 press(t,77.9e−3,19.1e−3,−0.035e6,0); %N/m2
57
58 ft = @(t) pt(t)*[Aint;Aint;Asup]; %N
59
60 % Definicion de Condiciones Iniciales y Tiempo Final
61 t0 = 0; %sec
62 u0 = [0;0;0]; %m
63 v0 = [0;0;0]; %m
64 ac0 = M\(ft(t0)−C*v0−K*u0); % de ec de movimiento Mu.. + Ku = ft
114 Capítulo 4. Introducción al Análisis Dinámico
65
66 tf = 2*max(T); % 2 periodos del modo fundamental
67
68 % Inicializacion Difrerencia Centrada
69 dt = dtcr/20;
70 a0 = 1/dt^2; a1=1/2/dt; a2=2*a0; a3=1/a2;
71
72 u(:,1) = u0 − dt*v0 + a3*ac0; % u(−dt)
73 u(:,2) = u0; % u(0)
74
75 Meff = a0*M+a1*C; Keff = (K−a2*M); M2 = a0*M−a1*C;
76
77 % Comienza Marcha en el Tiempo usando Diferencia Centrada
78 t(1) = −dt; t(2) = t0; k=2;
79 while t<tf
80 feff = ft(t(k)) − Keff*u(:,k) − M2*u(:,k−1);
81 u(:,k+1) = Meff\feff;
82 acc(:,k) = a0*(u(:,k+1)−2*u(:,k)+u(:,k−1));
83 vel(:,k) = a1*(u(:,k+1)−u(:,k−1));
84 k=k+1;
85 t(k)=t(k−1)+dt;
86 end
87
88 subplot(1,2,1), plot(t,pt(t)/1e3,'−b')
89 xlabel('t [sec]'), ylabel('presion [kPa]')
90 axis([t0,0.3*max(t),1.2*min(min(pt(t)/1e3)),1.2*max(max(pt(t)/1e3))])
91 title('Presion por Explosion')
92
93 subplot(1,2,2), plot(t,u(1,:),'−b',t,u(2,:),'−r',t,u(3,:),'−k')
94 legend('u1(t)','u2(t)','u3(t)'), xlabel('t [sec]'), ylabel('desplazamiento [m]')
95 axis([t0,1.1*max(t),1.5*min(min(u)),1.5*max(max(u))])
96 title('Respuesta de la Estructura')
Formulación y Pseudo-código
cidades se considera:
1 δ
K̂ = K + 2
M+ C, (4.30)
α∆t α∆t
Las Ecuaciones (4.29) a (4.33) son suficientes para poder llevar adelante el paso
integración en el tiempo de Newmark. En el Algoritmo 5 se presenta un Pseudo-
Código del Método de Newmark.
116 Capítulo 4. Introducción al Análisis Dinámico
Estabilidad Numérica
de la siguiente forma:
x̂t+∆t = Ax̂t . (4.35)
La matriz A, con β = (1/(ω ∆t ) + α) , tiene la siguiente forma:
2 2 −1
−β/∆t2
−(1/2 − α)β −β/∆t
A = ∆t[1 − δ − (1/2 − α)δβ] 1 − βδ −βδ/∆t . (4.36)
2
∆t [1/2 − α − (1/2 − α)αβ] ∆t(1 − αβ) 1 − αβ
4
Trapecio: α = 1/4, β = 1/2
Acc.Lineal: α = 1/6, β = 1/2
3.5
2.5
ρ(A)
1.5
0.5
0 0.5 1 1.5 2 2.5 3 3.5
∆t/T
Figura 4.5: Estabilidad de variantes del Método de Newmark (Trapecio y Acel. Li-
neal).
donde ẍG,t tiene el registro de aceleraciones de suelo medidos según una dirección
dada.
En el presente ejemplo, los datos de aceleración corresponden al sismo de Loma
Prieta, ocurrido el 17 de octubre de 1987 en California, U.S.A.
En la Figura 4.6 se muestran los resultados de la solución de la dinámica del
edificio para la aceleración de terreno elegida. Comparando los resultados obtenidos
con los de la solución de la dinámica del edificio sometido a una detonación, se
observa mayores desplazamientos en el caso del sismo.
0.01
0.4
0.005
0.2
u1 (t) [m]
accel [g]
0 0
-0.2 -0.005
-0.4 -0.01
0 10 20 30 40 0 10 20 30 40
t [s] t [s]
0.02
0.02
0.01 0.01
u2 (t) [m]
0 u3 (t) [m] 0
-0.01 -0.01
-0.02
-0.02
0 10 20 30 40 0 10 20 30 40
t [s] t [s]
Figura 4.6: Solución numérica del movimiento del edificio bajo acción de sismo.
La Ecuación (4.39) provee la regla mediante la cual se calcula ut+∆t . Notar que
en el caso de una estructura no lineal debemos en cada paso evaluar el vector de
fuerzas internas en el instante de tiempo actual. Es claro que el método continúa
siendo explícito.
Se destaca que sigue siendo necesario, tal como en la solución de problemas
lineales, que la matriz de masa y eventualmente la de amortiguamiento sean dia-
gonales de manera que la evaluación de ut+∆t a partir de la Ecuación (4.39) sea
computacionalmente económica.
A partir del desarrollo anterior, se puede actualizar el procedimiento de solución
del Método de Diferencia Centrada para el caso de estructuras no lineales tal como
se muestra en el Algoritmo 6.
En este momento es una buena idea comentar sobre la elección del paso temporal
(∆t) para el Método de Diferencia Centrada en problemas no lineales. Tal como se
mencionó anteriormente, la estabilidad numérica del método continúa restringiendo
el paso temporal máximo que se puede usar. La elección del paso se vuelve más difícil
en el caso de estructuras no lineales ya que el paso temporal crítico no es constante
durante la solución. La rigidez de la estructura puede variar a lo largo del tiempo
por efectos de rigidez geométrica o por no linealidad material. Si se espera que la
estructura solamente se flexibilice a medida que transcurre el tiempo, un criterio
razonable puede ser determinar el paso temporal en base a la rigidez inicial de la
estructura. Si este no es el caso, se debe prever el efecto de la rigidización de la
estructura en el paso temporal crítico y asegurar que el paso elegido no superará el
valor crítico mínimo previsto.
Al realizar análisis dinámicos de estructuras no lineales usando el Método de Di-
ferencia centrada, la elección de un paso temporal demasiado largo en comparación
con el mínimo paso temporal crítico previsto, puede llevar a acumular errores signi-
ficativos durante aquellos instantes de tiempo en los cuales el paso temporal elegido
superó al paso crítico. Esta acumulación acotada de errores es marcadamente dis-
tinta a la que se observa en análisis dinámicos lineales, en los cuales si se elije un
Sección 4.3. Dinámica No Lineal 121
paso temporal por encima del valor crítico, la solución diverge y el analista puede
reconocer fácilmente que seleccionó un paso temporal erróneo.
Se considera una cercha de tipo de Von Mises bajo la acción de una carga diná-
mica. La estructura consiste en dos columnas flexibles empotradas en la base y dos
bielas articuladas que trabajan como puntales de la cercha tal como se muestra en la
Figura 4.7.
Se coloca una masa m suspendida de la cumbre de la cercha. Todas las barras de
la cercha están formada por acero (E = 200 GPa) y sección rectangular (a = 3,2mm
y b = 25,4mm). Los parámetros que definen la geometría son: Lc = 240 mm, Lx =
187 mm y Lz = 84 mm. Esta estructura es la utilizada por el profesor A. Wadee en
la presentación referida al inicio del Capítulo 2 de este documento.
La forma de cargar la estructura consiste en suspender la masa en la configura-
ción indeformada y luego dejarla libre. Esto implica que el vector de deformación
inicial es nulo y también lo es el vector de velocidades inicial.
Dada la simetría axial, respecto al eje vertical, presente en la estructura y las
cargas aplicadas, se modela la mitad izquierda de la estructura. De esta forma solo
se debe analizar una biela de Green con condición de borde deslizante vertical en
el extremo derecho y deslizante horizontal con resorte en el extremo izquierdo tal
como se muestra en la Figura 4.8. La constante elástica del resorte está dado por la
rigidez flexional de los pilares.
El vector de desplazamientos es considerado ut = (u1,t , u2,t )T y por lo tanto la
122 Capítulo 4. Introducción al Análisis Dinámico
u2
Lz
u1 m
Lc
Lx
u2
u1
donde el segundo sumando del trabajo virtual interno corresponde al trabajo virtual
interno de la fuerza realizada por el resorte horizontal de constante kc .
lt2 − l02
εG = . (4.42)
2l02
con lo cual, el vector de fuerzas internas para la barra de Green y el resorte lineal
resulta:
EAl0 (u21,t + u22,t − 2Lx u1,t + 2Lz u2,t )
−Lx + u1,t
fint (ut ) =
2l04 Lz + u2,t
(4.48)
kc u1,t
+ .
0
Actividad
Discutir y estimar cuánto vale ∆Tcr para el ejemplo de la cercha de
Von Mises.
14 mb = l0*Ac*rho; %kg
15 m = 1.4; %kg Pandeo incipiente en 1.4
16 c = 2; %kg/s (amortiguamiento por friccion juntas y arrastre pesa)
17 g = 9.81; %m/s2
18
19 %% Defino Vector de fuerzas Internas: fint(u) − u = [u1 , u2]^T
20 Fint = @(u) E*Ac*l0*(u(1)^2+u(2)^2−2*Lx*u(1)+2*Lz*u(2))/2/l0^4*[−Lx+u(1);Lz+u(2)
]+[kc*u(1);0];
21
22 %% Defino Vector de fuerzas Externas: gravedad
23 ft = @(t) [0;−(m+mb)/2*g]; %N
24
25 %% Defino Matriz de Masa Concentrada
26 M = [mb 0 ; 0 (mb+m)/2];
27
28 %% Defino Matriz de Amortiguamiento
29 C = [c/10 0 ; 0 c];
30
31 %% Defino Condiciones Iniciales
32 t0 = 0;
33 u0 = [0;0];
34 v0 = [0;0];
35 ac0 = M\(ft(t0)−C*v0−Fint(u0)); % de ec de movimiento Mu.. + Fint(u) = ft
36
37 % Inicializacion Difrerencia Centrada
38 tf = 2.0; dt = .000025; % sec
39
40 a0 = 1/dt^2; a1=1/2/dt; a2=2*a0; a3=1/a2;
41
42 u(:,1) = u0 − dt*v0 + a3*ac0; % u(−dt)
43 u(:,2) = u0; % u(0)
44
45 Meff = a0*M+a1*C;
46 M2 = a0*M−a1*C;
47
48 % Comienza Marcha en el Tiempo usando Diferencia Centrada
49 t(1) = t0−dt; t(2) = t0; k=2;
50
51 epsg(k) = (u(1,k)^2+2*Lz*u(2,k)−2*Lx*u(1,k)+u(2,k)^2)/2/l0^2;
52
53 while t<tf
54 feff = ft(t(k)) −Fint(u(:,k)) +a2*M*u(:,k) − M2*u(:,k−1);
55 u(:,k+1) = Meff\feff;
56 acc(:,k) = a0*(u(:,k+1)−2*u(:,k)+u(:,k−1));
57 vel(:,k) = a1*(u(:,k+1)−u(:,k−1));
58 epsg(k+1) = (u(1,k+1)^2+2*Lz*u(2,k+1)−2*Lx*u(1,k+1)+u(2,k+1)^2)/2/l0^2;
59 k=k+1; t(k)=t(k−1)+dt;
60 end
61
62 subplot(3,1,1)
63 plot(t(1:10:end),1000*u(1,1:10:end))
64 xlabel('t [s]'); ylabel('u_1 [mm]');
65 axis([0 2 1e3*min(u(1,:))*1.1 1e3*max(u(1,:))*1.1]);
66 subplot(3,1,2)
126 Capítulo 4. Introducción al Análisis Dinámico
67 plot(t(1:10:end),1000*u(2,1:10:end))
68 xlabel('t [s]'); ylabel('u_2 [mm]');
69 axis([0 2 1e3*min(u(2,:))*1.1 1e3*max(u(2,:))*1.1]);
70 subplot(3,1,3)
71 plot(t(1:10:end),E*Ac*epsg(1:10:end))
72 xlabel('t [s]'); ylabel('Directa [N]')
73 axis([0 2 E*Ac*min(epsg)*1.1 E*Ac*max(epsg)*1.1]);
30
20
u1 [mm]
10
0
-10
0 0.5 1 1.5 2
t [s]
0
-50
u2 [mm]
-100
-150
-200
0 0.5 1 1.5 2
t [s]
100
Directa [N]
50
0
-50
0 0.5 1 1.5 2
t [s]
Actividad
Compare este valor con el valor de carga crítica correspondiente a
carga cuasi-estática.
Sección 4.3. Dinámica No Lineal 127
Actividad
Formular la solución de la ecuación no lineal definida en el paso de
Newmark aplicando soluciones iterativas de tipo Newton-Raphson
o incluso Newton-Raphson modificado.
presentan los resultados numéricos obtenidos al resolver dos ejemplos con solución
analítica conocida de inestabilidad en pórticos.
Módulo A Módulo B
Aco Aco
Adi Adi
h Ave Ave h Ave Ave
Aco Aco
` `
parámetros Aco , Ave y Adi representan las áreas de las secciones transversales de
las barras horizontales, verticales y diagonales respectivamente. Los parámetros h
y ` representan los distancias entre nodos en vertical y horizontal respectivamente,
donde la repetición de los módulos se produce en la dirección horizontal.
A continuación se realiza el desarrollo de las ecuaciones que permiten obtener
parámetros de los reticulados a partir de los parámetros de vigas para cada módulo.
En ambos casos se considerará que Aviga , Iviga y Aviga,s representan el área de
sección transversal, el segundo momento de inercia de sección transversal y el área
equivalente de cortante de la viga cuyo comportamiento se desea representar.
Se buscará representar el comportamiento de vigas bajo la hipótesis de Euler-
Bernoulli, es decir que no se producen deformaciones por cortante. Para esto, en
los modelos de reticulados, se buscará que la rigidez de deformación a cortante sea
elevada de manera que no hayan deformaciones de cortante significativas. El módulo
de Young de la viga E será idéntico al valor considerado para las barras del reticulado
equivalente.
A.1.1. Módulo A
obtiene:
El valor β podría estar dado por las propiedades del elemento de viga cuyo com-
portamiento se desea modelar. A pesar de esto y con el objetivo de simplificar el
desarrollo y obtener una expresión que permita su aplicación directa, se establece
como criterio que β debe ser elevado. Con esto se buscará reproducir el comporta-
miento de elementos de viga según la teoría de viga de Euler-Bernoulli en la que la
energía de deformación por cortante es despreciada.
Se debe determinar las áreas Ave y Adi con una condición dada por la Ecua-
ción (A.5) y la restricción de que ambos valores deben ser positivos:
Dado que el problema no tiene solución única se opta por buscar soluciones que
cumplan Ave = Adi . Sustituyendo se tiene:
`3d + h3 1 `2
2
Aco = − 2 Ave > 0 (A.8)
Ave `h 2β 3h
132 Apéndice A. Modelos equivalentes reticulados-vigas
lo que es equivalente a
3h2
Ave > 0 ⇐⇒ β < . (A.10)
2`2
Esto establece una cota superior para la máxima rigidez a deformación por corte que
el reticulado (con las hipótesis establecidas hasta el momento) es capaz de represen-
tar. Si se considera que el largo de la celda ` está dado por L/nc donde L es el largo
de la viga y nc el número de celdas. Se puede ver que si nc tiene a infinito se puede
reproducir el comportamiento de vigas de gran rigidez a deformación por cortante.
De esta forma el número de celdas será una herramienta de control de esta rigidez.
Entre los posibles valores de β se opta como criterio por el valor β = h2 /`2 .
Sustituyendo el valor de β en la ecuación anterior se tiene:
`3d + h3
Ave = 6Aco (A.11)
`3
Aviga `3d + h3
Aco = Adi = Ave = 6Aco h = 2ρviga (A.12)
2 `3
A.1.2. Módulo B
2EAve Adi `3
EAviga = 2EAco + (A.13)
(`3d Ave + 2Adi h3 )
1 `3d `2
= − (A.14)
GAviga,s 2EAdi `h2 12EIviga
Aco h2
EIviga = E (A.15)
2
nes se obtiene:
3
6Aviga Iviga β h2 + l2 2
Adi = (A.16)
h2 l (Aviga βl2 + 12Iviga )
2Iviga
Aco = (A.17)
h2
12Aviga Iviga βh Aviga h2 − 4Iviga
Ave = − (A.18)
l A2viga βh2 l2 − 16Aviga Iviga βl2 + 12Aviga Iviga h2 − 48Iviga
2
h2
Usando nuevamente el criterio β = `2 se tiene
3
6Aviga Iviga h2 + l2 2
Adi = (A.19)
l3 (Aviga h2 + 12Iviga )
2Iviga
Aco = (A.20)
h2
12Aviga Iviga h3 Aviga h2 − 4Iviga
Ave = − 3 2 2
(A.21)
l Aviga h4 − 4Aviga Iviga h2 − 48Iviga
(A.22)
Se debe garantizar nuevamente que todos los valores de las áreas sean positivos
para lo cual se debe estudiar Ave . Estudiando las raíces de los polinomios del nume-
rador y denominador se obtiene que un rango de valores de h que garantizan que
Ave sea positivo es:
q √
h ∈ (2ρviga , 2(1 + 13)ρviga ) ⊃ (2ρviga , 3ρviga ) (A.23)
π 2 EIviga
Pcr = ≈ 142 kN. (A.24)
L2viga
El factor de carga crítica del análisis lineal de pandeo del modelo del reticulado es
0,9968. En la Figura A.2 se muestran los gráficos obtenidos al utilizar el módulo A. En
1.4 1
1.2
1 0
Load factors
0.8
y
-1
0.6
0.4
-2
0.2
0 0 1 2 3 4 5
0 0.5 1 1.5 2
x
Displacements
Figura A.2: Gráficos obtenidos para ejemplo Pinned Buckling con módulo tipo A.
clase se analizarán los resultados obtenidos al variar los parámetros del reticulado,
particularmente el número de celdas de la discretización.
Para el módulo B se tiene un valor β ≈ 0,26. El factor de pandeo del LBA es
0,99708. En la Figura A.3 se muestran los gráficos obtenidos al utilizar el módulo
B. Dado que el módulo B es simétrico y que la carga es aplicada de forma simétrica
respecto al eje horizontal, la dirección en la cual se produce la deformación post-
pandeo estará dada por la imperfección. La imperfección está dada por el sentido
del desplazamiento asociado al vector propio del problema LBA. Por ser un vector
propio se cumple que el opuesto también lo es por lo que la imperfección puede ser
considerada en una u otra dirección en función de con qué criterio la función eig
devuelva los vectores.
Finalmente se desea obtener alguna verificación analítica de las curvas obtenidas
Sección A.2. Ejemplos de validación 135
1.4
1.2
2
1
Load factors
0.8
1
y
0.6
0.4
0
0.2
0 -1
0 0.5 1 1.5 2 0 1 2 3 4 5
Displacements x
Figura A.3: Gráficos obtenidos en ejemplo Pinned Buckling con módulo tipo B.
para las ramas post-críticas. Realizando el análisis con ambos módulos y comparan-
do con la curva dada por la solución analítica de segundo orden:
P (u) 1 π 2 u2
'1+ , (A.25)
Pcr 2 4 L2
se obtiene la gráfica mostrada en la Figura A.4.
1.3
1.2
Load factors
1.1
Módulo A
Módulo B
Aprox. orden 2
0.9
0 0.5 1 1.5
Displacements
Figura A.4: Gráficos de ramas post-críticas obtenidas usando los módulos reticulados
descritos y la expresión analítica.
π 2 EIviga
Pcr = ≈ 35,4 kN, (A.26)
4L2viga
1.4
1.2
3
1
Load factors
0.8
2
y
0.6
0.4 1
0.2
0 0
0 1 2 3 4 0 1 2 3 4 5
Displacements x
Figura A.5: Gráficos obtenidos en ejemplo Cantilever Buckling con módulo tipo A.
ONSAS: es el archivo principal del código, el cual debe ser ejecutado desde
octave para iniciar la ejecución de la herramienta. Este script es presentado
en el Código B.1. En la Figura B.2 se muestra un esquema de tipo diagrama
de flujo donde se muestran las etapas de ejecución. En el diagrama de flujo,
las elipses representan bloques de código m con una función específica, no
incluidos dentro de un archivo independiente. Los restantes bloques rectan-
gulares representan ejecuciones de archivos y los bloques con forma de rombo
representan estructuras de control de toma de decisión.
reading: script que realiza la lectura de datos a partir de los archivos de en-
trada (colocados en la carpeta input). El archivo se presenta en el Código B.2.
138 Apéndice B. ONSAS
ONSAS
ONSAS.m
sources
assembly.m
hyperElasModels.m
linearAnalysis.m
NLAPlots.m
nodes2dofs.m
NRALIter.m
NRIter.m
reading.m
setini.m
stressUpdate.m
varsVerifs.m
input
output
setini: script que crea matrices y vectores y asigna valores iniciales necesarios
para la ejecución del código. El archivo se presenta en el Código B.4.
reading.m y
varsVerifs.m
setini.m
¿LinBuckMode-
si Imperf > 0?
no
assembly.m
NRIter.m o
NRALIter
¿Criterio
no parada?
si
¿Carga objetivo no
alcanzada?
si
NLAPlots
88
89 if LinBuckModeImperf > 0
90
91 fprintf([' Imperfection is added. Max distance imperfection: %7.2e % % of the'
...
92 ' normalized %2i−th LBA mode.\n'], absoluteimperfection, LinBuckModeImperf )
93 % deforms the reference structure considering the chosen imperfection
94 Nodes = Nodes + reshape(UsModesLinBuck(:,LinBuckModeImperf),2,nnodes)' ...
95 * absoluteimperfection ;
96
97 % update geometry and related parameters
98 setini
99
100 fprintf(' Linear analysis with imperfection perfomed. \n')
101 linearAnalysis
102
103 end
104
105 % plots the reference structure
106 lw = 3;
107 if plotflag > 0
108 figure, hold on, grid on
109 for i=1:nelems
110 plot( xelems(i,:), yelems(i,:), 'b−o', 'linewidth', lw)
111 end
112 axis equal
113 end
114
115 quiver( Nodes(:,1), Nodes(:,2) , Fext(1:2:end)*visualloadfactor ...
116 , Fext(2:2:end)*visualloadfactor ,0,'g',filled,'linewidth',lw)
117 title( ['Reference configuration without imperfections (blue) and external ' ...
118 ' loads (green).'] ) ; xlabel('x'), ylabel('y')
119
120
121 %% Nonlinear analysis
122 % Core of the code: nonlinear incremental and iterative analysis. Newton−Raphson
and Cylindrical−Arc−Length variation.
123 %
124 %%
125
126 % initialize
127 reachedtargetLF = 0;
128 loadIter = 0 ; % counter of the load iterations
129 loadfactors = 0;
130
131 % parameters for the Arc−Length iterations
132 convDeltau = zeros(length(neumdofs),1) ;
133 currDeltau = zeros(length(neumdofs),1) ;
134 currLoadFactor = 0 ;
135
136 fprintf('−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− \n')
137 fprintf([ 'loadIncr & LoadFactr & dispits & maxStrain ( % %) & BucklingFac ' ...
138 ' & npos & nneg & \n'] )
139 prodvals = 1 ;
144 Apéndice B. ONSAS
140
141 filestressoutput = [outputdir '/output_stresses.txt' ];
142 fidstress = fopen(filestressoutput,'w');
143
144 STRESSESMAT = [] ;
145
146 tic;
147 % load step increment stopping criteria: target load
148 while ( reachedtargetLF == 0 )
149
150 loadIter += 1 ;
151 dispConverged = 0 ;
152 dispIter = 0 ;
153 currDeltau = zeros(length(neumdofs),1) ;
154
155 %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
156 % iteration in displacements (NR) or load−displacements (NR−AL)
157 % stopping criteria: convergence displacements
158 while ( dispConverged == 0 )
159 dispIter += 1 ;
160
161 assembly
162
163 if solutionMethod == 1
164 % performs one newton−raphson iteration
165 NRIter
166
167 elseif solutionMethod == 2
168 % performs one newton−raphson−arc−length iteration
169 NRALIter
170 end
171
172 stressUpdate
173
174 % stopping criteria
175 if ( ( normadeltau < ( normaUk * NRtoldeltau ) ) || dispIter > NRtolits )
176 if dispIter > NRtolits,
177 stopcritpar = 2;
178 fprintf('Warning: displacements iteration stopped by max iterations.\n');
179 else,
180 stopcritpar = 1;
181 end
182 dispConverged = 1 ;
183 Itsnums( loadIter ) = dispIter;
184 end
185 %−−−−−−−−−−−−−−−−−−−−−−−−−
186
187 end
188 %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
189
190 loadfactors(loadIter) = currLoadFactor ;
191 uks(loadIter) = sign(controldof) * Uk(abs(controldof)) ;
192 convDeltau = currDeltau ;
193
Sección B.1. Descripción de código 145
45
46 Uklocelem = R' * Uklin(dofselem) ;
47
48 epsk = ( Uklocelem(3) − Uklocelem(1) ) / largosini(elem) ;
49
50 stresslink(elem) = dsigdepsk(elem) * epsk ;
51
52 A = As(Conec(elem,4)) ;
53 le = largosini(elem) ;
54
55 Ksige = A * stresslink(elem) / le * Ge ;
56 KGlin (dofselem,dofselem) = KGlin (dofselem,dofselem) + Ksige ;
57
58 end
59
60 KGlinred = KGlin ( neumdofs,neumdofs) ;
61
62 % linear buckling analysis − LBA
63 [ ModesLinBuck, lambdas] = eig( KL0red , − KGlinred ) ;
64
65 [lambdas,perm] = sort(diag(lambdas));
66 ModesLinBuck=ModesLinBuck(:,perm);
67
68 factor_crit_lin = min( lambdas );
69
70 if factor_crit_lin < 0,
71
72 indsneg = find( lambdas <0) ;
73
74 fprintf(' ===== Warning: %3i negative linear buckling analysis factors! ===== \n
Removing negative eigenmodes.... \n',length(indsneg) );
75
76 lambdas(indsneg) = [] ;
77 ModesLinBuck(:,indsneg) = [] ;
78
79 fprintf(' =========================================================\n');
80 end
81
82 factor_crit_lin = min( lambdas ) ;
83
84 nmodesLinBuck = length(lambdas) ;
85 UsModesLinBuck = zeros( 2*nnodes, nmodesLinBuck ) ;
86
87 for i=1:nmodesLinBuck
88 uaux = ModesLinBuck(:,i) ;
89 UsModesLinBuck(neumdofs,i) = uaux / norm( uaux, 'inf' ) ;
90 end
91
92
93 % −−−−−−−−−−− plots −−−−−−−−−−−−−−−−−−−−−−−−−−−−
94
95 lw = 1.5;
96 ms = 3;
97
Sección B.1. Descripción de código 151
98 linfig = figure ;
99 title('Linear analysis results')
100
101 subplot(2,2,1)
102 hold on, grid on
103
104 displacmat = reshape(Uklin,2,nnodes)' / norm( Uklin, 'inf' ) * (0.1 * strucsize) ;
105 Nodesdefaux = Nodes + displacmat ;
106
107 for i=1:nelems
108 plot(xelems(i,:),yelems(i,:),'b−o','linewidth',lw,'markersize',ms)
109 end
110
111 for i=1:nelems
112 xelemsaux(i,:) = Nodesdefaux( Conec(i,1:2) , 1);
113 yelemsaux(i,:) = Nodesdefaux( Conec(i,1:2) , 2);
114 end
115
116 for i=1:nelems
117 plot(xelemsaux(i,:),yelemsaux(i,:),'r−s','linewidth',lw,'markersize',ms)
118 end
119
120 axis equal
121 title('Deformed structure (scaled)')
122
123 subplot(2,2,2); hold on, grid on
124
125 tolstress = 1e−7*max(abs(stresslink));
126
127 for i=1:nelems
128 if stresslink(i) > tolstress, colornormalforce='b';
129 elseif stresslink(i) <−tolstress, colornormalforce='r';
130 else colornormalforce='w'; end
131
132 plot(xelems(i,:),yelems(i,:),[colornormalforce '−o'],'linewidth',lw,'markersize'
,ms)
133 end
134
135 axis equal
136 title('Normal forces: >0 blue, <0 red, ==0 white.')
137
138 % −−−−−−−−−−− plots −−−−−−−−−−−−−−−−−−−−−−−−−−−−
139 subplot(2,2,3), hold on, grid on
140
141 modeaux = 1;
142
143 displacmat = reshape(UsModesLinBuck(:,modeaux),2,nnodes)' * (0.1*strucsize) ;
144 Nodesdefaux = Nodes + displacmat ;
145
146
147 for i=1:nelems
148 plot(xelems(i,:),yelems(i,:),'b−o','linewidth',lw,'markersize',ms)
149 end
150
152 Apéndice B. ONSAS
18
19 B1e = 1.0 / ( le^2 ) * Xe' * Ge' ; B2e = 1.0 / ( le^2 ) * Ue' * Ge' ;
20
21 KT1e = dsigdepsk(elem) * A * le * ( B1e' * B1e ) ;
22 KT2e = dsigdepsk(elem) * A * le * ( B2e' * B1e + B1e' * B2e + B2e' * B2e ) ;
23 Ksige = A * Stressk(elem) / le * Ge ;
24 Finte = (B1e+B2e)' * A * le * Stressk(elem) ;
25
26 % matrices assembly
27 KL0 (dofselem,dofselem) = KL0 (dofselem,dofselem) + KT1e ;
28 KM (dofselem,dofselem) = KM (dofselem,dofselem) + KT1e + KT2e ;
29 KG (dofselem,dofselem) = KG (dofselem,dofselem) + Ksige ;
30
31 % internal loads vector assembly
32 FintGk ( dofselem) = FintGk(dofselem) + Finte ;
33
34 end
35
36 KM = KM + KS ;
37 KL0 = KL0 + KS ;
38
39 FintGk = FintGk + KS * Uk ;
40
41 % boundary conditions are applied
42 KL0red = KL0 ( neumdofs, neumdofs );
43 KMred = KM ( neumdofs, neumdofs );
44 KGred = KG ( neumdofs, neumdofs );
45 KTred = KMred + KGred ;
10
11 Resred = FintGk(neumdofs) − FextG(neumdofs) ;
12
13 % incremental displacement
14 deltaured = KTred \ ( − Resred ) ;
15
16 normadeltau = norm( deltaured ) ;
17 normaUk = norm( Uk(neumdofs ) ) ;
18
19 Uk ( neumdofs ) = Uk(neumdofs ) + deltaured ;
42
43 Uk ( neumdofs ) = Uk(neumdofs ) + deltaured ;
44 currDeltau = currDeltau + deltaured ;
24 end
25
26 case 3
27 % Bi−modulus material model with pre−strain: param1: Tension young modulus,
28 % param2: Compression modulus, param3: epszero (pre−strain)
29 ET = paramsmodel(1) ;
30 EC = paramsmodel(2) ;
31 epszero = paramsmodel(3) ;
32
33 if ( epsilon − epszero ) >=(−1e−15)
34 sigma = ET * ( epsilon − epszero ) ;
35 dsigdeps = ET ;
36 else
37 sigma = EC * ( epsilon − epszero ) ;
38 dsigdeps = EC ;
39 end
40
41 case 4
42 % Example 5.11 from Reddy nonlinear sqrt stress−strain constitutive relation
43 E = paramsmodel(1) ;
44 sigma = sign(epsilon) * E * sqrt( abs( epsilon) ) ;
45 if epsilon == 0
46 dsigdeps = E^2 ;
47 else
48 dsigdeps = 0.5 * E / sqrt( abs(epsilon) )
49 end
50
51 end
24
25 axis equal
26
27 quiver( Nodesdef(:,1), Nodesdef(:,2) , Fext(1:2:end)*visualloadfactor ...
28 , Fext(2:2:end)*visualloadfactor ,0,'g',filled,'linewidth',lw2)
29
30 labx=xlabel('x'), laby=ylabel('y')
31 set(gca, 'linewidth', 1.2, 'fontsize', plotfontsize )
32 set(labx, FontSize, plotfontsize); set(laby, FontSize, plotfontsize);
33
34 if printflag == 1
35 cd(outputdir )
36 print( [ problemname '_deform' ] ,'−depslatex') ;
37 cd(currdir)
38 end
39
40 figure
41 plot(uks,LoadFactrs,'b.−','markersize',ms2,'linewidth',lw2)
42 grid on
43 labx=xlabel('Displacements'), laby=ylabel('Load factors')
44 set(gca, 'linewidth', 1, 'fontsize', plotfontsize )
45 set(labx, FontSize, plotfontsize); set(laby, FontSize, plotfontsize);
46
47 if printflag == 1
48 cd( outputdir )
49 print( problemname ,'−depslatex') ;
50 cd(currdir)
51 end
13 nelemsviga = 22 ;
14 nnodesviga = nelemsviga +1 ;
15 angle_alpha_arch = pi/6 ;
16 xnodsviga = linspace( −angle_alpha_arch+pi/2,angle_alpha_arch+pi/2,
nnodesviga )' ;
17
18 ri = 49; re = 50;
19
20 Rv = (ri+re)*.5 ;
21 hv = (re−ri) ;
22 Iv = As*hv^2/2 ;
23
24 PcritTeo = E*Iv / (Rv^3) * ( pi^2/angle_alpha_arch^2 −1 )
25
26 % −−−−−−−−− Node coordinates matrix −−−−−−−−
27 nnodes = 2*nnodesviga ;
28
29 Nodes = zeros(nnodes,2);
30
31 % polar coordinates matrix
32 Nodes(1:2:end,1) = xnodsviga; Nodes(2:2:end,1) = xnodsviga;
33 Nodes(1:2:end,2) = ones(nnodesviga,1)*ri; Nodes(2:2:end,2) = ones(nnodesviga,1)*re
;
34
35 radios = Nodes(:,2); thetas = Nodes(:,1);
36
37 % polar to cartesian transformation
38 Nodes = [ radios.*cos(thetas) radios.*sin(thetas) ] ;
39
40 # modify nodes to get pinned nodes
41 aux1 = Nodes (1 ,:)*0.5 + Nodes (2 ,:)*0.5 ;
42 aux2 = Nodes (end−1,:)*0.5 + Nodes (end,:)*0.5 ;
43
44 Nodes(1,:) = aux1 ; Nodes(2,:) = [] ;
45 Nodes(end−1,:) = aux2 ; Nodes(end ,:) = [] ;
46
47 nnodes = size(Nodes,1);
48 %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
49
50 ndia=5;
51 for i=2:(nelemsviga−1)
52 auxi = (i−1)*2+1−1 ;
53 Conec ( ((i−2)*ndia+1):((i−1)*ndia) ,: ) = [ ...
54 auxi auxi+2 1 1 ; ... % 1 3 cordon
55 auxi auxi+3 1 1 ; ... % 1 4 diagon
56 auxi+1 auxi+2 1 1 ; ... % 2 3 diagon
57 auxi+1 auxi+3 1 1 ; ... % 2 4 cordon
58 auxi+2 auxi+3 1 1 ] ; % 3 4 alma
59 end
60
61 Conec = [ 1 2 1 1; 1 3 1 1; 2 3 1 1; Conec ];
62 Conec = [ Conec; nnodes−2 nnodes−1 1 1; nnodes−2 nnodes 1 1; nnodes−1 nnodes 1 1
];
63
Sección B.2. Ejemplos de archivos de entrada 159
64 fixednodes = [ 1 nnodes ];
65 fixeddofs = nodes2dofs( fixednodes,2) ;
66
67 Lv = Rv*(2*angle_alpha_arch) ;
68 lelemviga = Lv / nelemsviga ;
69
70 nodesload = (2):2:(nnodes−2) ;
71
72 % external unitary loads vector
73 Fext = zeros(2*nnodes,1);
74 for i=1:length( nodesload )
75 coordaux = Nodes(nodesload(i),:)' ;
76 dofsnodesload = [ 2*nodesload(i)−1 2*nodesload(i) ]' ;
77 Fext( dofsnodesload ) = − coordaux / norm(coordaux) * lelemviga ;
78 end
79
80 refnodeload = (nnodesviga+1) ;
81
82 Springs = [ ];
83
84 controldof = −2*(nnodes/2);
85
86 targetLoadFactr = 2.5*PcritTeo ;
87
88 % params imperfeccion
89 %~ nLoadSteps = 25 ; incremarclen = 0.1 ; % limit point
90 %~ nLoadSteps = 200 ; incremarclen = 0.5 ; % post snap−through point
91 nLoadSteps = 150 ; incremarclen = 0.26 ; % zero load unstable point
92
93 LinBuckModeImperf = 0 ;
94 printflag = 1 ;
95
96 % Iterative method used: 1 incremental NR 2 Combined NR/arc−length method
97 solutionMethod = 2 ;
98
99 % nonlinear iteration parameters
100 NRtolits = 30 ;
101 NRtoldeltau = 1.0e−5 ;
102
103 diridofs = fixeddofs ;
Código B.14: Archivo de entrada para ejemplo de arco bi-articulado con imperfec-
ción.
1 %−−−−−−−−−−−−−−−−−
2 % A two−dimensional pinned circular arch submitted to radial distributed load.
3 %−−−−−−−−−−−−−−−−−
4
5 problemname = 'Pinned_Circular_Arch_imperf' ;
6
7 E = 210e9 ;
8 HyperElasParams = cell(1,1) ;
9 HyperElasParams{1} = [1 E] ;
160 Apéndice B. ONSAS
10
11 As = 0.4*.01 ;
12
13 nelemsviga = 22 ;
14 nnodesviga = nelemsviga +1 ;
15 angle_alpha_arch = pi/6 ;
16 xnodsviga = linspace( −angle_alpha_arch+pi/2,angle_alpha_arch+pi/2,
nnodesviga )' ;
17
18 ri = 49;
19 re = 50;
20
21 Rv = (ri+re)*.5 ;
22 hv = (re−ri) ;
23 Iv = As*hv^2/2 ;
24
25 PcritTeo = E*Iv / (Rv^3) * ( pi^2/angle_alpha_arch^2 −1 )
26
27
28 % −−−−−−−−− Node coordinates matrix −−−−−−−−
29 nnodes = 2*nnodesviga ;
30
31 Nodes = zeros(nnodes,2);
32
33 % polar coordinates matrix
34 Nodes(1:2:end,1) = xnodsviga; Nodes(2:2:end,1) = xnodsviga;
35 Nodes(1:2:end,2) = ones(nnodesviga,1)*ri; Nodes(2:2:end,2) = ones(nnodesviga,1)*re
;
36
37 radios = Nodes(:,2); thetas = Nodes(:,1);
38
39 % polar to cartesian transformation
40 Nodes = [ radios.*cos(thetas) radios.*sin(thetas) ] ;
41
42 # modify nodes to get pinned nodes
43 aux1 = Nodes (1 ,:)*0.5 + Nodes (2 ,:)*0.5 ;
44 aux2 = Nodes (end−1,:)*0.5 + Nodes (end,:)*0.5 ;
45
46 Nodes(1,:) = aux1 ; Nodes(2,:) = [] ;
47 Nodes(end−1,:) = aux2 ; Nodes(end ,:) = [] ;
48
49 nnodes = size(Nodes,1);
50 %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
51
52 ndia=5;
53 for i=2:(nelemsviga−1)
54 auxi = (i−1)*2+1−1 ;
55 Conec ( ((i−2)*ndia+1):((i−1)*ndia) ,: ) = [ ...
56 auxi auxi+2 1 1 ; ... % 1 3 cordon
57 auxi auxi+3 1 1 ; ... % 1 4 diagon
58 auxi+1 auxi+2 1 1 ; ... % 2 3 diagon
59 auxi+1 auxi+3 1 1 ; ... % 2 4 cordon
60 auxi+2 auxi+3 1 1 ] ; % 3 4 alma
61 end
Sección B.2. Ejemplos de archivos de entrada 161
62
63 Conec = [ 1 2 1 1; 1 3 1 1; 2 3 1 1; Conec ];
64 Conec = [ Conec; nnodes−2 nnodes−1 1 1; nnodes−2 nnodes 1 1; nnodes−1 nnodes 1 1
];
65
66 %~ fixednodes = [ 1 2 nnodes−1 nnodes ];
67 fixednodes = [ 1 nnodes ];
68 fixeddofs = nodes2dofs( fixednodes,2) ;
69
70 Lv = Rv*(2*angle_alpha_arch) ;
71 lelemviga = Lv / nelemsviga ;
72
73 nodesload = (2):2:(nnodes−2) ;
74
75 % external unitary loads vector
76 Fext = zeros(2*nnodes,1);
77 for i=1:length( nodesload )
78 coordaux = Nodes(nodesload(i),:)' ;
79 dofsnodesload = [ 2*nodesload(i)−1 2*nodesload(i) ]' ;
80 Fext( dofsnodesload ) = − coordaux / norm(coordaux) * lelemviga ;
81 end
82
83 refnodeload = (nnodesviga+1) ;
84
85 Springs = [ ];
86
87 controldof = −2*(nnodes/2);
88 targetLoadFactr = 2.5*PcritTeo ;
89
90 nLoadSteps = 400 ;
91 incremarclen = 0.08;
92 LinBuckModeImperf = 1 ;
93
94 printflag = 1 ;
95
96 % Iterative method used: 1 incremental NR 2 Combined NR/arc−length method
97 solutionMethod = 2 ;
98
99 % nonlinear iteration parameters
100 NRtolits = 30 ;
101 NRtoldeltau = 1.0e−5 ;
102
103 diridofs = fixeddofs ;
Código B.16: Archivo de entrada para ejemplo de cercha de Von Mises con imper-
fección.
1 %−−−−−−−−−−−−−−−−−−−
2 % Example 3.1 from journal article: Li and Khandelwal, 'Topology
Sección B.2. Ejemplos de archivos de entrada 163
Alhasawi, A., Heng, P., Hjiaj, M., Guezouli, S., and Battini, J.-M. (2017). Co-rotational
planar beam element with generalized elasto-plastic hinges. Engineering Struc-
tures, 151:188–205.
Belytschko, T., Liu, W. K., and Moran, B. (2014). Nonlinear Finite Elements for Conti-
nua and Structures. Wiley, 2 edition.
Bigoni, D., Bosi, F., Dal Corso, F., and Misseroni, D. (2014). Instability of a penetrating
blade. Journal of the Mechanics and Physics of Solids, 64:411–425.
Bonet, J., Barros, M., and Romero, M. (2006). Comparative study of analytical and
numerical algorithms for designing reinforced concrete sections under biaxial
bending. Computers & Structures, 84(31-32):2184–2193.
Bonet, J. and Wood, R. (2008). Nonlinear Continuum Mechanics for Finite Element
Analysis. Cambridge University Press.
Bridgens, B. and Birchall, M. (2012). Form and function: The significance of mate-
rial properties in the design of tensile fabric structures. Engineering Structures,
44:1–12.
Castrillo, P., Mondino, F., Pérez Zerpa, J. M., and Canelas, A. (2014). Desarrollo y
extensión de una herramienta numérica de elementos finitos para el dictado de
cursos de grado y de posgrado. In Mecánica Computacional, pages 2073–2086.
Asociacion Argentina de Mecanica Computacional.
Crisfield, M. A. (1997). Non-Linear Finite Element Analysis Solids and Structure, Vo-
lume 2, Advanced Topics. Wiley.
de Borst, R., Crisfield, M. A., Remmers, J. J. C., and Verhoosel, C. V. (2012). Non-
Linear Finite Element Analysis of Solids and Structures. John Wiley & Sons, Ltd,
Chichester, UK.
de Souza Neto, E. A., Peri, D., and Owen, D. R. J. (2008). Computational Methods for
Plasticity. John Wiley & Sons, Ltd, Chichester, UK.
Delfino, A., Stergiopulos, N., Moore, J. E., and Meister, J. J. (1997). Residual strain
effects on the stress field in a thick wall finite element model of the human
carotid bifurcation. Journal of biomechanics, 30(8):777–786.
Du, Z. and Guo, X. (2014). Variational principles and the related bounding theorems
for bi-modulus materials. Journal of the Mechanics and Physics of Solids, 73:183–
211.
Eaton, J. W., Bateman, D., Hauberg, S., and Wehbring, R. (2015). GNU Octave: a high-
level interactive language for numerical computations. CreateSpace Independent
Publishing Platform, 4 edition.
Feng, R.-Q., Ye, J., Yan, G., and Ge, J.-M. (2013). Dynamic Nonlinearity and Nonlinear
Single-Degree-of-Freedom Model for Cable Net Glazing. Journal of Engineering
Mechanics, 139(10):1446–1459.
Gilsanz, R., Hamburger, R., Barker, D., Smith, J. L., and Rahimian, A. (2013). De-
sign Guide 26: Design of Blast Resistant Structures. American Institute of Steel
Construction.
Goenezen, S., Dord, J.-F., Sink, Z., Barbone, P. E., Jiang, J., Hall, T. J., and Oberai,
A. A. (2012). Linear and Nonlinear Elastic Modulus Imaging: An Application
to Breast Cancer Diagnosis. IEEE Transactions on Medical Imaging, 31(8):1628–
1637.
Zaccaria, D., Bigoni, D., Noselli, G., and Misseroni, D. (2011). Structures buckling un-
der tensile dead load. Proceedings of the Royal Society A: Mathematical, Physical
and Engineering Sciences, 467(2130):1686–1700.
Zhang, L., Gao, Q., and Zhang, H. (2013). An efficient algorithm for mechanical
analysis of bimodular truss and tensegrity structures. International Journal of
Mechanical Sciences, 70:57–68.
Zhang, L., Zhang, H. W., Wu, J., and Yan, B. (2016). A stabilized complementarity
formulation for nonlinear analysis of 3D bimodular materials. Acta Mechanica
Sinica, 32(3):481–490.
Zienkiewicz, O. C. (1972). Introductory Lectures on the Finite Element Method. Sprin-
ger Vienna, Vienna.
Zienkiewicz, O. C., Taylor, R. L., and Fox, D. D. (2014). The Finite Element Method for
Solid and Structural Mechanics. Elsevier.