JP7133145B2

JP7133145B2 - シミュレーション方法およびシミュレーションプログラム

Info

Publication number: JP7133145B2
Application number: JP2018180716A
Authority: JP
Inventors: 孝介愛内
Original assignee: Nippon Electric Glass Co Ltd
Current assignee: Nippon Electric Glass Co Ltd
Priority date: 2018-09-26
Filing date: 2018-09-26
Publication date: 2022-09-08
Anticipated expiration: 2038-09-26
Also published as: JP2020050539A

Description

本発明は、ジュール熱の分布を求めるシミュレーション方法およびシミュレーションプログラムに関する。

交流電源により生じるジュール熱による加熱は多くの産業で利用される技術である。例えばガラス溶融工程においては、棒状電極や板状電極を溶融炉に設置し、溶融ガラス中に電流を流すことによってジュール熱を発生させ、ガラスを直接加熱する。高品位なガラス製品を生産するためには、溶融炉において溶融ガラスを適切な温度に加熱して、均質な溶融ガラスを生成することが求められる。そして、この要求を満たすには、電極の配置や電力条件等について最適な条件を選択することが求められ、これを予測するために数値シミュレーションが従来から使用されている。シミュレーションでは、対象とする領域をメッシュ、セル、コントロールボリューム等と呼ばれる微小領域に分割することで定まる計算点における電位を変数として、有限体積法、有限要素法、有限差分法等の方法により前記電位の関係を表す連立方程式を構築し、得られた電位の分布と対象領域を占める物質の導電率からジュール熱を求める。ガラス溶融工程では複数の電源が用いられることが多くあり、周波数の等しい複数の電源を考慮するシミュレーション方法が非特許文献１に記載されている。

P. Simons, K. Jochem, and K. Aiuchi, A power consistent mathematical formulation for Joulean heat release, Glass Technol.: Eur. J. Glass Sci. Technol. A, 49 (3), 109-118 (2008).

複数の電源が存在する場合、各電源において電流が保存される条件を課す必要がある。非特許文献１では、重ね合わせの理を適用することにより、各電源における電流保存条件を成立させる方法が記載されている。この方法では、（電源の数×２）の組の電位ポテンシャルの連立方程式を解く必要がある。ここで、前記の「×２」は交流の周期内において互いに１／４周期ずれた２つの時刻を表す。この方法により、周波数の等しい複数の電源を考慮したシミュレーションが可能になるが、電源が１つ追加されるごとに解かなければならない連立方程式が２組追加される。例えば電源の数が１０個では解かなければならない連立方程式は２０組となって、単一の電源の場合の２組に対して１０倍となる。すなわち、電源の数が多くなるほど計算時間が長くなり、大規模な溶融炉を対象とする場合は特にその負荷が大きくなる。

そこで、本発明の課題は、電源の数が増えても解くべき連立方程式の数の増大を抑制でき、これによって従来よりも計算時間の短いシミュレーション方法およびシミュレーションプログラムを提供することにある。

前記課題を解決するために創案された本発明に係るシミュレーション方法は、コンピュータが、周波数の等しい複数の交流電源により電流を流すことで２次元又は３次元の領域に生じるジュール熱の分布を求めるシミュレーション方法において、前記領域を分割して形成された複数の分割領域について、前記コンピュータが、前記分割領域のそれぞれの計算点における電位と、前記交流電源のそれぞれに定義した計算点における電位とを変数とし、前記分割領域ごとに成立する電流保存条件を表す方程式に、前記交流電源ごとに成立する電流保存条件を表す方程式を加えて連立方程式を構成し、前記コンピュータが、２つ以下の時刻について前記連立方程式を解くことを特徴とする。

ここで、分割領域とは、有限体積法、有限要素法、有限差分法等を用いたシミュレーション方法で、メッシュ、セル、コントロールボリューム等と呼ばれる微小領域のことである（以下、同様）。

