Luis Carlos Muñiz MenesesInstituto Politécnico Nacional
Programas Scilab Dos de los programas son de cálculos estructurales y el ultimo es del ciclo bryton Método de los 3 momentos (basado en el método de superposición) Este programa te calcula las reacciones y los momentos intermedios así como el área y el centroide de cada claro function pb(nc) //Se declaro una funcion en la cual se le va a indicar unicamente el numero de claros presente en toda la viga // Ciclo para determinar AREA, CENTROIDE, REACCIONES de cada claro arm=nc+1 for cc=1:nc printf("\n\nCual es la longitud del claro No %2.0f",cc); //Se pide al la longitud del claro a estudiar lc(cc)=input("Longitud=") printf("Cuantas cargas puntuales deceas en el claro No %2.0f",cc); //Aqui se ingresa el numero de cargas puntuales presentes en el claro p=input("Cargas puntuales=") byt=0 ayt=0 //Se inicializan determinado numero de variables para poder aplicar el metodo de superposicion att=0 xt=0 suax=0 pAX=0 for w=1:p // Se implementa otro ciclo para calcular y sumar las REACCIONES, AREAS y MOMENTOS que genera cada carga puntual en el claro printf("Cual es el valor de la carga No %2.0f",w) //Se le pide al que ingrese el valor especifico de X carga qp(w)=input("Valor=") if qp(w)==0 then //Se genera una condicionante en dado caso que quiera un claro con ninguna carga puntual a la hora de realizar by(w)=0 //Calculos no marque error por diviciones entre cero ay(w)=0 At(w)=0 X(w)=0 sax(w)=0 else printf("A que distancia se encuentra la carga puntual No %2.0f",w) //Si dio un valor distinto a cero se le pide al la distancia dp(w)=input("Distancia=") //a la que se encuentra dicha carga aplicando otra condicionante en dado caso de que la carga se if dp(w)>lc(cc) then //encuentre fuera del claro indicarle que esta equivocado printf("Hay una incongruencia la distancia se encuentra fuera del claro")
Luis Carlos Muñiz MenesesInstituto Politécnico Nacional il_fi ... Programas Scilab else by(w)=[qp(w)*dp(w)]/lc(cc) ay(w)=qp(w)-by(w); //Ecuaciones utilizadas para calcular AREAS, CENTROIDES Y REACCIONES Ai(w)=(((qp(w)-((qp(w)*dp(w))/lc(cc)))*dp(w))*dp(w))/2; //del claro Ad(w)=((((qp(w)-((qp(w)*dp(w))/lc(cc)))*dp(w))*(lc(cc)-dp(w)))/2); At(w)=Ai(w)+Ad(w); X(w)=(Ad(w)*2/3*(lc(cc)-dp(w))+Ai(w)*((lc(cc)-dp(w))+dp(w)/3))/At(w); //Metodo de superposicion para cargas puntuales byt=byt+by(w) ayt=ayt+ay(w) att=att+At(w) sax(w)=At(w)*X(w) suax=suax+sax(w) end end end LC(arm)=0 //Despues de las cargas puntuales se le indica al que si decea agregar una carga distribuida a todo el claro qdis=input("Deceas poner una carga distribuida por todo el claro Si[1], No[2]") if qdis==1 then //Se genera una condicionante para implementar o hacer caso printf("Cual es el valor de la carga distribuida del claro No %2.0f",cc) //omiso de la carga distribuida, en dado caso de decear la qd=input("Valor=") //carga se le indica al que ingrese el valor de la By=[qd*lc(cc)]/2 //carga para despues calcular de igual forma que en las cargas Ay=[qd*lc(cc)]-By //puntuales AREAS, CENTROIDES, REACCIONES e implementar at=2*[(2/3)*(lc(cc)/2)*((qd*((lc(cc))^2))/8)] //nuevamente el metodo de superpocision y obtener los valores Xo=lc(cc)/2 //reales del claro byt=byt+By ayt=ayt+Ay att=att+at pAX=at*Xo end if att==0 then //se genera un condicionante en dado caso de que el centroide total del claro sea igual a cero xt=0 //ya que se produciria una indeterminacion else xt=(suax+pAX)/att
Luis Carlos Muñiz MenesesInstituto Politécnico Nacional il_fi ... Programas Scilab end LC(cc)=lc(cc) printf("By total:%f",byt) printf("\nAy total:%f",ayt) printf("\nArea total:%f",att) printf("\nLocalizacion del centroide:%f\n",xt) ATT(cc)=att XT(cc)=xt Ayr(cc)=ayt Byr(cc)=byt end //Arreglo matricial de los 3 momentos Arm=(nc-1) //Se genera una matriz de ceros para guardar los escalares resultantes de las ecuaciones de los 3 momentos for t=1:Arm for y=1:arm x(t,y)=0 end end //Mediante este ciclo se genera la matriz de los escalares gracias a la ecuacion de los 3 momentos for i=1:Arm x(i,i)=lc(i) x(i,i+1)=2*(lc(i)+LC(i+1)) x(i,i+2)=LC(i+1) end //Se arregla a una matriz cuadrada la generada por las ecuaciones de los 3 momentos X=x(1:Arm,2:nc) //Se implementa un vector resultados para calcular los momentos distintos a cero m=zeros(Arm,1) for n=1:Arm m(n)=-(6*((ATT(n)*XT(n))/lc(n)))-6*((ATT(n+1)*XT(n+1))/lc(n+1)) end; //Se calculan los Momentos que son distintos a cero gracias al producto de la transpuesta de la matriz escalar producto con el //vector de los resultados M=inv(X)*m // Se agregan los momentos que se eliminaron para hacer el arreglo maticial Mo=zeros(arm,1) Mo(1,1)=0 Mo(arm,1) for z=2:nc Mo(z,1)=M(z-1,1) end
Luis Carlos Muñiz MenesesInstituto Politécnico Nacional il_fi ... Programas Scilab //cambio el verctor columna a vector fila mo=Mo' //Se implementa un siclo para obtener las cargas reales claro por claro para despues poder sumar todos los claros y convertirla en una viga continua ayr=zeros(1,nc) byr=zeros(1,nc) //Ayr,Byr Guardan las reacciones que generan todas las cargas en los apoyos para despues sumarle las reacciones que generan los momentos //en los apoyos de cada claro for o=1:nc ayr(o)=Ayr(o)+((-mo(o)+mo(o+1))/lc(o)) byr(o)=Byr(o)+((mo(o)-mo(o+1))/lc(o)) end //reacciones totales //Se imprimen los valores de las reacciones en los apoyos inicial y final ademas de implementar un ciclo para sumar //las reacciones de los claros que se encuentran en medio de la viga printf("\nApoyo inicial:%f",ayr(1)) for g=2:nc rreales(g)=byr(g-1)+ayr(g) qw=g-1 printf("\nApoyo%2.0f:%f\n",g,rreales(g)) end printf("Apoyo final:%f",byr(nc)) //Matriz inicial donde se obtiene los escalares de la ecuacion de los 3 momentos printf("\nMatriz escalar de los 3 momentos") disp (x) printf("Vector resultante") disp (m) printf("Matriz cuadrada de los 3 momentos eliminando Momento inicial y Momento final") disp (X) //disp (M) //Vector de los valores de los momentos sin los momentos inciales y finales //disp (Ayr) //Vector donde se guardan las reacciones Ay de cada claro //disp (Byr) //Vector donde se guardan las reacciones By de cada claro printf("Verctor de los momentos presentes en los apoyos recordando que el momento incial y final es cero") disp (mo) endfunction //El programa fue realizado/calibrado mediante el ejemplo de esta pagina http://www.docstoc.com/docs/17222644/Teorema-de-los-Tres-Momentos
Este es un método matricial para una viga (rigidez) sinceramente no recuerdo el tema como tal pero te pide que ingreses el número de elementos (e), número de nodos (n) y el número de variables conocidas (vc) cabe mencionar que las
Luis Carlos Muñiz MenesesInstituto Politécnico Nacional il_fi ... Programas Scilab variables conocidas son las condiciones de frontera las cuales te ayudaran a resolver el problema donde el vector de fuerzas es igual a la matriz de rigidez multiplicada por el vector de desplazamientos si mal no recuerdo. function mv(e, n, vc) ma=n*3; a=ma-vc; MAG=zeros(ma,ma) MFI=zeros(vc,vc) MF=zeros(vc,vc) MU=zeros(a,vc) for i=1:e Lx=0; Ly=0; A=0; E=0; I=0; fx1=input("numero de la variable fx1:"); fy1=input("numero de la variable fy1:"); M1=input("numero de la variable M1:"); fx2=input("numero de la variable fx2:"); fy2=input("numero de la variable fy2:"); M2=input("numero de la variable M2:"); x1=input("Ingresa la coordenada X del nodo 1:"); y1=input("Ingresa la coordenada Y del nodo 1:"); x2=input("Ingresa la coordenada X del nodo 2:"); y2=input("Ingresa la coordenada Y del nodo 2:"); L=sqrt(((x2-x1)^2)+((y2-y1)^2)) disp("La longitud de tu viga es de: "+string(L)); A=input("Area transversal del elemento (m^2):"); E=input("Modulo elastico del elemento:") I=input("Momento de inercia del elemento (m^4):") Lx=(x2-x1)/L; Ly=(y2-y1)/L; MAL=zeros(ma,ma) MAL(fx1,fx1)=((((A*E)/L)*Lx^2)+(((12*E*I)/L^3)*Ly^2)); MAL(fx1,fy1)=(((A*E)/L)-((12*E*I)/L^3))*Lx*Ly; MAL(fy1,fx1)=MAL(fx1,fy1); MAL(fx1,M1)=-((6*E*I)/L^2)*Ly; MAL(M1,fx1)=MAL(fx1,M1); MAL(fx1,fx2)=-MAL(fx1,fx1); MAL(fx2,fx1)=MAL(fx1,fx2); MAL(fx1,fy2)=-MAL(fx1,fy1); MAL(fy2,fx1)=MAL(fx1,fy2); MAL(fx1,M2)=-((6*E*I)/L^2)*Ly; MAL(M2,fx1)=MAL(fx1,M2); MAL(fy1,fy1)=((((A*E)/L)*Ly^2)+(((12*E*I)/L^3)*Lx^2));
Luis Carlos Muñiz MenesesInstituto Politécnico Nacional il_fi ... Programas Scilab MAL(fy1,M1)=((6*E*I)/L^2)*Lx; MAL(M1,fy1)=MAL(fy1,M1); MAL(fy1,fx2)=-(((A*E)/L)-((12*E*I)/L^3))*Lx*Ly; MAL(fx2,fy1)=MAL(fy1,fx2); MAL(fy1,fy2)=-((((A*E)/L)*Ly^2)+(((12*E*I)/L^3)*Lx^2)); MAL(fy2,fy1)=MAL(fy1,fy2); MAL(fy1,M2)=((6*E*I)/L^2)*Lx; MAL(M2,fy1)=MAL(fy1,M2); MAL(M1,M1)=((4*E*I)/L); MAL(M1,fx2)=((6*E*I)/L^2)*Ly; MAL(fx2,M1)=MAL(M1,fx2); MAL(M1,fy2)=-((6*E*I)/L^2)*Lx; MAL(fy2,M1)=MAL(M1,fy2); MAL(M1,M2)=((2*E*I)/L); MAL(M2,M1)=MAL(M1,M2); MAL(fx2,fx2)=MAL(fx1,fx1); MAL(fx2,fy2)=MAL(fx1,fy1); MAL(fy2,fx2)=MAL(fx2,fy2); MAL(fx2,M2)=((6*E*I)/L^2)*Ly; MAL(M2,fx2)=MAL(fx2,M2); MAL(fy2,fy2)=MAL(fy1,fy1); MAL(fy2,M2)=-((6*E*I)/L^2)*Lx; MAL(M2,fy2)=MAL(fy2,M2); MAL(M2,M2)=((4*E*I)/L); MAG=MAG+MAL disp("Lx:"+string(Lx)); disp("Ly:"+string(Ly)); disp(MAL); end VF=zeros(vc,1) for j=1:vc VF(j,1)=input("valores de fuerzas conocidas:"); end for k=1:vc for l=1:vc MF(k,l)=MAG(k,l) end end MFI=inv(MF) U=MFI*VF; for q=vc+1:ma for w=1:vc MU(q,w)=MAG(q,w)
Luis Carlos Muñiz MenesesInstituto Politécnico Nacional il_fi ... Programas Scilab end end F=MU*U; disp("Matriz Global:") disp(MAG) disp("Matriz inversa para el calculo de Desplazamientos") disp(MFI) disp("Matriz para el calculo de Fuerzas") disp(MU) disp("Desplazamientos: "+string(U)) disp("Fuerzas: "+string(F)) disp("SON FUERZAS SIN NINGUNA CARGA APLICADA SOLO FALTA AGREGARLE LAS CARGAS DEL VECTOR DE FUERZAS PLANTEADO AL PRINCIPIO (TU LIBRETA)") endfunction por ultimo aquí les dejo el programa del calculo termodinámico de un motor a reacción (ciclo bryton) las ecuaciones esta calibradas de 0 a 1023 ya que lo implementamos en un pic el cual tiene un adc de 10 bits, esta calibrado hasta 1300 K arrojando valores muy cercanos, existe una pequeña variación simplemente por las ecuaciones, ya que no describían las tablas termodinámicas del aire como un gas ideal al 100%, se supone que el valor máximo de temperatura es de 2250 k pero no estoy seguro de que tan exacto sea después de los 1300 K. function raiz(T1, rc, T3) //agregamos temperaturas en °C y la relacion de compresion como tal //Proceso 1 Ta1=T1+273.15 //convertimos los °C a K x1=(Ta1-200.000000000005)/2.00391006940562 h1=-6.74350407421917E-16*x1^6+2.45955412273339E-12*x1^53.44913088021204E-09*x1^4+0.0000021376018031205*x1^30.000179620062453978*x1^2+2.01421105823829*x1+199.910974626429 //Correccion presion relativa eficas para el rango de temperatura 258-488 if T1>=258 & T1<488 then prr1=1.0772706884655E-17*x1^6+7.65512300901147E13*x1^5+1.81621872923088E-09*x1^4+3.42023787069934E07*x1^3+0.000224193968733744*x1^2+0.00543167471187189*x1+0.4632799 08988625 r1=0.0000261978266133944*prr1^6 - 0.000769635601555637*prr1^5 + 0.00907080763507707*prr1^4 - 0.0550498738107573*prr1^3 + 0.182867580946214*prr1^2 - 0.309056603839044*prr1 + 0.152984217551028 pr1=prr1-r1 else pr1=1.0772706884655E-17*x1^6+7.65512300901147E13*x1^5+1.81621872923088E-09*x1^4+3.42023787069934E-
Luis Carlos Muñiz MenesesInstituto Politécnico Nacional il_fi ... Programas Scilab 07*x1^3+0.000224193968733744*x1^2+0.00543167471187189*x1+0.4632799 08988625 end //Proceso 2 pr2=rc*pr1 //Correccion de temperatura a la salida del compresor Ta2=287.318709393591*pr2^0.257288147654972 if Ta2<=464 then Td2=0.0000000000000298293123522627*Ta2^60.0000000000533987791502879*Ta2^5+0.0000000379673625363624*Ta2^40.0000131563359764458*Ta2^3+0.00197399763062981*Ta2^20.067985741104833*Ta2+11.7728341341899 else Td2=0.000000000000000031688404545543*Ta2^6+0.00000000000026515619372586 1*Ta2^50.000000000857868534910672*Ta2^4+0.00000138006438858453*Ta2^30.00107420057613197*Ta2^2+0.303620125127829*Ta2-13.1998416304612 end T2=Ta2-Td2 x2=(T2-200.000000000005)/2.00391006940562 h2=-6.74350407421917E-16*x2^6+2.45955412273339E-12*x2^53.44913088021204E-09*x2^4+ 0.0000021376018031205*x2^30.000179620062453978*x2^2+2.01421105823829*x2+199.910974626429 //Proceso 3 Ta3=T3+273.15 //convertimos los °C a K x3=(Ta3-200.000000000005)/2.00391006940562 h3=-6.74350407421917E-16*x3^6+2.45955412273339E-12*x3^53.44913088021204E-09*x3^4+0.0000021376018031205*x3^30.000179620062453978*x3^2+2.01421105823829*x3+199.910974626429 pr3=1.0772706884655E-17*x3^6+7.65512300901147E13*x3^5+1.81621872923088E-09*x3^4+3.42023787069934E07*x3^3+0.000224193968733744*x3^2+0.00543167471187189*x3+0.4632799 08988625 //Proceso 4 pr4=pr3/rc //Correccion de temperatura a la salida de la turbina Ta4=287.318709393591*pr4^0.257288147654972 Td4=0.000000000000000031688404545543*Ta4^6+0.00000000000026515619372586 1*Ta4^50.000000000857868534910672*Ta4^4+0.00000138006438858453*Ta4^30.00107420057613197*Ta4^2+0.303620125127829*Ta4-13.1998416304612 T4=Ta4-Td4 x4=(T4-200.000000000005)/2.00391006940562
Luis Carlos Muñiz MenesesInstituto Politécnico Nacional il_fi ... Programas Scilab h4=-6.74350407421917E-16*x4^6+2.45955412273339E-12*x4^53.44913088021204E-09*x4^4+0.0000021376018031205*x4^30.000179620062453978*x4^2+2.01421105823829*x4+199.910974626429 //Trabajo de entrado wc=h2-h1 //Trabajo de salida wt=h3-h4 //Trabajo neto wn=wt-wc //Calor de entrada qe=h3-h2 //Eficiencia termica n=(wn/qe)*100 printf("Temperatura 1:") disp(T1) printf("Entalpia 1:") disp(h1) printf("Presion relativa 1:") disp(pr1) printf("Temperatura 2:") disp(T2) printf("Entalpia 2:") disp(h2) printf("Presion relativa 2:") disp(pr2) printf("Temperatura 3:") disp(T3) printf("Entalpia 3:") disp(h3) printf("Presion relativa 3:") disp(pr3) printf("Temperatura 4:") disp(T4) printf("Entalpia 4:") disp(h4) printf("Presion relativa 4:") disp(pr4) printf("Trabajo de entrada") disp(wc) printf("Trabajo de salida") disp(wt) printf("Trabajo neto") disp(wn) printf("Calor de entrada") disp(qe) printf("Eficiencia termica") disp(n) endfunction
Luis Carlos Muñiz MenesesInstituto Politécnico Nacional il_fi ... Programas Scilab