Need help with this MATLAB HW
2 views (last 30 days)
Show older comments
%% Steady-state 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;
%% 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(:)=0; % right boundary
BC.top.a(:) = 1; BC.top.b(:)=0; BC.top.c(:)=bc(x,1); % top boundary
BC.bottom.a(:) = 1; BC.bottom.b(:)=0; BC.bottom.c(:)=bc(x,0); % bottom boundary
%% define the transfer coeffs
D_x = 1;
D_y = 1;
D = createFaceVariable(m, [D_x, D_y]);
Mdiff = diffusionTerm(D);
%% define source term
[X, Y] = ndgrid(x, y);
S = @(X,Y)(6*X.*Y.*(1-Y)-2*X.^3)
s1 = constantSourceTerm(createCellVariable(m,S(X,Y)));
%s1 = 0*s1;
[Mbc, RHSbc] = boundaryCondition(BC);
M = Mdiff+Mbc;
RHS = -s1+RHSbc;
c = solvePDE(m,M, RHS);
figure(1);visualizeCells(c);caxis([0,10]);drawnow;
% view 3D plot
% add ghost points to x and y
x = [-0.5*dx; x; L+0.5*dx];
y = [-0.5*dx; y; L+0.5*dx];
[X, Y] = ndgrid(x, y);
% 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);
figure(2); surf(X, Y, u);
%axis([0 L 0 L 0 10]);
function g = bc(x,y)
g = -sin(10*pi*x);
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!