FACULTAD DE INGENIERIA Y ARQUITECTURA


UNIVERSIDAD DE EL SALVADOR

MAESTRIA EN DISEÑO ESTRUCTURAL

SEGUNDA TAREA

Asignación de Proyecto: La tarea consiste en escribir un programa de elementos finitos para el analisis de esfuerzos y deformaciones de una estructura bidimensional, Programa escribir usando lenguaje MATLAB

INTRODUCCION:

El objetivo de esta tarea es presentar un modelo de Elementos Finitos que permita observar el comportamiento mecanico (patrón de deformación y esfuerzo) de una placa estructural y a su vez pueda ser el modelo utilizado en el analisis de elementos estructurales como losas, vigas, etc.
Se utilizó un programa Matlab de Elementos para desarrollar el calculo bajo el modelo de elementos finitos y de esta manera presentamos una herramienta que simule la condición de carga debida con respecto al peso propio de este tipo de estructuras

El problema se fundamenta en el enfoque plano debido a que la distribución de carga entre en el cuerpo es inicialmente desconocida.
Una solución típica para este tipo de sistemas requiere demasiadas iteraciones, volviéndose engorroso y tardada su solución, por tal motivo se formula dicho programa tomando en consideración las características abajo detalladas.
Los nodos y elementos del modelo se generaron basandonos en la forma CST.

OBJETIVO GENERAL:
Escribir un programa de elemento finito quepermita realizar analisis de esfuerzos y deformaciones de estructuras bidimensionales

OBJETIVO ESPECIFICO:
Programa en escritura MATLAB para calcular esfuerzos y deformaciones de estructuras bidimensionales.

ALCANCES:
* La discretización del continuo se limitara al uso de elementos triangulares (CST)
* No se considerara el efecto de la temperatura.
* Se contemplaran las fuerzas de tracción, así como también las fuerzas de campo y las puntuales.
* Se trataran los vectores de fuerza.
* Se consideraran las Matrices de rigidez .
* Se consideraran las condiciones de frontera.
* Los Esquemas de definición de los grados de libertad se podran retomar de los fundamentos vistos en cada ponencia diaria.

MANUAL DE USUARIO

Para el uso de este pequeño programa sera necesario que se tengan conocimientos basicos acerca del Método de Elementos Finitos, en relación a la información que debera ser extraída del problema a resolver.

PROBLEMA PROPUESTO DEL LIBRO DE CHANDRUPATLA 5.9 (DEL PARCIAL)

5.9- Resuelva el problema de esfuerzo plano en la figura P5.9 por el método CST

figura P5.9 : placa a resolver

Discretizacion:
Fig 1: coordenadas de los nodos

Para enumerar los desplazamientos Q se adopta la siguiente solucion:

Qy= 2i para un nudo i = 4 tenemos que:
Qy= 2(4) = 8
Qx=2(4) -1 = 7
Qx=2i-1
nudo i

Fig 2: grados de libertad de cada nudo.Luego se procede a extraer los datos para colocarlos en el fichero de entrada “PROB5.9CST.p2d”:
Donde:
NN= Número de nudos
NE= Número de elementos
NM= Número de materiales
NDIM= Número de coordenadas por nudo
NEN= Número de nudos por elemento
NDN= Número de grados de libertad por nudo
lx= longitud de la placa en dirección x
ly= longitud de la placa en dirección y
nelex= Número de triangulos en dirección x
neley= Número de triangulos en dirección y

ND= Número de Grados de Libertad con desplazamientos especificados (apoyos)
NL= Número de cargas aplicadas
NMPC= Número de constricciones multi-punto (no aplica)
FC= Número de fuerzas de cuerpo
TS= Número de tracciones superficiales

Quedando esta parte de la siguiente manera:


Luego se colocan las coordenadas de los nodos X y Y; estos los extraemos de la figura 1
Quedando esta parte de la siguiente manera:

Ahora procedemos a hacer el cuadro de conectividad de los CST; esto lo podemos extraer de la figura 1 o 2 Quedando esta parte de la siguiente manera:

En la primera columna se coloca el numero de elemento
N1, N2 y N3 es el numero de nodo que corresponde a cada vértice del elemento.
Mat es el tipo de material del elemento.
Espesor de cada elemento triangular
TempRise corresponde a la temperatura (en nuestro caso no aplica)

