program temp C==注意　この計算プログラムの結果は計算精度的に問題がある場合があります= C==論文等に利用する場合は適切に書き換えてください======================= C==プログラム製作者である上田はプログラミングに関しては全くの素人です=== C==プログラム製作者である上田はプログラミングに関しては全くの素人です=== C==プログラム製作者である上田はプログラミングに関しては全くの素人です=== C==必要があれば適切に書き換えてください================================= C==========使い方 C==========A．実験条件を入れる C==========B. 入口の圧力・流速を複素数で入れる c==========C. エンタルピーの値を入れる． C==========D. 出口温度と想定温度を比較 C==========ダメな場合はC. のエンタルピーの値を入れなおす． C======生成されるファイルは C Acoustic_filed.TXT 位置，圧力振幅，仕事流 C Temperature.TXT 位置，時間平均温度 C energy_flow.TXT 位置，エンタルピー流，仕事流，熱流 C======の3つ． IMPLICIT NONE integer b,bmax,m,mmax DOUBLE PRECISION cp,rhom,Mu,Nyu,Kappa,Alpha,Pr,Ss,Tm,Gam,Pm DOUBLE PRECISION R,PI,f DOUBLE PRECISION Ta,Tb,DelTm DOUBLE PRECISION Lstack,DelX,w DOUBLE PRECISION wtu,wta COMPLEX*16 Pin,Uin,Pout,Uout,P,U DOUBLE PRECISION Iin,Iout,In,Ha,H,C COMPLEX*16 Ya,Yu,Xa,Xu,I COMPLEX*16 ME(2,2),Ms(2,2), Mf(2,2),Mff(2,2) COMPLEX*16 DJ1J0 C================================================= mmax=1 !エンタルピーの値を探索する用 bmax=100 !スタック内を刻む用 C===========PiとIの定義===================================== PI=4D0*DATAN(1D0) !円周率の定義 I =CMPLX(0D0,1D0) !虚数項の定義 C===========単位行列の定義================================== ME(1,1)=CMPLX(1E0,0E0) ME(1,2)=CMPLX(0E0,0E0) ME(2,1)=CMPLX(0E0,0E0) ME(2,2)=CMPLX(1E0,0E0) C=====================A.実験条件==================== C=========================================================== R=0.1D-3 !蓄熱器の流路半径R粘性 Pm=1*101E3 !平均気圧 Ta=18+273 !蓄熱器一端の温度 Tb=278+273 !蓄熱器他端の温度 Lstack=35E-3 !蓄熱器の長さ f=41D0 !周波数 c====================入り口の計算====================================== w=2D0*f*Pi Tm=Ta CALL bussei(cp,rhom,Mu,Nyu,Kappa,Alpha,Pr,Ss,Tm,Gam,Pm) C============================================================== C==========B.入口の圧力と流速================================== C============================================================== Pin=cmplx(3.4D3,0D0) Uin=cmplx(3.4D3/(9D0*rhom*SS),0D0) Iin=0.5*DBLE(CONJG(Pin)*Uin) C============================================================== C==========C.エンタルピーの値================================== C============================================================== Ha=-0.21*Iin !エンタルピー流の値　右辺の最初の-0.21を変えると値が代わって，出口の温度が代わる．出口温度=想定温度となるような値を見つける． Mff(1,1)=ME(1,1) Mff(2,1)=ME(2,1) Mff(1,2)=ME(1,2) Mff(2,2)=ME(2,2) c========================================================== DelX=Lstack/bmax !管の刻み幅 Do 13 m=1,mmax !このdoを利用してニュートン法などを用いてエンタルピーの値を見つけ出すプルグラムを追加してください． OPEN(50,FILE='Acoustic_filed.TXT') OPEN(51,FILE='Temperature.TXT') OPEN(52,FILE='energy_flow.TXT') H=Ha Tm=Ta P=Pin U=Uin DO 3 b=1,bmax !蓄熱器内の計算スタート CALL bussei(cp,rhom,Mu,Nyu,Kappa,Alpha,Pr,Ss,Tm,Gam,Pm) !物性値の計算 wtu =R**2*w/(2D0*Nyu ) !このあたりの記号の使い方はは富永著の熱音響工学の基礎を参照 wta =R**2*w/(2D0*Alpha ) Ya =CMPLX(1D0,1D0)*CMPLX(DSQRT(wta),0D0) Yu =CMPLX(1D0,1D0)*CMPLX(DSQRT(wtu),0D0) Xa =2D0/I/Ya*DJ1J0(I*Ya) Xu =2D0/I/Yu*DJ1J0(I*Yu) C=rhom*cp*CDABS(U)*CDABS(U)/ &(2D0*w*(1D0-Pr*Pr)*CDABS(1D0-Xu)*CDABS(1D0-Xu) )* &AIMAG(Xa+Pr*CONJG(Xu)) !計算では熱伝導の項を無視した．実機の性能予測にはこの項は不可欠なので追加する必要あり DelTm=( &H-0.5D0*DBLE( &P*CONJG(U)*( 1D0-(Xa-CONJG(Xu) )/( (1+Pr)*(1-CONJG(Xu)) ) ) &) &)/C Ms(1,1)=0D0 Ms(1,2)=(-I*w*rhom)/(1D0-Xu) Ms(2,1)=-I*w*(1D0+(Gam-1D0)*Xa)/(Gam*Pm) Ms(2,2)=(Xa-Xu)*DelTm/(1D0-Xu)/(1D0-Pr)/Tm Mf(1,1)=ME(1,1)+delX*Ms(1,1) Mf(1,2)=ME(1,2)+delX*Ms(1,2) Mf(2,1)=ME(2,1)+delX*Ms(2,1) Mf(2,2)=ME(2,2)+delX*Ms(2,2) P=P*Mf(1,1)+U*Mf(1,2) U=P*Mf(2,1)+U*Mf(2,2) Tm=Tm+DelTm*DelX In=0.5*DBLE(P*CONJG(U)) write(50,*)DelX*b,CDABS(P),CDABS(U),In write(51,*)DelX*b,Tm write(52,*)DelX*b,H,In,H-In 3 continue close(50) close(51) close(52) 13 continue Pout=P Uout=U Iout=0.5*DBLE(CONJG(Pout)*Uout) write(*,*)'想定温度',Tb write(*,*)'計算の結果得られた出口温度',Tm write(*,*)'差',Tb-Tm write(*,*)'増幅量',In-Iin write(*,*)'出口における熱流',H-In write(*,*)'効率',-(In-Iin)/(H-In) stop End C=========物性値の計算====================================== SUBROUTINE bussei(cp,rhom,Mu,Nyu,Kappa,Alpha,Pr,Ss,Tm,Gam,Pm) DOUBLE PRECISION cp,rhom,Mu,Nyu,Kappa,Alpha,Pr,Ss,Tm,Gam,Pm rhom=(4.1461D0 & -0.019705*Tm & +4.7881e-05*Tm*Tm & -6.2842e-08*Tm*Tm*Tm & +4.2339e-11*Tm*Tm*Tm*Tm & -1.1472e-14*Tm*Tm*Tm*Tm*Tm)*Pm/101D3 Mu=-7.0031D-06 & +0.14018D-06*Tm & -0.00026899D-6*Tm*Tm & +3.5727D-13*Tm*Tm*Tm & -2.4676D-16*Tm*Tm*Tm*Tm & +6.8396D-20*Tm*Tm*Tm*Tm*Tm Kappa=4.9135D-3 & +0.053336D-3*Tm & +0.00013329D-3*Tm*Tm & -3.4467D-10*Tm*Tm*Tm & +3.4584D-13*Tm*Tm*Tm*Tm & -1.2551D-16*Tm*Tm*Tm*Tm*Tm Gam=1.3803D0 & +0.00020301*Tm & -5.1868D-7*Tm*Tm & +2.4927D-10*Tm*Tm*Tm & +1.0045D-13*Tm*Tm*Tm*Tm & -7.8378D-17*Tm*Tm*Tm*Tm*Tm C WRITE(*,*)Tm,Gam Cp=1057.9D0 & -0.3684*Tm & +0.00063128*Tm*Tm & +3.4653D-7*Tm*Tm*Tm & -8.7851D-10*Tm*Tm*Tm*Tm & +3.6441D-13*Tm*Tm*Tm*Tm*Tm C WRITE(*,*)Tm,rhom,Mu,Kappa Nyu =Mu /rhom Alpha = Kappa / rhom / Cp Pr =Mu / Kappa *Cp Ss =DSQRT(Gam * Pm /rhom ) c WRITE(*,*)Tm,Pr RETURN END C******************************************************************************* C=========行列の掛け算====================================== SUBROUTINE MATBY(Maa,Mbb,Mcc) COMPLEX*16 Maa(2,2),Mbb(2,2),Mcc(2,2) Mcc(1,1)=Maa(1,1)*Mbb(1,1)+Maa(1,2)*Mbb(2,1) Mcc(1,2)=Maa(1,1)*Mbb(1,2)+Maa(1,2)*Mbb(2,2) Mcc(2,1)=Maa(2,1)*Mbb(1,1)+Maa(2,2)*Mbb(2,1) Mcc(2,2)=Maa(2,1)*Mbb(1,2)+Maa(2,2)*Mbb(2,2) RETURN END C============================================================== C******************************************************************************* C FUNCTION OF J1(Z)/J0(Z) DEVELOPED BY SEIICHI WASHIO C J1(Z), J0(Z): 1ST AND 0-TH ORDER BESSEL FUNCTIONS OF THE FIRST KIND C DEVELOPED BY SEIICHI WASHIO C******************************************************************************* C COMPLEX*16 FUNCTION DJ1J0(Z) C IMPLICIT COMPLEX*16 (E,F,Z) IMPLICIT DOUBLE PRECISION (A,D,X,Y) COMPLEX*16 Z IF(CDABS(Z).GT.26.0) GO TO 9 M=1.581*CDABS(Z)+1.0 Z1=Z/2.0D+00 Z2=Z1*Z1 J1=0 J0=0 K=0 E1=1.0D+00 E0=1.0D+00 F1=1.0D+00 F0=1.0D+00 AK1=1.0D+00 1 K=K+1 AK=AK1 AK1=AK1+1.0D+00 F0=-F1*Z2/AK F1=F0/AK1 IF(K.GT.M)GO TO 2 E1=E1+F1 E0=E0+F0 GO TO 1 2 IF(J1.EQ.1)GO TO 5 S1=CDABS(F1) IF(ABS(S1/SNGL(DREAL(E1))).GE.1.0D-16)GO TO 3 IF(ABS(S1/SNGL(DIMAG(E1))).LT.1.0D-16)GO TO 4 3 E1=E1+F1 IF(J0.EQ.1)GO TO 1 GO TO 5 4 J1=1 IF(J0.EQ.1)GO TO 8 5 S0=CDABS(F0) IF(ABS(S0/SNGL(DREAL(E0))).GE.1.0D-16)GO TO 6 IF(ABS(S0/SNGL(DIMAG(E0))).LT.1.0D-16)GO TO 7 6 E0=E0+F0 GO TO 1 7 J0=1 IF(J1.EQ.1)GO TO 8 GO TO 1 8 DJ1J0=Z1*E1/E0 RETURN C APPROXIMATION BY HANKEL'S ASYMPTOTIC EXPANSION 9 X=DREAL(Z) Y=DIMAG(Z) I=1 IF(X.GE.0.0)GO TO 10 X=-X Y=-Y I=2 10 Z1=8.0D+00*DCMPLX(X,Y) ZE0=1.0D+00 E0=1.0D+00 ZE1=1.0D+00 E1=1.0D+00 ZF0=-1.0D+00/Z1 F0=-1.0D+00/Z1 ZF1=3.0D+00/Z1 F1=3.0D+00/Z1 Z2=Z1**2 IF(Y.LT.0.0)GO TO 11 D1=DEXP(-2.0D+00*Y) D0=1.0D+0-D1 D1=1.0D+0+D1 GO TO 12 11 D1=DEXP(2.0D+00*Y) D0=D1-1.0D+00 D1=1.0D+00+D1 12 AC=DCOS(X) AS=DSIN(X) DO 13 L=1,8 A1=DFLOAT((4*L-3)**2) A2=DFLOAT((4*L-1)**2) A3=DFLOAT(2*L*(2*L-1)) ZE0=-ZE0*A1*A2/Z2/A3 ZE1=-ZE1*(4.0D+00-A1)*(4.0D+00-A2)/Z2/A3 E0=E0+ZE0 13 E1=E1+ZE1 DO 14 L=1,7 A1=DFLOAT((4*L-1)**2) A2=DFLOAT((4*L+1)**2) A3=DFLOAT(2*L*(2*L+1)) ZF0=-ZF0*A1*A2/Z2/A3 ZF1=-ZF1*(4.0D+00-A1)*(4.0D+00-A2)/Z2/A3 F0=F0+ZF0 14 F1=F1+ZF1 ZA=DCMPLX(D1*AC,-D0*AS) ZB=DCMPLX(D1*AS,D0*AC) DJ1J0=((F1-E1)*ZA+(E1+F1)*ZB)/((E0+F0)*ZA+(E0-F0)*ZB) IF(I.EQ.1)RETURN DJ1J0=-DJ1J0 RETURN END