コンピュータが、流路の構成に関するデータと前記流路を流動する複数の粉体の物理量に関するデータとを少なくとも含む初期条件データを取得する取得ステップと、コンピュータが、前記初期条件データを用い、前記複数の粉体の混合比の前記流路内での分布の時間変化を数値計算する計算ステップと、コンピュータが、前記計算ステップの計算結果を出力装置に出力する出力ステップと、を含み、前記計算ステップでは、式1を用いて、各時刻での、前記流路内の各位置における前記複数の粉体の混合比が計算されることを特徴とする混合比計算方法。前記計算ステップは、各時刻で、前記流路内の各位置における前記粉体の流速分布および粉体存在密度を更新する更新ステップを含んでおり、前記更新された流速分布および粉体存在密度を用いて、前記複数の粉体の混合比が計算されることを特徴とする請求項1に記載の混合比計算方法。前記更新ステップでは、前記流路内の粉体の流動を表す連続の式と、前記流路内の粉体の流動性能特性とを用いて、前記粉体の流速分布および粉体存在密度が計算されることを特徴とする請求項2に記載の混合比計算方法。前記流路は、電子写真方式の画像形成装置における二成分現像器内の流路であり、前記複数の粉体は、二成分現像剤におけるトナーとキャリアであることを特徴とする請求項1~3のうちいずれか1項に記載の混合比計算方法。流路の構成に関するデータと前記流路を流動する複数の粉体の物理量に関するデータとを少なくとも含む初期条件データを取得する取得手段と、前記初期条件データを用い、前記複数の粉体の混合比の前記流路内での分布の時間変化を数値計算する計算手段と、前記計算手段の計算結果を出力する出力手段と、を含み、前記計算手段は、式1を用いて、各時刻での、前記流路内の各位置における前記複数の粉体の混合比を計算することを特徴とする混合比計算装置。前記計算手段は、各時刻で、前記流路内の各位置における前記粉体の流速分布および粉体存在密度を更新する更新手段を含んでおり、前記更新された流速分布および粉体存在密度を用いて、前記複数の粉体の混合比を計算することを特徴とする請求項5に記載の混合比計算装置。前記更新手段は、前記流路内の粉体の流動を表す連続の式と、前記流路内の粉体の流動性能特性とを用いて、前記粉体の流速分布および粉体存在密度を計算することを特徴とする請求項6に記載の混合比計算装置。前記流路は、電子写真方式の画像形成装置における二成分現像器内の流路であり、前記複数の粉体は、二成分現像剤におけるトナーとキャリアであることを特徴とする請求項5~7のうちいずれか1項に記載の混合比計算装置。請求項1~4のうちいずれか1項に記載の混合比計算方法の各ステップをコンピュータに実行させることを特徴とするプログラム。従来技術においては、現像器内の現像剤質量密度の空間分布がいたるところで一定と仮定をすることによって、移流拡散方程式の未知数をトナー濃度としている。しかしながら、実際の現像器においては、スクリューピッチや形状が異なる箇所が存在する場合がある。そのような構成においては現像剤質量密度の空間分布は一様ではなくなる。また、スリーブによって現像剤が汲み上げられる部分が存在する現像室と、そのような部分が存在しない撹拌室とでは、搬送に寄与する現像剤の空間分布が異なることは容易に想像できる。また、現像器内の現像剤質量密度分布が時間変化する場合もある。トナー補給や消費に伴う現像剤質量増減や現像剤排出機構がある場合などである。より詳細に現像器内のトナー濃度分布の経時変化を計算するためには、現像剤質量密度の空間分布が一定でない場合を考慮する必要がある。従来技術においては、そのような場合にはトナー濃度分布の経時変化を計算することができなかった。従って、本発明の目的は、粉体の存在密度の空間分布が一様ではない場合に、流路内の粉体の混合比分布の経時変化を計算するための粉体の混合比計算方法及び装置を提供することである。本発明の第1態様は、コンピュータにより実行される混合比計算方法であって、コンピュータが、流路の構成に関するデータと前記流路を流動する複数の粉体の物理量に関するデータとを少なくとも含む初期条件データを取得する取得ステップと、コンピュータが、前記初期条件データを用い、前記複数の粉体の混合比の前記流路内での分布の時間変化を数値計算する計算ステップと、コンピュータが、前記計算ステップの計算結果を出力装置に出力する出力ステップと、を含み、前記計算ステップでは、式1を用いて、各時刻での、前記流路内の各位置における前記複数の粉体の混合比が計算されることを特徴とする混合比計算方法である。また本発明の第2態様は、混合比計算装置であって、流路の構成に関するデータと前記流路を流動する複数の粉体の物理量に関するデータとを少なくとも含む初期条件データを取得する取得手段と、前記初期条件データを用い、前記複数の粉体の混合比の前記流路内での分布の時間変化を数値計算する計算手段と、前記計算手段の計算結果を出力する出力手段と、を含み、前記計算手段は、前記式1を用いて、各時刻での、前記流路内の各位置における前記複数の粉体の混合比を計算することを特徴とする混合比計算装置である。また本発明の第3態様は、混合比計算のプログラムであって、上述した混合比計算方法の各ステップをコンピュータに実行させることを特徴とするプログラムである。本発明によれば、粉体の存在密度の空間分布が一様でない場合に、流路内の粉体の混合比分布の経時変化を計算することが可能となる。また、粉体の存在密度の空間分布を計算で求めて、流路内の粉体の混合比分布の経時変化を計算することが可能となる。[実施例1]以下に、本発明を実施するための形態について必要に応じて図面を参照しつつ詳しく説明する。実施例1においては、電子写真方式の画像形成装置における二成分現像容器内を循環する現像剤の質量密度に既知の空間分布がある場合に、現像剤を連続体とみなして移流拡散方程式を解いてトナー濃度分布の経時変化を求める方法を代表例として説明する。本発明者らは、二成分現像器内ではキャリアの流れが支配的であり、トナーがキャリアとは別のふるまいをすることなく、キャリアに追従して流れていることを見出した。そこで、本実施例では、移流拡散方程式の未知数を式2のようにトナー濃度T[wt]と現像剤存在密度m[g/m]の積として解き、トナー濃度分布の経時変化を求める。なお、vは流速ベクトル、Dは拡散係数である。上記の計算手法は、二成分現像剤のトナー濃度分布の経時変化の計算だけでなく、ある流路を流動する複数二種類以上の粉体の混合比Cの空間時間変化の計算にも適用可能である。以下の式1は、式2を一般化したものであり、mは粉体存在密度、vは流速ベクトル、Dは拡散係数である。本実施例の混合比計算装置の構成を図1に示す。本装置は、CPU100、RAM101、ディスプレイやプリンタなどの出力装置102、キーボードやマウスなどの入力部103、ハードディスクなどの外部記憶装置104及びバス105を備えるコンピュータにより構成される。更に、上記CPU100内部において、符号100aは制御部であり、プログラム全体を制御する。符号100bは初期条件設定部であり、初期条件データを設定する。符号100cはトナー濃度計算部であり、移流拡散方程式を解き、トナー濃度を計算する。制御部100a、初期条件設定部100b、トナー濃度計算部100cは、プログラムを実行するCPU100によって実現される機能である。また、上記RAM101は、プログラム101a、現像器構成データ101b、計算条件データ101c、物理量データ101d、トナー濃度分布データ101eを格納する。これらのプログラムおよび初期設定として与えられるデータは、外部記憶装置104または不図示のROMなどの不揮発性の記憶媒体に格納されており、混合比計算を行う際にRAM101にロードされる。なお、本実施例では、初期条件設定部100b、トナー濃度計算部100c、出力装置102が、それぞれ、本発明の取得手段、計算手段、出力手段に対応する。ここで、各データの内容を説明する。現像器構成データ101aとは、寸法、スクリュー構成、スリーブ長さなど、現像器構成に関する値である。計算条件データ101cとは、初期トナー濃度、トナー補給消費のシーケンス、メッシュ数、タイムステップなど、計算条件に関する値である。物理量データ101dとは、流速、拡散係数、現像剤質量分布など、物理量の値である。トナー濃度分布データ101eとは、タイムステップ毎に計算されたトナー濃度分布を逐次更新していく値である。以下、実施例1における二成分現像剤のトナー濃度計算方法の一例について説明する。具体的にはトナーの補給や消費がない場合に、現像器中の現像室と撹拌室にそれぞれトナー濃度の異なる現像剤を初期設定として配置し、本発明を適用して現像剤のトナー濃度分布の経時変化を計算する場合について説明する。計算対象の現像器は、スクリューピッチの異なる領域が3つあり、搬送状態がそれぞれの領域で異なっている。現像剤質量密度の分布は定常状態として、実験値を用いている。図2は本発明の計算方法の流れを表すフローチャートである。同図を用いて現像剤のトナー濃度分布計算方法の流れを説明する。1初期条件設定部100bが、現像器の流路長やスクリュー部分の長さ、スリーブの長さなどの寸法を設定し、スクリューピッチやスクリューの種類などのスクリュー構成を設定して計算メッシュを定義する。また初期条件設定部100bは、初期トナー濃度、トナー補給消費シーケンス、タイムステップ、計算結果描画ステップ、計算終了時間などの計算条件を設定する。さらに初期条件設定部100bは、各メッシュに流速、拡散係数、現像剤質量密度分布を設定するステップS201。図3aに現像剤質量分布を示す。スクリューの違いに応じて現像剤の空間分布が変化している。図中、中心部分の点線は現像容器内の現像室と撹拌室とを仕切る線であり、線の左側は撹拌室であり、右側は現像室に対応している。図3bに初期トナー濃度を示す。撹拌室の初期トナー濃度を4wt、現像室の初期トナー濃度を12wtに設定している。図3cに現像剤質量とトナー濃度から計算されるトナー質量分布を示す。現像剤質量の空間分布とトナー濃度の空間分布に応じて、トナー質量も空間分布を持っている。2次に、トナー濃度計算部100cが、式2の移流拡散方程式を計算するステップS202。差分法には数値拡散を避けるため、CIPConstrainedInterpolationProfile法を使用する。3次に、制御部100aは、計算時間が1において設定した計算終了時間になっているか判定し、偽であれば時刻を更新し、2に戻るステップS203、S204。計算時間が計算終了時間を超えていれば、計算を終了するステップS205。計算結果は、ディスプレイやプリンタなどの出力装置102に出力される。図4a~図4dは、トナー濃度分布計算結果の一例である。縦軸はトナー濃度であり、横軸は位置を表している。スクリューは左から右方向へ現像剤を搬送しており、最右端の次のメッシュが最左端のメッシュと対応しており、循環するようになっている。図4a~図4dのグラフは、各々、初期状態0sec、10sec後、20sec後、30sec後のトナー濃度分布を示している。初期状態の分布から全体的にトナー濃度分布が変動していき、時間が経つほどトナー濃度の偏りがなくなり、トナー濃度がフラットに近づいているのがわかる。さらに、図5a~図5cにトナー濃度分布の時間変化が収束し、定常状態となったときの現像剤量分布、トナー濃度分布、トナー質量分布を示す。時間が十分に経過すると、トナー濃度は8.32に収束している。また、トナー質量分布は現像剤質量分布に基づいて空間分布を持っている。トナー質量の積算値は計算前と後で不変となっている。このように、初期状態から定常状態までの間で、トナー濃度がフラットになるということと、トナー総質量が保存されているという二点の物理現象を本実施例において表すことができた。また、トナーを補給あるいは消費させたければ、移流拡散計算のステップS202の後に、所定の箇所補給ならば補給領域、消費ならば現像領域に対応するメッシュに対して、所望のトナー質量を補給ならば加え、消費ならば差し引いてから、トナー濃度を計算しなおせばよい。このようにすることによって、トナー補給消費時においてもトナー濃度分布の経時変化を計算することが可能である。なお、本実施例においては、現像剤の存在密度として質量密度を用いているが、現像剤の個数密度などを用いてもよい。なお、本実施例においては、現像器内の流路が現像室と撹拌室を循環する1次元で記述される流路を例に説明したが、縦撹拌型の現像器のように、流路が分岐や合流する場合に対しても適用可能である。また、流路を2次元や3次元にしてもよい。これにより、より複雑な流路を持つ現像器に対しても容易に適用することができる。なお、本実施例においては、移流拡散方程式を解く際にCIP法を用いたが、数値計算の差分スキームは多々有り、風上差分、Leap-Frogスキーム、Lax-Wendoroffスキーム、流束制限法などの他の差分スキームを用いても構わない。また、本実施例においては、電子写真の二成分現像器を例に説明したが、流路が規定できる粉体搬送装置に対しても容易に適用可能であることは言うまでもない。[実施例2]実施例1では、現像剤質量密度分布を既知としてトナー濃度分布の経時変化を求める方法を説明した。実施例2においては、現像剤質量密度分布に時間変化がある場合に、トナー濃度分布の経時変化を計算する方法を説明する。本実施例の計算方法の特徴は、実施例1では現像剤質量密度分布を既知としてトナー濃度分布の経時変化を求めていたのに対して、非定常状態の現像剤質量密度分布とトナー濃度分布の経時変化を連成して計算することである。このことについて詳しく説明する。非定常状態の粉体の質量密度分布は、存在密度に対する拡散項を考慮した非定常の連続の式を用いる。粉体の存在密度に対する拡散項を考慮した連続の式は式4で表される。ここで、ρは質量密度を、vは流速ベクトルを、Dは拡散係数を表す。一方、粉体の流量性能を表す流動性能特性は、具体的には質量密度ρと流量ρvの関係を用いることとする。この関係は、質量密度が決まったときの流量が規定されることを意味する。そこで、式4と流動性能特性を連成して解くことにより、流路内の存在量分布ρと速度分布vの時間変化を求めることが可能となる。特に、連続の式に現像剤の質量密度に対する拡散項を考慮することで、非平衡な現像剤の質量密度分布が徐々に平衡状態に変化する現象を再現できる。また、その時定数は拡散係数Dで考慮可能である。この非定常状態の粉体の質量密度分布計算方法と、本発明のトナー濃度分布の経時変化計算方法を連成して解くことで、現像剤質量密度分布に時間変化がある場合でも、トナー濃度分布の経時変化を計算することが可能となる。本実施例の混合比計算装置の構成を図6に示す。図6中、符号100~105、100a~100c、及び101a~101eは図1と同じであることから、説明を省略する。図中、符号601は連続の式の計算部であり、式4を計算することで1タイムステップ後の質量密度を求める機能である。符号602は流速分布の修正部であり、流動性能特性を用いて流速分布を修正する機能である。RAM101に流動性能特性データ603、質量密度分布データ604、流速分布データ605を格納する。本実施例では、連続の式の計算部601および流速分布の修正部602が、本発明の更新手段に対応する。ここで各データの内容を説明する。流動性能特性データ603は、流速と質量密度の積で表される流量と質量密度の関係であり、質量密度分布データ604は流路内の現像剤の質量密度分布である。流速分布データ605は、流路内の現像剤の流速分布である。以下、実施例2における二成分現像剤のトナー濃度計算方法の一例について説明する。具体的には実施例1と同様のハードウェア構成、プロセス条件、計算条件であるが、現像剤質量密度分布とトナー濃度分布の経時変化を連成して計算する点が実施例1の方法と異なる。図7は計算方法の流れを表すフローチャートである。同図を用いて現像剤のトナー濃度分布計算方法の流れを説明する。1初期条件設定部100bが、現像器の流路長やスクリュー部分の長さ、スリーブの長さなどの寸法を設定し、スクリューピッチやスクリューの種類などのスクリュー構成を設定して計算メッシュを定義する。また初期条件設定部100bは、初期トナー濃度、トナー補給消費シーケンス、タイムステップ、計算結果描画ステップ、計算終了時間などの計算条件を設定する。さらに初期条件設定部100bは、各メッシュに拡散係数を設定する。さらに初期条件設定部100bは、所望の現像剤総質量、流速と質量の積で表される流量と質量密度の関係で表される流動性能特性をスクリューの種類ごとに設定する。また初期条件設定部100bは、質量密度と流量と流速の初期分布を設定するステップS201。2次に、連続の式の計算部601が、式4を数値計算することで、各位置での質量密度を求めるステップS701。数値計算手法としては、数値拡散を避けるため、CIP法を使用する。3次に、流速分布の修正部602が、質量密度分布と流動性能特性から、各位置での流速分布を求めるステップS702。4次に、トナー濃度計算部100cが、式2の移流拡散方程式を計算するステップS202。5次に、制御部100aは、計算時間が1において設定した計算終了時間になっているか判定し、偽であれば時刻を更新し、2に戻るステップS203、S204。計算時間が計算終了時間を超えていれば、計算を終了するステップS205。計算結果は、ディスプレイやプリンタなどの出力装置102に出力される。このように、現像剤質量密度分布とトナー濃度分布の経時変化を連成して求めることによって、トナーの補給又は消費や、現像剤排出機構などにより現像剤質量密度分布に時間変化がある場合でも、トナー濃度分布の経時変化を計算することが可能となる。なお、本実施例においては、連続の式と流動性能特性を連成することで現像剤の質量密度を求めているが、粉体を流体として考慮し、ナビエストークス方程式を使用して求めても良い。また、個別要素法で現像器全体の現像剤の挙動を解いて現像剤質量密度分布を求めても構わない。なお、本実施例においては、現像剤質量密度分布を非定常状態として求めたが、現像剤質量密度分布が定常状態とみなせる場合には、初期設定後に一度だけ定常状態の現像剤質量密度分布を求めればよい。なお、本実施例においては、現像剤の存在密度として質量密度を用いているが、現像剤の個数密度などを用いても良い。また、粉体の存在密度と流量の関係を表す流動性能特性は、実験や別の現像剤流動シミュレーションで取得した特性を用いればよい。また、流動性能特性は、存在密度と流速の2変数の関係を表すので、存在密度と流量の代わりに、存在密度と流速、流量と流速でもよい。なお、本実施例においては、現像器内の流路が現像室と撹拌室を循環する1次元で記述される流路を例に説明したが、縦撹拌型の現像器のように、流路が分岐や合流する場合に対しても適用可能である。また、流路を2次元や3次元にしてもよい。これにより、より複雑な流路を持つ現像器に対しても容易に適用することができる。なお、本実施例においては、移流拡散方程式を解く際にCIP法を用いたが、数値計算の差分スキームは多々有り、風上差分、Leap-Frogスキーム、Lax-Wendoroffスキーム、流束制限法などの他の差分スキームを用いても構わない。また、本実施例においては、電子写真の二成分現像器を例に説明したが、流路が規定できる粉体搬送装置に対しても容易に適用可能であることは言うまでもない。[実施例3]実施例3においては、食品の粉体混合工程において、本発明を適用する例を説明する。図8に食品の粉体混合装置の一例を示す。この装置は麺類を製造するために小麦粉とデンプンの混合物とアルギン酸を混合するための装置である。粉体供給ホッパー801に小麦粉とデンプンが充填されており搬送スクリュー802によって搬送される。また、粉体供給ホッパー803にはアルギン酸が充填されており、搬送スクリュー804によって搬送される。これら各粉体は主流路806に移動し、主流路806に移動した粉体は空気輸送によって搬送される。そして、混合機805で全ての粉体が混合される。このような装置において、本発明を適用することが可能である。具体的には、混合機805において所望の粉体混合比を得るためのハードウェア構成や運転条件を検討することができる。ハードウェア構成、運転条件とは、搬送スクリュー802や搬送スクリュー804の長さや搬送速度、搬送タイミング等である。その検討方法としては、流路を設定し、実施例1と同様、式1の移流拡散方程式を繰り返し解く。このように計算することで、各粉体の濃度の時間空間変化を求めることができる。これにより混合機805での粉体濃度を予測することができ、ハードウェア構成や運転条件を変えて計算し、所望の粉体濃度を得られる条件を探索することができる。主流路806は常に粉体質量密度が一様とは限らないため本発明は有効である。なお、実施例3においては、食品の粉体混合工程に関して説明したが、食品以外にも、薬品、化粧品、金属、化学、顔料、土石、石陶鉱、鋳物砂などの粉体を扱う工業において、類似の装置に対して本発明を適用することが可能である。