上記構成では、コンピュータが解くべき連立方程式は、分割領域ごとに成立する電流保存条件を表す方程式に、交流電源ごとに成立する電流保存条件を表す方程式を加えたものである。従って、連立方程式を構成する方程式の数は、分割領域の計算点の数と、交流電源の数の和となる。ガラス溶融炉内のシミュレーションでは、一般に数万～数百万の分割領域の計算点数が用いられており、この数に比べると追加する方程式の数（交流電源の数）は僅かであり、この方程式を追加することによる計算時間の増大はほぼ無視できる。

用いられる複数の交流電源が、単相電源であって、全ての電源の位相が揃っている場合は、例えば交流電圧がピークとなる時刻で代表される１つの時刻における電位を求めることにより、ジュール熱を計算することができる。この場合、前記１つの時刻について連立方程式を解く。

用いられる複数の交流電源が、電源ごとに位相が揃っていない単相電源を組み合わせる場合や、三相電源を含む場合、スコット結線で形成される単相電源を含む場合などでは、交流の周期内において互いに１／４周期ずれた２つの時刻について連立方程式を解く。

上記の構成において、前記交流電源が、単相電源又は三相電源を含んでもよい。

上記の構成において、前記交流電源が、スコット結線で形成された単相電源を含んでもよい。

上記の構成において、前記コンピュータが、前記交流電源のそれぞれから供給される実効電力又は皮相電力が事前に設定した値となるように、前記交流電源のそれぞれの電圧を調整してもよい。

この構成によれば、所定の実効電力または皮相電力が加えられた状態を効率良く計算することができる。また、本発明のシミュレーション方法においては、各電源から供給される電力として、実効電力か皮相電力のいずれかを選択する機能を有することがより好ましい。

上記の構成において、前記領域の一部が、ガラス溶融炉内であってもよい。

また、前記課題を解決するために創案された本発明に係るシミュレーションプログラムは、請求項１に記載のシミュレーション方法をコンピュータに実行させることを特徴とする。

この構成でも、冒頭で説明したシミュレーション方法と同様の作用及び効果を享受できる。

本発明によるシミュレーション方法およびシミュレーションプログラムでは、電源の数が増えても解くべき連立方程式の数の増大を抑制でき、これによって従来よりも短い計算時間でジュール熱の分布を求めることが可能となる。

図１は、三相星形結線の電源を含む複数の電源を用いる場合の概念図である。図２は、電源の端子を接続する電極の境界フェイスおよび周辺のコントロールボリュームを表す垂直断面図である。図３は、実施例１に使用した、３次元の解析領域の平面図である。図４は、実施例１における電源の数と計算時間との関係を示した図である。

以下、本発明の実施形態に係るシミュレーション方法およびシミュレーションプログラムについて、説明する。

本実施形態に係るシミュレーション方法は、コンピュータが、周波数の等しい複数の交流電源により電流を流すことで２次元又は３次元の領域に生じるジュール熱の分布を求めるものである。前記領域を分割して形成された複数の分割領域について、前記コンピュータが、前記分割領域のそれぞれの計算点における電位と、前記交流電源のそれぞれに定義した計算点における電位とを変数とする。そして、前記分割領域ごとに成立する電流保存条件を表す方程式に、前記交流電源ごとに成立する電流保存条件を表す方程式を加えて連立方程式を構成する。そして、前記コンピュータが、２つ以下の時刻について前記連立方程式を解く。

また、本実施形態に係るシミュレーションプログラムは、上述したシミュレーション方法をコンピュータに実行させるものである。

ここでは三相電源を含む複数の電源を用いる場合を例として、本発明の実施形態を説明する。方程式を離散化する方法として、ここでは有限体積法を採用する。計算対象とする領域をコントロールボリュームに分割し、各コントロールボリュームにおける計算点の電位を、連立方程式の変数とする。

単一の交流起電力を生じる電源gがあるとする。任意の各位置において、gによって生じる電位は次式で表される。

