Computational Fluid Dynamics is the Future
  • Main Page
    • Purpose of Website
    • About Me >
      • My PhD Thesis
      • My Teaching at the American University of the Middle East
      • My Teaching at the University of Sharjah
      • My Books & Codes
      • My CFD Projects
      • My SolidWorks Projects
      • My Family
      • In the Living Memory >
        • Family who contributed to My Personality
        • Lecturers Who Contributed to My Personality
      • Testimonials
    • CFD2012 Blog
    • معلومات عن الموقع
    • Page Contact >
      • Feedback Form
  • Research
    • C++ for Excel
    • Undergrad Stage Advice >
      • Cheat Sheet
      • Dealing with Dyslexia
      • Work/Research Placement
    • Masters Stage Advice >
      • PIV Lab
      • Prototype Modelling Lab
      • Field Trips 2006
      • Fuel Atomization Lab
      • Cardiff Airbus Seminar
      • Research Labs
      • GT Onsite Trips
      • On-Site Visits
    • PhD Stage Advice >
      • CFD Modelling >
        • Online CFD Codes
        • CFD Mandatory Reading List >
          • CFD Book Reviews >
            • Turbulence book Reviews
          • Turbulence Modelling Books
          • Finite Element Reading List
        • Eulerian and Lagrangian Descriptions
        • Multigrid Method
        • Finite Volume Method
        • Quantum Physics
        • Navier-Stokes Equations >
          • Atmospheric
          • CFD Simulation Validation
        • Numerical Methods >
          • Numerical Methods Book Reviews
          • Derivatives of Different Orders
          • Vector Calculus
          • Solvers >
            • Iterative Method
          • Data Structures
        • Grid Classification >
          • Mesh Geneation Book Reviews
          • Delaunay Trangulation
        • Reactive Flows >
          • Anaerobic Digestion
          • Combustion Theory Notes
          • Combustion Book Reviews
          • Swirl Flow and Combustion
          • Combustion Reading List
          • Working with Mixtures
          • Combustion Theory
      • Experimentation for CFD >
        • Diffusion Sensors
        • Experimental Wind Tunnels >
          • Wind Tunnel Walk Around
          • Wind Tunnels Books
      • Design of Experiment
      • Optimization
      • My Teaching Activities
      • Google Scholar Profile
      • Yearly Upgrade Report
      • Attending Regular Seminars >
        • Publication Reading
        • Making a Presentation
        • Research Collaboration Fundamentals >
          • Writting a Publication
          • Trip Planning
      • Supervisor Guidance >
        • Studies Budget >
          • Making a Pause for Your Studies
      • Thesis Writting >
        • Putting Together the Thesis
        • List of Symbols
        • Thesis Abstract
        • Thesis Rational and Finding the Gap
        • Thesis Literature Review
        • Thesis Methodolgy
        • Thesis Hypothesis
        • Thesis Conclusion
        • Thesis Check List
      • Referencing >
        • LaTeX
        • EndNote >
          • Setting Up EndNote with Google
          • EndNote Selecting Reference Method
      • Ideal Image of the VIVA >
        • Dealing with Correction Obstactles
        • Thesis Corrections
    • Postdoc >
      • Whitelee Windfarm
      • Dumbarton Scottish Maritime Museum
      • Meetings and Negotiating
      • Student Recommendations
      • Team Work
      • Writting a Research Proposal
      • Conference Organization
      • Research Networking
      • Supervising Students >
        • Types of PhD students >
          • Mind Mapping
    • Research Trends
    • Photoshop >
      • Inserting Text in Photoshop
      • Using Brush in Photoshop
      • Using Blur in Photoshop
    • Microsoft Office Skills >
      • Flow Chart
      • Microsoft Outlook
      • Making a Poster
      • Power Point
      • Making Gantt Chart
      • Mircosoft Word Thesis Layout >
        • Adding Rear Color in Word
        • Applying a Border Line in Word
        • Steps to Write a Thesis
        • Document Heading
        • Applying Chapter Headings
        • Document Footer
        • Using Text Box
        • Using Quick Parts
      • Microsoft Excel >
        • Reading Data into Excel
    • Jobs Search >
      • Jobs Requirments >
        • Requesting Copy of Reference
        • Work Email Formats
        • Regular Rejection Excuses >
          • Email Job Replies
      • CFD Job Sites >
        • Other Jobs Sites
      • Telephone Job Interview
      • Interview Clothing
  • ANSYS-Flow Modelling
    • ANSYS CFX Tutorials >
      • ANSYS CFX Introductory Tutorials >
        • Boundary Conditions
        • ANSYS CFX Introductory 2
      • ANSYS-CFX Porous Media >
        • ANSYS-CFX Porous Media Models
        • ANSYS-CFX Porous Media Bugs
      • ANSYS CFX Single Domain Wind Turbines >
        • ANSYS CFX Double Domain Wind Turbines
        • Rotating Wind Turbine
        • Wind Turbine Stress Analysis
      • ANSYS-CFX Turbo Machinery >
        • ANSYS-CFX Pump Simulation
        • ANSYS-CFX Turbine Cooling >
          • Steps to Model Gas Turbine Blades
      • ANSYS-CFX Formula One CAR >
        • Car Aerodynamics Books
      • ANSYS-CFX Heat Exchanger >
        • ANSYS-CFX Finned Heat Exchanger
        • Phase Change Heat Exchangers
        • Shell and Tube Heat Exchangers
        • Heat Exchangers Reading List
        • ANSYS CFX Heat Exchanger Tube Wear
        • ANSYS CFX Thermal Radiation
      • ANSYS CFX Combustion >
        • ANSYS-CFX Gas Turbine Combustor
        • ANSYS-CFX Multiphase Combustion Modelling
        • ANSYS CFX Flamelet
      • ANSYS-CFX Spary Modelling >
        • ANSYS-CFX Spray Modelling in Car Engines
        • ANSYS-CFX Resolving Multiphase Interface
        • ANSYS CFX Fluid/Solid Transport
        • ANSYS-CFX Air Assisted Sprays
        • Applying a Function of Time
        • ANSYS CFX Rosin Rammler
        • ANSYS CFX Nukiyama Tanasawa
        • SAUTER MEAN DIAMETER
      • ANSYS Geophysical Simualtions >
        • ANSYS Modelling Terrain
      • ANSYS-CFX Tank Sloshing
      • ANSYS Combustion Engines >
        • Setting up ICE Simulation
        • Applying Radiation in the Engine
      • ANSYS-CFX Pelton Turbine
      • ANSYS Flow Modelling Throttle Valve
      • ANSYS-CFX Immersed Solid
      • ANSYS CFX Changing Airfoil Para
      • ANSYS Flow around Buildings
      • ANSYS CFX Compressible Flows
      • ANSYS CFX Solid Particles
      • ANSYS Horizontal Francis Turbine
      • ANSYS-CFX Kaplan Turbine
      • ANSYS-CFX Hair Dryer
      • Types of Flaps >
        • Wings
    • ANSYS-FLUENT >
      • ANSYS FLUENT Simulation Setup >
        • ANSYS FLUENT Boundary Conditions
      • ANSYS FLUENT UDF
      • ANSYS-FLUENT Heat Exchanger Tutorial
    • ANSYS Design Modeller >
      • Blade Modelling >
        • Modelling a Turbine Blade
        • ANSYS Blade Modelling
        • ANSYS BladeGen Axial Compressor
        • ANSYS BladeGen Centrifugal Compressor
      • ANSYS Design Modeler Operations
      • ANSYS Design Modeler Boolean Operations
      • Design Modeller for Turbo Machinery
    • ANSYS CFX Meshing >
      • ICEM >
        • ICEM Introductory Tutorials
        • ICEM Surfacing
        • ICEM Parallel Meshing and Repair
        • ICEM Theory
      • Importing External Geometry to ANSYS
      • ANSYS CFX Types of Meshing >
        • Mesh (Refinement/Relevance)
        • Mesh Number of Cells
        • Mesh Types
        • Mesh Sizing
      • CFX Meshing Online Material
    • ANSYS Simulation Setup >
      • ANSYS CFX Lang CEL >
        • ANSYS-CFX Time Dependent Boundary Conditions
        • Applying a time dependent Velocity Profile
        • Inserting Equation into ANSYS
      • Workbech File Structure
      • Dealing With Memory Problems >
        • Calculating Resources
        • Calculation Guide Lines
      • Applying a Source Term
      • Applying a Velocity Profile to BC
      • Varabile Dependent Boundary Condition >
        • Time dependent Boundary Condition
        • Temperature Dependent Boundary Condition
      • Applying a Riged Body
      • CFX Data Transfer
      • CFX User Functions
      • ANSYS CFX Material Addition
    • ANSYS CFX Steady/Unsteady >
      • ANSYS CFX Time Stepping
      • Time Step Planning
      • Adaptive Time Stepping
      • Length Scale
      • Time Scale
      • Number of Iterations
      • CFL Condition
      • ANSYS CFX Transient Blade Flow
    • ANSYS CFX Data Analysis >
      • ANSYS CFX Moving Mesh >
        • ANSYS Dynamic Mesh
        • ANSYS CFX GGI Interface
      • CFD-Post >
        • CFX Point Parameter
        • ANSYS Parameter Analysis
        • Fatigue Life Optimization with ANSYS nCode DesignLife
        • ANSYS-CFX Probe Tool
        • Exporting Plane Data
        • CFD-Post Parameters
        • Using FFT for Data Analysis
        • Finding Paramters of a Close Surface
        • Histograms in CFD-Post
        • Extracting Data from a Stream Line
        • Transient or Seqence in CFD-Post
        • Linear Data Analysis in CFD-Post
        • CFD-Post Point Cloud
      • ANSYS Making Movies
      • ANSYS CFD-Post Data Loading
      • ANSYS-CFX Data Importing for Analysis
      • ANSYS Loading Simulation Data
      • ANSYS CFX Analysis Tools
      • ANSYS and Excel Data Analysis
      • Using Report Viewer
    • ANSYS Simulation Validation >
      • Flow Validation Around a Cylinder
      • Grid Sensitivity Analysis >
        • ANSYS-FLUENT Time Dependent Boundary Condition
      • Validation Steps
    • ANSYS Blogs >
      • ANSYS HPC
      • ANSYS Stress Analysis >
        • Wing Stress Analysis
        • ANSYS Mechanical APDL
        • Tail Fin Stress Analysis
        • ANSYS Stress Analysis Material >
          • Assigning Boundary Conditions
        • ANSYS Friction Modelling
        • ANSYS Static Structure
        • ANSYS Rigid Dynamics
        • ANSYS Explicit Dynamics
        • ANSYS Assembly Manager >
          • ANSYS Rigid Dynamics
  • MATLAB
    • MATLAB Control Circuits
    • MATLAB GUI
    • MATLAB Fourier Transform >
      • MATLAB Fourier Transform
    • MATLAB Numerical Analysis >
      • Newton Raphson Method
      • MATLAB Derivatives
      • Coding MATLAB EQUATIONS
    • MATLAB APPLICATION >
      • MATLAB Time Dependent Beams
      • MATLAB Hydrulic Circuts Losses
      • MATLAB Sensitivity Analysis
      • MATLAB Water Resources
      • MATLAB Reciprocating Engine
      • MATLAB Ready Polygon Data
      • MATLAB Aerofoil Lift Calculation
      • MATLAB Gas Diffusion
      • MATLAB Studying Drag
      • MATLAB Fuel Droplet Studies
      • MATLAB Atmospheric
      • MATLAB Gas Turbine Code
      • MATLAB Combustion
      • MATLAB Spray Modelling
      • MATLAB Moving Sets of Data
      • MATLAB Applying Non-Slip Conditions
      • MATLAB Fuel Gas Diffusion
      • MATLAB Landing Gear
      • MATLAB Beam Analysis >
        • MATLAB Bending of Plates
      • MATLAB Wind Analysis
      • MATLAB Code for Solar Radiation
      • MATLAB SIMULATION
    • MATLAB Data Analysis >
      • Adding Descriptive Text to Images
      • MATLAB Multiplying Two Functions
      • MATLAB Image Analysis >
        • MATLAB Image Simulation
      • MATLAB Movies >
        • MATLAB Cameras
      • MATLAB Plotting Functions
      • MATLAB PDF Methods
      • MATLAB Adding Two Functions
      • MATLAB Area Segmintation
      • MATLAB Reading Data In and Out
      • MATLAB Functions Written by User
      • COUPLING MATLAB WITH SOFTWARE
    • MATLAB FlOW MODELLING >
      • MATLAB Continuity Equation
      • MATLAB Navier Stokes Equations >
        • Navier Stokes U Velocity in 2D
      • MATLAB Flow Diffusion
      • Gauss Seidel Matlab
      • MATLAB Partical Motion
      • Matlab Gauss Elimination
      • MATLAB Ideal Gas Equation
      • MATLAB Fluid Properties
      • MATLAB Gauss-Seidel Method
      • MATLAB Boundary Layer
      • MATLAB Infinitesimal strain theory
      • MATLAB Stream Functions >
        • MATLAB Studying Vorticity
        • MATLAB 2D Heat Diffusion
      • MATLAB Atmospheric Analysis
      • MATLAB Crank Nicolson
      • Building Codes >
        • MATLAB Data Generation Algorithm
        • MATLAB DNS Subsonic Code
        • MATLAB Runge Kutta
        • MATLAB DNS Sonic Code
        • MATLAB DNS Incompressible Code
      • MATLAB Flow Applications
      • MATLAB Species Concentration
      • MATLAB Wind Flow Analysis
      • MATLAB Turbulence Modelling
      • MATLAB VECTOR FIELD PLOTS >
        • Vector Fields Sites
        • MATLAB Vector Arrow Function
        • MATLAB Vector Gradient
    • MATLAB Working with Different Coordinates >
      • MATLAB Cylindrical Coordinates
      • MATLAB Spherical Coordinates
    • MATLAB Algebric Operations >
      • MATLAB Diagonal Matrix Construction
      • MATLAB Applying Shear to a Box
      • MATLAB Rotating a Set of Points
      • MATLAB Translation
      • MATLAB Scaling a Box
    • MATLAB MESH GENERATION >
      • MATLAB GEOMETRICAL MODELLING
      • MATLAB Geometric Operations
      • MATLAB Mesh Simulation
      • Delaunay Trangulation >
        • Mesh Genration Code Trials
      • MATLAB Uniform Mesh
    • MATLAB PDE Problems >
      • MATLAB Vibrations Modelling >
        • MATLAB Harmonic Motion
        • Molecular vibration
      • MATLAB Solving ODEs
    • MATLAB Reading List
  • SolidWorks
    • AutoDesk 3ds Max
    • AutoCAD
    • Aircraft Design >
      • Aircraft Design Data Base 1
      • Aircraft Structures Books
      • Aircraft Cutaway Drawings
      • Aerodynamics Book Reviews
  • Programs
    • STAR-CCM+ Tutorials
    • FORTRAN90 >
      • The Netlib
      • Salome-Platform
    • CHEMKIN
    • OpenFoam >
      • OpenFOAM Installation
      • OpenFoam Links
    • C++ >
      • C++ Compiling Your First Code
      • C++ Delaunay Triangulation