Ahora procedemos a encontrar las condiciones de contorno (los desplazamientos son cero)
Como se puede observar los grados de libertad donde el desplazamiento es cero corresponde alos nodos que se encuentran en la parte empotrada siendo estos los siguientes:

Luego procedemos a colocar las cargas externas aplicadas según el numero de nodo (ver figura 2 carga de 10 kN en el nodo 12 y es negativa por su dirección hacia abajo)

A continuación se colocan las fuerzas de cuerpo y las tracciones superficiales (para este ejemplo se coloca cero debido a que el peso de la placa en estudio es despreciable o no ha sido considerado al igual que las tracciones superficiales)

A continuación se coloca el modulo de Young (E) el alpha que para nuestro caso es cero y el Nu (V) que es 0.33, la fila 34 es la construcción del multipunto.



Habiendo finalizado con la entrada de datos en el fichero “PROB5.9CST.p2d” se procede a la ejecución del programa en MATLAB de la siguiente forma:

Click en run y desplega la siguiente ventana:

Click en aceptar y desplega otra ventana:

Si aparece una ventana donde no se encuentra el archivo de entrada “PROB5.9CST.p2d” se procede a buscarlo en la direccion que a sido guardado; luego de encontrarlo se selecciona y se da click en abrir

Y el programa se ejecuta de forma instantanea y genera un archivo de salida llamado “EJEMPLO DE APLICACION CHANDRUPATLA 5.9(PARCIAL).sal”

Tambien desplega en pantalla los elementos triangulares:

A continuacion se presenta los datos guardados en el archivo “EJEMPLO DE APLICACION CHANDRUPATLA 5.9(PARCIAL).sal”

Discusión:
Este ejemplo se desarrollo discretizando el continuo en 4 triangulosaun siendo las dimensiones de esta placa muy pequeñas (60 mm x 30 mm) esto a manera de poder ver que se puede realizar para varios triangulos CST.
En cuanto a los resultados podemos observar que los desplazamientos en los nodos 1 y 4 son cero esto debido a las condiciones de frontera.
En los desplazamientos de los nodos 2 y 3 estos son negativos tanto en X como en Y esto era de esperarse porque estan siendo comprimidos no asi para los nodos 5 y 6 en la dirección X que presentan desplazamientos positivo y por lógica los desplazamientos en Y que son negativos.

Fig. 3. Bosquejo de como se flexiono la placa en estudio.

En cuanto a los esfuerzos: en los elementos 1 y 2 se observan que se encuentran a compresión, no asi los elementos 3 y 4 que se encuentran en tracción.

En cuanto a las reacciones en los apoyos se observa que la reacción 1 que se encuentra en la dirección x es de 20,000 y en el nudo 7 es -20,000 anulandose así; y quedando en equilibrio estatico en las componentes en X; en cuanto a las componentes en Y (reacción 2 y 8) estas se equilibran con la fuerza externa de 10,000 N.

Este trabajo podría mejorarse si se implementara una rutina que visualizara los desplazamientos en los nodos (por ejemplo la fig 3); también si fuese en Excel ya que es mas manipulable, ademas si se le hiciera un entorno mas vistoso tanto al archivo de entrada como al de salida; en cuanto a los resultados que arroja se debería de tener en cuenta presentar la matriz de rigidez ya ensamblada y lamatriz de rigidez reducida.

CÓDIGO DE PROGRAMA (DATOS DE ENTRADA EN FICHERO “PROB5.9CST.p2d”)

CÓDIGO DE PROGRAMA EN LENGUAJE MATLAB

%Analisis de Estructuras Sólidas bidimensionales utilizando elementos
%tipo CST (triangulos de deformación constante)

%presenta:
%WILSON DANILO CHINCHILLA LOPEZ
%YAMIL ABELARDO YAFFAR UMAÑA
%UNIVERSIDAD DE EL SALVADOR

clear
%% %%%% Se obtiene el nombre del archivo de entrada
button = questdlg('Se debe especificar un archivo de entrada de datos','CST 2D','Aceptar','Cancelar','Aceptar');
if strcmp(button,'Aceptar')
[archentrada,ruta]=uigetfile('*.p2d','Archivo de entrada - CST 2D');
else
archentrada=0;
end
if archentrada==0
disp('No se ha especificado ningún archivo de entrada.');
break
end

