Some help with a Fortran source code

[bigduke88]

Hello everyone
I would like to ask you for some help with a Fortran source code.
I have to make a program, which solves numerically an ordinary differential equation and after that an integral.
The equation, is going to be solved using the finite difference approximation.
This is equation "developed" in finite differences:
Ui+1 = 0.998*h - 0.7366*(Ui)^2*Cd*h + Ui
This is an equation, which describes the velocity of a moving spherical particle.
"Cd" is a drag coefficient and is function of U:
And this is the equation, which describes the drag coefficient:
Cd = (0.087108/U) + (0.622126/U^0.28)

I'm quite new in any use of an algorythmic "language", but I have to solve the equation and then an intergral in Fortran.
Could you help with some advice, how to do this?
Thank you!!!

 

[bigduke88]

Here is the code, I've tried to fix the things, but no success
The text file is still empty

IMPLICIT DOUBLE PRECISION (A-H,O-Z)
PARAMETER (NR=100)
DIMENSION U(NR)
open(UNIT=13,STATUS='new',file='dU_by_dT.txt')
NIT=2500
jprint=500
KPRINT=1
DT1=1D-8
DT=DT1
DTau=4.D0*DT
Dp=0.00271D0
Ro=998.2D0
Mu=0.001005D0
Diff=2.97271D-9
Tau=(4.D0*Diff*T)/(Dp**2.D0)
T=(Tau*(Dp**2.D0))/(Diff*4.D0)
imax=NR
U(i)=0.
T=0.
DO 100 i=1,imax
T=T+DT
Tau=Tau+DTau
Re=(Dp*U(i)*Ro)/(Mu)
Cd=(24/Re)*(1+0.125*(Re**0.72))
U(i+1)=U(i)*(1-0.7366*U(i)*Cd*DT)+0.998*DT
C DISTANCE
Z=Z+DT*(U(I)+U(I+1))

100 continue


write(UNIT=13,FMT=*) 'Z,T,U'
write(unit=13,FMT=*) Z,T,U
STOP
END

 

[salgerman]

bigduke:

You program has a few problems.

Line 15, you calculate Tau as funtion of T, which has not been initialized, so T=0 => Tau=0
Line 16, you calculate T as function of Tau, Tau=0 => T=0
Line 18, U(i) = 0.0; but i has not been initialized, so i=0, U(0)=0
Line 24, for i=1, Re is 0.0 because U(i=1) has not been initialized and so U(1)=0
Line 25, Cd calculation will divide by zero because of 24/Re

Also, Z and T are scalars, if you care to see their intermediate values, you may want to move the write statement into the loop and write Z,T,U(i)

You may want to move the header line to before the loop.

If you just wants to see the final Z and T and the entire U vector, you can leave the write statements where they are.

....and I am afraid I am running out of good mood, here.

 

[bigduke88]

Thank you for the help!!!
Please don't be running out fo good mood, obviously the "art of programming" is beyond my skills, at this moment, but hpefully I'll advance in it.... but the progress is really slow...
Best regards!!!

 

[gummibaer]

And...

In your do-loop (line 20) i runs to imax which is Nr which is 100
You reference u(i+1) in line 21, whereas u is diemnsioned to NR = 100. This will give an 'array bounds exceeded' error.

You should set

iMax = NR - 1

Norbert


The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.

 

[PaulWDent]

If you have never programmed before, it is good to start with solving the problem with an EXCEL spreadsheet. Your problem is very appropriate to solve on an Excel spread sheet. You jsut have columns labelled T, U, Cd, Z etc along the top, (columsn A, B, C, Din axecl) and then you start with 0.0 in box A2, put the initial velocity in box B2 and I guess the initial integral value Z is zero in box D2.
Then in boc A3 you write "=A1 + 1.0E-8" which adds dT (do you really want 10 nanoseconds!)
Then copy that box A3 to A4 to A101 and you have your T's.
Then in box C2 write the formula that calculates the drag coeffcient Cd from the velocity in box B2
Then in box B3 you write the equation that updates the initial velocity at T=0.0 to the
velociy at T+dT based on the drag coefficient in box C2. Then copy and paste from boc C2 to C3.

In box D3 you write ="D2+1.0E-8*(C2+C3)/2.0 which is doing trapezoidal rule integration for Z

Now just copy and paste from boxes A3,B3,C3,D3 into boxes A4,B4,C4,D4 to 101 and you
will have your answer.

You will find it a much easier step from there to FORTRAN

Now one of the big problems with your thinking is this: Drag reduces velocity, so the particle is slowing. But from what intial velocity? You have set the initial velocity to zero!!
You must know the initial velocity to solve this problem.
You must write, before the Do-loop, U(1) = something positive and probably large.

Alternatively, you could start with a very small velocity (not zero) and work in time-reversed order to find out what the initial velocity would have to be 10uS ago to have a terminal velocity of say 1m/sec 10uS later.
Another thing you must do for physical problems is to work in real engineering units, i.e. all quantities should be in the meter-kilo-second units, as otherwise you cannot write Acceleration = Force/Mass without some fudge factors. Also, with a Do-loop of 100 and dT = 1e-8, you do realize that your simulation is only simulating 10uS of real time, do you?

 

[gummibaer]

Another thing you must do for physical problems is to work in real engineering units

Depends on the point of view.
If you are dealing with physical scientific problems you should work with dimension-free numbers to do your calculations. I am not a specialist for two phase flow but I would guess Reynolds-Number, Grashof-Number and maybe a few others would be present in this problem when you did it.

If you are dealing with an engineering problem or this is a side-aspect of your task only, then you most probably will do your computations in engineering units.

Norbert




The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.