A Tutorial on Solving ODEs using MATLAB

Picture

Example  Using TDMA

Using TDMA, deriving the TDMA method using MATLAB, solving the problem in a step by step method:

%Forming TDMA matrix
clc
clear
N=5;
A=1;
K=1;
DX=1;
for i=1:N;
aw(i)=(K*A)/DX;
ae(i)=(K*A)/DX;
ap(i)=aw(i)+ae(i);
end

for i=1:N;
for j=1:N;
if ((i==1)&(j==1))

a(i,j)=ap(i);
elseif ((i==1)&(j==2))

a(i,j)=-ae(i);
elseif ((i==1)&(j>2))

a(i,j)=0;
end

if ((i==2)&(j==1))
a(i,j)=-ae(i);
elseif ((i==2)&(j==2))

a(i,j)=ap(i);
elseif ((i==2)&(j==3))

a(i,j)=-aw(i);
elseif ((i==2)&(j>3))

a(i,j)=0;
end

if ((i==3)&(j==1))
a(i,j)=0;
elseif ((i==3)&(j==2))

a(i,j)=-ae(i);
elseif ((i==3)&(j==3))

a(i,j)=ap(i);
elseif ((i==3)&(j==4))

a(i,j)=-aw(i);
elseif ((i==3)&(j>4))

