Modélisation de pompe centrifuge.

[tchouchedi]

Monsieur,
Je travaille sur un projet de modélisation de la montée en vitesse d'une pompe centrifuge.
J'ai trouvé un algorithme sur fortran mais je n'arrive pas à le comprendre surtout les variables.
Si quelqu'un peut m'aider et c'est urgent.
merci

****************************** l'algorithme************
Methode des caractéristiques
! Exemple 
integer N1,N2,N3
real R,rp,kL,la,kv,L1,L3,kk,delta,L 2,kvv

dimension V1(2001),Vp1(2001),H1(2001),Hp 1(2001)
dimension V2(2001),Vp2(2001),H2(2001),Hp 2(2001),kL(2001)
dimension X(2001),Dh(2001)
dimension V3(2001),Vp3(2001),H3(2001),Hp 3(2001)

open(1,file='C:\pression.dat')
open(2,file='C:\vitesse.dat')


!	 *************** données de la pompe *****************
pi=3.14
! N2=5
R=0.064
!	R=0.084
rp=0.019
Nm=1500
omo=Nm*pi/30
b=0.0045
c=900

kvv=8

teta=2*pi/kvv
g=9.81
tp=0.15
la=0.02 !6

L2=R-rp



! ************* données des conduites amont et aval**********************
D1=0.04
L1=0.5 !0.5
D3=0.04
L3=1.5 !0.5
!	write(*,*)omo,la,teta,1/2
!	stop

! ***************	calcul du diametre hydraulique pompe ****************************

dt=0.000005
N2=L2/dt/c

dx2=(R-rp)/N2
!	write(*,*)dx2,N2
!	stop


do 10 i=1,N2+1

kL(i)=(rp+(i-1)*dx2)*teta

Dh(i)=4*b*kL(i)/(2*kL(i)+2*b)
X(i)=rp+(i-1)*dx2

!	 write(*,*) Dh(i),x(i),dh(i)/dh(1)

10 continue

! ****************	calcul régime permanant pompe ***************************

!	stop

do 20 i=1,N2+1

H2(i)=0

V2(i)=0

20 continue




!	dt=dx2/c
N1=L1/(c*dt)
!	L1=c*dt*N1

N3=L3/(c*dt)
!	L3=c*dt*N3

!	write(*,*) N1,N3,L1,L3,N2

!	stop

do 21 i=1,N1+1

H1(i)=0

V1(i)=0

21 continue

do 22 i=1,N3+1

H3(i)=0

V3(i)=0

22 continue



! ****************************** *Calcul Du Régime Transitoire

t=0

!	write(1,*) t,H3(N3+1)
!	 write(2,*) t,V2(N2+1)

do 31 j=1,150000
t=t+dt

if (t.ge.tp) then 
om=omo 
else 
om=(omo/tp)*t
end if

!	 ******* Calcul de la conduite amont 1 ******************************

do 25 i=2, N1

Vp1(i)=0.5*((g/c)*(H1(i-1)-H1(i+1))+V1(i+1)+V1(i-1)+(g/c)*dt*
1(V1(i+1)-V1(i-1))-(la*dt/(2*D1))*(V1(i-1)*abs(V1(i-1))+V1(i+1)*
1abs(V1(i+1))))


Hp1(i)=(c/(2*g))*(V1(i-1)-V1(i+1)-(g/c)*dt*(V1(i-1)+
1V1(i+1))+(la*dt/(2*D1))*(-V1(i-1)
1*abs(V1(i-1))+V1(i+1)*abs(V1(i+1))))+(H1 (i-1)+H1(i+1))/2

25	continue



!	 ******** calcul au niveau de la pompe 2	*****************************

do 26 i=2, N2


!	Vp2(i)=0.5*((g/c)*(H2(i-1)-H2(i+1))+V2(i+1)+V2(i-1)+
!	1(g/c)*dt*(V2(i+1)-V2(i-1))-(la*dt/2)*
!	1(V2(i-1)*abs(V2(i-1))/Dh(i-1)+V2(i+1)*abs(V2(i+1))/Dh(i+1))
! 1+(om**2)*dt*(X(i-1)+X(i+1)))

!	vp2(i)=vp2(i-1)*kl(i-1)/kl(i)

!	Hp2(i)=(c/(2*g))*(V2(i-1)-V2(i+1)-(g/c)*dt*
!	1(V2(i+1)+V2(i-1))+(la*dt/2)*(-V2(i-1)*abs(V2(i-1))
!	1/Dh(i-1)+V2(i+1)*abs(V2(i+1))/Dh(i+1))+(om**2)*dt*(X(i-1)-X(i+1)))
! 1+(H2(i-1)+H2(i+1))/2

C1=g/c

xl1=C1*H2(i-1)+V2(i-1)-(la*dt/2)*V2(i-1)*abs(V2(i-1))/Dh(i-1)+
1(om**2)*dt*X(i)-C1*dt*V2(i-1)

xl2=-C1*H2(i+1)+V2(i+1)-(la*dt/2)*V2(i+1)*abs(V2(i+1))/Dh(i+1)+
1(om**2)*dt*X(i)+C1*dt*V2(i+1)

Hp2(i)=(xl2-xl1)/(-2*C1)

vp2(i)=xl1-C1*Hp2(i)



!	write(*,*) t,xl1

!	stop

26	continue

! stop


!	 ************ calcul de la conduite aval	3 ****************************** *