ここで、fは交流の周波数、tは時刻、α_gは基準の位相に対するgのずれを表す角度、φ_g0はその位置におけるピーク電位である。交流起電力を生じる電源が複数ある場合、重ね合わせの理により、各位置による電位は次式で表される。

電流密度についても数１および数２と同様の式が成立するため、以下の数３、数４が得られる。

ここで、i_g0は単独のgの場合に生じる電流密度ベクトルのピークである。この時、i_Reおよびi_Imを以下のように定義する。

数５、数６を数４に代入すると次式を得る。

数７より、i_Reはt=(n+1/4)/f（ただしnは任意の整数）の時刻における各位置の電流密度、i_Imはt=n/f（ただしnは任意の整数）の時刻における各位置の電流密度をそれぞれ表す。

解析領域内の任意の体積領域について、体積領域の表面積を通過する電流の合計は、電流保存条件により０となる。これを数７に適用すると次式が得られる。

ここでΓは前記体積領域を囲む面であり、dΓは体積領域を囲む面のうちの微小領域に対し、外向きの法線方向で大きさが前記微小領域の面積のベクトル（以降面積ベクトルと呼ぶ）である。

全ての時刻tについて数８が成立するため、この条件を課すことにより、次の数９、数１０が得られる。

ここで数５、６におけるi_g0は、領域を占める物質の導電率γを用いて次式で表される。

数１１を数５、数６に代入し、さらにこれらを数９、数１０にそれぞれ代入して整理すると、数１２、数１３が得られる。

ここで、φ_Re、φ_Imを以下のように定義する。

数１４、数１５をそれぞれ数１２、数１３に代入すると、φ_Reおよびφ_Imに関する以下の方程式が得られる。

また、数１４および数１５を数２に適用すると次式を得る。

数１８より、φ_Reはt=(n+1/4)/f（ただしnは任意の整数）の時刻における各位置の電位、φ_Imはt=n/f（ただしnは任意の整数）の時刻における各位置の電位をそれぞれ表す。したがって、φ_Reとφ_Imは２つの時刻における電位を表す。

ここで、数５、数６、数１１、数１４、数１５を整理すると以下の関係が得られ、各時刻におけるi_Reおよびi_Imは、φ_Reおよびφ_Imから求めることができる。

解析領域を分割した全てのコントロールボリュームについて、数１６および数１７が成立するようにφ_Reとφ_Imの分布を決定することにより、全てのコントロールボリュームにおいて電流の保存を成立させることができる。しかしながら、複数の電源を用いる場合、全てのコントロールボリュームにおいて電流の保存が成立するだけでは不十分であり、これに加えて各電源において電流が保存される条件を課す必要がある。ここでは図１に示した概念図を用いて本発明の方法の一形態を説明する。

図１にはガラス溶融炉を模擬した３次元の解析領域の平面図１が示されている。解析領域には溶融ガラス２、溶融ガラスを囲む耐火物３、電極４a～４fからなる複数の電極４が含まれ、電極は底面から溶融ガラス内に所定の長さにて挿入されている。これらの領域は微小なコントロールボリュームに分割され、これにより前記計算点が定義される。この解析領域１に対して、ここでは複数の三相星形結線の電源を前記電極に接続する状態を考える。図１には第一の三相星形結線の電源５と第二の三相星形結線の電源６が示されているが、ここでこれらの電源および接続するための導線は前記３次元の解析領域には含まれない。前記三相星形結線の電源５は、中性点をP1として、端子R1の位相が0°、端子S1の位相が120°、端子T1の位相が240°となるよう構成されている。前記電源６についても中性点をP2として、端子R2、S2、T2が同様の関係になるよう構成されている。電極４a、４b、４cには電源５の端子R1、S1、T1をそれぞれ接続し、電極４d、４e、４fには電源６の端子R2、S2、T2をそれぞれ接続する。より具体的には、電極４aの場合、例えば電極の下端に位置する複数のコントロールボリュームの下側の境界フェイスに対して、端子R1が接続された状態を想定し、この扱いは他の電極でも同様である。

