Pompe centrifuge

*
**
***
****
*****
******
*******
********
*********
**********
***********
************
*************
**************
***************
****************
*****************
******************
*******************
********************
*********************
**********************
***********************
************************
*************************
**************************
***************************
****************************
*****************************
****************************** 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