Autor Tema: MEF - sistema acoplado

0 Usuarios y 1 Visitante están viendo este tema.

23 Octubre, 2019, 11:45 am
Leído 501 veces

adhemir

  • Aprendiz
  • Mensajes: 283
  • Karma: +0/-0
  • Sexo: Masculino
Hola,  estudiando por youtube me encontre com el siguiente codigo em matlab ,  sigue el hilo
https://www.youtube.com/watch?v=44vDk7VeoNU&list=PLccOt6x5p76ch1l5vydzhG5NgR0-DCMkk&index=8
solo que no compila pues el programa esta imcompleto ya que falta definir en la funcion linearelement los parametros
 k311, k312,k312, k313;k314,f11 e f12 , ??? ??? ??? ???  alguna ayuda como definir esos parametros.

23 Octubre, 2019, 12:53 pm
Respuesta #1

Luis Fuentes

  • el_manco
  • Administrador
  • Mensajes: 47,035
  • País: es
  • Karma: +1/-0
  • Sexo: Masculino
Hola

Hola,  estudiando por youtube me encontre com el siguiente codigo em matlab ,  sigue el hilo
https://www.youtube.com/watch?v=44vDk7VeoNU&list=PLccOt6x5p76ch1l5vydzhG5NgR0-DCMkk&index=8
solo que no compila pues el programa esta imcompleto ya que falta definir en la funcion linearelement los parametros
 k311, k312,k312, k313;k314,f11 e f12 , ??? ??? ??? ???  alguna ayuda como definir esos parametros.


Sería deseable que planteases la cuestión de una forma más manejable que no exija "capturar" un código de un programa  de un video.

Saludos.

23 Octubre, 2019, 02:40 pm
Respuesta #2

adhemir

  • Aprendiz
  • Mensajes: 283
  • Karma: +0/-0
  • Sexo: Masculino
Hola, es verdad tienes razón,  el código es el siguiente,


function FEM_code__Coupled()
clear all;close all; clc
n = 5;
nn = n+1;
lgth = 1;
he = lgth/n;
x = (0:he:lgth);
AC = 0.00005;
U= zeros(nn,1);
V=zeros(nn,1);
U(1)=0 ; V(1)=1;
U(nn)=1; V(nn)=exp(1);
c = 1.0;
count = 0;
tic
     while(c>0)
         
         [U1,V1] = assembly(U,V,n,he);
         c =0.0;
         for i=1:nn
             if (abs(U(i)-U1(i))>AC)
                 c= c+1;
                 break;
             end
             if (abs(V(i)-V1(i))>AC)
                 c= c+1;
                 break ;
             end
         end
         U =U1;
         V = V1;
         
         count = count +1;
     end
     
     disp('hence solution:');
     
     % saida de variaveis primaria e sequandarias.
     diff = abs(U-(x.^2)');
     diff1 = abs(V-exp(x)');
     
     fprintf(' No de elementos = %d\n',n);
     disp('     x       u(FEM)       u(Exact)    error    v(FEM)      v(Exact)     Error')
     disp([x',U,(x.^2)',diff,V,exp(x)',diff1])
     fprintf('No de iteracoes=%d\n ',count);
     
     %grafica variaveis primarias
     figure(1)
     plot(x,U,'--rs','Linewidth',2)
     xlabel('x');
     ylabel('u(x)');
     title('Solucao grafica para o pvc')
     figure(2)
     plot(x,V,'--rs','Linewidth',2)
     xlabel('x')
     ylabel('v(x)')
     title('Solucao grafica para o PVC')
     toc
     
end

%derivacao de elemntos da matriz e emsamble
function [U1,V1] = assembly(U,V,n,he)
nn=n+1;
L11=zeros(nn,nn); L12=zeros(nn,nn); R1=zeros(nn,1);
L21=zeros(nn,nn); L22=zeros(nn,nn); R2=zeros(nn,1);

lmm=[];
for i=1:n
    lmm=[lmm;[i,i+1]];
end
 for i=1:n
     lm=lmm(i,:);
     
  %geracao de elemnto da matriz k11 e RMS matriz (f1)
    % Integracao usando quadratura gaussiana
   
[k111,k112,k113,k311,k312,k313,k314,f11,f12] = linearelement(he,i);

% Para primeira equacao
k11 = -k111;
k12 = k313*V(lm(1))+k314*V(lm(2))+(k311*U(lm(1))+k312*U(lm(2)));

% para sequnda equcao
k21 = k112;
k22 = -k111-k113;
f1 = f11;
f2 = f12;

% matriz de conetividade assembly
L11(lm,lm) = L11(lm,lm)+k11;
L12(lm,lm) = L12(lm,lm)+k12;
L21(lm,lm) = L21(lm,lm)+k21;
L22(lm,lm) = L22(lm,lm)+k22;
R1(lm) = R1(lm)+f1;
R2(lm) = R2(lm) +f2;

end

% imposicao das condicoes de fronteira

L11(1,:)=0.0;   L12(1,:)=0.0;
L11(nn,:)=0.0;  L12(nn,:)=0.0;
L21(1,:)=0.0  ;  L22(1,:)=0.0;
L21(n,:)=0.0  ; L22(n,:)=0.0;
L11(1,1)=1;   L11(nn,nn)=1;
L22(1,1)=1;   L22(nn,nn)=1;


R1(1,1)=U(1);   R1(nn,1)=U(nn);
R2(1,1)=V(1);   R2(nn,1)=V(nn);


% solucao do sistema de equacoes F1
K= [L11,L12;L21,L22];

R=[R1;R2];
d=K\R;
U1= d(1:nn);
V1= d(nn+1:2*nn);

end



function [k111,k112,k113,k311,k312,k313,k314,f11,f12] = linearelement(he,i);
  gL = [-sqrt(3/5);0;sqrt(3/5)];
  gW = [5/9,8/9,5/9];
 
 
 
 
k111=zeros(2); k112=zeros(2); k113=zeros(2); k221=zeros(2); k222=zeros(2);k311=zeros(2); k312=zeros(2);
k313=zeros(2); k314=zeros(2); f11=zeros(2,1); f12=zeros(2,1);

for k=1:length(gW)
    s = gL(k); w = gW(k);
    n = [(1-s)/2, (1+s)/2];
    dns = [-1/2, 1/2];
    coord = [(i-1)*he;i*he];
    x = n*coord(:,1);
   
    J=he/2;
    dx=(1/J)*dns;
   
 k111=k111+J*w*dx'*dx;
  k112=k112+J*w*n'*dx;
  k113=k113+J*w*n'*n;
   
  k221=k221+J*w*n'*dx*n(1);
   k222=k222+J*w*n'*dx*n(2);
   

% Cuando digito los siguientes parametros me sale matriz no singular

%    % falta definir
%    k311=k311+J*w*n'*dx*n(1);
%    k312=k312+J*w*n'*dx*n(2);
%    k313=k313+J*w*n'*dx*n(1);
%    k314=k314+J*w*n'*dx*n(2);
%   
%    f11=f11+J*w*n'*(x.^2)';
%    f12=f12+J*w*n'*exp(x)';
%   
   
   
end
end