%% %%%% SE INICIA LA LECTURA DEL ARCHIVO

%% Lectura de datos generales del problema
fid=fopen(archentrada);
saltar=fgets(fid);
titulo=fgets(fid);
saltar=fgets(fid);
NN=fscanf(fid,'%f',1); %Número de nudos
NE=fscanf(fid,'%f',1); %Número de elementos
NM=fscanf(fid,'%f',1); %Número de materiales
NDIM=fscanf(fid,'%f',1); %Número de coordenadas por nudo
NEN=fscanf(fid,'%f',1); %Número de nudos por elemento
NDN=fscanf(fid,'%f',1); %Número de grados de libertad por nudo
lx=fscanf(fid,'%f',1); %Número de grados de libertad por nudo
ly=fscanf(fid,'%f',1); %Número de grados de libertad por nudo
nelex=fscanf(fid,'%f',1); %Número de grados de libertad por nudoneley=fscanf(fid,'%f',1); %Número de grados de libertad por nudo
NQ=NDN*NN; %Calcula el número total de grados de libertad de la estructura
saltar=fgetl(fid);
saltar=fgetl(fid);
ND=fscanf(fid,'%f',1); %Número de Grados de Libertad con desplazamientos especificados (apoyos)
NL=fscanf(fid,'%f',1); %Número de cargas aplicadas
NMPC=fscanf(fid,'%f',1); %Número de constricciones multi-punto (no aplica)
FC=fscanf(fid,'%f',1); %Número de fuerzas de cuerpo
TS=fscanf(fid,'%f',1); %Número de tracciones superficiales
saltar=fgetl(fid);

%% Lectura de coordenadas de nudos
coordnud=zeros(NN,NDIM);
saltar=fgetl(fid);
for i=1:NN
N=fscanf(fid,'%f',1); %Número de nudos
coordnud(N,1)=fscanf(fid,'%f',1); %Coordenada X de nudo N
coordnud(N,2)=fscanf(fid,'%f',1); %Coordenada Y de nudo N
end
saltar=fgetl(fid);

%% Lectura de conectividad, material, area y temperatura de elementos
conecelem=zeros(NE,NEN);
MAT=zeros(NE,1);
espesor=zeros(NE,1);
DT=zeros(NE,1);

saltar=fgetl(fid);
for i=1:NE
N=fscanf(fid,'%f',1); %Número de elemento
conecelem(N,1)=fscanf(fid,'%f',1); %Número local 1
conecelem(N,2)=fscanf(fid,'%f',1); %Número local 2
conecelem(N,3)=fscanf(fid,'%f',1); %Número local 3
MAT(N)=fscanf(fid,'%f',1); %Tipo de material del elemento
espesor(N)=fscanf(fid,'%f',1); %Espesor del elemento
DT(N)=fscanf(fid,'%f',1); %Temperatura del elemento(no aplica)
end
saltar=fgetl(fid);

%% Lectura de desplazamientos especificados (condiciones frontera)
apoyos=zeros(ND,2);
saltar=fgetl(fid);
for i=1:ND
apoyos(i,1)=fscanf(fid,'%f',1); %Número de Grados de Libertad frontera (apoyo)
apoyos(i,2)=fscanf(fid,'%f',1); %Desplazamiento especificado
end
saltar=fgetl(fid);

%% Lectua de cargas aplicadas a cada DOF
F=zeros(NQ,1);
saltar=fgetl(fid);
for i=1:NL
N=fscanf(fid,'%f',1); %Número de grado de libertad
F(N)=fscanf(fid,'%f',1); %Carga aplicada en el libertad
end
saltar=fgetl(fid);

%% Lectura de Fuerzas de cuerpo
Fcuerpo=zeros(FC,3);
saltar=fgetl(fid);
for i=1:FC
Fcuerpo(i,1)=fscanf(fid,'%f',1); %Número de elemento
Fcuerpo(i,2)=fscanf(fid,'%f',1); %Fuerza de cuerpo X
Fcuerpo(i,3)=fscanf(fid,'%f',1); %Fuerza de cuerpo Y
end
saltar=fgetl(fid);