ここで、電源５を接続した電極４a、４b、４cに対して電流保存条件を課す方法の一形態を説明する。本発明を適用するため、ここでは電源５における中性点P1を計算点とし、前記２つの時刻における中性点P1の電位に対応する２つの変数をφ_P1__Re、φ_{P1_Im}として定義する。この時、端子R1、S1、T1における２つの時刻の電位は、以下の式で表される。

ここで、V₁は第一の三相星形結線の電源５における相電圧のピーク値を表す。α_R、α_S、α_Tはそれぞれ0°、120°、240°である。

数１６および数１７は各コントロールボリュームにおいて成立する式であるが、電極端子が接続されるコントロールボリュームにおいてこれらの式を適用する方法を以下に示す。ここで、数１６および数１７は、有限体積法により離散化した表現として、次の数２７および数２８となる。

ここで、ΔΓは前記各コントロールボリュームを囲むフェイスΓの面積ベクトルである。

図２は、電源の端子を接続する電極の境界フェイスおよび周辺のコントロールボリュームを表す垂直断面図である。最下部に電源の端子を接続する電極の境界フェイスがある。図２には前記境界フェイスの中心点８と、この境界フェイスに隣接するコントロールボリュームの中心点９が示されている。連立方程式の変数となる電位は、このコントロールボリューム中心点（計算点）において定義される。なお、電源の端子を接続する電極の境界フェイスは一つ以上存在するが、ここではその一つに注目して示している。この境界フェイスに隣接するコントロールボリュームについて成立する数２７および数２８を、連立方程式を構成する変数を用いて表現する必要がある。この時、数２７および数２８における各フェイスΓの寄与を表現するのに、電源の端子を接続するフェイス以外の場合は、従来の方法を用いれば良く、例えばフェイスの両側にコントロールボリュームが存在するフェイスの場合は前記両側のコントロールボリュームの中心点におけるφ_Reまたはφ_Imとγおよび中心点間距離等のコントロールボリューム形状に関わる定数を用いて、数２７または数２８内のγ▽φ_Re・ΔΓまたはγ▽φ_Im・ΔΓを表すことができる。

数２７および数２８においてフェイスΓが電源の端子を接続するフェイスである場合として、ここでは端子R1が接続され、このフェイスをΓ_R1_iと呼ぶこととする。ここでのiは端子R1に接続される複数のフェイスのうちのi番目であることを意味している。数２７および数２８内のΓ_R1_iの寄与を表す(γ▽φ_Re・ΔΓ)Γ__{R1_i}および(γ▽φ_Im・ΔΓ)Γ__{R1_i}は、以下の数２９および数３０で表される。

ここで、ΔΓ_iはフェイスΓ_R1_iの面積である。Δd_iはコントロールボリューム中心点９からフェイスΓ_R1_iまでの距離であり、直交格子の場合は前記コントロールボリューム中心点９からフェイス中心点８までの距離に等しい。φ_{R1_i_nb_Re}およびφ_{R1_i_nb_Im}はコントロールボリューム中心点９における電位を表す。数２９および数３０に数２１および数２２をそれぞれ代入すると、次の数３１、数３２になる。

数３１および数３２を数２７および数２８にそれぞれ代入することにより、追加された変数φ_{P1_Re}およびφ_{P1_Im}を含む変数φ_Reおよびφ_Imの組み合わせにより、端子R1が接続されたコントロールボリュームにおける方程式を得ることができる。端子S1、T1等の他の端子が接続されたコントロールボリュームについても、同様の操作によりそれぞれのコントロールボリュームで成立する方程式を得ることができる。

ここで、数３１および数３２の左辺は、それぞれの時刻においてフェイスΓ_R1_iを通過する電流に-1を乗じたものを表している。

このように、数２７と数２８により、コントロールボリュームごとに成立する電流保存条件を表す方程式が得られる。

