• Author / Uploaded
• Lucas Nieva

Citation preview

33 4.

CÁLCULO DE LAS DEFORMACIONES ESTÁTICAS EN SISTEMAS CONTINUOS POR FEM 4.1

Matriz de rigidez de una barra (o cuerda)

Para obtener la matriz de rigidez del elemento “barra” supongamos que una barra elástica uniforme (longitud L, sección transversal A, el módulo de elasticidad E) está sometida a las dos fuerzas F1 y F2 , aplicadas a los extremos, y se encuentra en el estado de equilibrio: k F F 2

1

FA

FA

u2

u1 Fig. 4.1

Entonces

(4.1)

F1 + Aσ = 0, F2 − Aσ = 0 siendo F la tensión dentro de la barra (F = F/A). Las fuerzas F1,2 se han dirigido a la derecha formalmente, en realidad las fuerzas estarán opuestas para mantener el equilibrio, las fuerzas AF corresponden a la barra estirada. Si recordamos la ley de Hook en forma diferencial: u −u σ = E ε, ε= 2 1 (4.2) L donde g es la deformación relativa de la barra, u1 y u2 son los desplazamientos del extremos (nodos), llegamos a las dos ecuaciones u −u u −u F1 + AE 2 1 = 0, F2 − AE 2 1 = 0 (4.3) L L que pueden ser escritas en forma matricial: (4.4) A E  1 − 1  u1   F1  = .       L − 1 1  u 2  F2  La matriz AE  1 − 1 (4.5) k= L − 1 1  y es la matriz de rigidez del elemento barra. Es la misma que en el caso de resorte (1.1).

4.2

Barra cónica (método Kattan)

Una barra cónica está fija por un extremo y libre por el otro (fig.4.2). Al extremo libre se aplica la fuerza P = 18 kN. Otros datos de la barra son: E = 210 GPa, las secciones transversales Aizq = 0.012 m2 y Adrch = 0.002 m2 . Determinar el desplazamiento del extremo libre.

P

Calcular el desplazamiento del extremo de la barra cónica .

x=0

x=3 Fig. 4.2

4. Cálculo de las deformaciones estáticas en sistemas continuos por MEF

34

Solución exacta: Destacamos en la barra una capa fina de espesor dx, sie ndo “dy” el alargamiento de esta capa. De acuerdo con la ley de Hook (ver 4.3): dy P dx P = A( x) E → dy = (4.6) dx E A( x) Para la función A(x) tenemos la expresión: x (4.7) A ( x ) = A izq − (A izq − A drch ) 3 Integrando desde x= 0 hasta x=3 y sustituyendo los valores numéricos, obtenemos el alargamiento total de la barra: 3 3 P dx 18000 dx y= ∫ = = 4.607·10 - 5 m (4.8) 9 ∫ E 0 A(x ) 210·10 0 0.012 − 0.01· x 3 Solución por MEF (método Kattan) (ver Probl. 3.2 en Soluciones) FUNCIONES UTILIZADAS LinearBarElementStiffness(E,A,L) La misma función que la función SpringAssemble(K,k,i,j)) sustituyendo la constante k por E*A/L. LinearBarAssemble Es exactamente la misma que SpringAssemble(K,k,i,j)

•

Discretización del dominio

1 1

2 2

3 3

4 4

elementos

5 5

6

nodos

Fig. 4.3 Dividimos la barra en cinco elementos iguales de longitud (fig. 4.3). Consideramos que cada elemento tiene la sección transversal constante, igual a su sección en el punto medio. Por tanto: P = 18·103 , E = 210·109 , L = 3/5, A1 = 0.011, A2 = 0.009, A3 = 0.007, A4 = 0.005, A5 = 0.003 (los elementos están numerados de 1 a 5 desde el extremo izquierdo hasta el derecho). La tabla de conectividad: Tabla 2 Número del elemento

Nodo i (izq)

Nodo j (drch)

1

1

2

2

2

3

3

3

4

35

4. Cálculo de las deformaciones estáticas en sistemas continuos por MEF

‚

4

4

5

5

5

6

Creación de la matriz de rigidez elemental

Aplicamos cinco veces la función LinearBarElementStiffness, cada vez con la sección transversal correspondiente al elemento: ki=LinearBarElementStiffness(E,Ai,L)

ƒ

i=1,2,3,4,5

Ensamblaje de la matriz de rigidez global

El sistema tiene seis nodos , por tanto las dimensiones de la matriz de rigidez global son 6 x 6. » K=zeros(6,6) K = 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

Llamamos la función LinearBarAssemble seis veces, de acuerdo con el número de los elementos: » » » » »

K= K= K= K= K=