a(i,j)=0;
end

if ((i==4)&(j==1))
a(i,j)=0;
elseif ((i==4)&(j==2))

a(i,j)=0;
elseif ((i==4)&(j==3))

a(i,j)=-ae(i);
elseif ((i==4)&(j==4))

a(i,j)=ap(i);
elseif ((i==4)&(j==5))

a(i,j)=-aw(i);
end

if ((i==5)&(j==1))
a(i,j)=0;
elseif ((i==5)&(j==2))

a(i,j)=0;
elseif ((i==5)&(j==3))

a(i,j)=-0;
elseif ((i==5)&(j==4))

a(i,j)=-aw(i);
elseif ((i==5)&(j==5))

a(i,j)=ap(i);
end
end
end
a


Picture

b-The provided code is a TDMA function. This following code constructs the diagonal matrix in a different method.

function x = TDMAsolver(a,b,c,d)
%a, b, c, and d are the column vectors for the compressed tridiagonal matrix
n = length(b);
% n is the number of rows
% Modify the first-row coefficients
c(1) = c(1) / b(1);
% Division by zero risk.

d(1) = d(1) / b(1);
% Division by zero would imply a singular matrix.
for i = 2:n id = 1 / (b(i) - c(i-1) * a(i));
% Division by zero risk.