%% Lectura de Tracciones Superficiales
Tsup=zeros(TS,4);
saltar=fgetl(fid);
for i=1:TS
Tsup(i,1)=fscanf(fid,'%f',1); %Número de elemento
Tsup(i,2)=fscanf(fid,'%f',1); %Nudo inicial de la carga
Tsup(i,3)=fscanf(fid,'%f',1); %Nudo final de la carga
Tsup(i,4)=fscanf(fid,'%f',1); %Tracción superficial en nudo inicial
Tsup(i,5)=fscanf(fid,'%f',1); %Tracción superficial en nudo final
end
saltar=fgetl(fid);

%% Lectura de Propiedades de materiales
PM=zeros(NM,3);
saltar=fgetl(fid);
for i=1:NM
N=fscanf(fid,'%f',1); %Número de materialPM(N,1)=fscanf(fid,'%f',1); %Módulo de elasticidad
PM(N,2)=fscanf(fid,'%f',1); %Alpha del material (no aplica)
PM(N,3)=fscanf(fid,'%f',1); %NU del material
end
saltar=fgetl(fid);

%% Lectura de constricciones multi-punto B1*Qi+B2*Qj=B0
if NMPC>0
BT=zeros(NMPC,3);
MPC=zeros(NMPC,2);
saltar=fgetl(fid);
for i=1:NMPC
BT(i,1)=fscanf(fid,'%f',1);
MPC(i,1)=fscanf(fid,'%f',1);
BT(i,2)=fscanf(fid,'%f',1);
MPC(i,2)=fscanf(fid,'%f',1);
BT(i,3)=fscanf(fid,'%f',1);
end
end

fclose(fid); %Fin de lectura de datos y cierre del archivo de entrada

%% %%%% EVALUACIÓN DEL ANCHO DE BANDA

NBW=0;
for N=1:NE
NMIN=conecelem(N,1);
NMAX=conecelem(N,1);
for j=2:NEN
if NMIN>conecelem(N,j)
NMIN=conecelem(N,j);
end
if NMAX<conecelem(N,j)
NMAX=conecelem(N,j);
end
end
NTMP=NDN*(NMAX-NMIN+1);
if NBW<NTMP
NBW=NTMP;
end
end
for I=1:NMPC
NABS=abs(MPC(I,1)-MPC(I,2))+1;
if NBW<NABS
NBW=NABS;
end
end

%% %%%% ENSAMBLAJE DE MATRIZ DE RIGIDEZ GLOBAL BANDEADA

TLe=zeros(2*NDN);
S=zeros(NQ,NBW+1);
Bs=cell(1,NE);
eos=cell(1,NE);
Ds=cell(1,NE);

for N=1:NE
%Calculo de cosenos directores
I1=conecelem(N,1); %I1 nudo local 1
I2=conecelem(N,2); %I2 nudo local 2
I3=conecelem(N,3); %I3 nudolocal 3
E=PM(MAT(N),1); %Módulo de elasticidad del elemento
nu=PM(MAT(N),3); %nu del elemento
eo=[PM(MAT(N),2)*DT(N);PM(MAT(N),2)*DT(N);0];
eos=eo;
X13=coordnud(I1,1)-coordnud(I3,1);
X23=coordnud(I2,1)-coordnud(I3,1);
X21=coordnud(I2,1)-coordnud(I1,1);
X32=coordnud(I3,1)-coordnud(I2,1);
Y12=coordnud(I1,2)-coordnud(I2,2);
Y13=coordnud(I1,2)-coordnud(I3,2);
Y23=coordnud(I2,2)-coordnud(I3,2);
Y31=coordnud(I3,2)-coordnud(I1,2);

D=E/(1-nu^2)*[1,nu,0;nu,1,0;0,0,(1-nu)/2]; %Matriz de material para esfuerzo plano
Ds=D;
J=[X13,Y13;X23,Y23];
B=1/det(J)*[Y23,0,Y31,0,Y12,0;0,X32,0,X13,0,X21;X32,Y23,X13,Y31,X21,Y12];
Bs=B;
Ae=1/2*det(J);

SE=espesor(N)*Ae*B'*D*B; %Matriz de rigidez del elemento N

%Grados de libertad a los que contribuye el elemento N
gle(1)=NDN*I1-1;
gle(2)=NDN*I1;
gle(3)=NDN*I2-1;
gle(4)=NDN*I2;
gle(5)=NDN*I3-1;
gle(6)=NDN*I3;

%Carga de temperatura en el elemento N (no aplica)
%TLe=espesor(N)*Ae*B'*D*eo;