LinearBarAssemble(K,k1,1,2); LinearBarAssemble(K,k2,2,3); LinearBarAssemble(K,k3,3,4); LinearBarAssemble(K,k4,4,5); LinearBarAssemble(K,k5,5,6);

Al final la matriz de rigidez global será: K = 1.0e+006 * 3.8500 -3.8500 -3.8500 7.0000 0 -3.1500 0 0 0 0 0 0

0 -3.1500 5.6000 -2.4500 0 0

0 0 -2.4500 4.2000 -1.7500 0

0 0 0 -1.7500 2.8000 -1.0500

0 0 0 0 -1.0500 1.0500

(4.9)

„ Aplicación de las condiciones frontera La ecuación matricial de nuestro problema tiene el aspecto:

106

3.85 -3.85 0 0 0 0 -3.85 7.00 -3.15 0 0 0 0 -3.15 5.60 -2.45 0 0 0 0 -2.45 4.20 -1.75 0 0 0 0 -1.75 2.80 -1.05 0 0 0 0 -1.05 1.05

u1 u2 u3 u4 u5 u6

=

f1 f2 f3 f4 f5 f6

(4.10)

Las condiciones frontera son: u1 = 0 , f2 = f3 = f4 = f5 = 0, f6 = 18. Las insertamos en (4.10):

106

3.85 -3.85 0 0 0 0 -3.85 7.00 -3.15 0 0 0 0 -3.15 5.60 -2.45 0 0 0 0 -2.45 4.20 -1.75 0 0 0 0 -1.75 2.80 -1.05 0 0 0 0 -1.05 1.05

0 u2 u3 u4 u5 u6

=

f1 0 0 0 0 18

(4.11)

36

4. Cálculo de las deformaciones estáticas en sistemas continuos por MEF

… Resolución de la ecuación matricial Para resolver el sistema (4.11) extraemos submatrices de las filas y columnas 2,3,4,5,6 (ver Apéndice_submatriz)):

106

7.00 -3.15 0 0 0 -3.15 5.60 -2.45 0 0 0 -2.45 4.20 -1.75 0 0 0 -1.75 2.80 -1.05 0 0 0 -1.05 1.05

u2 u3 u4 u5 u6

=

0 0 0 0 18

(4.12)

Los comandos de MATLAB necesarios para resolución serán: k f u u

= K(2:6,2:6) = [0;0;0;0;18] = k\f = 1.0e-004 * 0.0468 0.1039 0.1774 0.2802 0.4517

Por tanto, el desplazamiento del extremo libre es igual a -4.517·10-5 m, unos 2% menos que el valor exacto.

4.3

Barra cónica (método Kwon)

% DATOS DE ENTRADA % Datos de control nel=5; nnel=2; ndof=1; nnode=6; sdof=nnode*ndof;

% % % % %

número número número número número

de elementos de nodos por elemento de los grados de libertad por nodo total de nodos total de los grados de libertad

% módulo de elasticidad E=210e9; % área del extremo gordo A=0.012; % área del extremo fino a=0.002; % Coordenadas de los nodos Introducimos las coordenadas globales del nodo: gcoord(1)=0.0; gcoord(5)=2.4;

gcoord(2)=0.6; gcoord(6)=3.0;

gcoord(3)=1.2;

gcoord(4)=1.8;

El array gcoord(n) almacena los valores de las coordenadas del nodo n (en nuestro caso solo una coordenada, puesto que el dominio es 1D). % Conectividad. Se trata de la correspondencia entre el número local (dentro del elemento) y número global de cada nodo. (No confundir con la relación entre los GDL locales y globales)

En nuestro caso cada elemento tiene dos nodos, es decir, hay dos nodos locales: izquierdo (1) y dere-

4. Cálculo de las deformaciones estáticas en sistemas continuos por MEF

37

cho (2). La sentencia “nodes(i,j)= k” significa: el nodo local “j” del elemento “i” se coloca en el nodo global “k”. nodes(1,1)=1; nodes(3,1)=3; nodes(5,1)=5;

nodes(1,2)=2; nodes(3,2)=4; nodes(5,2)=6;

nodes(2,1)=2; nodes(4,1)=4;

nodes(2,2)=3; nodes(4,2)=5;

% Condiciones frontera (CF)

Las condiciones frontera incluyen las restricciones impuestas para las coordenadas de los grados de libertad por el enunciado del problema. Todos los grados de libertad están numerados continuamente. Cuando los nodos tienen más de un grado de libertad (por ejemplo, un extremo de una viga), la numeración continua global de los grados de libertad lógicamente no coincide con la numeración continua global de los nodos. Las condiciones frontera también están numeradas globalmente e independientemente. Una condición frontera para el grado de libertad “i” se describe con dos sentencias: bcdof(n) = i; bcval(n) = val. Aquí “n” es el número de la condición frontera, “i” es el número global del GDL condicionado, “val” es el valor impuesto de la coordenada correspondiente. bcdof(1)=1; bcval(1)=0;