c(i) = c(i)* id;
% Last value calculated is redundant.
d(i) = (d(i) - d(i-1) * a(i)) * id;
end
% Now back substitute.
x(n) = d(n);
for i = n-1:-1:1

x(i) = d(i) - c(i) * x(i + 1);
end
end


Other Proposed Methods

Some provided functions in MATLAB  used to solve a set of algebraic equations

1-Cholesky Factorization. Example R=chol(A).

2-LU: Factorization.Example [L,U]=lu(A)

3-QR:Factorization.

Example Coding the Standard Deviation Method for a Set of 1D Velocity and then comparing the Output with the Built in Function in MATLAB

clc
clear
N=10;
for i=1:N+1
mb=randn(1,1);
if mb>1
V(i)=-rand(1,1);
else
V(i)=rand(1,1);
end
end
%Finding the standard deviation direct;y from matlab
s=std(V(:));
%Programing the standard deviation
jh=size(V(:));
mu=sum(V(:))/jh(1);
for i=1:N+1
vl(i)=(V(i)-mu)^2;
end
ss=(sum(vl(:))/(jh(1)-1))^0.5;
%Standard Deviation finished code
%Next finding the Skewness
for i=1:N+1
VVV(i)=(V(i))^3;
end
VV=mean(sum(VVV(:)));
skewness(V(:));
SKEWN=VV/(s)^3;
%Finding the variance
f=var(V(:));

