• jon98
    了解作者
  • matlab
    开发工具
  • 21KB
    文件大小
  • rar
    文件格式
  • 0
    收藏次数
  • 1 积分
    下载积分
  • 6
    下载次数
  • 2018-09-03 20:46
    上传日期
Finite element model for poroelastic problems.
Poroelasticity.rar
  • FintQ4.m
    2.1KB
  • DataMesh.m
    1.7KB
  • Parabolic_solver_vForm.m
    1.8KB
  • mainroutine.m
    9.2KB
  • PRQ4.m
    2KB
  • MatGlob.m
    2.8KB
  • KeQ4.m
    3.2KB
  • DEFQ4.m
    2.1KB
  • drawingMesh.p
    4.1KB
  • drawingField.p
    4KB
  • BCDirichlet.m
    1.2KB
  • DataMaterial.m
    752B
  • calc_pres.m
    1KB
  • Parabolic_solver_dForm.m
    1.8KB
  • ElastMESHQ4.m
    536B
  • Parabolic_solver_vForm_Spilker.m
    758B
  • DataBC.m
    1.8KB
  • rectangularMesh.p
    1.9KB
内容介绍
%************************************************************************* % ************************************************************************* % Programa de Elementos Finitos para Poroviscoelasticidade % ************************************************************************* %------------------------------------------------------------------------- % Limpeza da tela e da memoria antes da inicialização %------------------------------------------------------------------------- clear all close all clc; format long; % Precisao de 16 digitos %------------------------------------------------------------------------- % Leitura dos Dados DataMesh; DataMaterial; DataBC; %------------------------------------------------------------------------- % Processamento %------------------------------------------------------------------------- % Gera plot da malha ( acho que vou usar o drwing mesh no lugar) Mesh = ElastMESHQ4(NumNos, Incid); gplot(Mesh,Coord); % using gplot: Mesh é a matriz de adjacencia, Mesh(i,j) é não zero se i é % conectado com j. Coord(i,:) é a coord do nó i. % ------------------------------------------------------------------------ % Construção das matrizes Globais % ------------------------------------------------------------------------ [K,C] = MatGlob(NumElem, NumNosEl,Incid,Coord,NumGLTot,PropMat,M_C); % matriz K e C calcula uma só vez, não se atualiza elas % ------------------------------------------------------------------------ % Parâmetros de Inicialização % ------------------------------------------------------------------------ Uan=zeros(NumGLTot*2,1); %Uan é o deslocamento d_m Van=zeros(NumGLTot*2,1); %Van é a velocidade v_m DisplacementSolid=zeros(length(time),NumGLTot); DisplacementFluid=zeros(length(time),NumGLTot); VeloSolid=zeros(length(time),NumGLTot); VeloFluid=zeros(length(time),NumGLTot); Press=zeros(length(time),NumNos); Force=zeros(length(time),2*NumGLTot); % ------------------------------------------------------------------------ % Inicio do Loop % ------------------------------------------------------------------------ for i=2:length(time) deltat = time(i)-time(i-1); ddesloc = desloc(i)-desloc(i-1); presc_velo=ddesloc/deltat; Us=zeros(NumGLTot,1); Uf=zeros(NumGLTot,1); Vs=zeros(NumGLTot,1); Vf=zeros(NumGLTot,1); %% BPE [U,V] = Parabolic_solver_vForm_Spilker(C,K,F,deltat,Uan,Van,presc_velo); % Alternativas % [U,V] = Parabolic_solver_dForm(C,K,F,deltat,Uan,Van,ID_solid,ID_fluid,desloc(i));% Original do algoritmo % [U,V] = Parabolic_solver_vForm(C,K,F,deltat,Uan,Van,ID_solid,ID_fluid,presc_velo); % Armazenando cada parcela do sólido e fluido em vetores Vs=V(1:NumGLTot); Us=U(1:NumGLTot); Vf=V(NumGLTot+1:NumGLTot*2); % ngltot2 é o nº de gdl vezes 2 Uf=U(NumGLTot+1:NumGLTot*2); Vel=(fif*Vf)+(PropMat(1,1)*Vs); % Para o calculo da pressão % armazena resp de cada instante de tempo em uma linha DisplacementSolid(i,:)=Us; DisplacementFluid(i,:)=Uf; VeloSolid(i,:)=Vs; VeloFluid(i,:)=Vf; Force(i,:)=C*V+K*U; %% pressão Pressure=calc_pres(Coord,Incid,Vel,PropMat(1,5)); Press(i,:)=Pressure; %% update Van=V; % o Van e Uan de t=0 é zero, pois o spikler diz isso, também tenho que % fazer algo similar em t = to Uan=U; end %% Pos-Processamento % Separa a força da fase fluida e sólida F_s=Force(:,1:NumGLTot); F_f=Force(:,NumGLTot+1:end); % Plotar a força no nó 36 em y em funçõ do tempo: plot(time/t0,-F_s(:,72),'-') xlabel('Tempo [t/t_0]','FontSize',12); ylabel('Força [N]','FontSize',12); title({'Força y no nó 36 em função do tempo'},'FontSize',14); legend('Força no nó') % Problema, a força no nó 35 em y devia ser igual ao do nó 36 em y %já que ambas tem o mesmo deslocamento aplicado, entretanto isso não ocorre plot(time/t0,-F_s(:,70),'-') xlabel('Tempo [t/t_0]','FontSize',12); ylabel('Força [N]','FontSize',12); title({'Força y no nó 35 em função do tempo'},'FontSize',14); legend('Força no nó') % Calculo da Velocidade Relativa VeloRelative=VeloFluid-VeloSolid; %X e Y VeloRelative_Y=VeloRelative(:,2:2:end); %Apenas Y % Calculo da Deslocamento Relativa DisplacementRelative=DisplacementFluid-DisplacementSolid; %X e Y DisplacementRelative_Y=DisplacementRelative(:,2:2:end); %Apenas Y DisplacementSolid_Y=DisplacementSolid(:,2:2:end); %Apenas Y % gerar o depthu automaticamente no futuro, a partir das coordenadas % gera prof realtiva aux=1; for l=1:2:length(Coord) DepthU(aux)=Coord(l,2)/altura*100; DepthU(aux)=100-DepthU(aux); aux=aux+1; end %% sólido for i=1: NumNos Ux_s(i,1)= U(2*(i-1)+1); Uy_s(i,1)= U(2*(i-1)+2); end %% fluido for i=1: NumNos Ux_f(i,1)= U(2*(i-1)+1+NumGLTot); Uy_f(i,1)= U(2*(i-1)+2+NumGLTot); end %% contorno % contorno deslocamento y no solido em tf scaleFactor=-7.5; figure drawingField(Coord+scaleFactor*[Ux_s Uy_s],Incid,'Q4',Uy_s); hold on drawingMesh(Coord+scaleFactor*[Ux_s Uy_s],Incid,'Q4','k-'); c = colorbar; c.Label.String = 'm'; title('uy sólido no ultimo del_t') % contorno deslocamento y no fluído em tf scaleFactor=-10; figure drawingField(Coord+scaleFactor*[Ux_s Uy_s],Incid,'Q4',Uy_f); hold on drawingMesh(Coord+scaleFactor*[Ux_s Uy_s],Incid,'Q4','k-'); c = colorbar; c.Label.String = 'm'; title('uy fluído no ultimo del_t') %obs: achar o pocisao com a maior pressão (valor absoluto) [A,ind] = max(abs(Press(:))) [m,n] = ind2sub(size(Press),ind) %m-> diz o t em que ocorre %n-> diz o nó % se quiser fazer o contorno pega toda a linha % contorno pressão tf scaleFactor=-10; figure drawingField(Coord+scaleFactor*[Ux_s Uy_s],Incid,'Q4',1e-3*Press(m,:)'); hold on drawingMesh(Coord+scaleFactor*[Ux_s Uy_s],Incid,'Q4','k-'); c = colorbar; c.Label.String = 'KPa'; title('contorno no del_t de maior pressão') %% Comparações % Soluções analiticas fig4_3 = csvread('vel_ana_fig_4.3.csv'); fig4_4 = csvread('pres_ana_fig_4.4.csv'); fig4_5 = csvread('vel_ana_fig_4.5.csv'); % %% BPE SPILKER 1990 % % Plot Velocidade Relativa at 4% depth (nó 28) a=altura-0.04*altura; [~,row]=min(abs(Coord-a)); figure(1) plot(time/t0,1e6*VeloRelative_Y(:,row(2)),'o') hold on plot(fig4_3(:,1),fig4_3(:,2),'-') xlabel('Tempo [t/t_0]','FontSize',12); ylabel('Velocidade [\mum/s]','FontSize',12); title({'Velocidade Relativa em 4% de Profundidade'},'FontSize',14); set(gca,'XTick',0:0.05:0.1,'FontSize',12,'XMinorTick','on') xlim([0 0.1]) ylim([0 0.25]) legend('Solução Numérica','Solução Analitica') grid on % Plot variação da velo relativa do fluido em t=200 com a profundidade figure(4) plot(DepthU,1e6*VeloRelative_Y(401,1:2:NumNos),'ro') hold on plot(fig4_5(:,1),fig4_5(:,2),'-') ylabel('Velocidade Relativa [\mum/s]','FontSize',12); xlabel('Profundidade [%]','FontSize',12); title({'Velocidade Relativa para t=200'},'FontSize',14); set(gca,'XTick',0:20:100,'FontSize',12,'XMinorTick','on') legend('Solução Numérica','Solução exata') grid on % % Plot Pressão figure(3) plot(DepthU,Press(601,1:2:NumNos)*1e-3,'m+') hold on plot(fig4_4(:,1),fig4_4(:,2),'-') ylabel('Pressão (KPa)','FontSize',12); xlabel('Profundidade (%)','FontSize',12); title({'Pressão para t=300seg'},'FontSize',14); set(gca,'XTick',0:20:100,'FontSize',12,'XMinorTick','on') xlim([0 100]) legend('Solução Numérica','Solução Numérica') grid on hold on % % % % % % figure2 = figure(2); % set(figure2,'Color',[1 1 1]); % axes1 = axes('Parent',figure2); % box(axes1,'on'); % grid(axes1,'on'); % hold(axes1,'on'); % plot(axes1,time/t0,abs(DisplacementRelative_Y(:,28)),'LineWidth',1,... % 'Marker','o','LineWidth',1,'Color','b'); % xlim([0 0.1]) % ylim([0 max(abs(DisplacementRelative_Y(:,28)))]) % legend('Numerical','Locat
评论
    相关推荐
    • matlabcnhelp.rar
      matlab中文帮助很难找的,快速下载
    • MobilePolice.rar
      移动警察,车牌识别,车牌定位系统源代码,已经运用在移动车载稽查系统中。
    • SVM(matlab).rar
      支持向量机(SVM)实现的分类算法源码[matlab]
    • svm.zip
      用MATLAB编写的svm源程序,可以实现支持向量机,用于特征分类或提取
    • Classification-MatLab-Toolbox.rar
      模式识别matlab工具箱,包括SVM,ICA,PCA,NN等等模式识别算法,很有参考价值
    • VC++人脸定位实例.rar
      一个经典的人脸识别算法实例,提供人脸五官定位具体算法及两种实现流程.
    • QPSK_Simulink.rar
      QPSK的Matlab/Simulink的调制解调仿真系统,给出接收信号眼图及系统仿真误码率,包含载波恢复,匹配滤波,定时恢复等重要模块,帮助理解QPSK的系统
    • LPRBPDemo2009KV.rar
      车牌识别,神经网络算法,识别率高达95%,识别时间低于80ms。
    • MODULATION.RAR
      这个源程序代码包提供了通信系统中BPSK,QPSK,OQPSK,MSK,MSK2,GMSK,QAM,QAM16等调制解调方式 用matlab的实现,以及它们在AWGN和Rayleigh信道下的通信系统实现及误码率性能
    • algorithms.rar
      十大算法论文,包括遗传算法,模拟退火,蒙特卡罗法等等,对于初学者很有帮助!!