% primera CF es la del primer GDL … % …que está inmovilizado

Las fuerzas externas se aplican aparte. %

Puesta a cero de las matrices y vectores

kk=zeros(sdof,sdof); % matriz global de rigidez index=zeros(nnel*ndof,1); % vector de los GDL de uno de los elementos. % Mas explicaciones sobre “index” ver abajo. ff=zeros(sdof,1); % vector global de carga % vector de cargas externas % aplicamos la fuerza externa en el GDL número 6 ff(6)=18e3; %

CÁLCULO DE LA MATRIZ GLOBAL

for iel=1:nel

% para todos los elementos

%¿cuales son los nodos del elemento ‘iel‘? nl=nodes(iel,1); nr=nodes(iel,2); % el nodo izquierdo y el nodo derecho %¿dónde están estos nodos? xl=gcoord(nl); xr=gcoord(nr); % ¿cuál es la longitud del elemento? (por si los elementos son de % diferentes longitudes) eleng=xr-xl; % ¿cuáles son los GDL del elemento ‘iel’? index=feeldof1(iel,nnel,ndof); % function incorporada, ver abajo % En el vector “index” se colocan los GDL globales de todos los GDL % del elemento; podemos llamar este vector “el vector de % conectividad”; % El vector “index” cambia del contenido al pasar de un elemento al

4. Cálculo de las deformaciones estáticas en sistemas continuos por MEF

%

38

otro, es un vector temporal, se utiliza como un “intermediario”

% área del elemento area=A+(A-a)*(0.5-iel)/5; % matriz de rigidez del elemento (matriz del elemento, ver (4.5)) k= (area*E/eleng)*[ 1 -1; -1 1]; % ensamblaje final [kk]=feasmbl1(kk,k,index);

% function incorporada

end % APLICACIÓN DE LAS CONDICIONES FRONTERA

% Antes de resolver el sistema de ecuaciones para los valores nodales es necesario aplicar las condiciones frontera. Sin aplicar las condiciones frontera el número de ecuaciones puede no ser igual al número de incógnitas (es decir, al número de los grados de libertad). [kk,ff]=feaplyc2(kk,ff,bcdof,bcval); % function incorporada

% RESOLVEMOS EL SISTEMA DE ECUACIONES u=kk\ff El resultado es el siguiente (el mismo que por el método de Catan): u = 1.0e-004 * -0.0000 0.0468 0.1039 0.1774 0.2802 0.4517 %--------------------------------------------------------------FUNCIONES INCORPORADAS function [index]=feeldof1(iel,nnel,ndof) % Objetivo: % Especificar los números de los GDL del elemento “iel”.Caso 1D % % Variables % index - vector de los GDL del elemento ‘iel’ % iel – número del elemento % nnel – número de los nodos del elemento % ndof – número de los GDL por nodo % edof – número de los GDL por elemento % start - número del GDL que está inmediatamente “detrás” del primer % GDL del elemento “iel” en la lista global del GDL del sistema. edof = nnel*ndof; start = (iel-1)*(nnel-1)*ndof; % % % % % %

(iel-1): restamos 1 porque si se trata, por ejemplo, del elemento 8, entonces detrás están 7 elementos. (nnel-1): restamos 1 porque en todos los (iel-1) elementos “traseros” el último nodo del elemento anterior es al mismo tiempo el primer nodo del elemento siguiente y se cuenta solamente una vez, como el primero del elemento siguiente.

4. Cálculo de las deformaciones estáticas en sistemas continuos por MEF

39

for i=1:edof index(i)=start+i; end …………………………………………………………………………………………………………. function [kk]=feasmbl1(kk,k,index) % % % % % % %

Objetivo Ensamblar las matrices [k] elementales en globales [kk] Variables kk – matriz global k - matriz elemental index – vector de los números de GDL de un elemento

edof = length(index); for i=1:edof ii=index(i); for j=1:edof jj=index(j);

% cantidad de GDL en el elemento % para todo GDL del elemento % encontramos su fila global % y su columna global

A cada GDL local (1,2,3…i,… edof) le corresponde su GDL global index(i)=ii. A cada elemento de la matriz elemental k(i,j) le corresponden dos GDL locales (el de la fila: “i” y el de la columna: “j”) y, por tanto, dos GDL globales: “ii” y “jj”. Así se determina el lugar del k(i,j) en la matriz global del sistema: fila “ii” y columna “jj”. kk(ii,jj)=kk(ii,jj)+k(i,j); end end ... ... ... ... ... ... ... ... ... ... ... ... ... ...