これに対し、各電源における電流保存条件を成立させるには、電源の各端子に接続された全てのフェイスを通過する電流の合計が0になる必要がある。これにより、本発明において、新たに加える電流保存条件を表す方程式は、例えば電源５について以下の式となる。

ここで、i_R1、i_S1、i_T1はそれぞれ端子R1、S1、T1に接続されたフェイスに隣接するコントロールボリュームiについて積算することを意味する。数３３または数３４は、それぞれ各時刻における連立方程式（コントロールボリュームごとに成立する電流保存条件を表す方程式で構成された連立方程式）に追加され、新たな連立方程式を構成する。この新たな連立方程式を解くことにより、電源における電流保存条件を満たす各時刻の電位分布を得ることができる。

さらに電源を追加する場合も、同様の扱いにより、各電源における電流保存条件を成立させることが可能である。なお、上記の説明において、三相星形結線における中性点を計算点としたが、本発明はこれに限定されるものではない。例えば電源５において、端子R1またはS1またはT1を計算点として、端子間の電位差の関係を用いて新たな方程式を構築することも可能である。また、本発明の方法は、三相三角結線の電源や、単相電源、三相電源からスコット結線により形成された単相電源にも適用することができる。

また、本発明の方法によれば、連立方程式の係数行列を対称行列とすることができる。このため、ICCG法などの対称行列専用の効率の良い連立方程式解法を適用することができる。

得られた２つの時刻における電位分布から、各コントロールボリュームにおけるジュール熱Qは、交流周期の時間平均のジュール熱として次式で表される。

ここで、i_Reおよびi_imは、既に得られているφ_Reおよびφ_Imを用いて、数１９および数２０により計算できる。これを数３５に適用すれば、各位置におけるジュール熱が得られる。

上記の方法では、各電源におけるピーク電圧（例えば電源５におけるV₁）を事前に設定し、この条件におけるジュール熱の分布が得られる。

以上のように構成された本実施形態のシミュレーション方法では、コンピュータが解くべき連立方程式は、分割領域ごとに成立する電流保存条件を表す方程式（数２７と数２８）に、交流電源ごとに成立する電流保存条件を表す方程式（数３３と数３４）を加えたものである。従って、連立方程式を構成する方程式の数は、分割領域の計算点の数と、交流電源の数の和となる。ガラス溶融炉内のシミュレーションでは一般に数万～数百万の分割領域の計算点数が用いられており、この数に比べると追加する方程式の数（交流電源の数）は僅かであり、この方程式を追加することによる計算時間の増大はほぼ無視できる。

次に本発明の別の実施形態として、各電源から供給される実効電力または皮相電力が事前に設定した値となる状態を求める方法について説明する。このためには、上述した電位分布計算の実施後、各電源から供給される実効電力または皮相電力を計算する。供給される実効電力または皮相電力が事前に設定した値より小さければ、各電源におけるピーク電圧を上昇させる調整を行い、逆に大きければ、各電源におけるピーク電圧を低下させる調整を行う。この調整を、十分誤差が小さくなるまで繰り返し実行すればよい。

各電源に供給される実効電力および皮相電力の計算は、例えば前記三相星形結線の電源５の場合は次のように計算できる。まず、数７を以下のように変形する。

ただし、βは次式で表される。

数３６より、実効電流は次式により求めることができる。

各端子R1、S1、T1における実効電流を数３８で求め、各端子の相電圧の実効値であるV1/(√2)をこれに乗じることにより、各端子の皮相電力が得られ、全端子について合計すれば、電源５の皮相電力が得られる。

また、電流と電圧の位相差は数３７と各端子の相電圧の位相との差により求めることができる。この位相差をδとすると、皮相電力にcosδを乗じたものが実効電力となり、これを全端子について合計すれば、電源５の実効電力が得られる。

