有限元法模拟牛顿流体matlab程序

BE2=[[Nx:Nx:Ny*Nx]',2*ones(size([1:Ny]')),BBB1,BBB2];
BE3=[[Nx*(Ny-1)+1:Nx*Ny]',3*ones(size([1:Nx]')),CCC1,CCC2];
BE4=[[1:Nx:(Ny-1)*Nx+1]',4*ones(size([1:Ny]')),DDD1,DDD2];
JXYV(i,1)=R*cos(theta/180*pi+theta1);
JXYV(i,2)=R*sin(theta/180*pi+theta1);
end
%%%%%%% 网格平面旋转
%%%%%%% 四边形二次单元JMV生成
k=0;
for i=1:Ny
for j=1:Nx
Dx=L/Nx/2; %水平方向网格间距
Dy=H/Ny/2; %竖直方向网格间距
%%%%%%% 单元间距
%%%%%%% 结点分布拓扑
AAA=zeros(Ny*2+1,Nx*2+1);
for i=1:2:2*Nx+1
AAA(1,i)=(i+1)/2;
end
for i=1:Nx 3
JBV=[JBV1;JBV3];
P2=0;
P4=1000;
JBP2=[BE2,ones(size(BE2(:,1)))*P2,ones(size(BE2(:,1)))*P2];
JBP4=[BE4,ones(size(BE4(:,1)))*P4,ones(size(BE4(:,1)))*P4];
k=k+1;
JMV(k,:)=[AAA(2*i-1,2*j-1),AAA(2*i-1,2*j),...
AAA(2*i-1,2*j+1),AAA(2*i,2*j-1),AAA(2*i,2*j),...
AAA(2*i,2*j+1),AAA(2*i+1,2*j-1),AAA(2*i+1,2*j),...
%%%%%%%%%%%%% 迭代步骤A结束 %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% 迭代步骤B开始 %%%%%%%%%%%%
thetay1=pi-theta/180*pi; % 1号边界外发现方向与Y轴夹角
thetax2=theta/180*pi; % 2号边界外发现方向与X轴夹角
thetay2=pi/2-thetax2; % 2号边界外发现方向与Y轴夹角
thetax3=pi-pi/2+theta/180*pi; % 3号边界外发现方向与X轴夹角
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% 迭代步骤E开始 %%%%%%%%%%%%
%%%%%%%%% 计算单元方程系数矩阵子块并组装 %%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% BE数据生成
%%%%%%% 调用四边形网格绘制程序
rectangle_grid(JMP,JXYV);
%%%%%%% 调用四边形网格绘制程序
%%%%%%% 清除多余变量，存储网格数据
clear Dx Dy H L Nx Ny i j k
clear theta theta1 R AAA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u1=0;v1=0;
u3=0.01;v3=0;
JBV1=[BP1',u1*ones(size(BP1))',v1*ones(size(BP1))'];
JBV3=[BP3',u3*ones(size(BP3))',v3*ones(size(BP3))'];
CCC2=ones(Nx,1)*cos(thetay3); % 上边方向余弦
DDD1=ones(Ny,1)*cos(thetax4); % 左侧边方向余弦
DDD2=ones(Ny,1)*cos(thetay4); % 左侧边方向余弦
BE1=[[1:Nx]',ones(size([1:Nx]')),AAA1,AAA2];
%%%%%%% 区域几何尺寸及网格划分参数
%%%%%%% 总单元数和结点数
E=Nx*Ny; %总单元数
Nz=(2*Nx+1)*(2*Ny+1); % 二次单元结点总数
Nd=(Nx+1)*(Ny+1); % 线性单元结点总数
%%%%%%% 总单元数和结点数
%%%%%%% 单元间距
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B1=zeros(Nd,Nz);
B2=zeros(Nd,Nz);,Nz);
D21=zeros(Nz,Nz);
D22=zeros(Nz,Nz);
C1=zeros(Nz,Nd);
AAA(1,2*i)=(Nx+1)*(Ny+1)+i;
end
for i=1:2*Nx+1
AAA(2,i)=(Nx+1)*(Ny+1)+Nx+i;
end
for j=3:2:2*Ny+1
for i=1:2:2*Nx+1
AAA(j,i)=(i+1)/2+(Nx+1)*(j-1)/2;
%%%%%%%%%%%%% 设定物料黏度 %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
niandu=1000;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% 迭代步骤B结束 %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% 迭代步骤D开始 %%%%%%%%%%%%
%%%%%%%%%% 初始化总体方程各个子块 %%%%%%%%%%
end
end
for j=3:2:2*Ny+1
for i=1:Nx
AAA(j,2*i)=(Nx+1)*(Ny+1)+(Nx+2*Nx+1)*(j-1)/2+i;
end
end
for j=4:2:2*Ny
for i=1:2*Nx+1
AAA(j,i)=(Nx+1)*(Ny+1)+(j/2)*Nx+(2*Nx+1)*(j/2-1)+i;
%%%%%%% 四边形线性单元JMP和JXYP生成
%%%%%%% BP数据生成
BP1=AAA(1,:);
BP2=AAA(:,2*Nx+1);
BP3=AAA(2*Ny+1,:);
BP4=AAA(:,1);
%%%%%%% BP数据生成
%%%%%%% BE数据生成
thetax1=pi/2-theta/180*pi; % 1号边界外发现方向与X轴夹角
thetay3=pi-theta/180*pi+pi; % 3号边界外发现方向与Y轴夹角
thetax4=(180+theta)/180*pi; % 4号边界外发现方向与X轴夹角
thetay4=pi/2+theta/180*pi; % 4号边界外发现方向与Y轴夹角
AAA1=ones(Nx,1)*cos(thetax1); % 底边方向余弦
for i=1:E
for ie=1:9;
clear thetax1 thetax2 thetax3 thetax4
clear thetay1 thetay2 thetay3 thetay4
save msh
%%%%%%% 清除多余变量，存储网格数据
第一部分：速度压力有限元求解
（2）主程序
clc
clear
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% 迭代步骤A开始 %%%%%%%%%%%%
%%%%%%%%%%%%% 读取网格数据 %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load msh
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
JBP=[JBP2;JBP4];
clear JBV1 JBV3 BP1 BP3 BP4
clear JBP2 JBP4 P2 P4
clear BE1 BE2 BE3 BE4
clear u1 v1 u3 v3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% 迭代步骤C结束 %%%%%%%%%%%%
end
end
%%%%%%% 结点分布拓扑