CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > General Forums > Main CFD Forum

Block Tridiagonal Matlab code

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   April 25, 2020, 20:58
Default Block Tridiagonal Matlab code
  #1
New Member
 
RJ
Join Date: Apr 2020
Posts: 1
Rep Power: 0
mohsen1985 is on a distinguished road
Hello

I am trying to write a simple code to solve a block tridiagonal matrix with Tomas method
I developed the code same as a simple tridiagonal matrix but I don't know why my code is going wrong
the Matlab code is attached here
%-------------------------------------------------------------
clear all
clc;
y_0 = 0;
h = 40e-3;
y_h = y_0 + h;
dy = 1e-2;
U0 = 40;
nu = 217e-6;
y = y_0:dy:y_h;
Ny = length(y);
dt_B = 0.01;
t_end = 1.08;
t_B = 0:dt_B:t_end;
Nt_B = length(t_B);
d_B = (nu*dt_B/dy);
v_B1 = ones(1,2*Ny);
for j=1:Ny-1
v_B1(2*j)= -dy/2;
end
v_B1(2*Ny)=0;

v_B2 = ones(1,2*Ny-1);v_B2(1) = 0;
for j=2:Ny
v_B2(2*j-1) = -d_B;
end
v_B3 = -1*ones(1,2*Ny-1);v_B3(2*Ny-1)= 1;
for j=2:Ny
v_B3(2*(j-1)) = d_B;
end
v_B4 = zeros(1,2*Ny-2);
for j=2:Ny
v_B4(2*(j-1)) = -dy/2;
end
v_B5 = ones(1,2*Ny-2);
for j=2:Ny
v_B5(2*(j-1)) = 0;
end

v_B = diag(v_B1) + diag(v_B2,1) + diag(v_B3,-1)+diag(v_B4,2)+diag(v_B5,-2);
v_B_new = diag(ones(2*Ny,1));

uv_B = zeros(2*Ny,Nt_B);
%B.C
uv_B(2*Ny,= 0;
uv_B(1,= U0;
R_B = zeros(2*Ny,1);
R_B(1) = U0;
R_B(2*Ny)=0;
R_B_new = zeros(2*Ny,1);


k=1;
i=1;
A_old(1,1) = v_B(2*i-1,2*i-1);
A_old(1,2) = v_B(2*i-1,2*i);
A_old(2,1) = v_B(2*i,2*i-1);
A_old(2,2) = v_B(2*i,2*i);
%--------------------------------
R_old(1,1) = U0;
R_old(2,1) = 0;
%--------------------------------
R_new = A_old\R_old;
%--------------------------------
R_B_new(2*i-1) = R_new(1,1);
R_B_new(2*i,1) = R_new(2,1);
%--------------------------------
for k = 1:Nt_B-1
for j=2:Ny-1
R_B(2*j-1)= uv_B(2*j-1,k) + uv_B(2*(j-1)-1,k) + d_B*(uv_B(2*j,k) - uv_B(2*(j-1),k));
end
for i=2:Ny
%i = 2;
%get B
B_old(1,1) = v_B(2*(i-1)-1,2*i-1);
B_old(1,2) = v_B(2*(i-1)-1,2*i);
B_old(2,1) = v_B(2*(i-1),2*i-1);
B_old(2,2) = v_B(2*(i-1),2*i);
%get C
C_old(1,1) = v_B(2*i-1,2*(i-1)-1);
C_old(1,2) = v_B(2*i-1,2*(i-1));
C_old(2,1) = v_B(2*i,2*(i-1)-1);
C_old(2,2) = v_B(2*i,2*(i-1));
%get R
R_old(1,1) = uv_B(2*i-1,k) + uv_B(2*(i-1)-1,k) + d_B*(uv_B(2*i,k) - uv_B(2*(i-1),k));
R_old(2,1) = 0;
A_new(1,1) = v_B(2*i-1,2*i-1);
A_new(1,2) = v_B(2*i-1,2*i);
A_new(2,1) = v_B(2*i,2*i-1);
A_new(2,2) = v_B(2*i,2*i);
%----------------------------------
B_new = A_old\B_old;
%----------------------------------
v_B_new(2*(i-1)-1,2*i-1) = B_new(1,1);
v_B_new(2*(i-1)-1,2*i) = B_new(1,2);
v_B_new(2*(i-1),2*i-1) = B_new(2,1);
v_B_new(2*(i-1),2*i) = B_new(2,2);
%----------------------------------
A_old = A_new - C_old*B_new;
R_new = A_old\(R_old - C_old*R_new);
%--------------------------------
R_B_new(2*i-1) = R_new(1,1);
R_B_new(2*i,1) = R_new(2,1);
end
i = Ny;
uv_B(2*i-1,k+1)= R_B_new(2*i-1);
uv_B(2*i,k+1)= R_B_new(2*i);
for i = Ny-1:-1:1
uv_B(2*i-1,k+1)= R_B_new(2*i-1) - v_B_new(2*i-1,2*(i+1)-1)*R_B_new(2*(i+1)-1)-...
v_B_new(2*i-1,2*(i+1))*R_B_new(2*(i+1));
uv_B(2*i,k+1)= R_B_new(2*i) - v_B_new(2*i,2*(i+1)-1)*R_B_new(2*(i+1)-1)-...
v_B_new(2*i,2*(i+1))*R_B_new(2*(i+1));
end
end
mohsen1985 is offline   Reply With Quote

Old   April 26, 2020, 05:55
Default
  #2
Senior Member
 
sbaffini's Avatar
 
Paolo Lampitella
Join Date: Mar 2009
Location: Italy
Posts: 1,850
Blog Entries: 29
Rep Power: 36
sbaffini will become famous soon enoughsbaffini will become famous soon enough
Send a message via Skype™ to sbaffini
I haven't tried it, but you may want to give a look here:

https://it.mathworks.com/matlabcentr...iagonal-solver
sbaffini is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
MATLAB fractional step code Darren Main CFD Forum 7 September 25, 2018 14:58
matlab code for fluid flow shilpa sood Main CFD Forum 1 February 24, 2015 09:09
Fluent FFT functionality equivalent MATLAB code Mojtaba.a FLUENT 0 January 27, 2015 07:56
Matlab code for the numerical solution of a Prandtl-Meyer expansion wave flow field ivan_CFD Main CFD Forum 4 March 8, 2012 15:17
code for solving block tridiagonal system cfd Main CFD Forum 22 December 24, 2010 17:15


All times are GMT -4. The time now is 20:10.