%Modificación del vector de cargas por la temperatura (no aplica)
%for ii=1:NEN*NDN
% F(gle(ii))=F(gle(ii))+TLe(ii);
%end

%Fuerza de cuerpo en el elemento N
for z=1:FC
if Fcuerpo(z,1)==N
%Vector de fuerzas de cuerpo
fx=Fcuerpo(z,2);
fy=Fcuerpo(z,3);fce=espesor(N)*Ae/3*[fx;fy;fx;fy;fx;fy];
%Modificación de vector de cargas por fuerzas de cuerpo
for ii=1:NEN*NDN
F(gle(ii))=F(gle(ii))+fce(ii);
end
end
end


%Tracción superficial en el elemento N
for z=1:TS
if Tsup(z,1)==N
%Vector de fuerzas de tracción
ni=Tsup(i,2);
nf=Tsup(i,3);
p1=Tsup(i,4);
p2=Tsup(i,5);
Xt1=coordnud(ni,1);
Yt1=coordnud(ni,2);
Xt2=coordnud(nf,1);
Yt2=coordnud(nf,2);
lt12=sqrt((Xt2-Xt1)^2+(Yt2-Yt1)^2);
c=(Yt2-Yt1)/lt12;
s=(Xt1-Xt2)/lt12;
Tx1=-c*p1;
Tx2=-c*p2;
Ty1=-s*p1;
Ty2=-s*p2;
TSe=espesor(N)*lt12/6*[2*Tx1+Tx2;2*Ty1+Ty2;Tx1+2*Tx2;Ty1+2*Ty2];
%Modificación de vector de cargas por fuerzas de tracción
F(2*ni-1)=F(2*ni-1)+TSe(1);
F(2*ni)=F(2*ni)+TSe(2);
F(2*nf-1)=F(2*nf-1)+TSe(3);
F(2*nf)=F(2*nf)+TSe(4);
end
end


%Matriz de rigidez bandeada
for ii=1:NEN*NDN
for jj=1:NEN*NDN
if gle(jj)>=gle(ii)
S(gle(ii),gle(jj)-gle(ii)+1)=S(gle(ii),gle(jj)-gle(ii)+1)+SE(ii,jj);
end
end
end

end

%% %%%% MODIFICACIÓN PARA LAS CONDICIONES FRONTERA

%Valor de C para afectar matriz de rigidez
CNST=max(max(S))*10000;%Modificaciones por desplazamientos impuestos en rigidez y cargas
for i=1:ND
F(apoyos(i,1))=F(apoyos(i,1))+CNST*apoyos(i,2);
S(apoyos(i,1),1)=S(apoyos(i,1),1)+CNST;
end

for I=1:NMPC
%Modificaciones de rigidez por constricciones multi-punto
S(MPC(I,1),MPC(I,1)-MPC(I,1)+1)=S(MPC(I,1),MPC(I,1)-MPC(I,1)+1)+CNST*BT(I,1)^2;
S(MPC(I,1),MPC(I,2)-MPC(I,1)+1)=S(MPC(I,1),MPC(I,2)-MPC(I,1)+1)+CNST*BT(I,1)*BT(I,2);
S(MPC(I,2),MPC(I,2)-MPC(I,2)+1)=S(MPC(I,2),MPC(I,2)-MPC(I,2)+1)+CNST*BT(I,2)^2;

%Modificaciones de cargas por constricciones multi-punto
F(MPC(I,1))=F(MPC(I,1))+CNST*BT(I,3)*BT(I,1);
F(MPC(I,2))=F(MPC(I,2))+CNST*BT(I,3)*BT(I,2);
end

%% %%%% SOLUCIÓN DE ECUACIONES USANDO GAUSS MODIFICADO PARA MATRIZ DE RIGIDEZ BANDEADA

%inicialización de variables de proceso
n=NQ;
Q=zeros(1,n);
salir=0;
g=NBW;
b=F;

% Eliminacion hacia adelante
for k=1:n-1
fin=min(n,k+g);
for i=k+1:fin
%Verificación de singularidad de la matriz de rigidez
if S(k,1)==0
salir=1;
disp('Detención prematura de la solución del sistema');
disp('Sistema mal condicionado. Verifique condiciones frontera.');
break
end
factor=S(k,i-k+1)/S(k,k-k+1);
for j=i:fin
S(i,j-i+1)=S(i,j-i+1)-factor*S(k,j-k+1);
end
b(i)=b(i)-factor*b(k);
end
if salir==1
k=n-1;
end
end