Example Conducting a study on a Heat Transfer Case in a Solid:

Finite Volume 1D Steady Heat Diffusion Studied Case, the code uses TDMA.

clc
clear
K=100;
A=0.01;
TA=100;
TB=500;
L=1;
N=500;
DX=L/N;
for I=1:N
X(I)=I*DX;
end
for I=1:N;
if (I==1)
SU(I)=(2*K*A*TA)/DX;
SP(I)=(-2*K*A)/DX;
AW(I)=0; AE(I)=(K*A)/DX;
AP(I)=AE(I)+AW(I)-SP(I);
AA(I)=-1*AW(I);
BB(I)=AP(I);
CC(I)=-1*AE(I);
DD(I)=SU(I);
elseif (I==N)
SU(I)=(2*K*A*TB)/DX;
SP(I)=(-2*K*A)/DX;
AW(I)=(K*A)/DX;
AE(I)=0;
AP(I)=AE(I)+AW(I)-SP(I);
AA(I)=-1*AW(I);
BB(I)=AP(I);
CC(I)=-1*AE(I);
DD(I)=SU(I);
elseif ((I~=1)|(I~=N))
SU(I)=0;
SP(I)=0;
AW(I)=(K*A)/DX;
AE(I)=(K*A)/DX;
AP(I)=AE(I)+AW(I);
AA(I)=-1*AW(I);
BB(I)=AP(I);
CC(I)=-1*AE(I);
DD(I)=0;
end
end
CC(1)=CC(1)/BB(1);
DD(1)=DD(1)/BB(1);
for I=2:N;
ID(I)=1/(BB(I)-CC(I-1)*AA(I));
CC(I)=CC(I)*ID(I);
DD(I)=(DD(I)-DD(I-1)*AA(I))*ID(I);
end

T(N)=DD(N);
I=N;
for J=1:N-1;
I=I-1;
T(I)=DD(I)-CC(I)*T(I+1);
end
%CHECKING OUTPUT
% AW' % AE' % SU' % SP' % AP'
% T' plot(X,T);
grid on;
xlabel('X');
Ylabel('T');



Finite Volume 1D Steady Heat Diffusion with Heat Source,the code uses TDMA

clc 
clear
K=0.5;
A=1;
TA=100;
TB=200;
L=0.02;
N=5;
DX=L/N;
q=1000000;
for I=1:N;
X(I)=I*DX;
end
for I=1:N;
if (I==1)

SU(I)=(2*K*A*TA)/DX+q*A*DX;
SP(I)=(-2*K*A)/DX;
AW(I)=0;
AE(I)=(K*A)/DX;
AP(I)=AE(I)+AW(I)-SP(I);

AA(I)=-1*AW(I);
BB(I)=AP(I);

CC(I)=-1*AE(I);
DD(I)=SU(I);
elseif (I==N) SU(I)=(2*K*A*TB)/DX+q*A*DX;
SP(I)=(-2*K*A)/DX;
AW(I)=(K*A)/DX;
AE(I)=0;
AP(I)=AE(I)+AW(I)-SP(I);
AA(I)=-1*AW(I);
BB(I)=AP(I);
CC(I)=-1*AE(I);
DD(I)=SU(I);
elseif ((I~=1)|(I~=N))
SU(I)=+q*A*DX;
SP(I)=0;
AW(I)=(K*A)/DX;
AE(I)=(K*A)/DX;
AP(I)=AE(I)+AW(I);
AA(I)=-1*AW(I);
BB(I)=AP(I);
CC(I)=-1*AE(I);
DD(I)=SU(I);
end
end
CC(1)=CC(1)/BB(1);
DD(1)=DD(1)/BB(1);
for I=2:N;
ID(I)=1/(BB(I)-CC(I-1)*AA(I));
CC(I)=CC(I)*ID(I);
DD(I)=(DD(I)-DD(I-1)*AA(I))*ID(I);
end

T(N)=DD(N);
I=N;
for J=1:N-1;

I=I-1;
T(I)=DD(I)-CC(I)*T(I+1);

end
%CHECKING OUTPUT % AW' %AE' %SU' %SP' %AP' %T' % EXACT solution
for i=1:N;
TT(i)=((TB-TA)/L)*X(i)+(0.5*q)*((L-X(i))*X(i))/K+TA;
end
plot(X,T);
grid on;
xlabel('X');
Ylabel('T');
hold on;
plot(X,TT,'-r')


