!= Module ChemCalc ! ! Authors:: SUGIYAMA Ko-ichiro ! Version:: $Id: chemcalc.f90,v 1.12 2006/09/30 05:53:42 odakker Exp $ ! Tag Name:: $Name: arare4-20061224 $ ! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved. ! License:: See COPYRIGHT[link:../../COPYRIGHT] ! !== Overview ! !化学関連の諸量を計算するためのモジュール. !AMP と Antoine の飽和蒸気圧式を用いて以下を求める. !デフォルトでは AMP 式を使うようにしてある. ! * 飽和蒸気圧 ! * 飽和蒸気圧の温度微分 ! * 潜熱 ! !== Error Handling ! !== Bugs ! !== Note ! !== Future Plans ! ! module ChemCalc ! !化学関連の諸量を計算するためのモジュール. !AMP と Antoine の飽和蒸気圧式を用いて以下を求める. !デフォルトでは AMP 式を使うようにしてある. ! * 飽和蒸気圧 ! * 飽和蒸気圧の温度微分 ! * 潜熱 ! !モジュール呼び出し use dc_message, only: MessageNotify use ChemData, only: ChemData_SpcNum, & !データベース上の化学種数 & ChemData_MolWt, & !分子量 & ChemData_SvapPress_AntoineA, & !Antoine 式の A 係数 & ChemData_SvapPress_AntoineB, & !Antoine 式の B 係数 & ChemData_SvapPress_AntoineC, & !Antoine 式の C 係数 & ChemData_SvapPress_AntoineUnit,& !単位変換用係数 & ChemData_SvapPress_AMPA, & !AMP 式の A 係数 & ChemData_SvapPress_AMPB, & !AMP 式の B 係数 & ChemData_SvapPress_AMPC, & !AMP 式の C 係数 & ChemData_SvapPress_AMPD, & !AMP 式の D 係数 & ChemData_SvapPress_AMPE !AMP 式の E 係数 !暗黙の型宣言禁止 implicit none !全体に対して private 属性を付ける. private !公開するものには public 属性を付ける public ChemCalc_Init !初期化ルーチン public ChemCalcSpc_Init, ChemCalcDim_Init !初期化ルーチン public MolWt !分子量 public GasR !気体定数 public Phase !相 public CpRef, CpPerMolRef, CvRef !定圧比熱, 定積比熱 public SvapPress, xz_SvapPress !飽和蒸気圧 [Pa] public xz_DSvapPressDTemp !飽和蒸気圧の温度微分 [Pa/K] public xz_LatentHeat !潜熱 [J/K kg] public LatentHeatPerMol !潜熱 [J/K mol] public ReactHeatNH4SH, ReactHeatNH4SHPerMol !NH4SH 反応熱 real(8) :: ReactHeatNH4SH !NH4SH 生成反応熱 [J/K kg] real(8) :: ReactHeatNH4SHPerMol !NH4SH 生成反応熱 [J/K mol] integer :: XMin = 10 !X 座標範囲(下限). 値はダミー integer :: XMax = 10 !X 座標範囲(上限). 値はダミー integer :: ZMin = 10 !Z 座標範囲(下限). 値はダミー integer :: ZMax = 10 !Z 座標範囲(上限). 値はダミー integer :: SpcNum = 1 !化学種の数. 値はダミー real(8) :: a_antA(ChemData_SpcNum) !Antoine の蒸気圧式の A 係数 real(8) :: a_antB(ChemData_SpcNum) !Antoine の蒸気圧式の B 係数 real(8) :: a_antC(ChemData_SpcNum) !Antoine の蒸気圧式の C 係数 real(8) :: a_antU(ChemData_SpcNum) !Antoine の蒸気圧式の単位換算のための係数 real(8) :: a_ampA(ChemData_SpcNum) !AMP 式の蒸気圧式の A 係数 real(8) :: a_ampB(ChemData_SpcNum) !AMP 式の蒸気圧式の B 係数 real(8) :: a_ampC(ChemData_SpcNum) !AMP 式の蒸気圧式の C 係数 real(8) :: a_ampD(ChemData_SpcNum) !AMP 式の蒸気圧式の D 係数 real(8) :: a_ampE(ChemData_SpcNum) !AMP 式の蒸気圧式の E 係数 real(8) :: a_MolWt(ChemData_SpcNum) !分子量 integer, allocatable :: SpcID(:) !系に含まれる化学種. 数字で判別 character(3), allocatable :: Phase(:) !相 save ReactHeatNH4SH, ReactHeatNH4SHPerMol save XMin, XMax, ZMin, ZMax, SpcNum save SpcID, Phase save a_antA, a_antB, a_antC, a_antU, a_MolWt save a_ampA, a_ampB, a_ampC, a_ampD, a_ampE contains !!! !!! 初期化ルーチン !!! !!!========================================================================== subroutine ChemCalc_Init(IMin, IMax, KMin, KMax, SNum, SID) ! !初期化ルーチンのラッパー. ! ChemCalcSpc_init と ChemCalcDim_init を実行 ! !モジュール呼び出し use ChemData, only: GasRUniv, & !気体定数 & ChemData_OneSpcID !化学種の ID 検索 !暗黙の型宣言禁止 implicit none !入出力変数 integer,intent(in) :: IMin !X 方向の配列の下限 integer,intent(in) :: IMax !X 方向の配列の上限 integer,intent(in) :: KMin !Z 方向の配列の下限 integer,intent(in) :: KMax !Z 方向の配列の上限 integer,intent(in) :: SNum !化学種の個数 integer,intent(in) :: SID(SNum) !系に存在する化学種の数 character(20) :: Name integer :: id !初期化ルーチンを 2 つ呼び出すだけ call ChemCalcSpc_Init(SNum, SID) call ChemCalcDim_Init(IMin, IMax, KMin, KMax) !NH4SH の反応熱の初期化 ! NH4SH 1kg に対する反応熱にする. Name = 'NH4SH-s' id = ChemData_OneSpcID( Name ) ReactHeatNH4SHPerMol = GasRUniv * 10834.0d0 ReactHeatNH4SH = GasRUniv * 10834.0d0 / MolWt( id ) ! write(*,*) ReactHeatNH4SH, MolWt( id ) call MessageNotify( "M", & & "ChemCalc_Init", "ReactHeatNH4SH = %f", d=(/ReactHeatNH4SH/) ) call MessageNotify( "M", & & "ChemCalc_Init", "NH4SH MolWt = %f", d=(/MolWt(id)/) ) end subroutine ChemCalc_Init !!!========================================================================== subroutine ChemCalcSpc_Init(SNum, SID) ! !初期化ルーチン(その 1). ! 利用する化学種の確定と, それぞれの物性値を決定する. ! 配列計算をしない場合には, この初期化をすれば十分である. ! !モジュール呼び出し use ChemData, only: ChemData_Phase !相 !暗黙の型宣言禁止 implicit none !入出力変数 integer,intent(in) :: SNum !化学種の個数 integer,intent(in) :: SID(SNum) !系に存在する化学種の数 !内部変数 integer :: s !----------------------------------------------------------- ! 初期化 !----------------------------------------------------------- SpcNum = SNum call MessageNotify( "M", & & "ChemCalcSpc_Init", "SpcNum = %d", i=(/SpcNum/) ) !配列割り当て if (allocated(SpcID)) deallocate( SpcID, Phase) allocate( SpcID(SpcNum), Phase(SpcNum) ) !----------------------------------------------------------- ! データベースから必要なデータだけ取り出す !----------------------------------------------------------- !系に存在する化学種. 入力された値を取り込む SpcID = SID !相 do s = 1, SpcNum Phase(s) = ChemData_Phase(SpcID(s)) end do !Antoine の飽和蒸気圧式の係数 a_antA = ChemData_SvapPress_AntoineA a_antB = ChemData_SvapPress_AntoineB a_antC = ChemData_SvapPress_AntoineC a_antU = ChemData_SvapPress_AntoineUnit !AMP 式の飽和蒸気圧式の係数 a_ampA = ChemData_SvapPress_AMPA a_ampB = ChemData_SvapPress_AMPB a_ampC = ChemData_SvapPress_AMPC a_ampD = ChemData_SvapPress_AMPD a_ampE = ChemData_SvapPress_AMPE !分子量 a_MolWt = ChemData_MolWt end subroutine ChemCalcSpc_Init !!!========================================================================== subroutine ChemCalcDim_Init(IMin, IMax, KMin, KMax) ! !初期化ルーチン(その 2) ! 配列サイズの確定を行う. 配列関数を利用する場合には, ! このサブルーチンで初期化する必要がある. ! 先に ChemCalcSpc_init を実行しておく必要がある. ! !暗黙の型宣言禁止 implicit none !入出力変数 integer,intent(in) :: IMin !X 方向の配列の下限 integer,intent(in) :: IMax !X 方向の配列の上限 integer,intent(in) :: KMin !Z 方向の配列の下限 integer,intent(in) :: KMax !Z 方向の配列の上限 !----------------------------------------------------------- ! 初期化 !----------------------------------------------------------- XMin = IMin; XMax = IMax ZMin = KMin; ZMax = KMax end subroutine ChemCalcDim_Init !!! !!! 飽和蒸気圧, 潜熱, etc. の基本関数. !!! 化学種の ID と温度に対して値を返す !!! !!!========================================================================== function CpRef(ID) ! !引数で与えられた化学種に対して, 標準状態での単位質量当たりの定圧比熱を計算 ! !モジュール呼び出し use ChemData, only: ChemData_CpRef !標準状態での単位質量当たりの比熱 !暗黙の型宣言禁止 implicit none !入出力変数 real(8) :: CpRef !標準状態での単位質量当たりの比熱 integer, intent(in) :: ID !化学種の ID !データベースから情報取得 CpRef = ChemData_CpRef(ID) end function CpRef !!!========================================================================== function CpPerMolRef(ID) ! !引数で与えられた化学種に対して, 標準状態での単位モル当たりの定圧比熱を計算 ! !モジュール呼び出し use ChemData, only: ChemData_CpPerMolRef !標準状態での単位モル当たりの比熱 !暗黙の型宣言禁止 implicit none !入出力変数 real(8) :: CpPerMolRef !標準状態での単位モル当たりの比熱 integer, intent(in) :: ID !化学種の ID !データベースから情報取得 CpPerMolRef = ChemData_CpPerMolRef(ID) end function CpPerMolRef !!!========================================================================== function CvRef(ID) ! !引数で与えられた化学種に対して, 標準状態での単位質量当たりの定圧比熱を計算 ! !モジュール呼び出し use ChemData, only: ChemData_CvRef !標準状態での単位質量当たりの比熱 !暗黙の型宣言禁止 implicit none !入出力変数 real(8) :: CvRef !標準状態での単位質量当たりの比熱 integer, intent(in) :: ID !化学種の ID !データベースから情報取得 CvRef = ChemData_CvRef(ID) end function CvRef !!!========================================================================== function MolWt(ID) ! !引数で与えられた化学種に対して, 分子量を計算 ! !モジュール呼び出し use ChemData, only: ChemData_MolWt !分子量 !暗黙の型宣言禁止 implicit none !入出力変数 real(8) :: MolWt !分子量 integer, intent(in) :: ID !化学種の ID !データベースから情報取得 MolWt = ChemData_MolWt(ID) end function MolWt !!!========================================================================== function GasR(ID) ! !引数で与えられた化学種に対して, 気体定数を計算 ! !モジュール呼び出し use ChemData, only: ChemData_GasR !気体定数 [J/K kg] !暗黙の型宣言禁止 implicit none !入出力変数 real(8) :: GasR !分子量 integer, intent(in) :: ID !化学種の ID !データベースから情報取得 GasR = ChemData_GasR(ID) end function GasR !!! !!! 空間 2 次元の関数群 !!! !!!========================================================================== function xz_SvapPressAnt( ID, xz_Temp ) ! !飽和蒸気圧の計算. Antoine の式を利用 ! !暗黙の型宣言禁止 implicit none ! 入出力変数 real(8) :: xz_SvapPressAnt(XMin:XMax, ZMin:ZMax) !飽和蒸気圧 real(8), intent(in) :: xz_Temp(XMin:XMax, ZMin:ZMax) !温度 integer, intent(in) :: ID !化学種の ID !内部変数 real(8) :: xz_LogSvapPress(XMin:XMax, ZMin:ZMax) real(8), parameter :: Temp0C = 273.15d0 !飽和蒸気圧の log を計算 !対数が大きくなりすぎないようにする. ! Fujitsu Fortran Compiler では 700 より大きい数の exp を取ると警告が出る. xz_LogSvapPress = & & min( & & ( & & a_antA(ID) & & - a_antB(ID) / (a_antC(ID) + xz_Temp - Temp0C) & & ) * dlog(10.0d0) & & + a_antU(ID), & & 7.0d2 & & ) !飽和蒸気圧を計算 xz_SvapPressAnt = dexp( xz_LogSvapPress ) end function xz_SvapPressAnt !!!========================================================================== function xz_DSvapPressDTempAnt(ID, xz_Temp) ! !飽和蒸気圧の温度微分の計算. Antoine の式を利用 ! !暗黙の型宣言禁止 implicit none !入出力変数 real(8) :: xz_DSvapPressDTempAnt(XMin:XMax, ZMin:ZMax) !飽和蒸気圧の温度微分 [Pa/K] real(8), intent(in) :: xz_Temp(XMin:XMax, ZMin:ZMax) !温度 [K] integer, intent(in) :: ID !化学種の ID !内部変数 real(8) :: xz_LogSvapPress(XMin:XMax, ZMin:ZMax) real(8) :: xz_DLogSvapPressDTemp(XMin:XMax, ZMin:ZMax) real(8), parameter :: Temp0C = 273.15d0 !飽和蒸気圧の log を計算 !対数が大きくなりすぎないようにする. ! Fujitsu Fortran Compiler では 700 より大きい数の exp を取ると警告が出る. xz_LogSvapPress = & & min( & & ( & & a_antA(ID) & & - a_antB(ID) / (a_antC(ID) + xz_Temp - Temp0C) & & ) * dlog(10.0d0) & & + a_antU(ID), & & 7.0d2 & & ) xz_DLogSvapPressDTemp = & & a_antB(ID) * dlog(10.0d0) & & / ( (a_antC(ID) + xz_Temp - Temp0C) ** 2.0d0 ) !飽和蒸気圧の温度微分を計算 xz_DSvapPressDTempAnt = xz_DLogSvapPressDTemp * dexp( xz_LogSvapPress ) end function xz_DSvapPressDTempAnt !!!========================================================================== function xz_LatentHeatAnt(ID, xz_Temp) ! !飽和蒸気圧より潜熱を計算する. Antoine の式を利用 ! !暗黙の型宣言禁止 implicit none !入出力変数 real(8) :: xz_LatentHeatAnt(XMin:XMax, ZMin:ZMax) !潜熱[J/K kg] real(8), intent(in) :: xz_Temp(XMin:XMax, ZMin:ZMax) !温度[K] integer, intent(in) :: ID !化学種の ID !内部関数 real(8) :: xz_DLogSvapPressDTemp(XMin:XMax, ZMin:ZMax) real(8), parameter :: GasRUniv = 8.314d0 real(8), parameter :: Temp0C = 273.15d0 !飽和蒸気圧の log を計算した後, 潜熱を計算する. xz_DLogSvapPressDTemp = & & a_antB(ID) * dlog(10.0d0) & & / ( (a_antC(ID) + xz_Temp - Temp0C) ** 2.0d0 ) xz_LatentHeatAnt = & & xz_DLogSvapPressDTemp * GasRUniv * (xz_Temp ** 2.0d0) & & / a_MolWt(ID) end function xz_LatentHeatAnt !!!========================================================================== function xz_SvapPress(ID, xz_Temp) ! !引数で与えられた化学種と温度に対して, 飽和蒸気圧を計算. AMP 式を利用 ! !暗黙の型宣言禁止 implicit none !入出力変数 real(8) :: xz_SvapPress(XMin:XMax, ZMin:ZMax) !飽和蒸気圧 real(8), intent(in) :: xz_Temp(XMin:XMax, ZMin:ZMax) !温度 [K] integer, intent(in) :: ID !化学種の ID !内部変数 real(8) :: xz_LogSvapPress(XMin:XMax, ZMin:ZMax) !飽和蒸気圧の対数を計算 !対数が大きくなりすぎないようにする. ! Fujitsu Fortran Compiler では 700 より大きい数の exp を取ると警告が出る. xz_LogSvapPress = & & min( & & ( & & a_ampA(ID) / xz_Temp & & + a_ampB(ID) & & + a_ampC(ID) * dlog( xz_Temp ) & & + a_ampD(ID) * xz_Temp & & + a_ampE(ID) * ( xz_temp ** 2 ) & & + dlog(1.0d-1) & & ), & & 7.0d2 & & ) !飽和蒸気圧を計算 xz_SvapPress = dexp( xz_LogSvapPress ) end function xz_SvapPress !!!========================================================================== function xz_DSvapPressDTemp(ID, xz_Temp) ! !引数で与えられた化学種と温度に対して, 飽和蒸気圧の温度微分を計算 ! !暗黙の型宣言禁止 implicit none !入出力変数 real(8) :: xz_DSvapPressDTemp(XMin:XMax, ZMin:ZMax) !飽和蒸気圧の温度微分 [Pa/K] real(8), intent(in) :: xz_Temp(XMin:XMax, ZMin:ZMax) !温度 [K] integer, intent(in) :: ID !化学種名 !内部変数 real(8) :: xz_LogSvapPress(XMin:XMax, ZMin:ZMax) real(8) :: xz_DLogSvapPressDTemp(XMin:XMax, ZMin:ZMax) !飽和蒸気圧の対数を計算 !対数が大きくなりすぎないようにする. ! Fujitsu Fortran Compiler では 700 より大きい数の exp を取ると警告が出る. xz_LogSvapPress = & & min( & & ( & & a_ampA(ID) / xz_Temp & & + a_ampB(ID) & & + a_ampC(ID) * dlog( xz_Temp ) & & + a_ampD(ID) * xz_Temp & & + a_ampE(ID) * ( xz_temp ** 2 ) & & + dlog(1.0d-1) & & ), & & 7.0d2 & & ) !飽和蒸気圧の温度微分 xz_DLogSvapPressDTemp = & & - a_ampA(ID) / (xz_Temp ** 2.0d0) & & + a_ampC(ID) / xz_Temp & & + a_ampD(ID) & & + a_ampE(ID) * 2.0d0 * xz_Temp !飽和蒸気圧の温度微分 xz_DSvapPressDTemp = & & xz_DLogSvapPressDTemp * dexp( xz_LogSvapPress ) end function xz_DSvapPressDTemp !!!========================================================================== function xz_LatentHeat(ID, xz_Temp) ! !引数で与えられた化学種と温度に対して, 潜熱を計算 ! !暗黙の型宣言禁止 implicit none !入出力変数 real(8) :: xz_LatentHeat(XMin:XMax, ZMin:ZMax) !潜熱 real(8), intent(in) :: xz_Temp(XMin:XMax, ZMin:ZMax) !温度 integer, intent(in) :: ID !化学種名 !内部変数 real(8) :: xz_DLogSvapPressDTemp(XMin:XMax, ZMin:ZMax) real(8), parameter :: GasRUniv = 8.314d0 !普遍気体定数 !飽和蒸気圧の温度微分 xz_DLogSvapPressDTemp = & & - a_ampA(ID) / (xz_Temp ** 2.0d0) & & + a_ampC(ID) / xz_Temp & & + a_ampD(ID) & & + a_ampE(ID) * 2.0d0 * xz_Temp !潜熱の計算 xz_LatentHeat = & & xz_DLogSvapPressDTemp * GasRUniv * (xz_Temp ** 2.0d0) & & / a_MolWt(ID) end function xz_LatentHeat !!!========================================================================== function SvapPress(ID, Temp) ! !引数で与えられた化学種と温度に対して, 飽和蒸気圧を計算. AMP 式を利用 ! !暗黙の型宣言禁止 implicit none !入出力変数 real(8) :: SvapPress !飽和蒸気圧 real(8), intent(in) :: Temp !温度 [K] integer, intent(in) :: ID !化学種の ID !内部変数 real(8) :: LogSvapPress !飽和蒸気圧の対数を計算 !対数が大きくなりすぎないようにする. ! Fujitsu Fortran Compiler では 700 より大きい数の exp を取ると警告が出る. LogSvapPress = & & min( & & ( & & a_ampA(ID) / Temp & & + a_ampB(ID) & & + a_ampC(ID) * dlog( Temp ) & & + a_ampD(ID) * Temp & & + a_ampE(ID) * ( temp ** 2 ) & & + dlog(1.0d-1) & & ), & & 7.0d2 & & ) !飽和蒸気圧を計算 SvapPress = dexp( LogSvapPress ) end function SvapPress !!!========================================================================== function LatentHeatPerMol(ID, Temp) ! !引数で与えられた化学種と温度に対して, 潜熱を計算 ! !暗黙の型宣言禁止 implicit none !入出力変数 real(8) :: LatentHeatPerMol !潜熱 real(8), intent(in) :: Temp !温度 integer, intent(in) :: ID !化学種名 !内部変数 real(8) :: DLogSvapPressDTemp real(8), parameter :: GasRUniv = 8.314d0 !普遍気体定数 !飽和蒸気圧の温度微分 DLogSvapPressDTemp = & & - a_ampA(ID) / (Temp ** 2.0d0) & & + a_ampC(ID) / Temp & & + a_ampD(ID) & & + a_ampE(ID) * 2.0d0 * Temp !潜熱の計算 LatentHeatPerMol = & & DLogSvapPressDTemp * GasRUniv * (Temp ** 2.0d0) end function LatentHeatPerMol end module ChemCalc