Fem Pract - 2
Fem Pract - 2
Fem Pract - 2
EA
L
EA
L
_
(1)
x
i E, A
L
j
Figura 1: Elemento de barra lineal.
Es evidente que la matriz de rigidez elemental para el elemento de barra lineal es similar
a la del resorte elemental con la rigidez reemplazada por EA/L. Es claro que el elemento de
barra elemental tiene s olo dos grados de libertad - uno en cada nodo. Consecuentemente para
una estructura con n nodos, la matriz de rigidez global K tendr a dimension n n (ya que
tenemos un grado de libertad en cada nodo). La matriz de rigidez global K es ensamblada por
la funcion de MatLab LinearBarAssemble que esta escrita especcamente para este prop osito.
Este proceso se ilustra en detalle en los ejemplos.
Una vez que ha sido obtenida la matriz de rigidez global K, tendremos la siguiente
ecuaci on matricial:
[K] {U} = {F} (2)
EAP IH - UNC MTODOS NUMRICOS APLICADOS 2012-1
1
donde U es el vector global de desplazamientos nodales y F es el vector de fuerzas global nodal.
En este paso las condiciones de contorno se aplican manualmente a los vectores U y F.
A continuaci on, el sistema (2) es resuelto por particionamiento seguido de eliminacion Gaussiana.
Por ultimo una vez que los desplazamientos y reacciones desconocidas han sido encontradas, las
fuerzas de los elementos se obtienen para cada elemento de la siguiente manera:
{f} = [k] {u} (3)
en donde f es un 2 1 vector elemental de fuerzas y u es un 2 1 vector elemental de despla-
zamientos. La tension elemental es obtenida dividiendo las fuerzas elementales por el area de la
secci on transversal A.
2. Uso de las funciones de MatLab.
Las cuatro funciones de MatLab usadas para el caso del elemento de barra lineal son:
LinearBarElementStiness(E,A,L) Esta funcion calcula la matriz de rigidez elemental
por cada barra elemental con m odulo de elasticidad E, area de secci on transversal A y
longitud L. Retorna una matriz de rigidez elemental de 2 2.
LinearBarAssemble(K, k, i, j) Estas funciones ensamblan la matriz de rigidez elemental
k de la barra lineal, uniendo los nodos i (del extremo izquierdo) y j (del extremo derecho)
en la matriz global de rigidez K. Se retorna una matriz global de rigidez K nn cada vez
que es ensamblado un elemento.
LinearBarElementForces(k, u) Esta funcion calcula el vector elemental de fuerzas
usando la matriz de rigidez elemental k y el vector elemental de desplazamiento u. Se
retorna un 2 1 vector elemental de fuerzas f.
LinearBarElementStresses(k, u, A) Esta funcion calcula el vector elemental de tension
usando la matriz elemental de rigidez k, el vector elemental de desplazamiento u y el area
de secci on transversal A. Retorna un 2 1 vector elemental de tension sigma o s.
EAP IH - UNC MTODOS NUMRICOS APLICADOS 2012-1
2
La siguiente es una lista de codigos fuente en MatLab para cada funcion:
function y = LinearBarElementStiffness(E,A,L)
% LinearBarElementStiffness Esta funci on retorna la matriz de rigidez
% elemental para una barra lineal con m odulo de elasticidad E, area de
% secci on transversas A y longitud L.
% La dimensi on de la matriz de rigidez elemental es 2 x 2.
y = [E*A/L -E*A/L ; -E*A/L E*A/L];
function y = LinearBarAssemble(K,k,i,j)
% LinearBarAssemble Esta funci on ensambla la matriz de rigidez elemental k
% de una barra lineal con nodos i y j en la matriz de rigidez global K.
% Esta funci on retorna una matriz de rigidez global K luego que hayan sido
% ensambladas las matrices de rigidez elemental.
K(i,i) = K(i,i) + k(1,1) ;
K(i,j) = K(i,j) + k(1,2) ;
K(j,i) = K(j,i) + k(2,1) ;
K(j,j) = K(j,j) + k(2,2) ;
y = K;
function y = LinearBarElementForces(k,u)
% LinearBarElementForces Esta funci on retorna un vector de fuerza
% nodal elemental dada la matriz de rigidez elemental k y
% el vector elemental de desplazamiento nodal u.
y = k * u;
function y = LinearBarElementStresses(k, u, A)
% LinearBarElementStresses Esta funci on retorna el vector de tensi on nodal
% elemental dada la matriz de rigidez elemental k, el vector elemental de
% desplazamiento nodal u y el area de secci on transversal A.
y = k * u/A;
EAP IH - UNC MTODOS NUMRICOS APLICADOS 2012-1
3
Ejemplo 1 Considerar la estructura conformada por dos barras lineales como se muestra en la
Fig. 2. Sea E = 210 GPa, A = 0,003 m
2
, P = 10 kN y el nodo 3 es desplazado hacia la derecha
en 0,002 m. Determinar:
1. la matriz de rigidez global para la estructura.
2. el desplazamiento en el nodo 2.
3. las reacciones en los nodos 1 y 3.
4. la tensi on en cada barra.
1
1.5 m 1 m
3 2
P
Figura 2: Estructura con dos barras para el ejemplo 1.
SOLUCI
ON
Paso 1. Discretizaci on del dominio:
Este problema ya est a discretizado. El dominio ha sido dividido en dos elementos y tres no-
dos. Las unidades usadas en los c alculos del MatLab son kN y metros. La Tabla 1 muestra la
conectividad elemental para este ejemplo:
TABLA 1. Conectividad elemental para el ejemplo 1.
N umero de elemento Nodo i Nodo j
1 1 2
2 2 3
Paso 2. Escritura de las matrices elementales de rigidez:
Las dos matrices de rigidez elemental k
1
y k
2
son obtenidas por medio de la funci on de MatLab
LinearBarElementStiness. Cada matriz tiene dimension 2 2.
EAP IH - UNC MTODOS NUMRICOS APLICADOS 2012-1
4
>> E=210e6
E =
210000000
>> A=0.003
A =
0.0030
>> L1=1.5
L1 =
1.5000
>> L2=1
L2 =
1
>> k1=LinearBarElementStiffness(E,A,L1)
k1 =
420000 -420000
-420000 420000
>> k2=LinearBarElementStiffness(E,A,L2)
k2 =
630000 -630000
-630000 630000
Paso 3. Ensamble de la matriz de rigidez global.
Dado que la estructura tiene tres nodos, el tama no de la matriz de rigidez global es de 33.
Por lo tanto para obtener K primeramente debemos construir una matriz nula de 3 3 luego
hacer dos llamados a la funci on de MatLab LinearBarAssemble dado que tenemos dos barras
lineales elementales en la estructura. Cada llamada a la funci on ensamblara un elemento,
EAP IH - UNC MTODOS NUMRICOS APLICADOS 2012-1
5
>> K=zeros(3,3)
K =
0 0 0
0 0 0
0 0 0
>> K=LinearBarAssemble(K,k1,1,2)
K =
420000 -420000 0
-420000 420000 0
0 0 0
>> K=LinearBarAssemble(K,k2,2,3)
K =
420000 -420000 0
-420000 105000 -630000
0 -630000 630000
Paso 4. Aplicacion de las condiciones de contorno.
El sistema (2) para este caso es obtenido utilizando la matriz de rigidez global del paso anterior:
_
_
420000 420000 0
420000 1050000 630000
0 630000 630000
_
_
_
_
U
1
U
2
U
3
_
_
=
_
_
F
1
F
2
F
3
_
_
(4)
Las condiciones de contorno para este problema son dadas por:
U
1
= 0, F
2
= 10, U
3
= 0,002 (5)
Insertando ambas condiciones en (4), obtenemos:
_
_
420000 420000 0
420000 1050000 630000
0 630000 630000
_
_
_
_
0
U
2
0,002
_
_
=
_
_
F
1
10
F
3
_
_
(6)
EAP IH - UNC MTODOS NUMRICOS APLICADOS 2012-1
6
Paso 5. Resolviendo las ecuaciones.
La soluci on del sistema de ecuaciones (6) se llevar a a cabo mediante un particionamiento
(manual) y posterior eliminaci on Gaussiana (con MatLab). Primeramente particionamos (6)
extrayendo la submatriz de la la 2 y columna 2 que resulta ser una matriz 1 1. Debido al
desplazamiento de 0.002 m efectuado por el nodo 3, tenemos que extraer la submatriz de la la
2 y columna 3, que tambien resulta ser una matriz 1 1 Luego obtenemos:
[1050000] U
2
+ [630000] (0,002) = {10} (7)
La soluci on de ambos sistemas es obtenido usando MatLab del siguiente modo. Notar que el
operador \ (backslash) es usado para la eliminaci on Gaussiana.
>> k=K(2,2)
k =
1050000
>> k0=K(2,3)
k0 =
-630000
>> u0=0.002
u0 =
0.0020
>> f=[-10]
f =
-10
>> f0=f-k0*u0
f0 =
1250
>> u=k\f0
EAP IH - UNC MTODOS NUMRICOS APLICADOS 2012-1
7
u =
0.0012
Ahora est a claro que el desplazamiento en el nodo 2 es 0.0012 m.
Paso 6. Post-procesamiento.
En este paso obtenemos las reacciones en los nodos 1 y 3, y la tensi on en cada barra usando
MatLab del siguiente modo. Primeramente creamos el vector global de desplazamientos nodales
U y luego calculamos el vector global de fuerzas nodales F.
>> U=[0;u;u0]
U =
0
0.0012
0.0020
>> F=K*U
F =
-500.0000
-10.0000
510.0000
As, las reacciones en los nodos 1 y 3 son fuerzas de 500 kN (dirigida hacia la izquierda)
y 510 kN (dirigida hacia la derecha), respectivamente. Queda claro que la fuerza de equilibrio
se cumple. Luego establecemos los vectores elementales de desplazamiento nodal u
1
y u
2
; luego
calculamos el vector elemental de fuerzas f
1
y f
2
por medio de la funci on de MATLAB Li-
nearBarElementForces. Finalmente dividimos cada fuerza elemental entre el area de la secci on
transversal del elemento para obtener las tensiones elementales,
>> u1=[0;U(2)]
u1 =
0
0.0012
EAP IH - UNC MTODOS NUMRICOS APLICADOS 2012-1
8
>> f1=LinearBarElementForces(k1,u1)
f1 =
-500.0000
500.0000
>> sigma1=f1/A
sigma1 =
1.0e+005 *
-1.6667
1.6667
>> u2=[U(2) ; U(3)]
u2 =
0.0012
0.0020
>> f2=LinearBarElementForces(k2,u2)
f2 =
-510.0000
510.0000
>> sigma2=f2/A
sigma2 =
1.0e+005 *
-1.7000
1.7000
La tensi on en el elemento 1 es 1,667 10
5
kN/m
2
(o 166,7 MPa de tensi on); y la tensi on en el
elemento 2 es 1,7 10
5
kN/m
2
(o 170 MPa de tensi on). Alternativamente, podemos obtener la
tensi on elemental directamente por medio de la funci on de MatLab LinearBarElementStresses.
Esto se realiza de la siguiente manera, obtieniendose los mismos resultados.
EAP IH - UNC MTODOS NUMRICOS APLICADOS 2012-1
9
>> s1=LinearBarElementStresses(k1,u1,A)
s1 =
1.0e+005 *
-1.6667
1.6667
>> s2=LinearBarElementStresses(k2,u2,A)
s2 =
1.0e+005 *
-1.7000
1.7000
Indice
1. Ecuaciones basicas. 1
2. Uso de las funciones de MatLab. 2
EAP IH - UNC MTODOS NUMRICOS APLICADOS 2012-1
10