Crank Nicolson scheme using Finite Volume for 1D Transient Problem.

1 次查看（过去 30 天）
Shahid Hasnain
Shahid Hasnain 2021-9-11
编辑： darova 2021-9-13
Dear Community,
I have attached discretized scheme (Problem_1_1 and Problem_1_2). According to complete MATLAB code which is working. I am getting suspicious about results when I compare them with books that are related to the Fully Implicit scheme. A few days earlier I posted my question. Any best idea which can do to satisfy?
Your help will be highly appreciated.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
clc
%% Sort out Inputs
Theta = 0.5; % Crank - Nicolson scheme
% Number of Control Volumes
N = 5;
% Domain length
L = 0.02; %[m]
% Grid Spacing (Wall thickness)
h = L/(N);
% Diffusivity (Thermal conductivity)
k = 10; %[W/m-K]
% Specific Heat Capacity
Row_c = 10e6; %[J/m^(3)-K]
% The center of the first cell is at dx/2 & the last cell is at L-dx/2.
x = h/2 : h : L-(h/2);
% Diffusivity (Heat)
alpha = k/Row_c;
% Simulation Time
t_Final = 40;
t_steps = 8;
% Discrete Time Steps
dt = t_Final/(t_steps);
% Left Surface Temperature
dTdt = 0; %[\circ c]
% Right Surface Temperature
T_b = 0; %[\circ c]
% Parameteric Setup
lambda = (alpha.*dt)/(h^2)
ap_0 = Row_c * h/(dt);
dt_Check = Row_c *(h^2)/k;
% Initializing Variable
% C*T^(n+1) = BB*T^(n)
%
% At Node # 1 dpdx=0 is adjusted during simplification of scheme
%
% (1+lamda/2)T_P = (1/2)*lamda*T_E +(T_P)_0
%
% At Node # 5
%
% (1+2*lamda)T_P = lamda*T_W + lamda*T_b+(T_P)_0
%
% At Node Interior Nodes
%
% -(lamda/2)*T_W + (1+lamda)T_P -(lamda/2)*T_E = (lamda/2)*T_W_0 + (1-lamda)T_P_0 +(lamda/2)*T_E_0
%
% Unknowns at time level n
T_Old = zeros(N,1);
% Unknowns at time level n+1
T_New = zeros(N,1);
% Initial Value (Initial Condition)
T_Old(:) = 200;
C = zeros(N,N);
D = zeros(N,1);
for t=1:dt:t_Final
for i = 1: N
if i > 1 && i < N
C(i,i) = (1 + lambda);
BB(i,i) = (1-lambda);
C(i,i+1) = -0.5*lambda;
BB(i,i+1) = 0.5*(lambda);
C(i,i-1) = -0.5*lambda;
BB(i,i-1) = 0.5*(lambda);
D(i)=T_Old(i);
elseif i == 1
% For 1st boundary node
C(i,i) = (1 + 0.5*lambda);
C(i,i+1) = -0.5*lambda;
D(i) = T_Old(i);
else
% For Last boundary node
C(i,i) = (1 + 2*lambda);
C(i,i-1) = -lambda;
BB(i,i) = 0;
D(i) = T_b*lambda + T_Old(i);
end
end
Temp = inv(C)*BB;
T_New = Temp*T_Old;
T_Old = T_New;
end
T_New

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Gas Dynamics 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by