求助,function POINTEHL
syms ENDA,A1,A2,A3,Z,HM0,DH
pi=3.14159265 Z=0.68
[N,PH,E1,ENDA0,RX,US,X0,XE]=deal(33,0.8E9,2.21E11,0.05,0.02,1.0,-2.5,1.5)
MM=N-1
A1=log(EDA0)+9.67
A2=5.1E-9*PH
A3=0.59/(PH*1.E-9)
U=EDA0*US/(2.*E1*RX)
B=pi*PH*RX/E1
W0=2.*PAI*PH/(3.*E1)*(B/RX)^2
ALFA=Z*5.1E-9*A1
G=ALFA*E1
HM0=3.63*(RX/B)^2*G^0.49*U^0.68*W0^(-0.073)
ENDA=12.*U*(E1/PH)*(RX/B)^3
fprintf('%d,%d,%d,%d,%d,%d,%d,%d,%d'\n,N,X0,XE,W0,PH,E1,EDA0,RX,US)
fprintf('wait please!')
function SUBAK(MM)
function EHL(N,X0,XE)
stop
end
function EHL(N,X0,XE)
X=cell(65),Y=cell(65),H=cell(4500),RO=cell(4500),EPS=cell(4500),EDA=cell(4500),P=cell(4500),POLD=cell(4500),V=cell(4500)
syms ENDA,A1,A2,A3,Z,HM0,DH
[ MK,G0]=deal(1,2.0943951)
function INITI(N,DX,X0,XE,X,Y,P,POLD)
KK=0
function HREE(N,DX,KK,H00,G0,X,Y,H,RO,EPS,EDA,P,V)
function ITER(N,KK,DX,H00,G0,X,Y,H,RO,EPS,EDA,P,V)
MK=MK+1
function ERP(N,ER,P,POLD)
fprintf('ER=%d',ER)
if(ER.>1.E-5&DH.>1.E-6)
elseif(MK>=20)
MK=1
DH=0.5*DH
end
KK=15
end
if (DH.L<=1.E-6)fprintf('Pressure are not convergent!!!')
function OUTPUT
syms N,DX,X,Y,H,P
return
end
function ERP
syms N,ER,P,POLD
P=cell(N,N),POLD=cell(N,N)
ER=0.0
SUM=0.0
for I=1:N; J=1:N
ER=ER+ABS(P(I,J)-POLD(I,J))
POLD(I,J)=P(I,J)
SUM=SUM+P(I,J)
continue
ER=ER/SUM
return
end
function INITI
syms N,DX,X0,XE,X,Y,P,POLD
X=cell(N),Y=cell(N),P=cell(N,N),POLD=cell(N,N)
NN=(N+1)/2
DX=(XE-X0)/(N-1.)
Y0=-0.5*(XE-X0)
for I=1:N
X(I)=X0+(I-1)*DX
Y(I)=Y0+(I-1)*DX
continue
for I=1:N
D=1.-X(I)*X(I)
for J=1:NN
C=D-Y(J)*Y(J)
if(C<=0.0) P(I,J)=0.0
elseif(C>=0.0) P(I,J)=sqrt(C)
POLD(I,J)=P(I,J)
end
end
return
end
function HREE
syms N,DX,KK,H00,G0,X,Y,H,R0,EPS,EDA,P,V
X=cell(N),Y=cell(N),P=cell(N,N),H=cell(N,N),R0=cell(N,N),EPS=cell(N,N),EDA=cell(N,N),V=cell(N,N)
syms ENDA,A1,A2,A3,Z,HM0,DH
AK=(0:65,0:65)
DATA PAI,PAI1/3.14159265,0.2026423/
NN=(N+1)/2
fonction VI(N,DX,P,V)
HMIN=1.E3
for I=1:N;J=1:NN
RAD=X(I)*X(I)+Y(J)*Y(J)
W1=0.5*RAD
H0=W1+V(I,J)
if (H0<HMIN) HMIN=H0
H(I,J)=H0
if(KK==0)
KK=1
DH=0.01*HM0
H00=-HMIN+HM0
end
W1=0.0
for I=1:N; J=1:N
W1=W1+P(I,J)
W1=DX*DX*W1/G0
DW=1.-W1
if(DW<0.0) H00=H00+DH
if(DW>0.0)H00=H00-DH
end
end
for I=1:N; J=1:NN
H(I,J)=H00+H(I,J)
ENDA1=EXP(A1*(-1.+(1.+A2*P(I,J))^Z))
EDA(I,J)=EDA1
R0(I,J)=(A3+1.34*P(I,J))/(A3+P(I,J))
EPS(I,J)=R0(I,J)*H(I,J)^3/(ENDA*EDA1)
for J=NN+1:N
JJ=N-J+1
for I=1:N
H(I,J)=H(I,JJ)
R0(I,J)=R0(I,JJ)
EDA(I,J)=EDA(I,JJ)
EPS(I,J)=EPS(I,JJ)
return
end
function ITER
syms N,KK,DX,H00,G0,X,Y,H,RO,EPS,EDA,P,V
X=cell(N),Y=cell(N),P=cell(N,N),H=cell(N,N),RO=cell(N,N),EPS=cell(N,N),EDA=cell(N,N),V=cell(N,N)
AK=cell(0:65,0:65)
pi=0,3.14159265
if KG1~=0 NN=(N+1)/2
KG1=1
AK00=AK(0,0)
AK10=AK(1,0)
end
for K=1:KK
for J=2:NN
J0=J-1
J1=J+1
D2=0.5*(EPS(1,J)+EPS(2,J))
for I=2:N-1
I0=I-1
I1=I+1
D1=D2
D2=0.5*(EPS(I1,J)+EPS(I,J))
D4=0.5*(EPS(I,J0)+EPS(I,J))
D5=0.5*(EPS(I,J1)+EPS(I,J))
D8=2.0*RO(I,J)*AK00/PAI^2
D9=2.0*RO(I0,J)*AK10/PAI^2
D10=D1+D2+D4+D5+D8*DX-D9*DX
D11=D1*P(I0,J)+D2*P(I1,J)+D4*P(I,J0)+D5*P(I,J1)
D12=(RO(I,J)*H(I,J)-D8*P(I,J)-RO(I0,J)*H(I0,J)+D9*P(I,J))*DX
P(I,J)=(D11-D12)/D10
if(P(I,J)<=0.0) P(I,J)=0.0
continue
for J=1:NN
JJ=N+1-J P(I,JJ)=P(I,J)
for I=1:N
P(I,JJ)=P(I,J)
function HREE
continue
return
end
end
end
function VI
syms N,DX,P,V
P=cell(N,N),V=cell(N,N)
AK=cell(0:65,0:65)
PAI1=0.2026423
for I=1:N J=1:N
H0=0.0 V(I,J)=H0*DX*PAI1
end
for K=1:N
IK=abs(I-K) H0=H0+AK(IK,JL)*P(K,L)
for L=1:N
JL=abs(J-L) H0=H0+AK(IK,JL)*P(K,L)
return
end
function SUBAK
syms MM
AK=cell(0:65,0:65)
S(X,Y)=X+sqrt(X^2+Y^2)
for I=0:MM
XP=I+0.5
XM=I-0.5
for J=0:I
YP=J+0.5
YM=J-0.5
A1=S(YP,XP)/S(YM,XP)
A2=S(XM,YM)/S(XP,YM)
A3=S(YM,XM)/S(YP,XM)
A4=S(XP,YP)/S(XM,YP)
AK(I,J)=XP*ALOG(A1)+YM*ALOG(A2)+XM*ALOG(A3)+YP*ALOG(A4)
AK(J,I)=AK(I,J)
return
end
function OUTPUT
syms N,DX,X,Y,H,P
X=cell(N),Y=cell(N),H=cell(N,N),P=cell(N,N)
A=0.0
fprintf('1X=%d' I,Y(I))
for I=1:N
fprintf ('X(I)',I,J,H(I,J))
end
fprintf('A=%d' I,Y(I))
for I=1:N
fprintf('X(I)=%d',I,J,(P(I,J) )
end
return
end |