program temp C==注意　この計算プログラムの結果は計算精度的に問題がある場合があります= C==論文等に利用する場合は適切に書き換えてください======================= C==プログラム製作者である上田はプログラミングに関しては全くの素人です=== C==プログラム製作者である上田はプログラミングに関しては全くの素人です=== C==プログラム製作者である上田はプログラミングに関しては全くの素人です=== C==必要があれば適切に書き換えてください================================= C======================================================================= C======================================================================= C=========使い方======================================================== C42-47行目の条件を入力する C55,56行目の圧力と流速を複素数で入力する． IMPLICIT NONE integer b,bmax DOUBLE PRECISION cp,rhom,Mu,Nyu,Kappa,Alpha,Pr,Ss,Tm,Gam,Pm DOUBLE PRECISION R,PI,det,f DOUBLE PRECISION Tc,TH,dTm DOUBLE PRECISION Lstack,DelX,w DOUBLE PRECISION wtu,wta COMPLEX*16 Pin,Uin,Pout,Uout,P,U DOUBLE PRECISION Iin,Iout,In COMPLEX*16 Ya,Yu,Xa,Xu,I COMPLEX*16 ME(2,2),Ms(2,2), Mf(2,2),Mff(2,2),Msta(2,2) COMPLEX*16 DJ1J0 C================================================= OPEN(54,FILE='G.TXT') 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=====================実験条件==================== C=========================================================== R=0.34D-3 !蓄熱器の流路半径R粘性 Pm=1*101E3 !平均気圧 Tc=293 !蓄熱器低温端の温度 TH=293*2.3 !蓄熱器高温端の温度 Lstack=20E-3 !蓄熱器の長さ f=100D0 !周波数 c====================入り口の計算====================================== w=2D0*f*Pi Tm=TC CALL bussei(cp,rhom,Mu,Nyu,Kappa,Alpha,Pr,Ss,Tm,Gam,Pm) Pin=cmplx(1.08371789371,-0.19130178936) Uin=cmplx(1.777853001082D-03,8.733481219974D-05) Iin=0.5*DBLE(CONJG(Pin)*Uin) 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 !管の刻み幅 dTm=(TH-TC)/Lstack !温度勾配 DO 3 b=1,bmax !蓄熱器内の計算スタート Tm=TC+dTm*DelX*b 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) 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)*dTm/(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) CALL MATBY(Mf,Mff,Msta) Mff(1,1)=Msta(1,1) Mff(2,1)=Msta(2,1) Mff(1,2)=Msta(1,2) Mff(2,2)=Msta(2,2) P=Pin*Mff(1,1)+Uin*Mff(1,2) U=Pin*Mff(2,1)+Uin*Mff(2,2) In=0.5*DBLE(P*CONJG(U)) write(54,*)CDABS(P),CDABS(U),In/Iin 3 continue Pout=Mff(1,1)*Pin+Mff(1,2)*Uin Uout=Mff(2,1)*Pin+Mff(2,2)*Uin Iout=0.5*DBLE(CONJG(Pout)*Uout) write(*,*)'増幅率',Iout/Iin,'周波数',f write(*,*)'出口の音響インピーダンス',Pout/Uout write(*,*)'出口での特性音響インピーダンス',rhom*Ss 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