Simulating Contact Force in a Moving Dashpot System Using FORTRAN

[grecko94]

Hi there I am really new into programiing thing, and I am trying to make program of this problem usin FORTRAN

I want the output to be Time,poisiton,velocity and contact force

and i already know

= amplitude

U = velocity of the car

d = rail distance

and so far this is my program <Moderator's note: please use code tags>

program dashpot

implicit none

integer,parameter :: N=10001,IN=300

integer :: i

real :: pi,dt,k,m,c,U,d,wb,am

real,dimension (N)::T,V,X

!initial condition= position,velocity,time, pi

pi = 3.14159265359

dt = 2*pi/FLOAT(N-1)

x(1) = 0

t(1) = 0

v(1) = 1

k = 3533.66

m = 130

c = 0

U = 8

d = 10

wb = 2*pi*U/D

am =0.5
!looping

do I=1, N-1

T(i+1) = dt*i

X(i+1) = x(i) + v(i)*dt

V(i+1) = v(i) + (-k/m*(x(i)-am*sin(wb*u*dt))-c/m*(v(i)-wb*U*am*cos(wb*dt)))*dt

end do

print*, '      T(I)        X(I)         V(I)'

WRITE (*,"(3F16.8)") (T(I),X(I),V(I),I=1,N,IN)

END PROGRAM dashpot

But I am still have no idea how do i make contact force loop
Can someone help me please...

 

[BvU]

But i still have no idea how do i make contact force loop

What makes you say that ? You have the equation of motion in there, isn't it ?

First thing I notice is the different arguments for the sine and the cosine -- one of the two is wrong.

So what was it you expected to get ?

 

[grecko94]

since F = -K(Y-Yo) -C(y°-yo°)
then a = -k/m(Y-yo)-C(y°-yo°)

isn't the code supposed to be

-k/m*(x(i)-am*sin(wb*u*dt))-c/m*(v(i)-wb*U*am*cos(wb*dt))

and maybe i should make the initial condition velocity of spring to zero

 

[grecko94]

for the contact force I'm thinking about something like this

F = -k*(x(i)-am*sin(wb*u*dt))-c*(v(i)-wb*U*am*cos(wb*dt))

but i know something is off with this code

 

[grecko94]

