C使い方 C34-41行目の装置の条件に値を入れる C46-47行目の境界条件を入れる Cコンパイル実行 Cacoustic_field.txt と work_flow.txtが生成 前者には位置，圧力振幅，流速振幅，圧力の位相，後者には位置，仕事流が書き込まれる C5 プログラムを書くことに関してはど素人のため，スマートな命令分ではない C===================================================================================================== program tube IMPLICIT NONE integer N,Nmax COMPLEX*16 P0,U0,P,U DOUBLE PRECISION Pamp,Uamp,Phi,W,Pphi,Uphi COMPLEX*16 I,Ya,Yu,Xa,Xu,k,Mtube(2,2),Mtube1(2,2) DOUBLE PRECISION Pi,r,L,A,TH,k0,omega,f,cp,rhom DOUBLE PRECISION Mu,Nyu,Kappa,Alpha,Pr,Ss,Tm,Gam,Pm DOUBLE PRECISION wta,wtu,X,dx COMPLEX*16 DJ1J0 C============================================================ OPEN(10,FILE='acoustic_field.TXT') OPEN(11,FILE='work_flow.TXT') C============================================================ Pi=4D0*DATAN(1D0) !円周率の定義 I =CMPLX(0D0,1D0) !複素数iの定義 Nmax=100 !管の分割数 C============================================================ C=========管の形状など======================================= C============================================================ Pm=101D3 !mean pressure r=10.5D-3 !radius of tube L=1.00D0 !length of tube A=r*r*Pi !cross-sectional area of tube TH=300D0 !mean temperature f=148.5D0 !frequency omega=2D0*Pi*f !angular frequency dx=-L/real(Nmax) !ここの-をなくせば，正の方向に進む．***重要事項*** C============================================================ C=========境界条件======================================= C============================================================ P0=CMPLX(0.78D3,0D0) !pressure at tube-end 閉端の場合は有限値，開端の場合はゼロ U0=CMPLX(0D0,0D0) !velocity at tube-end 開端の場合はゼロ , 開端の場合は有限値 Tm=TH CALL bussei(cp,rhom,Mu,Nyu,Kappa,Alpha,Pr,Ss,Tm,Gam,Pm) !物性値の読み込み，初期設定は空気 wta =(r*r)*omega/(2D0*Alpha ) !おめが　たう　アルファ wtu =(r*r)*omega/(2D0*Nyu ) !おめが　たう　ニュー Ya =CMPLX(1D0,1D0)*CMPLX(DSQRT(wta),0D0) Yu =CMPLX(1D0,1D0)*CMPLX(DSQRT(wtu),0D0) Xa =2D0/I/Ya*DJ1J0(I*Ya) !かい　アルファー　DJ1J0はベッセル関数の比 Xu =2D0/I/Yu*DJ1J0(I*Yu) !かい　ニュー　DJ1J0はベッセル関数の比 K0=omega/Ss !自由空間中の波数 K=k0*CSQRT( ( 1D0+(Gam-1D0)*Xa ) /(1D0-Xu) ) !半径rの管中の波数 Do 100 N=1,Nmax !長さLの管をNmax分割して音場を計算 X=(real(N)*dx) Mtube(1,1)=CDCOS(K*X) !閉端からXの位置までの管の伝達行列1 1成分 Mtube(1,2)=-(I*omega*rhom)/(K*(1-Xu))*CDSIN(K*X) !閉端からXの位置までの管の伝達行列1 2成分 Mtube(2,1)=(K*(1-Xu))/(i*omega*rhom)*CDSIN(K*X) !閉端からXの位置までの管の伝達行列2 1成分 Mtube(2,2)=CDCOS(K*X) !閉端からXの位置までの管の伝達行列2 2成分 P=P0*Mtube(1,1)+U0*Mtube(1,2) !位置xでの圧力(複素数） U=P0*Mtube(2,1)+U0*Mtube(2,2) !位置xでの流速(複素数) Pamp=CDABS(P) !位置xでの圧力振幅(実数) Uamp=CDABS(U) !位置xでの流速振幅(実数) Pphi=DATAN2(AIMAG(P),DBLE(P)) !位置xでの圧力の位相(実数，ラジアン)，基準はP0 Uphi=DATAN2(AIMAG(U),DBLE(U)) !位置xでの流速の位相(実数，ラジアン)，基準はP0 Phi=(Uphi-Pphi) !位置xでの圧力から見た流速の位相（実数，ラジアン） W=0.5*DBLE(P*CONJG(U)) !位置xでの仕事流，断面積Aを掛ければpower flowになる. write(10,*)x,Pamp,Uamp,Pphi !ファイル10（acoustic_field.txtに（位置，圧力振幅，流速振幅，圧力の位相）を書きこみ write(11,*)x,W !ファイル10（work_flow.txtに(位置，仕事流(W/m^2))を書きこみ 100 CONTINUE close(10) close(11) 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 C密度 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 C粘性係数 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 C熱伝導度 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 C比熱比 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定圧比熱 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動粘性係数 Nyu =Mu /rhom C熱拡散係数 Alpha = Kappa / rhom / Cp Cプラントル数 Pr =Mu / Kappa *Cp C音速 Ss =DSQRT(Gam * Pm /rhom ) RETURN END 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 このFunctionはベッセル関数の比を計算するのに使います．少々，tanhを使用する場合（平行平板）に比べて計算負荷が増える 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