図３は３次元の解析領域の平面図であり、縦４ｍ、横１１ｍ、紙面奥行き方向を表す高さ１ｍの寸法となっている。領域内には電極１１が縦方向に３本、横方向に10本並び、１ｍ間隔の等間隔で合計３０本設置されている。ここで、電極は高さ方向全体を占める円柱形状であって、導電率が3×10⁶S/mの導体である。電極を除く領域は1300℃のソーダ石灰ガラスに相当する溶融ガラス１０で満たされ、この溶融ガラスの導電率は16S/mである。

この時、電極（Ra、Sa、Ta）の３本を１番目の組とし、同様に２番目、３番目、・・・、１０番目の組を電極（Rb、Sb、Tb）、（Rc、Sc、Tc）、・・・、（Rj、Sj、Tj）とする。各組にそれぞれ独立した三相星形結線の電源を割り当て、各電源における３つの端子に１本ずつ電極を接続することを考える。このとき電源の数は１０個となるが、１０個の電源のうち１つだけ有効にした場合、２つ有効にした場合、５つ有効にした場合、１０個全て有効にした場合それぞれについて、本発明を適用してジュール熱の計算を行い、計算時間を比較した。いずれのケースでも一つの三相星形結線の電源の電力は10kWである。ここでは、解析領域を51.7万のコントロールボリュームに分割して有限体積法による離散化を行った。連立方程式の解法には並列計算に対応したICCG法を用い、２０並列により計算を実施した。計算時間は３回の平均により評価した。

評価結果を図４と表１に示す。図４は、有効にした電源の数と計算時間の関係を示したものである。また、表１は、図４の数値を示す。

電源の数によって若干計算時間にばらつきがあるのは、電源の数によって電位分布やジュール熱の分布が異なるため、連立方程式計算における収束に要する反復回数に差が生じることが影響したものと考えられるが、この結果より電源の数を増やしても概ね計算時間は変わらないことが分かる。従来法によれば電源の数を増やすごとに解くべき連立方程式の組が増えるため、計算時間は電源の数に概ね比例して増大することになるが、本発明では電源の数を増やしても計算時間はほぼ変わらず、多数の電源が用いられる場合でも効率良く計算を行うことができると言える。

１ガラス溶融炉を模擬した３次元の解析領域の平面図
２溶融ガラス
３耐火物
４電極
４a～４f 各電極
５第一の三相星形結線の電源
６第二の三相星形結線の電源
７コントロールボリューム
８電源の端子を接続する電極の境界フェイス中心点
９境界に接するコントロールボリュームの中心点
１０溶融ガラス
１１電極

Claims

コンピュータが、周波数の等しい複数の交流電源により電流を流すことで２次元又は３次元の領域に生じるジュール熱の分布を求めるシミュレーション方法において、
前記領域を分割して形成された複数の分割領域について、前記コンピュータが、前記分割領域のそれぞれの計算点における電位と、前記交流電源のそれぞれに定義した計算点における電位とを変数とし、前記分割領域ごとに成立する電流保存条件を表す方程式に、前記交流電源ごとに成立する電流保存条件を表す方程式を加えて連立方程式を構成し、
前記コンピュータが、２つ以下の時刻について前記連立方程式を解くことを特徴とするシミュレーション方法。
前記交流電源が、単相電源又は三相電源を含むことを特徴とする請求項１に記載のシミュレーション方法。
前記交流電源が、スコット結線で形成された単相電源を含むことを特徴とする請求項１に記載のシミュレーション方法。
前記コンピュータが、前記交流電源のそれぞれから供給される実効電力又は皮相電力が事前に設定した値となるように、前記交流電源のそれぞれの電圧を調整することを特徴とする請求項１～３の何れか１項に記載のシミュレーション方法。
前記領域の一部が、ガラス溶融炉内であることを特徴とする請求項１～４の何れか１項に記載のシミュレーション方法。
請求項１に記載のシミュレーション方法をコンピュータに実行させることを特徴とするシミュレーションプログラム。