function [kk,ff]=feaplyc2(kk,ff,bcdof,bcval) Para comprender mejor la actuación de esta función veamos el siguiente comentario. Al aplicar una condición frontera a uno de los GDL, por ejemplo, GDL número “id”, tenemos que sustituir la ecuación número “id” del sistema por una simple igualdad xld = val (siendo “val” el valor impuesto de la coordenada correspondiente. Para hacerlo automáticamente sin cambiar las dimensiones de las matrices podemos realizar lo siguiente: En la matriz de rigidez: igualar a cero todos los componentes de la fila número “id” igualar a 1 el componente diagonal de la fila número “id” En el vector de carga: igualar a “val” el componente número “id”. Recordamos que esta manipulación del vector de carga es independiente de la imposición de las fuerzas externas, que se hace aparte como datos de entrada. %---------------------------------------------------------% Objetivo: % Aplicar las condiciones frontera a la ecuación matricial [kk]{x}={ff} % sin romper la simetría de la matriz global [kk]. % % Variables: % kk – matriz global antes de aplicar CF % ff - matriz global antes de aplicar CF % bcdof - vector con los GDL que se limitan % bcval - vector con los valores de GDL que se imponen

4. Cálculo de las deformaciones estáticas en sistemas continuos por MEF

40

% Por ejemplo, los GDL limitados son 2 y 10, cuyos valores % son, respectivamente, 0.0 y 2.5: % bcdof(1)=2, bcdof(2)=10; bcval(1)=1.0, bcval(2)=2.5. %----------------------------------------------------------n=length(bcdof); sdof=size(kk); for i=1:n c=bcdof(i); % el número de GDL limitado en la CF “i” for j=1:sdof % los componentes de toda la fila “c”de la matriz kk ... kk(c,j)=0; % ...que corresponde al GDL señalado, serán ceros... end kk(c,c)=1; % ff(c)=bcval(i); end

...salvo el de la diagonal % en la misma fila del vector de carga ponemos... % ... el valor exigido por la CF

%----------------------------------------------------------Si además no queremos perder la simetría de la matriz de rigidez, entonces hay que igualar a cero también los componentes de toda la columna “id” (salvo el componente diagonal). Para que esta operación no cambie el sistema tenemos que al mismo tiempo restar de todos los elementos del vector de carga (salvo el número “id”)el producto “val* (elemento de la columna “id” de la matriz de rigidez)” (ver el Apéndice_cond_front). FINAL FUNCIONES INCORPORADAS

4.4

Barra con resorte (método Kwon)

Veamos ahora el mismo problema con una modificación: supongamos que entre el extremo izquierdo de la misma barra y la pared rígida se coloca un resorte con la constante de elasticidad “q” (fig. 4.5). Para introducir este resorte en nuestro esquema del MEF lo podemos ver como un elemento más. Sea el eleme nto número 1. El resorte consideramos coP mo un elemento “concentrado” (igual que en los sistemas discretos), sin dimensión propia. Para evitar un elemento con la lonx=3 x=-10-6 x=0 gitud cero podemos asignarle una longitud Fig. 4.5 e, no nula, pero tan pequeña que no produzca ningún efecto numérico(por ejemplo,10-6). Los cambios que tenemos que introducir en el fichero “bar_con_estat_kw.m” serán: (ver el fichero “bar_cil_res_kw.m”) % Datos de control . . . . . . nel=6; %% número de elementos . . . . . . nnode=7; %% número total de nodos . . . . . . . . % área del extremo gordo A=0.007; % área del extremo fino a=0.007;

4. Cálculo de las deformaciones estáticas en sistemas continuos por MEF

%% Coordenadas de los nodos gcoord(1)=-1e-6; gcoord(2)=0.0; gcoord(3)=0.6; gcoord(4)=1.2; gcoord(5)=1.8; gcoord(6)=2.4; gcoord(7)=3.0; % Conectividad . . . . . . nodes(6,1)=6; nodes(6,2)=7; %% aplicamos la fuerza externa en el GDL número 7 ff(7)=18e3; %

CÁLCULO DE LA MATRIZ GLOBAL

%la colocación de las matrices elementales de los elementos 2 – 7 hacemos % como antes, pero empezando por el iel=2 for iel=2:nel % para todos los elementos salvo el primero . . . . . % la matriz del primer elemento podemos colocar antes o después q=1e10; % elegimos el valor de q de tal manera que la elasticidad del primer ele % mento sea de orden de magnitud de la de los demás: A·E/eleng = 2.45·109 kk(1,1)=kk(1,1)+q; kk(1,2)=kk(1,2)-q; kk(2,1)=kk(2,1)-q; kk(2,2)=kk(2,2)+q; El resultado es: u = 1.0e-004 * 0.0000 0.0180 0.0915 0.1649 0.2384 0.3119 0.3853.

41