okay it`s seem, I am able make contact force output

program dashpot

implicit none

integer,parameter :: N=10001,IN=300

integer :: i

real :: pi,dt,k,m,c,U,d,wb,am

real,dimension (N)::T,V,X,f

!initial condition= position,velocity,time, pi

pi = 3.14159265359

dt = 2*pi/FLOAT(N-1)

x(1) = 0

t(1) = 0

v(1) = 0

k = 3281.27

m = 130

c = 0

U = 8

d = 10

wb = 2*pi*U/D

am =0.5
!looping

do I=1, N-1

T(i+1) = dt*i

X(i+1) = x(i) + v(i)*dt

V(i+1) = v(i) + (-k/m*(x(i)-am*sin(wb*u*dt))-c/m*(v(i)-wb*U*am*cos(wb*u*dt)))*dt
f(i) = -k*(x(i)-am*sin(wb*u*dt))-c*(v(i)-wb*U*am*cos(wb*u*dt))end do

print*, '      T(I)        X(I)         V(I)'

WRITE (*,"(4F16.8)") (T(I),X(I),V(I),f(I),I=1,N,IN)

END PROGRAM dashpot

what do you think about that? i know its not enough, but i still can't figured out when you said my argument for sin and consine are wrong.

First thing I notice is the different arguments for the sine and the cosine -- one of the two is wrong.

 

[BvU]

I see, you want an extra printout of the force.
Why do you multiply wd dt with u ?
(perhaps you should explain the symbols you use?)

 

[grecko94]

because
F = -K(Y-Yo) -C(y°-yo°)

Y is the displacement of the spring which in my program is x
Yo is the wavy surface profile
and since Yo = A sin(Ωt)
then the code

am*sin(wb*dt)

and yo° = ΩA cos(Ωt)
and the code

wb*am*cos(wb*dt))

Ω= wb
U= car velocity
dt = time

Why do you multiply wd dt with u ?

you are right, i don't need u

 

[BvU]

So you are happy now ? Any other problems ? Then it's time to make $c\ne 0$ and see if the damping works.

maybe i should make the initial condition velocity of spring to zero

Meaning you assume the road is flat before the start of the wobbling at $x=0$. Reasonable enough.

 

[grecko94]

The damping work fine by following damping ratio rule
ζ=C/Cc, Cc = 2√km

program dashpot

implicit none

integer,parameter :: N=10001,IN=300

integer :: i

real :: pi,dt,k,m,c,U,d,wb,am

real,dimension (N)::T,V,X,f

!initial condition= position,velocity,time, pi

pi = 3.14159265359

dt = 2*pi/FLOAT(N-1)

x(1) = 0

t(1) = 0

v(1) = 0

k = 625
m = 25

c = 0

U = 5

d = 2*pi

wb = 2*pi*U/D

am = 1
!looping

do I=1, N-1

T(i+1) = dt*i

X(i+1) = x(i) + v(i)*dt

V(i+1) = v(i) + (-k/m*(x(i)-(am*sin(wb*dt)))-c/m*(v(i)-(wb*am*cos(wb*dt))))*dt
f(i) = ABS(-k*(x(i)-am*sin(wb*dt))-c*(v(i)-wb*U*am*cos(wb*dt)))end do

print*, '      T(I)          X(I)            V(I)        F'

WRITE (*,"(4F16.8)") (T(I),X(I),V(I),f(I),I=1,N,IN)

END PROGRAM dashpot

however there is one thing bugging me
when ω=Ω isn't the value of spring displacement (X) supposed to be higher than ω<Ω and ω>Ω

ω=√(k/m) Ω= 2*velocity of the vehicle(u)/ Distance(d)

 

[BvU]

Ω= 2* π * velocity of the vehicle(u)/ Distance(d)

however there is one thing bugging me

And your program doesn't demonstrate that ? What does it show ?

 

[BvU]

Ah! found a big bug ...

Yo = A sin(Ωt)

And what is the argument of the sine in your program ? Oops !

 

[BvU]

There's a small bug and a remark as well:

f(i) = ABS(-k*(x(i)-am*sin(wb*dt))-c*(v(i)-wb* U *am*cos(wb*dt))

I don't think you want that U there when $c\ne0$
And re the ABS: I wonder how the ground can exert a negative force ?

 

[grecko94]

And what is the argument of the sine in your program ? Oops !

Sorry, i try to look at my program and still don't get what you said. can you explain a little bit about it
damn,i feel kinda stupid about it now

 

[BvU]

Post #11: Look at the expression you want to integrate too:

you want $\ \ \sin (\Omega t)\ \$ and not $\ \$ sin(wb*dt)​

Post #12: we'll get to that once you have some nice results using c = 0 as you have now.

 

[grecko94]

alright, I'm going to take a look after i finish my part time job.

but just one thing, in my program the expression for Ω is wb
so Ω = 2* π * velocity of the vehicle(u)/ Distance(d)
and i wrote wb = 2*π*U/d

isn't the same thing?

 

[BvU]

Let me give you a hint : The $\Omega$ part in $\ \ \sin (\Omega t)\ \$ is correct in the program.
So yes, it's the same thing in the formulas and in the Fortran.

 

[grecko94]

Sooo i assumed the time function is not correct ?
not dt bu T(i) ?

 

[BvU]

YESSSS ! dt is a constant ! T(i) increases from loop pass to loop pass

 

[grecko94]

hmm... it does seem getting better, with or without damp.
but hey, what do you think about the force? the ground exert a negative force at the hill and big positive force on the valley.
should i add another condition?, so when the F value goes to negative it equal to zero

 

[BvU]

Not so fast.

There is a force that has been neglected so far: gravity. A baseload that compresses the spring, so the actual F doesn't go below zero so quickly any more. But it still can (if you drive fast enough, one or more of your wheels can leave the ground).

And there also is a mass that has been neglected so far: the parts of your contraption other than M. Once the little wheel is off the ground, the spring and dashpot can't press it down infinitely fast.
Furthermore, the spring can't be compressed further than a certain length; at that point you no longer have F = -k(x-x₀ .

Of course, all these nice complications can be translated to something you can incorporate in your model. However:

I propose you first look into $c>0$ if you haven't done that already.

Then:
The problem you are working on is known as a damped, driven harmonic oscillator. Didactically a gold mine: it occurs all over physics and a good understandoing is worth a great deal. I don't know where you are in your curriculum, but the solution of your differential equation consists of a transient response (dependent on intital conditions) and a long-term response that is proportional to the driving force. See if you can make those two visible explicitly.

 

[grecko94]

YES SIR, i already did c>0 , and it does fulfill the condition for underdamped , overdamp and critical.

There is a force that has been neglected so far: gravity

And there also is a mass that has been neglected so far:

If that the case, then i need to add something to my equation of motion. I will try again on Monday at school.

By the way, I'm 4th year undergraduate student in japan and on my way to master. So basically, this semester its pretty much only self study for me. Every week, i have to report what i learned to my teacher. and since my research field going to be related to rail corrugation wear, fortran is compulsory for me.

and oh , if you would be so kind as to recommend me some fortran textbook i would be so grateful . cause, so far youtube is the only place i learned about fortran.

 

[BvU]

I greatly benefited from studying the Digital Equipment Corporation "Fortran User Guide" in the 1980's . Not much use to you, I'm afraid . So I wouldn't have a better advice than google "learning fortran" and find places like https://www.fortrantutorial.com/ or https://www.tutorialspoint.com/fortran/

(perhaps Monday I'll find a few more resources on my PC at work... )

 

[BvU]

Well, that went wrong altogether; I was thinking of a link in my browser favorites but it took me a long time to come up with this one:
http://www.fortran.com/the-fortran-company-homepage/fortran-tools-libraries-and-application-software/

Perhaps these might help you (but you can of course google yet another thousand of their kind ) :
https://www.nsc.liu.se/~boein/f77to90/
http://www.fortran.bcs.org/resources.php
http://www.fortranlib.com/freesoft.htm#Application
http://www.star.le.ac.uk/~cgp/fortran.html

 

[grecko94]

Thank you it going to help me a lot ! Btw just want to inform you, that i already fix the program and it really have beautiful output

 

[grecko94]

Hey there, sorry when i said i already finished my program. But my teacher told me to find the error in this data. i`m just wondering do you know any method to calculate error in my program. Thank you in advance