Picture

Finite Volume 1D Steady Convection Diffusion ,the code uses TDMA

clc 
clear
N=5;
FFA=1;
FFB=0;
DX=0.2;
for I=1:N;
X(I)=I*DX;
if (I==1)
DENSITYw(I)=0;
DENSITYe(I)=1;
Uw(I)=0;
Ue(I)=0.1;
Rw(I)=0;
Re(I)=0.1;
DELTAXWP(I)=0.2;
DELTAXPE(I)=0.2;
Fw(I)=DENSITYw(I)*Uw(I);
Fe(I)=DENSITYe(I)*Ue(I);
Dw(I)=Rw(I)/DELTAXWP(I);
De(I)=Re(I)/DELTAXPE(I);
AW(I)=0;
AE(I)=De(I)-0.5*Fe(I);
SP(I)=-(2*De(I)+Fe(I));
AP(I)=AW(I)+AE(I)-SP(I);
Su(I)=(2*De(I)+Fe(I))*FFA;
AA(I)=-1*AW(I);
BB(I)=AP(I);
CC(I)=-1*AE(I);
DD(I)=Su(I);
elseif (I==N) DENSITYw(I)=1;
DENSITYe(I)=0;
Uw(I)=0.1;
Ue(I)=0;
Rw(I)=0.1;
Re(I)=0;
DELTAXWP(I)=0.2;
DELTAXPE(I)=0.2;
Fw(I)=DENSITYw(I)*Uw(I);
Fe(I)=DENSITYe(I)*Ue(I);
Dw(I)=Rw(I)/DELTAXWP(I);
De(I)=Re(I)/DELTAXPE(I);
AW(I)=Dw(I)+0.5*Fw(I);
AE(I)=0;
SP(I)=-(2*Dw(I)-Fw(I));
AP(I)=AW(I)+AE(I)-SP(I);
Su(I)=(2*Dw(I)+Fw(I))*FFB; AA(I)=-1*AW(I); BB(I)=AP(I);
CC(I)=-1*AE(I); DD(I)=Su(I); elseif ((I~=I)|(I~=N)) DENSITYw(I)=1; DENSITYe(I)=1; Uw(I)=0.1; Ue(I)=0.1; Rw(I)=0.1; Re(I)=0.1; DELTAXWP(I)=0.2; DELTAXPE(I)=0.2; Fw(I)=DENSITYw(I)*Uw(I); Fe(I)=DENSITYe(I)*Ue(I); Dw(I)=Rw(I)/DELTAXWP(I); De(I)=Re(I)/DELTAXPE(I); AW(I)=Dw(I)+0.5*Fw(I); AE(I)=De(I)-0.5*Fe(I); AP(I)=AW(I)+AE(I)+(Fe(I)-Fw(I)); SP(I)=0; Su(I)=0; AA(I)=-1*AW(I);
BB(I)=AP(I); CC(I)=-1*AE(I);
DD(I)=Su(I);
end
end
%checking output
% AW' % AE' % Su' % SP' % AP'
CC(1)=CC(1)/BB(1);
DD(1)=DD(1)/BB(1);
for I=2:N;
ID(I)=1/(BB(I)-CC(I-1)*AA(I));
CC(I)=CC(I)*ID(I);
DD(I)=(DD(I)-DD(I-1)*AA(I))*ID(I);
end
FF(N)=DD(N);
I=N;
for J=1:N-1; I=I-1;
FF(I)=DD(I)-CC(I)*FF(I+1);
end
for I=1:N;
t(I)=I*DX;
FFF(I)=(2.7183-exp(t(I)))/1.7183;
end
%THIS IS JUST TO VALIDATE OUTPUT
% PP=[1.55 -0.45 0 0 0; -0.55 1.0 -0.45 0 0; 0 -0.55 1.0 -0.45 0; 0 0 -0.55 1.0 -0.45; 0 0 0 -0.55 1.45];
% PPP=[1.1 0 0 0 0]';
% PPPP=(PP^(-1))*(PPP);
%checking output
% ID= ID'
% CC=CC'
% BB=BB'
% DD=DD'
% AA=AA'
% FF=FF'
plot(X,FF,'--B')
grid on;
xlabel('Distance')
ylabel('FF')
hold on
plot(X,FFF,'-r')
%THIS IS JUST TO VALIDATE OUTPUT
% hold on
% plot(X,PPPP,'--r')


Picture

Finite Volume 1D Heat Diffusion Studied Case, that offers the option to show different heat profiles for a changing temperature boundary the code uses TDMA.

