基于有限元的流体模拟仿真源码程序 - matlab模式识别

下载频道> 资源分类> matlab源码> 模式识别> 基于有限元的流体模拟仿真源码程序

标题：基于有限元的流体模拟仿真源码程序
分享到：

所属分类: 模式识别资源类型：程序源码文件大小: 140.45 KB 上传时间: 2016-01-24 22:47:13 下载次数: 12 资源积分：1分提供者: matlab源码 20160124104659645

内容：
基于有限元的流体模拟仿真源码程序，程序员在编程的过程中可以参考学习使用，希望对IT程序员有用，此源码程序简单易懂、方便阅读，有很好的学习价值！

部分代码如下：

% See also http://www.particleincell.com/blog/2012/matlab-fem/

close all

% Load mesh data, coordinates of element 'p' and connectivity 't'==========

load MeshVessel;

% Set nodes of Dirichlet- and Neumann bounaries============================

Neumann=[1 14;14 2;7 25;25 6];

Dirichlet=[4;24;5];

% Inlet velocity===========================================================

V0=5;

% =========================================================================

TotalNumberOfNodes=size(p,1);

NumberOfElements=size(t,1);

% Initialisation of K and F================================================

K=zeros(TotalNumberOfNodes,TotalNumberOfNodes);

F=zeros(TotalNumberOfNodes,1);

% Computation of Ke & assemblage of K======================================

for e=1:NumberOfElements

NodesOfElement=t(e,:);

P=[ones(1,3);p(NodesOfElement,:)'];

AreaOfElement=abs(det(P))/2;

Coeffs=inv(P);

diffPhi=Coeffs(:,2:end);

B{e}=[];

for i=1:3

Bsmall=[diffPhi(i,1);

diffPhi(i,2)];

B{e}=[B{e},Bsmall];

end

Ke=diffPhi*diffPhi'*AreaOfElement;

K(NodesOfElement,NodesOfElement)=K(NodesOfElement,NodesOfElement)+Ke;

end

% F_Gamma==================================================================

for e=1:size(Neumann,1)

knooppunten=Neumann(e,:);

P=p(knooppunten,:);

LengthOfElement=norm(diff(P));

% Value of directional derivative, perpendicular to the edge, and

% pointed inside, resulting in a negative sign

FeGamma=(-V0)*1/2*[1;1]*LengthOfElement;

F(knooppunten)=F(knooppunten)+FeGamma;

end

%==========================================================================

K(Dirichlet,:)=0;% put zeros in boundary rows of K and F

K(Dirichlet,Dirichlet)=eye(numel(Dirichlet),numel(Dirichlet));

% Set Dirichlet boundaries (value of fixed phi)============================

F(Dirichlet)=[0;0;0];

%==========================================================================

phi=K\F;

% Initalize matrix with centroids of each element==========================

Pmean=[];

% Computation of element velocities and magnitude of nodal velocities=====

V=[];

MagnitudesOfNodalVelocities=zeros(TotalNumberOfNodes,1);

for e=1:NumberOfElements

NodesOfElement=t(e,:);

P=[ones(1,3);p(NodesOfElement,:)'];

% Storage of centroids of each element

Pmean=[Pmean;mean(p(NodesOfElement,:))];

Coeffs=inv(P);

diffPhi=Coeffs(:,2:end);

phi_e=phi(NodesOfElement);

% Storage of element velocities [vxe,vye] in V

ve=(diffPhi'*phi_e)'; % velocity in element e

V=[V;ve];

MagnitudeOfVelocity=sqrt(ve(1)^2+ve(2)^2);

MagnitudesOfNodalVelocities(NodesOfElement)=MagnitudesOfNodalVelocities(NodesOfElement)+MagnitudeOfVelocity;

end

% Determine node frequencies (occurrences)=================================

NodeFrequencies=[];

for i=1:TotalNumberOfNodes

NodeFrequencies=[NodeFrequencies;numel(find(t==i))];

end

MagnitudesOfNodalVelocities=MagnitudesOfNodalVelocities./NodeFrequencies;

% Plotting solution========================================================

set(gcf,'color','w');

hold on;

axis equal;

axis off;

for e=1:NumberOfElements

NodesOfElement=t(e,:);

C=p(NodesOfElement,:);

fill(C(:,1),C(:,2),MagnitudesOfNodalVelocities(t(e,:)),'EdgeColor',[0.5,0.5,0.5]);

end

% Plotting velocity vecotrs================================================

quiver(Pmean(:,1),Pmean(:,2),V(:,1),V(:,2),'color',[1 1 0],'LineWidth',1.5)

% Show node numbers========================================================

text(p(:,1),p(:,2),phi,num2str([1:size(p,1)]'),'FontSize',10','Color','k',...

'FontName','Cambria','HorizontalAlignment','left','VerticalAlignment','bottom');

% Computation of streamlines===============================================

Xp=p(:,1);

Yp=p(:,2);

% Timestep

dt=1e-2;

% Locate inital position of particles, add some, loose some================

Particles=[0 ,0.09;

0 ,0.25;

0 ,0.55;

0.0582 ,1.4658;

0.3073 ,1.8348];

% Loop over all elements===================================================

for j=1:size(Particles,1);

counter=0;

particle=Particles(j,:);

for i=1:60

for e=1:NumberOfElements

% Locate the element where the particle is to be found

if inpolygon(particle(1),particle(2),Xp(t(e,:)),Yp(t(e,:)))==1

%==========================================================

POS=[particle(1),particle(2),V(e,1),V(e,2)];

H(4*counter+1:4*(counter+1),j)=[POS];

particle=particle+[V(e,1)*dt,V(e,2)*dt];

counter=counter+1;

break

end

e=NumberOfElements;

end

% Plot the trajectories====================================================

for k=0:floor(size(H,1)/4)-2

% Plot the final position of the particles=============================

for i=1:size(H,2);

plot([H(4*k+1,i),H(4*(k+1)+1,i)],[H(4*k+2,i),H(4*(k+1)+2,i)],'color',[1 1 1],'LineWidth',1.5);

end

%======================================================================

drawnow

end

文件列表（点击上边下载按钮，如果是垃圾文件请在下面评价差评或者投诉）:

基于有限元的流体模拟仿真源码程序/
基于有限元的流体模拟仿真源码程序/Flow.jpg
基于有限元的流体模拟仿真源码程序/MeshVessel.mat
基于有限元的流体模拟仿真源码程序/PotentialFlow2D.m

关键词：有限元 流体 源码

matlab源码下载排行

回到顶部

友情链接

matlab源码网站源码编程语言网络编程手机通讯编程 KUNAG 官网 KUNAG 电商更多

: 联系方式| 版权声明| 招聘信息| 广告服务| 银行汇款| 法律顾问| 兼职技术| 付款方式| 关于我们|; 网站客服程序员兼职招聘 411241577@qq.com; 沪ICP备19040327号-3; 公安备案号：沪公网安备 31011802003874号; 库纳格流体控制系统（上海）有限公司版权所有; Copyright © 1999-2014, GUSUCODE.COM, All Rights Reserved

关键词： 有限元 流体 源码

相关推荐

matlab源码下载排行

友情链接

关键词：有限元流体源码