Translation of Fortran code to related Matlab code in Internal Fluid Dynamics
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
Hi,
I am a kind of new to Matlab. My supervisor gave me a task in Fortran. I am supposed to transfer the code into a Matlab code. The matlab code that I have written is:
clear all
close all
clc
dt=0.001;
t = 0:dt:3;
delta_x=0.01; %step size
x=0:delta_x:100;
rmass=0.5;
rmomi=1.0;
% dx=1/(ii-1);
pi=4.*atan(1);
fik=0.4;
H1=0.5;
H1D=0;
A=0;
AD=0.1;
ADinit=AD;
c1b=0.5;
c2b=1-c1b;
ub = cell(2, 1);
ini_cond = [0,0];
for i = 1:2;
ub{i} = zeros(1,length(x));
ub{i}(:, length(x)) = ini_cond(i) * rand(1, 1);
end
for i=1:length (x);
% x=i*delta_x ;
fikness=fik*sin(pi.*x);
ub{1}=(c1b-H1D*(x-0.5)+AD/2.*(x-0.5).^2)./(H1-0.5*fikness-A*(x-0.5));
ub{2}=(c2b+H1D*(x-0.5)-AD/2.*(x-0.5).^2)./(1.-H1+0.5*fikness+A*(x-0.5));
end
c1=c1b;
c2=1-c1;
u = cell(2, 1);
ini_cond = [0,0];
for i = 1:2;
u{i} = zeros(1,length(x));
u{i}(:, length(x)) = ini_cond(i) * rand(1, 1);
end
for i=1:length (x);
% x=i*delta_x ;
fikness=fik*sin(pi.*x);
u{1}=(c1-H1D*(x-0.5)+AD/2.*(x-0.5).^2)./(H1-0.5*fikness-A*(x-0.5));
u{2}=(c2+H1D*(x-0.5)-AD/2.*(x-0.5).^2)./(1.-H1+0.5*fikness+A*(x-0.5));
end
% p = cell(2, 1);
% q = cell(2, 1);
% % ini_cond = [0,0];
%
% for i = 1:2;
% p{i} = zeros(1,length(t));
% p{i}(:, length(t)) = ini_cond(i) * rand(1, 1);
% q{i} = zeros(1,length(t));
% q{i}(:, length(t)) = ini_cond(i) * rand(1, 1);
% end
p{1}(1)=0.5*(1.-u{1}(1).^2);
q{1}(1)=0;
for i=2:length(x)
q{1}(i)=q{1}(i-1)-delta_x*(u{1}(i-1)-ub{1}(i-1))./dt;
p{1}(i)=0.5*(1.-u{1}(i).^2)+q{1}(i);
end
p{2}(1)=0.5*(1.-u{2}(1).^2);
q{2}(1)=0;
for i=2:length(x)
q{2}(i)=q{2}(i-1)-delta_x*(u{2}(i-1)-ub{2}(i-1))./dt;
p{2}(i)=0.5*(1.-u{2}(i).^2)+q{2}(i);
end
sum = cell(2, 1);
st = cell(2, 1);
ini_cond = [0,0];
for i = 1:2;
sum{i} = zeros(1,length(t));
sum{i}(:, length(t)) = ini_cond(i) * rand(1, 1);
st{i} = zeros(1,length(t));
st{i}(:, length(t)) = ini_cond(i) * rand(1, 1);
end
% for j=1:length (t);
% st(j)=p{1}(j)-p{2}(j);
% end
dc=0.001;
c1=c1+dc;
c1=(c1*st(1)-(c1-dc)*st(2))/(st(1)-st(2));
sum{1}(1)=0.5*(p{2}(1)-p{1}(1));
sum{2}(1)=0.5*(p{2}(1)-p{1}(1)).*(-1/2);
ADDOT(1)=sum{2}(1)*delta_x./rmomi;
H1DDOT(1)=-sum{1}(1).*delta_x./rmass;
H1(i)=H1(i-1)+dt*H1D(i-1);
ADDOT(i)=sum{2}(i)*delta_x./rmomi;
AD(i)=AD(i-1)+dt*ADDOT(i-1);
A(i)=A(i-1)+dt*AD(i-1);
for i=2:length(t)
% x=(i-1)*dx;
sum{1}(i)=sum{1}(i-1)+(p{2}(i)-p{1}(i));
sum{2}(i)=sum{2}(i-1)+(p{2}(i)-p{1}(i)).*(x-1/2);
H1DDOT(i)=-sum{1}(i).*delta_x./rmass;
H1D(i)=H1D(i-1)+dt*H1DDOT;
H1(i)=H1(i-1)+dt*H1D(i-1);
ADDOT(i)=sum{2}(i)*delta_x./rmomi;
AD(i)=AD(i-1)+dt*ADDOT(i-1);
A(i)=A(i-1)+dt*AD(i-1);
end
H1L=H1+A*1/2;
H1R=H1-A*1/2;
H2=1.-H1;
rat1=AD/ADinit;
rat2=0;
rat2=ADDOT/AD;
for i=1:length(x)
ub{1}(i)=u{1}(i);
ub{2}(i)=u{2}(i);
c1b=c1;
end
figure (1);
subplot(1,2,1);
plot(t,H1D,'.:',t,H1DDOT,'.:',t,c2,'*');
axis([0 1 0 2.5]);
xlabel('time');
ylabel('body leading and trailing edges');
hold on
subplot(1,2,2);
plot(t,A,t,AD,t,ADDOT);
xlabel('time');
ylabel('the angles');
figure (2);
subplot(1,2,1);
plot(x,u{1},x,u{2},':', 'LineWidth',1.5);
axis([0 10 0 2.5]);
xlabel('x lateral position');
ylabel('velocity profiles u1, u2 in the gaps');
hold on
subplot(1,2,2);
plot(x,p{1},':',x, p{2},':m', 'LineWidth', 2);
axis([0 10 -3 1]);
xlabel('x lateral position');
ylabel('pressure p1, p2 in the gaps');
hold on
suptitle('1 grain in a fluid (N=2)');
But 'sum' and 'st' parts (Newton-Raphson Method) do not work.
The fortran code is this:
program grab3
!!from grair2.f90 on 09.01.2009
!!1 grain in air between 2 walls: inviscid
!!iterates more than grair1
!!Here A = - theta
implicit double precision (a-h,o-z)
common u1(1001),u2(1001),u1b(1001),u2b(1001)&
,p1(1001),q1(1001),p2(1001),q2(1001),st(9)
ii=101
b=1.
dx=1./(ii-1.)
dt=0.0001
tmax=3.000000001
tmax=2.3999999999
!!tmax=4.500000001
!!tmax=dt*1.000001
nt=0
ntp=1000
ntp1=200
tstart=1.8999999999999
!!tstart=4.2999999999999
rmass=1.
rmomi=2.
!!rmomi=1.
!!rmass=10.
!!rmomi=20.
pi=4.*atan(1.)
!!thickness-effect parameter --actually camber in this prog
fik=0.4
fik=-0.4
t=0.
H1=0.5
H1D=0.
A=0.
AD=0.1
!!AD=0.001
ADinit=AD
c1b=0.5
c2b=1.-c1b
do 1 i=1,ii
x=(i-1.)*dx
fikness=fik*sin(pi*x)
u1b(i)=(c1b-H1D*(x-b/2.)+AD/2.*(x-b/2.)**2)/(H1-0.5*fikness-A*(x-b/2.))
u2b(i)=(c2b+H1D*(x-b/2.)-AD/2.*(x-b/2.)**2)/(1.-H1+0.5*fikness+A*(x-b/2.))
!!note sign above
1 continue
dc=0.001
88 t=t+dt
nt=nt+1
mmm=1
c1=c1b
666 m=1
77 c2=1.-c1
do 3 i=1,ii
x=(i-1.)*dx
fikness=fik*sin(pi*x)
u1(i)=(c1-H1D*(x-b/2.)+AD/2.*(x-b/2.)**2)/(H1-0.5*fikness-A*(x-b/2.))
u2(i)=(c2+H1D*(x-b/2.)-AD/2.*(x-b/2.)**2)/(1.-H1+0.5*fikness+A*(x-b/2.))
3 continue
p1(1)=0.5*(1.-u1(1)**2)
q1(1)=0.
do 4 i=2,ii
q1(i)=q1(i-1)-dx*(u1(i-1)-u1b(i-1))/dt
p1(i)=0.5*(1.-u1(i)**2)+q1(i)
4 continue
p2(1)=0.5*(1.-u2(1)**2)
q2(1)=0.
do 44 i=2,ii
q2(i)=q2(i-1)-dx*(u2(i-1)-u2b(i-1))/dt
p2(i)=0.5*(1.-u2(i)**2)+q2(i)
44 continue
!!write(6,*)m,m,m,m
st(m)=p1(ii)-p2(ii)
!!write(6,*)m,m
m=m+1
if(m.eq.3)go to 51
if(m.eq.4)go to 52
c1=c1+dc
go to 77
51 c1=(c1*st(1)-(c1-dc)*st(2))/(st(1)-st(2))
go to 77
52 continue
mmm=mmm+1
if(mmm.eq.2)go to 666
!!write(6,*)m,m,m,m
!!really need to iterate more
!!write(6,*)c1,st(3)
sum1=0.5*(p2(1)-p1(1))
sum2=0.5*(p2(1)-p1(1))*(-b/2.)
!!write(6,*)m,m,m,m
do 93 i=2,ii-1
x=(i-1.)*dx
sum1=sum1+(p2(i)-p1(i))
93 sum2=sum2+(p2(i)-p1(i))*(x-b/2.)
!!write(6,*)m,m,m,m
H1DDOT=-sum1*dx/rmass
H1D=H1D+dt*H1DDOT
H1=H1+dt*H1D
ADDOT=sum2*dx/rmomi
AD=AD+dt*ADDOT
A=A+dt*AD
!!write(6,*)m,m,m,m,m
!!if(((nt/ntp)*ntp).ne.nt)go to 5499
if(((nt/ntp)*ntp).ne.nt)go to 7021
write(6,*)t,H1,A,st(3)
H1L=H1+A*b/2.
H1R=H1-A*b/2.
H2=1.-H1
rat1=AD/ADinit
rat2=0.
if(abs(AD).lt.0.0000001)go to 9753
rat2=ADDOT/AD
9753 continue
write(5,*)t,H1L,H1R,c1
close(5)
open(5,file='g3bD.dat',status='unknown',position='append')
write(5,*)t,A,AD,ADDOT
close(5)
7021 continue
!!if(ii.gt.0)go to 5499
if(t.lt.tstart)go to 7023
if(((nt/ntp1)*ntp1).ne.nt)go to 7023
do 7022 i=1,ii
x=(i-1.)*dx
fikness=fik*sin(pi*x)
gap1=(H1-0.5*fikness-A*(x-b/2.))
gap2=(1.-H1+0.5*fikness+A*(x-b/2.))
open(5,file='g3cD.dat',status='unknown',position='append')
write(5,*)x,gap1,gap2
close(5)
7022 continue
7023 continue
if(ii.gt.0)go to 5499
if(t.lt.tstart)go to 5499
if(((nt/ntp1)*ntp1).ne.nt)go to 5499
do 702 i=1,ii
x=(i-1.)*dx
open(5,file='gibaS.dat',status='unknown',position='append')
write(5,*)x,u1(i),u2(i)
close(5)
open(5,file='gibbS.dat',status='unknown',position='append')
write(5,*)x,p1(i),p2(i)
close(5)
702 continue
5499 if(t.gt.tmax)go to 549
do 801 i=1,ii
u1b(i)=u1(i)
u2b(i)=u2(i)
801 continue
c1b=c1
go to 88
549 continue
stop
end
Could you please please please help me..
0 comentarios
Respuestas (1)
Walter Roberson
el 3 de En. de 2014
naming a variable "sum" almost always leads to trouble. Avoid naming variables the same as commonly-used MATLAB library routines.
"clear all" should only be used by experts.
Ver también
Categorías
Más información sobre Quantum Mechanics en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!