clc
clear
K=100;
A=0.01;
TIME=1;
TIMEA=1;
TIMEB=1;
TA(TIMEA)=100;
TB(TIMEB)=500;
L=1; N=5;
DX=L/N;
for I=1:N
X(I)=I*DX;
end
TIMEA=0;
for TIME=1:10; TIMEA=TIMEA+1;
TA(TIMEA+1)=TA(TIMEA)+130;
for I=1:N;
if (I==1) SU(I,TIME)=(2*K*A*TA(TIMEA))/DX;
SP(I,TIME)=(-2*K*A)/DX; AW(I,TIME)=0;
AE(I,TIME)=(K*A)/DX; AP(I,TIME)=AE(I,TIME)+AW(I,TIME)-SP(I,TIME);
AA(I)=-1*AW(I,TIME); BB(I)=AP(I,TIME);
CC(I)=-1*AE(I,TIME); DD(I)=SU(I,TIME);
elseif (I==N) SU(I,TIME)=(2*K*A*TB(TIMEB))/DX;
SP(I,TIME)=(-2*K*A)/DX; AW(I,TIME)=(K*A)/DX;
AE(I,TIME)=0; AP(I,TIME)=AE(I,TIME)+AW(I,TIME)-SP(I,TIME);
AA(I)=-1*AW(I,TIME);
BB(I)=AP(I,TIME);
CC(I)=-1*AE(I,TIME);
DD(I)=SU(I,TIME);
elseif ((I~=1)|(I~=N))
SU(I,TIME)=0; SP(I,TIME)=0;
AW(I,TIME)=(K*A)/DX;
AE(I,TIME)=(K*A)/DX;
AP(I,TIME)=AE(I,TIME)+AW(I,TIME);
AA(I)=-1*AW(I,TIME);
BB(I)=AP(I,TIME);
CC(I)=-1*AE(I,TIME);
DD(I)=0; end
end
CC(1)=CC(1)/BB(1);
DD(1)=DD(1)/BB(1);
for I=2:N;
ID(I)=1/(BB(I)-CC(I-1)*AA(I));
CC(I)=CC(I)*ID(I);
DD(I)=(DD(I)-DD(I-1)*AA(I))*ID(I);
end
T(N,TIME)=DD(N);
I=N;
for J=1:N-1;
I=I-1;
T(I,TIME)=DD(I)-CC(I)*T(I+1,TIME);
end
%CHECKING OUTPUT
% AW'
% AE'
% SU'
% SP'
% AP'
% T'
plot(X,T);
grid on;
xlabel('X');
Ylabel('TA');
hold on
end


Picture

Finite Volume 1D Unsteady Heat Diffusion Studied Case Using Crank Nicholson, the Code Uses TDMA.

NOTE: The Method will be Validated Soon.
clc
clear
K=10;
A=0.01;
TIME=1;
TIMEA=1;
TIMEB=1;
TA(TIMEA)=200;
TB(TIMEB)=-100*TIMEB+500;
L=0.02;
N=3;
DX=L/N;
THETA=0.5;
HEAT_CAPACITY=1000000;
DENSITY=10;

for I=1:N
X(I)=I*DX;
end
%-------------------------------------------------------------------------- %-INITILIZATION------------------------------------------------------------
TIME=1;
for I=1:N;
if (I==1)
SU(I,TIME)=(2*K*A*TA(TIMEA))/DX;
SP(I,TIME)=(-2*K*A)/DX;
AW(I,TIME)=0;
AE(I,TIME)=(K*A)/DX;
AP(I,TIME)=AE(I,TIME)+AW(I,TIME)-SP(I,TIME);
AA(I)=-1*AW(I,TIME);
BB(I)=AP(I,TIME);
CC(I)=-1*AE(I,TIME);
DD(I)=SU(I,TIME);
elseif (I==N)
SU(I,TIME)=(2*K*A*TB(TIMEB))/DX;
SP(I,TIME)=(-2*K*A)/DX;
AW(I,TIME)=(K*A)/DX; AE(I,TIME)=0;
AP(I,TIME)=AE(I,TIME)+AW(I,TIME)-SP(I,TIME);
AA(I)=-1*AW(I,TIME);
BB(I)=AP(I,TIME);
CC(I)=-1*AE(I,TIME);
DD(I)=SU(I,TIME);
elseif ((I~=1)|(I~=N))