%Verificación si salir es 0(no) secalcula sustitución hacia atras
%pero si salir es 1(si) no se realizan dichos calculos
if salir==0
% Sustitución hacia atras
Q(n)=b(n)/S(n,n-n+1);
for i=n-1:-1:1
fin=min(n,i+g);
sum=0;
for j=i+1:fin
sum=sum+S(i,j-i+1)*Q(j);
end
Q(i)=(b(i)-sum)/S(i,i-i+1);
end
end

%% %%%% CALCULO DE ESFUERZOS Y FUERZAS INTERNAS EN ELEMENTOS

%esfuerzos=zeros(NE,3);
%fuerzas=zeros(NE,2);
esfuerzos=cell(1,NE);
for I=1:NE
I1=conecelem(I,1); %Nudo local 1
I2=conecelem(I,2); %Nudo local 2
I3=conecelem(I,3); %Nudo local 3
q=[Q(NDN*I1-1);Q(NDN*I1);Q(NDN*I2-1);Q(NDN*I2);Q(NDN*I3-1);Q(NDN*I3)];
%esfuerzos(I,1)=I;
esfuerzos=(Ds*Bs)*q-Ds*eos;
%esfuerzos(I,2)=PM(me,1)/le*lm*q-PM(me,1)*PM(me,2)*DT(I);
%fuerzas(I,1)=I;
%fuerzas(I,2)=esfuerzos(I,2)*espesor(I);
end

%% %%%% CALCULO DE REACCIONES

reacciones=zeros(ND+NMPC*NDN,2);

%Calculo de reacciones en apoyos especificados
for i=1:ND
N=apoyos(i,1);
reacciones(i,1)=N;
reacciones(i,2)=CNST*(apoyos(i,2)-Q(N));
end

%Calculo de reacciones en constricciones multi-punto
for i=1:2:NMPC
reacciones(ND+i,1)=MPC(i,1);
reacciones(ND+i,2)=-CNST*BT(i,1)*(BT(i,1)*Q(MPC(i,1))+BT(i,2)*Q(MPC(i,2)-BT(i,3)));
reacciones(ND+i+1,1)=MPC(i,2);
reacciones(ND+i+1,2)=-CNST*BT(i,2)*(BT(i,1)*Q(MPC(i,1))+BT(i,2)*Q(MPC(i,2)-BT(i,3)));
end

%% %%%% SALIDA DE RESULTADOS

%Guardandoresultados en archivo de salida
archsalida=strcat(ruta,titulo,'.sal');
fid=fopen(archsalida,'w');
fprintf(fid,' rn');
fprintf(fid,'Salida de Datos - %srn',titulo);
fprintf(fid,'nttt ‘ ; ----‘-------‘-------‘-------‘-------‘------------ ');
fprintf(fid,'nttt ‘;| UNIVERSIDAD DE EL SALVADOR |')
fprintf(fid,'nttt ‘;');
fprintf(fid,'nttt ‘;| MAESTRIA EN INGENIERIA ESTRUCTURAL |')
fprintf(fid,'nttt ‘;| METODO DE ELEMENTOS FINITOS |')
fprintf(fid,'nttt ‘;');
fprintf(fid,'nttt ‘;| PROGRAMA PARA RESOLVER: |')
fprintf(fid,'nttt ‘;| ESTRUCTURAS SOLIDAS BIDIMENSIONALES |')
fprintf(fid,'nttt ‘;| UTILIZANDO ELEMENTOS TIPO CST (TRIANGULOS |')
fprintf(fid,'nttt ‘;| DE DEFORMACION CONSTANTE |')
fprintf(fid,'nttt ‘;| YAMIL ABELARDO YAFFAR U. |')
fprintf(fid,'nttt ‘;| WILSON DANILO CHINCHILLA L. |')
fprintf(fid,'nttt ‘;');
fprintf(fid,'nttt ‘;| JULIO DE 2012 |');
fprintf(fid,'nttt ‘;-----‘-------‘-------‘--------‘-------‘-----------');
fprintf(fid,'nnn');
fprintf(fid,'Analisis de Estructuras Sólidas bidimensionalesrn');
fprintf(fid,'utilizando elementos tipo CSTrn');
fprintf(fid,'(triangulos dedeformación constante)rn');
fprintf(fid,'DESPLAZAMIENTOS EN LOS NUDOSrn');
fprintf(fid,' Nudo X Yrn');
for I=1:NN
fprintf(fid,' %g %f %frn',I,Q(I*2-1),Q(I*2));
end

