请教一下 有大大会写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