do 27 i=2,N3

Vp3(i)=0.5*((g/c)*(H3(i-1)-H3(i+1))+V3(i+1)+V3(i-1)+(g/c)*dt*
1(V3(i+1)-V3(i-1))-(la*dt/(2*D3))
1*(V3(i-1)*abs(V3(i-1))+V3(i+1)*abs(V3(i+1))))


Hp3(i)=(c/(2*g))*(V3(i-1)-V3(i+1)-(g/c)*dt*(V3(i+1)+
1V3(i-1))+(la*dt/(2*D3))*(-V3(i-1)
1*abs(V3(i-1))+V3(i+1)*abs(V3(i+1))))+(H3 (i-1)+H3(i+1))/2

27	continue


!	 ********** condition au limite à l'extrémité amont de la conduite 1 *********

Hp1(1)=0
C1=(g/c)
C2=C1*H1(2)-V1(2)+(la*dt/(2*D1))*V1(2)*abs(V1(2))-C1*dt*V1(2)

!	A=C1*(2.14E-8)*D1**4
!	B=1-C1*(5.26E-6)*D1**2
!	C11=C2-C1*9.2384
!	delta1=B*B-4*A*C11


!	write(*,*) t,H3(N3+1)

!	stop	 -
!	Vp1(1)=(-B+sqrt(delta1))/2/A

Vp1(1)=C1*Hp1(1)-C2
!	Hp1(1)=(-Vp1(1)-C2)/C1
!	if (Hp1(1).lt.-10.33) goto 123

! goto 124
! 123 Hp1(1)=-10.33

!124	 continue 





!	*********** jonction entre conduite 1 aval et pompe extrémité amont ***************

s1=(pi*D1**2)/4
s3=(pi*D3**2)/4

!	sp=2*pi*rp*b ! surface latérale aspiration
!	spp=2*pi*R*b	 ! ,, ,, refoulement


Dx4=4*b*(kvv*kL(1))/(2*kvv*kL(1)+2*b)

Dx5=4*b*(kvv*kL(N2+1))/(2*kvv*kL(N2+1)+2*b)

sp=(pi*dx4**2)/4 
spp=(pi*dx5**2)/4

!	write(*,*) dx2,dx4,dx5
!	stop

kk=(sp+s1)/s1

Vp2(1)=(1/kk)*((g/c)*(H1(N1)-H2(2))+V2(2)+V1(N1)+(om**2)*dt *X(1)-
1(la*dt/2)*(V2(2)*abs(V2(2))/Dh(2)+V1(N1)*abs(V1(N1))/D1)+
1(g/c)*dt*(-V1(N1)+V2(2)))

Vp1(N1+1)=Vp2(1)*sp/s1

Hp1(N1+1)=(c/g)*(V1(N1)-Vp1(N1+1)-(la*dt/(2*D1))*
1V1(N1)*abs(V1(N1)))+H1(N1)-V1(N1)*dt

Hp2(1)=Hp1(N1+1)

!	write(*,*) 1/kk,Vp2(1),(spp+s3)/s3,(spp+s3)/spp

!	stop	

!	*********** jonction entre conduite 3 amont et pompe extrémité aval ***************

kk=(spp+s3)/s3

Vp2(N2+1)=(1/kk)*((g/c)*(H2(N2)-H3(2))+V2(N2)+V3(2)+(om**2)*dt *
1X(N2+1)+(g/c)*dt*(V3(2)-V2(N2))-(la*dt/2)*(V2(N2)*abs(V2(N2))/
1Dh(N2)+V3(2)*abs(V3(2))/D3))

Vp3(1)=Vp2(N2+1)*spp/s3

Hp3(1)=(c/g)*(Vp3(1)-V3(2)+(la*dt/(2*D3))*
1V3(2)*abs(V3(2)))+H3(2)+V3(2) *dt

Hp2(N2+1)=Hp3(1)


!	*********** jonction entre conduite 3 aval et vanne ***************

kv=0.15 !1000	 ! perte de charge singulière vanne
cv=1/sqrt(1+kv)
cc=0.8 !0.5	 ! coefficient de contaction
ag=1*s3	 ! ouverture de la vanne

bv=cv*cc*ag	 ! coefficient cd=cv*cc

ck=((s3/bv)**2)*(1./(2*g))	

bb=g/c

c3=(la*dt/(2*D3))*V3(N3)*abs(V3(N3))-bb*H3(N3)-V3(N3)+bb*V3(N3)*dt

delta=1-4*bb*ck*c3

Vp3(N3+1)=(-1+sqrt(delta))/(2*ck*bb)

!	Vp3(N3+1)=0

Hp3(N3+1)=(1/bb)*(-c3-Vp3(N3+1))


!	write(*,*) t,Hp3(N3+1),1/bb,delta


!	stop

do 28 i=1, N1+1

H1(i)=Hp1(i)
V1(i)=Vp1(i)

28	continue


do 29 i=1, N2+1

H2(i)=Hp2(i)
V2(i)=Vp2(i)



29	continue



do 30 i=1, N3+1

H3(i)=Hp3(i)
V3(i)=Vp3(i)



30	continue

!	if (IU/Kkk*Kkk-IU) 31,300,31

if (j/30*30-j) 39,300,39	

300	write(1,*) t,H3(1)-H1(N1+1)
write(2,*) t,V2(N2+1)*spp*1000

39 continue


31 continue
end

-----