SU(I,TIME)=0;
SP(I,TIME)=0;
AW(I,TIME)=(K*A)/DX;
AE(I,TIME)=(K*A)/DX;
AP(I,TIME)=AE(I,TIME)+AW(I,TIME);
AA(I)=-1*AW(I,TIME);
BB(I)=AP(I,TIME);
CC(I)=-1*AE(I,TIME);
DD(I)=0;
end
end
CC(1)=CC(1)/BB(1);
DD(1)=DD(1)/BB(1);
for I=2:N;
ID(I)=1/(BB(I)-CC(I-1)*AA(I));
CC(I)=CC(I)*ID(I);
DD(I)=(DD(I)-DD(I-1)*AA(I))*ID(I);
end
T(N,TIME)=DD(N);
I=N;
for J=1:N-1;
I=I-1;
T(I,TIME)=DD(I)-CC(I)*T(I+1,TIME);
end
%CHECKING OUTPUT
% AW'
% AE'
% SU'
% SP'
% AP'
% T'
plot(X,T);
grid on;
xlabel('X');
Ylabel('TA');
hold on
%-------------------------------------------------------------------------- %-END_INITILIZATION--------------------------------------------------------
I=1;
TIME=1;
DELTA_TIME=((DX)^2)*(DENSITY*HEAT_CAPACITY)/(2*K);
%STARTING_TIME_STEPPING----------------------------------------------------
TIMEB=0
FINAL_TIME=4;
for TIME=1:FINAL_TIME;
TIMEB=TIMEB+1; TB(TIMEB)=-100*TIMEB+500;
for I=1:N;
if (I==1) SU(I,TIME+1)=(2*K*A*TA(TIMEA))/DX;
SP(I,TIME+1)=(-2*K*A)/DX; AW(I,TIME+1)=0;
AE(I,TIME+1)=(K*A)/DX;
AP0(I,TIME)=(DX*DENSITY*HEAT_CAPACITY)/DELTA_TIME; AP(I,TIME+1)=THETA*(AW(I,TIME+1)+AE(I,TIME+1))+AP0(I,TIME); AA(I)=-1*THETA*AW(I,TIME+1);
BB(I)=AP(I,TIME+1); CC(I)=-1*THETA*AE(I,TIME+1);
DD(I)=SU(I,TIME+1)+(1-THETA)*(T(I+1,TIME))+AP0(I,TIME)*T(I,TIME)-(1-THETA)*((AW(I,TIME+1))*T(I,TIME)+ (AE(I,TIME+1))*(T(I,TIME))); elseif (I==N)
SU(I,TIME+1)=(2*K*A*TB(TIMEB))/DX;
SP(I,TIME+1)=(-2*K*A)/DX; AW(I,TIME+1)=(K*A)/DX; AE(I,TIME+1)=0; AP0(I,TIME)=(DX*DENSITY*HEAT_CAPACITY)/DELTA_TIME;
AP(I,TIME+1)=THETA*(AW(I,TIME+1)+AE(I,TIME+1))+AP0(I,TIME);
AA(I)=-1*THETA*AW(I,TIME); BB(I)=AP(I,TIME);
CC(I)=-1*THETA*AE(I,TIME);
DD(I)=SU(I,TIME+1)+(1-THETA)*(T(I-1,TIME))+AP0(I,TIME)*T(I,TIME)-(1-THETA)*((AW(I,TIME+1))*T(I,TIME)+ (AE(I,TIME+1))*(T(I,TIME)));
elseif ((I~=1)|(I~=N))
SU(I,TIME+1)=0;
SP(I,TIME+1)=0;
AW(I,TIME+1)=(K*A)/DX;
AE(I,TIME+1)=(K*A)/DX;
AP0(I,TIME)=(DX*DENSITY*HEAT_CAPACITY)/DELTA_TIME; AP(I,TIME+1)=THETA*(AW(I,TIME+1)+AE(I,TIME+1))+AP0(I,TIME); AA(I)=-1*AW(I,TIME+1);
BB(I)=AP(I,TIME+1); CC(I)=-1*AE(I,TIME+1); DD(I)=(1-THETA)*(T(I+1,TIME)+T(I-1,TIME))+AP0(I,TIME)*T(I,TIME)-(1-THETA)*((AW(I,TIME+1))*T(I,TIME)+ (AE(I,TIME+1))*(T(I,TIME)));
end
end
CC(1)=CC(1)/BB(1);
DD(1)=DD(1)/BB(1);
for I=2:N;
ID(I)=1/(BB(I)-CC(I-1)*AA(I));
CC(I)=CC(I)*ID(I);
DD(I)=(DD(I)-DD(I-1)*AA(I))*ID(I);
end
T(N,TIME+1)=DD(N);
I=N;
for J=1:N-1;
I=I-1;
T(I,TIME+1)=DD(I)-CC(I)*T(I+1,TIME+1);
end
plot(X,T);
end
Picture

Using MATLAB Functions to Solve ODEs

clc
clear
eqn2= 'D3y+3*D2y+3*Dy+2*y=x^3';
init2= 'y(0)=0,Dy(0)=1,D2y(0)=1';
y=dsolve(eqn2,init2,'x')
x=linspace(-1,1,10);
z=eval(vectorize(y));
plot(x,z,'-*')
grid on
Picture

Unless otherwise noted, all content on this site is @Copyright by Ahmed Al Makky 2012-2013 - http://cfd2012.com
Web Hosting by Just Host