program dashpot
implicit none
integer,parameter :: N=30001,IN=300
integer :: i
real :: pi,dt,k,m,c,U,d,wb,am
real,dimension (N)::T,V,X,f
!initial condition= position,velocity,time, pi
pi = 3.14159265359
dt = 2*pi/FLOAT(N-1)
x(1) = 0
t(1) = 0
v(1) = 0
k = 543
m = 17
c = 0
U = 3
d = 10
wb = 2*pi*U/D
am =0.1
!looping
do I=1, N-1
T(i+1) = dt*i
X(i+1) = x(i) + v(i)*dt
V(i+1) = v(i) + (-k/m*(x(i)-am*sin(wb*u*t(i)))-c/m*(v(i)-wb*U*am*cos(wb*u*t(i)))-am*(wb**2)*(U**2)*sin(t(i)*U*wb))*dt
f(i) = -k*(x(i)-am*sin(wb*u*t(i)))-c*(v(i)-wb*U*am*cos(wb*u*t(i)))
end do
print*, '   T(I)             X(I)         V(I)          F'
WRITE (*,"(4F16.8)") (T(I),X(I),V(I),f(I),I=1,N,IN)
END PROGRAM dashpot

 

[BvU]

Sorry I completely missed this extra. Did you find a decent method ? With this kind of simulation, simple forward Euler is 'always' too late and thereby accumulates a considerable error. Reducing step size is a start to get an idea of the error. Improving on Euler is a second way. Recently some threads revealed a very simple improvement: euler-cromer for a damped pendulum - not unlike a dashpot ! ) Some googling shows other possibilities.

Another approach: simplify the problem until you can get an analytical answer and compare that with the numerical result. There's a lot to learn that can come in useful in later practice.

 

[grecko94]

its fine, i already finished this quite long time ago