Need help with this MATLAB HW

2 次查看（过去 30 天）

Dhinesh Nagarajan 2021-4-20

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/807362-need-help-with-this-matlab-hw

评论： Walter Roberson 2021-4-20

% 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