Error using / Matrix dimensions must agree.
显示 更早的评论
Hi there, I'm trying to piece together this 2-D implicit heat code. However I am having trouble creating kappa which can be used in further loops when implicitly solving the equation. Firstly this does not work, if I create a matrix (Km) and run the code it says it cannot run as the sizes do not match on either side.
I have attached the code below:
% Solution of 2D temperature equation on a staggered grid
% for a non-moving medium with constant conductivity;
% conversion to implicit formulation.
% Clean all variables
clear all;
% Clear all figures
clf;
% Numerical model parameters
% Boundary conditions = constant temperature
% setup correspond to a squre temperature wave
% Model size, m
xsize=100000; % Horizontal
ysize=100000; % Vertical
% Numbers of nodes
Nx=51; % Horizontal
Ny=31; % Vertical
% Grid step
dx=xsize/(Nx-1); % Horizontal
dy=ysize/(Ny-1); % Vertical
% Defining coordinates of staggered points
% P-Points
xpr=-dx/2:dx:xsize+dx/2; % Horizontal
ypr=-dy/2:dy:ysize+dy/2; % Vertical
% Defining the number of timesteps in that the loop iterates through
ntimesteps=20; % number of timesteps
% Material properties
K=3; % Thermal conductivity, W/m/K
% Defining density in xy-points
RHOCP=zeros(Ny+1,Nx+1); % Volumetric heat capacity, J/K/m^3
% Making vectors for nodal points positions (basic nodes)
x=0:dx:xsize; % Horizontal
y=0:dy:ysize; % Vertical
% Timestep
dt=min(dx,dy)^2/(4*K/min(min(RHOCP))); % Timestep, s
% Creating temperature arrays (profiles) for iteration of temperature changes through time
TDT=zeros(Ny+1,Nx+1); % New temperature, K
T0=zeros(Ny+1,Nx+1); % Old temperature, K
% Populating T0 using assigned temperatures and RHOCP at eulerian nodes
for j=1:1:Nx+1
for i=1:1:Ny+1
d=((xpr(j)-xsize/2)^2+(ypr(i)-ysize/2)^2)^0.5; % distance to the centre from the given nodal point
if(d<20000)
% Plume
RHOCP(i,j)=3.2e+6;
T0(i,j)=1873;
else
% Medium
RHOCP(i,j)=3.3e+6;
T0(i,j)=1573;
end
% Sticky air layer
if(ypr(i)<ysize*0.2)
RHOCP(i,j)=3.5e+6;
T0(i,j)=273;
end
end
end
% Formulating kappa (κ) to simplify the implicit equation:
% Equation: κ = k/RHO*Cp
% Thereby formulated as:
kappa(i,j)=K/RHOCP;
% Beginning Time-Cycle for implicit solution
timesum=0; % Elapsed time
for t=1:1:ntimesteps
% Matrix of coefficients initialization for implicit solving
L=sparse(Nx*Ny,Nx*Ny);
% Vector of right part initialization for implicit solving
R=zeros(Nx*Ny,1);
% Implicitly solving for the 2-D temperature equation using 5-point cross:
% RHOCP*dT/dt=K(d2T/dx^2+d2T/dy^2)
% | ∆x
% _____|________________|___________________
% | -
% | To2
% | -
% | -
% ∆y---| --To1-------To3-------To5--
% | -
% | -
% | To4
% | -
% |
%
% Composing matrix of coefficients L()
% and vector (column) of right parts R()
% Process all grid points
for i=1:1:Ny
for j=1:1:Nx
% Global index
k=(j-1)*Ny+i;
% Boundary nodes
if(i==1 || i==Ny || j==1 || j==Nx)
% Upper boundary
if(i==1)
L(k,k)=1;
R(k,1)=273;
end
% Lower boundary
if(i==Ny)
L(k,k)=1;
R(k,1)=1573;
end
% Left boundary
if(j==1 && i>1 && i<Ny)
L(k,k)=1;
R(k,1)=0;
end
% Right boundary
if(j==Nx && i>1 && i<Ny)
L(k,k)=1;
R(k,1)=0;
end
% Internal nodes
else
% dT/dt=kappa*(d2T/dx2+d2T/dy2)
% TDT(i)/dt-kappa*(TDT(i,j-1)-2*TDT(i,j)+TDT(i,j+1))/dx^2-kappa*(TDT(i-1,j)-2*TDT(i,j)+TDT(i+1,j))/dy^2=T0(i,j)/dt
L(k,k-Ny)=-kappa/dx^2; % coefficient for T(i,j-1)
L(k,k-1)=-kappa/dy^2; % coefficient for T(i-1,j)
L(k,k)=1/dt+2*kappa/dx^2+2*kappa/dx^2; % coefficient for T(i,j)
L(k,k+1)=-kappa/dy^2; % coefficient for T(i+1,j)
L(k,k+Ny)=-kappa/dx^2; % coefficient for T(i,j+1)
R(k,1)=T0(i,j)/dt;
end
end
end
% Obtaining solution
S=L\R;
% Reloading solution to the new temperature array
for i=1:1:Ny
for j=1:1:Nx
% Global index
k=(j-1)*Ny+i;
% Reload solution
TDT(i,j)=S(k);
end
end
% Reloading TDT to T0 for the interation of the next timestep
T0=TDT;
% Openning figure
figure(1);
% Visualization of geometrical arrays OMEGA(), PSI(), vx(), vy()
figure(1);clf;colormap('Jet')
% Drawing RHOCP()
subplot(1,2,1)
pcolor(xpr,ypr,RHOCP)
shading flat
colorbar
axis ij image;
title('RHOCP, J/m^3/K')
% Plotting implicit solution (TDT) for every timestep
subplot(1,2,2);
pcolor(xpr,ypr,TDT);
shading interp;
colorbar;
axis ij image;
title(['Implicit: step=',num2str(t),' time,Myr=',num2str(timesum/(1e+6*365.25*24*3600))]);
% Stop for .1 second
pause(0.2);
% Add time counter
timesum=timesum+dt;
% Reassigning temperature profile (TDT) for the next step as base
% temperature profile (T0)
T0=TDT;
end
The error code is as follows:
Error using /
Matrix dimensions must agree.
Any help or a prod in a direction will be greatly appreciated.
Regards,
Aleks
采纳的回答
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Heat Transfer 的更多信息
Translated by 
