請教一下 有大大會寫fortran程式嗎 小弟目前有一個內燃引擎otto cycle的程式
用fortran寫的 但是有一些數值跑不出來不合現象 請問有大大可以幫忙解決嗎
ps:因為此程式關係到畢業論文 希望有大大可以伸出援手
道時順利畢業再另外答謝 程式碼的部份 附在下面 希望可以幫忙解決
謝謝 感恩
問題:為何第一狀態跟第二狀態的溫度跟壓力已經變了 CP值卻還是固定呢
問題出在哪裡?
!***************************************
program test
real y1(6),y2(10)
data ao/47870./,fs/0.06548/
r= 10.
phi= 1.2
f= 0.1
p1= 1.0
t1= 350.
call farg(p1,t1,phi,f,h1,u1,v1,s1,y1,cp,dlvlt,dlvlp)
write(*,*) 'enthalpy 1',h1
write(*,*) 'internal energy 1',u1
write(*,*) 'specific volume 1',v1
write(*,*) 'entropy 1',s1
write(*,*) 'mole fraction',y1
write(*,*) 'constant pressure',cp
write(*,*) 'DLVLT',DLVLT
write(*,*) 'dlvlp',dlvlp
v2= v1/r
s2= s1
gamma= cp/(cp+p1*v1/t1/10.*dlvlt**2./dlvlp)
t2= t1*(v1/v2)**(gamma-1.0)
p2= p1*(v1/v2)**(gamma)
call cmprss(v2,s2,phi,r,f,t2,p2,u2)
write(*,*) '==================================================='
write(*,*) 'internal energy 2',u2
write(*,*) 'specific volume 2',v2
write(*,*) 'entropoy 2',s2
write(*,*) 'mole fraction',y2
write(*,*) 'constant pressure',cp
write(*,*) 'DLVLT',DLVLT
write(*,*) 'dlvlp',dlvlp
v3= v2
u3= u2
p3= 200.
t3= 4000.
call combst(v3,u3,phi,p3,t3,s3)
stop
end
!====================================================================
subroutine combst(v3,u3,phi,p3,t3,s3)
real y2(10)
data maxits/650/,tol/.0001/,dx3/1000./,dx4/100./
u2= u3
v2= v3
do i= 1,maxits
call ecp(p3,t3,phi,h3,u3,v3,s3,y2,cp,dlvlt,dlvlp,ier)
if(ier.ne.0)then
write(6,10)ier
endif
f1= u2-u3
f2= v2-v3
det= v3/p3*dlvlp*(cp-p3*v3/t3*dlvlt)+v3**2/t3*dlvlt*(dlvlt+dlvlp)
dx3= (f2*v3*(dlvlt+dlvlp)+f1*v3/p3*dlvlp)/det
dx4= (-v3/t3*dlvlt*f1+f2*(cp-p3*v3/t3*dlvlt))/det
t3= t3+dx3
p3= p3+dx4
if(abs(dx3)/t3 < tol.and.abs(dx4)/p3 < tol)then
return
endif
print*,t3,p3,s3
pause
enddo
write(6,20)
return
10 format('ECP RETURNED IER=',i3)
20 format('CONVERGENCE FAILURE')
end
!=======================================================
subroutine cmprss(v2,s2,phi,r,f,t2,p2,u2)
! IMPLICIT REAL (A-H,O-Z)
! IMPLICIT INTEGER (I-N)
! REAL IMEP
real y1(6)
data maxits/550/,tol/.0001/,dx1/2000./,dx2/100./
s1= s2
v1= v2
! DO ITERATION
do i= 1,maxits
call farg(p2,t2,phi,f,h2,u2,v2,s2,y1,cp,dlvlt,dlvlp)
if(abs(dx1)/t2 < tol.and.abs(dx2)/p2 < tol)then
return
endif
f1= s1-s2
f2= v1-v2
det= (cp/p2*dlvlp+v2/10./t2*dlvlt**2)
dx1= (t2/p2*dlvlp*f1+dlvlt*f2/10.)/det
dx2= (-dlvlt*f1+cp/v2*f2)/det
t2= t2+dx1
p2= p2+dx2
print*,t2,p2
pause
enddo
write(6,30)
return
30 format('CONVERGENCE FAILURE')
end
!=======================================================
subroutine farg(p,t,phi,f,h,u,v,s,y,cp,dlvlt,dlvlp)
!
! PURPOSE :
! COMPUTE THE PROPERTIES OF A FUEL, AIR, RESIDUAL
! GAS MIXTURE
!
! GIVEN :
! P - PRESSURE (BARS)
! T - TEMPERATURE (K)
! PHI - FUEL AIR EQUIVALENCE RATIO
! F - RESIDUAL MASS FRACTION
!
! RETURNS :
! H - ENTHALPY (J/G)
! U - INTERNAL ENERGY (J/G)
! V - SPECIFIC VOLUNE (CM**3/G)
! Y - A SIX DIMENSIONAL COMPOSITION VECTOR OF
! MOLE FRACTION
! 1=CO2, 2=H2O, 3=N2, 4=O2, 5=CO, 6=H2
! CP - SPECIFIC HEAT AT CONSTANT PRESSURE (J/G/K)
! DLVLT - DERIVATIVE OF LOG VOLUME WITH RESPECT TO LOG
! TEMPERATURE AT CONSTANT PRESSURE
! DLVLP - DERIVATIVE OF LOG VOLUME WITH RESPECT TO LOG
! PRESSUREE AT CONSTANT TEMPERATURE
!
! REMARKS :
! 1. VALID FOR 300 < T < 1000
! 2. ASSUMES THE FUEL IS GASOLINE(C7H17). TO CHANGE FUELS
! ONLY THE FUEL DATA STATEMENT NEEDS MODIFICATION
! 3. ENTHALIES OF O2, H2, N2 AND C(S) ARE SET TO ZERO
! AT 298K
! 4. THE FUEL MOLE FRACTION IS UNITY MINUS THE SUM OF Y(I)
!
!***************************************
logical rich,lean
real a(7,6),table(6),y(6),cp0(6),h0(6),s0(6),mres,n2,mw,mfa,mfuel,k
real m(6),nu(6)
!..... FUEL DATA
data
alpha/7.0/,beta/17.0/,gamma/0./,delta/0./,a0/4.0652/,b0/6.0977e-02/,c0/-1.8801e-5/,d0/-3.588e+04/,e0/15.45/
!..... TABLE 3.1, I=1,6
data
a/.24007797e+01,.87350957e-02,-.6607878e-05,.20021861e-08,.63274039e-15,-.48377527e+05,.96951457e+01,.40701275e+01,-.11084499e-02,.41521180e-05,-.29637404e-08,.80702103e-12,-.30279722e+05,-.32270046e+00,.36748261e+01,-.12081500e-02,.23240102e-05,-.63217559e-09,-.22577253e-12,-.10611588e+04,.23580424e+01,.36255985e+01,-.18782184e-02,.70554544e-05,-.67635137e-08,.21555993e-11,-.10475226e+04,.43052778e+01,.37100928e+01,-.16190964e-02,.36923594e-05,-.20319674e-08,.23953344e-12,-.14356310e+05,.2955535e+01,.305
74451e+01,.26765200e-02,-.58099162e-05,.55210391e-08,-.18122739e-11,-.98890474e+03,-.22997056e+01/
!..... OTHER DATA
data
ru/8.315/,table/-1.,1.,0.,0.,1.,-1/,m/44.01,18.02,28.008,32.000,28.01,2.018/
!..... COMPUTE RESIDUAL GAS COMPOSITION ACCORDING TO TABLE 3.3
rich= phi > 1.0
lean= .not.rich
dlvlt= 1.0
dlvlp= -1.0
eps= .21/(alpha+0.25*beta-0.5*gamma)
if(.not.rich)then
nu(1)= alpha*phi*eps
nu(2)= beta*phi*eps/2.
nu(3)= 0.79+delta*phi*eps/2.
nu(4)= 0.21*(1.0-phi)
nu(5)= 0.
nu(6)= 0.
dcdt= 0.
else
z= 1000./t
k= exp(2.743+z*(-1.761+z*(-1.611+.2803)))
dkdt= -k*(-1.761+z*(-3.222+z*.8409))/1000.
a1= 1.0-k
b= .42-phi*eps*(2.*alpha-gamma)+k*(.42*(phi-1.)+alpha*phi*eps)
c= -.42*alpha*phi*eps*(phi-1.0)*k
nu(5)= (-b+sqrt(b*b-4.*a1*c))/2./a1
!write(*,*) 'nu5',nu(5)
dcdt=
dkdt*(nu(5)**2-nu(5)*(.42*(phi-1.)+alpha*phi*eps)+.42*alpha*phi*eps*(phi-1.))/(2.*nu(5)*a1+b)
nu(1)= alpha*phi*eps-nu(5)
nu(2)= .42-phi*eps*(2.*alpha-gamma)+nu(5)
nu(3)= .79+delta*phi*eps/2.
nu(4)= 0.
nu(6)= .42*(phi-1.)-nu(5)
!..... COMPUTE MOLE FRACTIONS AND MOLECULAR WEIGHT OF RESIDUAL
endif
tmoles= 0.
do i= 1,6
tmoles= tmoles+nu(i)
enddo
mres= 0.
do i= 1,6
y(i)= nu(i)/tmoles
mres= mres+y(i)*m(i)
!..... COMPUTE MOLE FRACTION AND MOLECULAR WEIGHT OF FUEL-AIR
enddo
fuel= eps*phi/(1.+eps*phi)
o2= .21/(1.+eps*phi)
n2= .79/(1.+eps*phi)
mfa= fuel*(12.01*alpha+1.008*beta+16.*gamma+14.01*delta)+32.*o2+28.02*n2
!..... COMPUTE MOLE FRACTION OF FUEL-AIR-RESIDUAL GAS
yres= f/(f+mres/mfa*(1.-f))
do i= 1,6
y(i)= y(i)*yres
enddo
yfuel= fuel*(1.-yres)
y(3)= y(3)+n2*(1.-yres)
y(4)= y(4)+o2*(1.-yres)
!..... COMPUTE COMPONENT PROPERITIES
do i= 1,6
cp0(i)= a(1,i)+a(2,i)*t+a(3,i)*t**2+a(4,i)*t**3+a(5,i)*t**4
h0(i)=
a(1,i)+a(2,i)/2.*t+a(3,i)/3.*t**2+a(4,i)/4.*t**3+a(5,i)/5.*t**4+a(6,i)/t
s0(i)=
a(1,i)*log(t)+a(2,i)*t+a(3,i)/2.*t**2+a(4,i)/3.*t**3+a(5,i)/4.*t**4+a(7,i)
enddo
mfuel= 12.01*alpha+1.008*beta+16.000*gamma+14.01*delta
cpfuel= a0+b0*t+c0*t**2
hfuel= a0+b0/2.*t+c0/3.*t**2+d0/t
s0fuel= a0*log(t)+b0*t+c0/2.*t**2+e0
!..... COMPUTE PROPERTIES OF MIXTURE
h= hfuel*yfuel
if(yfuel == 0.)then
s= 0.
else
s= (s0fuel-log(yfuel))*yfuel
endif
cp= cpfuel*yfuel
mw= mfuel*yfuel
do i= 1,6
h= h+h0(i)*y(i)
if(y(i).ne.0.)then
s= s+y(i)*(s0(i)-log(y(i)))
endif
cp= cp+cp0(i)*y(i)+h0(i)*t*table(i)*dcdt*yres/tmoles
mw= mw+y(i)*m(i)
enddo
r= ru/mw
h= r*t*h
u= h-r*t
v= 10.*r*t/p
s= r*(-log(.9869*p)+s)
cp= r*cp
return
end
程式尚未結束 因為要先解決第一狀態和第二狀態cp為何相同的問題才能往下繼續做
麻煩各位大大了
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 140.112.14.201