fprintf(fid,' rn');
fprintf(fid,'ESFUERZOS EN ELEMENTOSrn');
fprintf(fid,'Elemento Efuerzo X Esfuerzo Y Esfuerzo XYrn');
for i=1:NE
fprintf(fid,' %g %10.3f %10.3f %10.3frn',i,esfuerzos);
end

fprintf(fid,' rn');
fprintf(fid,'REACCIONES EN APOYOSrn');
fprintf(fid,' DOF# Reacciónrn');
for i=1:ND+NMPC*NDN
fprintf(fid,' %g %10.3frn',reacciones(i,1),reacciones(i,2));
end

fclose(fid); %Cierre del archivo de salida

%Desplegando resultados en pantalla
type(archsalida);

fprintf('Resultados guardados en %sn',archsalida);
fprintf('n');
fprintf('%s nn', '==========GENERACION DEL TRIANGULO========');
%Calculations
%No. of elements and increment in x & y direction
nnodex = nelex + 1;
nnodey = neley + 1;
dx = lx / nelex;
dy = ly / neley;
% Generate Coordinates
nelexy = nelex * neley*2; % total no. of Triangular elements
nnodexy = nnodex * nnodey; % total no. of nodes
% zero matrices and vectors
xc=zeros(nnodexy,1);
yc=zeros(nnodexy,1);
nodedata=zeros(nnodexy,3);
for row = 1: nnodey
for col = 1: nnodex
nt = (row-1)* nnodex + col;
xc(nt)= (col-1)*dx;
yc(nt)= (row-1)*dy;
end
end
% combining node data into a matrix nodedata=[no., xc, yc]
% and elementdata=[no., 3 nodes]
format shortfor i= 1: nnodexy
nodedata(i,1)= i;
nodedata(i,2)= xc(i);
nodedata(i,3)= yc(i);
end
for elerow = 1: neley
for elecol = 1: nelex
eleno = (elerow-1)* (nnodex-1) + elecol;
node1 = eleno + (elerow-1);
node2 = eleno + (elerow);
node3 = node2 + nnodex;
node4 = node1 + nnodex;
elementdata(eleno,1)=eleno;
elementdata(eleno,2)=node1;
elementdata(eleno,3)=node2;
elementdata(eleno,4)=node3;
elementdata(eleno+nelex * neley,1)=eleno+nelex * neley;
elementdata(eleno+nelex * neley,2)=node1;
elementdata(eleno+nelex * neley,3)=node3;
elementdata(eleno+nelex * neley,4)=node4;
end
end
% Plot of the structure
figure (1)
for n=1:nelexy
x(1)= xc(elementdata(n,2));
x(2)= xc(elementdata(n,3));
x(3)= xc(elementdata(n,4));
x(4)= x(1);
y(1)= yc(elementdata(n,2));
y(2)= yc(elementdata(n,3));
y(3)= yc(elementdata(n,4));
y(4)= y(1);
plot (x, y, '-o', 'linewidth', 0.05, 'markersize', 1.5,'markerfacecolor', 'b')
axis equal
hold on
end
leng=length(xc);
for i=1:leng
istr=num2str(i);
text(xc(i)+0.1*dx, yc(i)+0.1*dy, istr )
end
for i=1:length (elementdata)
elementno=elementdata(i,1);
nd1=elementdata((i),2);
x1=xc(nd1);
y1=yc(nd1);
nd2=elementdata((i),3);
x2=xc(nd2);
y2=yc(nd2);
nd3=elementdata((i),4);
x3=xc(nd3);
y3=yc(nd3);
xx= [x1; x2; x3];
yy= [y1; y2; y3];
xaverage=mean(xx);
yaverage=mean(yy);
elnStr = num2str(elementno);
text(xaverage, yaverage, elnStr )
end
nodedata
elementdata
NoOfElements = nelexy
NoOfNodes= nnodexy
%%%%%%% FIN DEL PROGRAMA