Need help with this MATLAB HW
1 view (last 30 days)
Show older comments
% Use FVTool to solve a 2D heat equation on \Omega=(0,1)x(0,1) and time
% domain [0, T]
%% Transient diffusion equation
%% PDE and boundary conditions
% The transient diffusion equation reads
%
% $$\alpha\frac{\partial c}{\partial t}+\nabla.\left(-D\nabla c\right)=0,$$
%
% where $c$ is the independent variable (concentration, temperature, etc)
% , $D$ is the diffusion coefficient, and $\alpha$ is a constant.
clc; clear;
%% add paths
addpath('FVTool');
FVToolStartUp
%% Define the domain and create a mesh structure
L = 1; % domain length
Nx = 20; % number of cells
dx = L/Nx;
m = createMesh2D(Nx,Nx, L,L);
x = m.cellcenters.x;
y = m.cellcenters.y;
[X, Y] = ndgrid(x, y);
%% Create the boundary condition structure
BC = createBC(m); % all Neumann boundary condition structure
BC.left.a(:) = 0; BC.left.b(:)=1; BC.left.c(:)=10; % left boundary
BC.right.a(:) = 0; BC.right.b(:)=1; BC.right.c(:)=10; % right boundary
BC.top.a(:) = 0; BC.top.b(:)=1; BC.top.c(:)=10; % top boundary
BC.bottom.a(:) = 0; BC.bottom.b(:)=1; BC.bottom.c(:)=10; % bottom boundary
%% define the transfer coeffs
D_x = 1;
D_y = 0.01;
D = createFaceVariable(m, [D_x, D_y]);
alfa = createCellVariable(m, 1);
%% define initial values
c_init = uinitial(X, Y);
c_old = createCellVariable(m, c_init,BC); % initial values
c = c_old; % assign the old value of the cells to the current values
S = @(X,Y)(10*exp(-50*((X-0.5).^2+(Y-0.5).^2)));
s1 = constantSourceTerm(createCellVariable(m,S(X,Y)));
s1 = 0*s1;
% add ghost points
x = [-0.5*dx; x; L+0.5*dx];
y = [-0.5*dx; y; L+0.5*dx];
[X, Y] = ndgrid(x, y);
%% loop
dt = 0.01; % time step
final_t = 0.2;
for t=dt:dt:final_t
[M_trans, RHS_trans] = transientTerm(c_old, dt, alfa);
%Dave = harmonicMean(D);
Dave = D;
Mdiff = diffusionTerm(Dave);
[Mbc, RHSbc] = boundaryCondition(BC);
M = M_trans-Mdiff+Mbc;
RHS = RHS_trans+RHSbc;
RHS = -s1 + RHS;
c = solvePDE(m,M, RHS);
c_old = c;
%figure(1);visualizeCells(c);drawnow;
% fix corner points as ghost point
u = c.value;
u(1,1) = u(1,2);
u(1,end) = u(1,end-1);
u(end,1) = u(end,2);
u(end,end) = u(end,end-1);
% view 3D plot
figure(1);
surf(X, Y, u);
axis([0 L 0 L 0 30]);
%axis([-0.5*dx L+0.5*dx -0.5*dx L+0.5*dx 0 30]);
pause(0.1)
end
function u = uinitial(x, y)
u = 0*x;
u(y<=0.5) = 30.0;
end
1 Comment
Answers (0)
See Also
Categories
Find more on Mathematics and Optimization in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!