WO2022210325A1

WO2022210325A1 - 推定装置、推定方法、及びコンピュータプログラム

Info

Publication number: WO2022210325A1
Application number: PCT/JP2022/014251
Authority: WO
Inventors: 洋輔岡部; 元信村上; 翔太山内; 茂樹山手
Original assignee: 株式会社Ｇｓユアサ
Priority date: 2021-03-31
Filing date: 2022-03-25
Publication date: 2022-10-06
Also published as: DE112022001905T5

Abstract

推定装置、推定方法、及びコンピュータプログラムの提供。蓄電素子に発生する歪みに係るデータを取得する取得部と、蓄電素子の内部の力学的状態を表すシミュレーションモデルを用いて、取得部が取得したデータに基づき、蓄電素子の内部応力を推定する推定部とを備える。

Description

推定装置、推定方法、及びコンピュータプログラム

　本発明は、推定装置、推定方法、及びコンピュータプログラムに関する。

　近年、リチウムイオン電池などの蓄電素子は、ノート型パーソナルコンピュータ、スマートフォンなどの携帯端末の電源、再生可能エネルギ蓄電システム、ＩｏＴデバイス電源など、幅広い分野において使用されている。

　リチウムイオン電池の開発は高容量、高エネルギ密度化を目指して進められており、新規の電極材料の探索が進められている。

特開２０１６－２０７３１８号公報特開２０１９－０９１６１５号公報

　多くの場合、高容量、高エネルギ密度を達成できる電極材料は、充放電や劣化にともない、体積膨張をすることが知られている（例えば、特許文献１、２を参照）。通常、電極は何らかの筐体内に置かれ、力学的な拘束を受けている。そのため、体積膨張によって電池の内部応力が発生する。

　内部応力は、電池の内部抵抗や、反応生成物の析出などの電池の特性に影響を及ぼすことが知られている。しかし、内部応力と電池の特性を関連付けて電池の性能評価や状態監視をおこなう方法については、提案されていない。

　本発明は、斯かる事情に鑑みてなされたものであり、電池内部の挙動の１つとして電池内部に発生する内部応力を推定する推定装置、推定方法、及びコンピュータプログラムを提供することを目的とする。

　推定装置は、蓄電素子に発生する歪みに係るデータを取得する取得部と、前記蓄電素子の内部の力学的状態を表すシミュレーションモデルを用いて、前記取得部が取得したデータに基づき、前記蓄電素子の内部応力を推定する推定部とを備える。

　推定方法は、蓄電素子に発生する歪みに係るデータを取得し、前記蓄電素子の内部の力学的状態を表すシミュレーションモデルを用いて、前記取得部が取得したデータに基づき、前記蓄電素子の内部応力を推定する処理をコンピュータによって実行する。

　上記構成によれば、電池内部の挙動の１つとして電池内部に発生する内部応力を推定できる。

実施の形態１に係る推定システムの全体構成を示す模式図である。蓄電素子の構成を説明する説明図である。固体電解質層の内部構成を説明する説明図である。推定装置の内部構成を示すブロック図である。実施の形態１における内部応力の推定手順を説明するフローチャートである。蓄電素子の内部応力とオーム抵抗との関係を示すグラフである。等価回路モデルの一例を示す回路図である。実施の形態５に係る推定装置の内部構成を示すブロック図である。有効拡散係数の計算手法を説明する説明図である。全固体電池における有効拡散係数の計算手法を説明する説明図である。全固体電池における有効拡散係数の計算結果を示す図である。全固体電池における有効拡散係数の計算結果を示す図である。全固体電池における有効拡散係数の計算結果を示す図である。接触面積比率と有効拡散係数との関係を示すグラフである。実施の形態５に係る推定装置が実行する演算処理の実行手順を説明するフローチャートである。実施の形態６に係る推定装置が実行する演算処理の実行手順を説明するフローチャートである。蓄電素子の構成を説明する説明図である。析出物の生成量と固有歪みとの関係の一例を示すグラフである。析出物の生成量と固有歪みの関係の他の例を示すグラフである。応力・歪み分布の計算手順を示すフローチャートである。等価回路モデルの一例を示す回路図である。

　推定装置は、蓄電素子に発生する歪みに係るデータを取得する取得部と、前記蓄電素子の内部の力学的状態を表すシミュレーションモデルを用いて、前記取得部が取得したデータに基づき、前記蓄電素子の内部応力を推定する推定部とを備える。
　歪みに係るデータは、歪みセンサによる計測データであってもよい。
　この構成によれば、蓄電素子に発生する歪みのデータから、直接的に観測し得ない蓄電素子の内部応力をシミュレーションによって推定できる。

　推定装置において、前記シミュレーションモデルは、前記蓄電素子の固有歪みと、前記蓄電素子に対する拘束力とをパラメータに含み、前記歪みに係るデータの入力に応じて、前記蓄電素子の内部応力に係るデータを出力するよう構成されてもよい。この構成によれば、例えば、蓄電素子に対する拘束力、蓄電素子の固有歪み、及び蓄電素子の内部応力に基づき、力のつり合いを考慮することによって、蓄電素子の内部応力を推定できる。

　推定装置において、前記固有歪みは、活物質粒子の孤立化、析出物の成長、及び熱膨張の少なくとも１つに起因して生じる前記蓄電素子の歪みであってもよい。この構成によれば、活物質粒子の孤立化、析出物の成長、及び熱膨張の少なくとも１つに起因した蓄電素子の固有歪みを考慮して、蓄電素子の内部応力を推定できる。

　推定装置において、前記推定部は、非線形フィルタを用いた状態推定器を備えてもよい。この構成によれば、アンサンブルカルマンフィルタ、粒子フィルタ、拡張カルマンフィルタ、無香料カルマンフィルタなどの非線形フィルタを用いるので、固有歪みと内部応力との間に線形性を仮定しない場合であっても、蓄電素子の内部応力を良好に推定できる。

　推定装置において、前記推定部は、前記蓄電素子の内部抵抗を前記内部応力の関数として推定してもよい。この構成によれば、内部応力の値を基に蓄電素子の内部抵抗を推定でき、内部応力を反映した蓄電素子の電気化学現象を推定できる。

　推定装置において、前記蓄電素子は、電解質が固体である全固体電池であってもよい。この構成によれば、全固体電池の性能に決定的な影響を及ぼす内部応力の値を推定できる。

　推定装置において、記蓄電素子は、負極に金属リチウムを用いた電池であってもよい。この構成によれば、析出物が生成されやすい電池について、析出物の成長に起因した内部応力を推定できる。
　ここで、正極材料や電解質の種類は問わない。蓄電素子は、負極に金属リチウムを用いた全固体電池であってもよいし、代替的に、正極に硫黄を用いたリチウム硫黄電池（ＬｉＳ電池）であってもよい。全固体電池や、負極に金属リチウムを用いた電池でなくとも、充放電や劣化にともない電極の体積膨張が生じる蓄電素子において、同様の論理が成立する。

　推定装置は、活物質粒子と固体電解質との接触面積をパラメータに含むシミュレーションモデルを用いて、前記固体電解質を含む全固体電池の電気化学現象をシミュレートする演算部を備える。
　この構成によれば、全固体電池の特徴的なパラメータの１つである接触面積を基に、全固体電池の電気化学現象をシミュレートできる。

　推定装置は、充放電に応じて析出物が生じる蓄電素子に関して、前記析出物の生成量及び析出形態に基づき、前記蓄電素子の固有歪みを算出し、算出した固有歪みに基づき、前記蓄電素子に生じる応力又は歪みの分布を算出する演算部を備える。この構成によれば、蓄電素子の固有歪みを算出し、算出した固有歪みに基づき、応力又は歪みの分布を算出できる。

　推定方法は、蓄電素子に発生する歪みに係るデータを取得し、前記蓄電素子の内部の力学的状態を表すシミュレーションモデルを用いて、前記取得部が取得したデータに基づき、前記蓄電素子の内部応力を推定する処理をコンピュータによって実行する。
　この構成によれば、蓄電素子に発生する歪みのデータから、直接的に観測し得ない蓄電素子の内部応力を推定できる。

　推定方法は、活物質粒子と固体電解質との接触面積をパラメータに含むシミュレーションモデルを用いて、前記固体電解質を含む全固体電池の電気化学現象をシミュレートする処理をコンピュータにより実行する。
　この構成によれば、全固体電池の特徴的なパラメータの１つである接触面積を基に、全固体電池の電気化学現象をシミュレートできる。

　推定方法は、充放電に応じて析出物が生じる蓄電素子に関して、前記析出物の生成量及び析出形態に基づき、前記蓄電素子の固有歪みを算出し、算出した固有歪みに基づき、前記蓄電素子に生じる応力又は歪みの分布を算出する処理をコンピュータにより実行する。この構成によれば、蓄電素子の固有歪みを算出し、算出した固有歪みに基づき、応力又は歪みの分布を算出できる。

　コンピュータプログラムは、コンピュータに、活物質粒子と固体電解質との接触面積をパラメータに含むシミュレーションモデルを用いて、前記固体電解質を含む全固体電池の電気化学現象をシミュレートする処理を実行させる。
　全固体電池は、電解液系リチウムイオン電池と異なり、活物質粒子と固体電解質との接触面積を積極的に大きくしようとしない限り、接触抵抗が大きい。すなわち、全固体電池は、外部から圧力を加え、活物質粒子と固体電解質との接触面積を増大させることによって、電気が流れやすいという特徴を有する。
　上記の構成によれば、全固体電池の特徴的なパラメータの１つである接触面積を基に、全固体電池の電気化学現象をシミュレートできる。

　コンピュータプログラムにおいて、前記シミュレーションモデルは、前記接触面積と前記活物質粒子の有効拡散係数との関係を規定してあり、前記接触面積の値に基づき、前記活物質粒子の有効拡散係数を推定する処理を前記コンピュータに実行させてもよい。この構成によれば、接触面積の値と活物質粒子の有効拡散係数とを紐付けることにより、接触面積の値から全固体電池の電気化学現象をシミュレートできる。

　コンピュータプログラムにおいて、前記シミュレーションモデルは、前記接触面積と前記活物質粒子の有効イオン伝導度との関係を規定してあり、前記接触面積の値に基づき、前記活物質粒子の有効イオン伝導度を推定する処理を前記コンピュータに実行させてもよい。この構成によれば、接触面積の値と活物質粒子の有効イオン伝導度とを紐付けることにより、接触面積の値から全固体電池の電気化学現象をシミュレートできる。

　コンピュータプログラムにおいて、前記接触面積は、前記全固体電池の内部応力の関数であり、前記内部応力の値に基づき、前記全固体電池の電気化学現象をシミュレートする処理を前記コンピュータに実行させてもよい。この構成によれば、接触面積の値と内部応力の値とを紐付けることにより、内部応力の値から全固体電池の電気化学現象をシミュレートできる。

　コンピュータプログラムにおいて、前記全固体電池に発生する歪みを計測する歪みセンサから、前記歪みに係る計測データを取得し、前記全固体電池の内部の力学的状態を表すモデルを用いて、取得した計測データに基づき、前記全固体電池の内部応力を推定する処理を前記コンピュータに実行させてもよい。この構成によれば、センサによって計測できる歪みのデータに基づき、内部応力を推定し、推定した内部応力の値から全固体電池の電気化学現象をシミュレートできる。

　コンピュータプログラムにおいて、前記全固体電池の内部抵抗は、前記内部応力の関数であり、前記内部応力の値に基づき、前記内部抵抗の値を推定する処理を前記コンピュータに実行させてもよい。この構成によれば、内部応力の値と内部抵抗の値とを紐付けることにより、全固体電池の内部抵抗を推定できる。

　コンピュータプログラムは、コンピュータに、充放電に応じて析出物が生じる蓄電素子に関して、前記析出物の生成量及び析出形態に基づき、前記蓄電素子の固有歪みを算出し、算出した固有歪みに基づき、前記蓄電素子に生じる応力又は歪みの分布を算出する処理を実行させる。
　蓄電素子の充放電を繰り返した場合、蓄電素子の内部には析出物が析出する。例えば、高エネルギー密度化が期待されるリチウム金属電池は、負極にリチウム金属を使用した電池であるが、充放電の繰り返しにより負極表面に疎な析出物（デンドライトなど）が析出する。リチウム金属電池に限らず、リチウムイオン電池、全固体電池など様々な電池において、様々な形態で金属などの析出物が析出することが知られている。特に、電池の内部応力や歪みにバラツキがあると、局所的に析出物の成長が加速しやすい。析出物の成長は、電池特性の変化、電池の膨張などを引き起こす要因となり得る。
　析出物の成長を抑制するために、電池内部に均一な応力が生じるように電池の外部から拘束力を加えることが効果的であるが、電池内部の応力又は歪みの分布を推定する手法は現在のところ提案されていない。
　上記の構成によれば、蓄電素子の固有歪みを算出し、算出した固有歪みに基づき、応力又は歪みの分布を算出できる。

　コンピュータプログラムにおいて、前記析出物の生成速度を、生成反応の場に生じる応力の関数として記述してあり、前記関数により算出される前記析出物の生成速度に基づき、前記生成量を算出する処理を前記コンピュータに実行させてもよい。この構成によれば、析出物の生成速度に基づき生成量を算出し、算出した析出物の生成量を基に固有歪みを算出できる。

　コンピュータプログラムにおいて、前記生成量及び応力場に基づき、前記蓄電素子の電気化学現象をシミュレートする処理を前記コンピュータに実行させてもよい。この構成によれば、析出物の生成量及び蓄電素子の応力場に基づき、蓄電素子の内部抵抗変化などを含む電気化学現象をシミュレートできる。

　コンピュータプログラムにおいて、前記生成量及び析出形態に基づき、前記蓄電素子における熱現象をシミュレートする処理を前記コンピュータに実行させてもよい。この構成によれば、析出物の生成量及び析出形態に基づき、蓄電素子の温度挙動を含む熱現象をシミュレートできる。

　コンピュータプログラムにおいて、前記蓄電素子内部のガス発生量を算出し、算出したガス発生量に基づき、前記蓄電素子に生じる応力又は歪みの分布を算出する処理を実行させてもよい。この構成によれば、蓄電素子の内部に発生するガスの量を基に蓄電素子の応力又は歪みの分布を算出できる。

　以下、本発明をその実施の形態を示す図面に基づいて具体的に説明する。
（実施の形態１）
　図１は実施の形態１に係る推定システムの全体構成を示す模式図である。実施の形態１に係る推定システムは、推定装置１及び蓄電素子２を備える。推定装置１は、例えば、ＢＭＵ（Battery Management Unit）などの装置であり、後述する手法にて蓄電素子２の内部応力を推定し、推定した内部応力に係る情報を出力する。図１の例では、便宜的に推定装置１と蓄電素子２とを別体として記載したが、推定装置１と蓄電素子２とは一体となった構成であってもよい。更に、推定装置１は、蓄電素子２を含むバッテリシステムと通信可能に接続されるコンピュータ、サーバ装置などの情報処理装置であってもよい。推定装置１は蓄電素子２に近接している必要はなく、別の建屋のサーバ室に設置されてもよいし、国内あるいは国外の遠隔地に設置されてもよい。さらには、蓄電素子２が大気圏や宇宙空間にあって推定装置１が地球上に存在してもよいし、蓄電素子２及び推定装置１の両方が大気圏や宇宙空間にあってもよい。

　実施の形態１における蓄電素子２は、例えば全固体電池である。放電時には、蓄電素子２は負荷７に接続される。蓄電素子２は、接続された負荷７に対して直流電力を供給する。充電時には、蓄電素子２は充電装置（不図示）に接続される。蓄電素子２は、接続された充電装置から供給される電力によって蓄電する。蓄電素子２は全固体電池に限定されず、電極の膨張がおこる電池であればよい。

　推定システムは、蓄電素子２の状態を計測する各種センサを備える。推定システムが備えるセンサの一例は、歪みセンサＳ１である。歪みセンサＳ１は、蓄電素子２に発生する歪みを時系列的に計測し、計測結果を示すデータを推定装置１へ出力する。

　推定システムは、蓄電素子２の温度を計測する温度センサＳ２を備えてもよい。温度センサＳ２は、蓄電素子２の温度を時系列的に計測し、計測結果を示すデータを推定装置１へ出力する。更に、推定システムは、蓄電素子２の環境温度を推定する温度センサＳ３を備えてもよい。温度センサＳ３は、蓄電素子２が設置されている周囲環境の温度を計測し、計測結果を示すデータを推定装置１へ出力する。

　推定システムは、蓄電素子２に流れる電流を計測する電流計Ｓ４を備えてもよい。電流計Ｓ４は、蓄電素子２に流れる電流を時系列的に計測し、計測結果を示すデータを推定装置１へ出力する。推定装置システムは、更に、蓄電素子２の電圧を計測する電圧計Ｓ５を備えてもよい。電圧計Ｓ５は、蓄電素子２の電圧を時系列的に計測し、計測結果を示すデータを推定装置１へ出力する。

　推定装置１は、各種センサによって計測された計測データを取得し、取得した計測データに基づき、蓄電素子２の内部応力を推定する。

　以下、蓄電素子２の構成について説明する。
　図２は蓄電素子２の構成を説明する説明図である。蓄電素子２は、例えば、正極集電体層２１、正極活物質層２２、固体電解質層２３、負極活物質層２４、及び負極集電体層２５からなる積層体を含む全固体電池である。

　正極集電体層２１は、金属箔、金属メッシュ等により構成される。正極集電体層２１を構成する金属は、アルミニウム、ニッケル、チタン、ステンレス鋼等の良好な導電性を有する金属である。正極集電体層２１の表面には、接触抵抗を調整するためのコート層が形成されていてもよい。コート層の一例は炭素コートである。正極集電体層２１の厚みは特に限定されるものではなく、例えば０．１μｍ以上１ｍｍ以下である。

　正極活物質層２２は、少なくとも正極活物質を含む層である。正極活物質層２２は、正極活物質の他に、固体電解質、導電助剤、及びバインダなどを含んでもよい。正極活物質層２２は、例えば０．１μｍ以上１ｍｍ以下の厚みを有する。

　正極活物質には、全固体電池に使用できる適宜の正極活物質が用いられる。例えば、コバルト酸リチウム、ニッケル酸リチウム、マンガン酸リチウム、スピネル系リチウム化合物等の各種のリチウム含有複合酸化物が、正極活物質として用いられる。正極活物質は、例えば、平均粒径（Ｄ５０）が０．５μｍ以上２０μｍ以下の粒子である。正極活物質を構成する粒子は、一次粒子であってもよく、二次粒子であってもよい。正極活物質は、粒子に限らず、薄膜状であってもよい。正極活物質層２２に含まれる固体電解質には、イオン伝導度が比較的高く、耐熱性に優れた無機固体電解質が用いられる。このような無機固体電解質として、ランタンジルコン酸リチウム等の酸化物固体電解質やＬｉ₂ Ｓ－Ｐ₂ Ｓ₅ 等の硫化物固体電解質が使用できる。導電助剤には、アセチレンブラックやケッチェンブラック等の炭素材料やニッケル、アルミニウム、ステンレス鋼等の金属材料が用いられる。バインダには、ブタジエンゴム（ＢＲ）、アクリレートブタジエンゴム（ＡＢＲ）、ポリフッ化ビニリデン（ＰＶｄＦ）等の材料が用いられる。

　固体電解質層２３は、少なくとも固体電解質を含む層である。固体電解質層２３は、固体電解質の他に、バインダなどを含んでもよい。固体電解質層２３は、例えば０．１μｍ以上１ｍｍ以下の厚みを有する。固体電解質層２３に含まれる固体電解質には、上述した酸化物固体電解質や硫化物固体電解質などの無機固体電解質が用いられる。バインダには、正極活物質層２２に用いられるバインダと同様のものが用いられる。

　負極活物質層２４は、少なくとも負極活物質を含む層である。負極活物質層２４は、負極活物質の他に、固体電解質、導電助剤、及びバインダなどを含んでもよい。負極活物質層２４は、例えば０．１μｍ以上１ｍｍ以下の厚みを有する。

　負極活物質には、全固体電池に使用できる適宜の負極活物質が用いられる。例えば、金属活物質及びカーボン活物質が負極活物質として用いられる。金属活物質として、Ｌｉ、Ｉｎ、Ａｌ、Ｓｉ、Ｓｎ等が挙げられる。金属活物質は、金属単体に限らず、金属複合酸化物であってもよい。カーボン活物質として、例えば、メソカーボンマイクロビーズ（ＭＣＭＢ）、高配向性グラファイト（ＨＯＰＧ）、ハードカーボン、ソフトカーボン等が挙げられる。負極活物質は、例えば、平均粒径（Ｄ５０）が０．５μｍ以上２０μｍ以下の粒子である。負極活物質を構成する粒子は、一次粒子であってもよく、二次粒子であってもよい。負極活物質は、粒子に限らず、薄膜状であってもよい。負極活物質層２４に用いられる固体電解質、導電助剤、及びバインダには、正極活物質層２２に用いられる固体電解質、導電助剤、及びバインダと同様のものが適宜用いられる。

　負極集電体層２５は、金属箔、金属メッシュ等により構成される。負極集電体層２５を構成する金属は、銅、ニッケル、チタン、ステンレス鋼等の良好な導電性を有する金属である。負極集電体層２５の表面には、接触抵抗を調整するためのコート層が形成されていてもよい。コート層の一例は炭素コートである。負極集電体層２５の厚みは特に限定されるものではなく、例えば０．１μｍ以上１ｍｍ以下である。

　蓄電素子２は、拘束部材３によって拘束される。拘束部材３は、例えば、蓄電素子２を収容するケース３１と、ケース３１内に圧縮された状態にて配置される弾性部材３２とを備える。ケース３１は、例えば直方体状の容器であり、底面部３１１と側面部３１２とにより構成されるケース本体３１０と、ケース本体３１０の開口を閉塞する蓋体３１３とを備える。ケース本体３１０（底面部３１１及び側面部３１２）並びに蓋体３１３は、例えば、ステンレス鋼、アルミニウム、アルミニウム合金など溶接可能な金属によって形成される。代替的に、ケース本体３１０（底面部３１１及び側面部３１２）並びに蓋体３１３は、樹脂によって形成されてもよい。ケース本体３１０は、蓄電素子２が収容された後、蓋体３１３によって密閉される。

　弾性部材３２は、蓄電素子２の最下層（図２の例では正極集電体層２１）と底面部３１１との間、及び、蓄電素子２の最上層（図２の例では負極集電体層２５）と蓋体３１３との間に、圧縮された状態にて配置される。弾性部材３２は例えばゴム状のシートである。弾性部材３２は、その弾性力によって、蓄電素子２に積層方向（図の上下方向）の拘束力を付与する。

　図２の例では、ケース３１の内部に弾性部材３２を配置することによって、蓄電素子２に対して拘束力を付与する構成とした。代替的に、ケース３１内に高圧の流体を充填することによって蓄電素子２に拘束力を付与してもよい。この場合、流体としては、電池材料に対して不要な反応を生じさせないものが好ましい。例えば、窒素等の不活性ガスや乾燥空気等が用いられる。代替的に、蓄電素子２を積層方向の両側から板部材によって挟み込み、蓄電素子２に拘束力が付与された状態にて板部材同士を連結することによって、蓄電素子２に拘束力を付与する構成としてもよい。

　蓄電素子２の歪みを計測する歪みセンサＳ１は、蓄電素子２の内部応力に応じて生じる歪みを計測し得る場所に取り付けられる。図２の例では、蓄電素子２の内部応力に応じた歪みはケース３１の側面部３１２に現れるので、歪みを計測する歪みセンサＳ１は、この側面部３１２の適宜箇所に取り付けられることが好ましい。代替的に、歪みセンサＳ１は、ケース３１の底面部３１１又は蓋体３１３に取り付けられてもよい。更に、歪みセンサＳ１は、蓄電素子２に取り付けられてもよい。

　図３は固体電解質層２３の内部構成を説明する説明図である。図３の例では、活物質粒子はハッチングを付した球として示され、固体電解質はハッチングを付していない球として示されている。簡略化のために、図３では、導電助剤及びバインダは省略されている。従来の電解液系リチウムイオン電池では、活物質粒子の周囲が電解液によって満たされており、活物質の全表面が電解液と接触した状態にある。これに対し、固体電解質を用いた全固体電池では、図中の黒丸によって示しているように、固体電解質と活物質粒子とが微小な接触面積（点）で接触した状態にある。固体電解質と活物質粒子との接触面積は、蓄電素子を拘束する拘束力や内部応力によって変化する。

　全固体電池では、拘束力や内部応力に応じて固体電解質及び活物質粒子の接触面積が変化し、これに伴って電池特性は大きく変化する。全固体電池の電池特性（充放電特性など）を正確に推定するためには、内部応力の推定が不可欠である。負極に金属リチウムを用いた電池では、内部応力によって析出物の生成速度が変化するので、内部応力の推定が不可欠である。

　以下、推定装置１の構成について説明する。
　図４は推定装置１の内部構成を示すブロック図である。推定装置１は、例えば、演算部（推定部）１１、記憶部１２、入力部１３、及び出力部１４を備える。

　演算部１１は、ＣＰＵ（Central Processing Unit）、ＲＯＭ（Read Only Memory）、ＲＡＭ（Random Access Memory）などを備える演算回路である。演算部１１が備えるＣＰＵは、ＲＯＭや記憶部１２に格納された各種コンピュータプログラムを実行し、上述したハードウェア各部の動作を制御することによって、装置全体を、蓄電素子２の内部応力を推定するための状態推定器（オブザーバともいう）として機能させる。具体的には、演算部１１は、蓄電素子２の内部の力学的状態をシミュレートするシミュレーションモデルＭＤ１を用いて、入力部１３を通じて入力される歪みの計測データに基づき、蓄電素子２の内部応力を推定するための演算を実行する。代替的に、演算部１１は、ユーザの手により作成される歪みの仮想的なデータを用いて、蓄電素子２の内部応力を推定するための演算を実行してよい。演算部１１は、推定装置１又は外部のコンピュータにより生成される歪みの仮想的なデータを用いて、蓄電素子２の内部応力を推定するための演算を実行してもよい。

　代替的に、演算部１１は、複数のＣＰＵ、マルチコアＣＰＵ、ＧＰＵ（Graphics Processing Unit）、マイコン、揮発性または不揮発性のメモリ等を備える任意の処理回路または演算回路であってもよい。演算部１１は、計測開始指示を与えてから計測終了指示を与えるまでの経過時間を計測するタイマ、数をカウントするカウンタ、日時情報を出力するクロック等の機能を備えてもよい。

　記憶部１２は、フラッシュメモリ、ハードディスクなどの記憶装置を備える。記憶部１２には、各種のコンピュータプログラム及びデータが記憶される。記憶部１２に記憶されるコンピュータプログラムは、シミュレーションモデルＭＤ１を用いて、蓄電素子２の内部応力を推定する処理をコンピュータに実行させるための推定プログラムＰＧ１を含む。シミュレーションモデルＭＤ１は、推定プログラムＰＧ１の中において記述されてもよい。記憶部１２に記憶されるデータには、シミュレーションモデルＭＤ１において用いられるパラメータ、推定プログラムＰＧ１において用いられるパラメータ、演算部１１によって生成されるデータなどが含まれる。

　推定プログラムＰＧ１は、ＭＡＴＬＡＢ（登録商標）、Ａｍｅｓｉｍ（登録商標）、Ｔｗｉｎ　Ｂｕｉｌｄｅｒ（登録商標）、ＭＡＴＬＡＢ＆Ｓｉｍｕｌｉｎｋ（登録商標）、Ｓｉｍｐｌｏｒｅｒ（登録商標）、ＡＮＳＹＳ（登録商標）、Ａｂａｑｕｓ（登録商標）、Ｍｏｄｅｌｉｃａ（登録商標）、ＶＨＤＬ－ＡＭＳ（登録商標）、Ｃ言語、Ｃ＋＋、Ｊａｖａ（登録商標）などの市販の数値解析ソフトウェア又はプログラミング言語によって記述されてもよい。数値解析ソフトウェアは、１Ｄ－ＣＡＥと称される回路シミュレータであってもよく、３Ｄ形状で行う有限要素法や有限体積法などのシミュレータであってもよい。代替的に、これらに基づいた縮退モデル（ＲＯＭ : Reduced-Order Model）を用いてもよい。

　推定プログラムＰＧ１を含むコンピュータプログラムは、当該コンピュータプログラムを読み取り可能に記録した非一時的な記録媒体Ｍにより提供される。記録媒体Ｍは、ＣＤ－ＲＯＭ、ＵＳＢメモリ、ＳＤ（Secure Digital）カードなどの可搬型メモリである。演算部１１は、図に示していない読取装置を用いて、記録媒体Ｍから所望のコンピュータプログラムを読み取り、読み取ったコンピュータプログラムを記憶部１２に記憶させる。代替的に、上記コンピュータプログラムは通信により提供されてもよい。

　入力部１３は、各種センサを接続するためのインタフェースを備える。入力部１３には、蓄電素子２に発生する歪みを計測する歪みセンサＳ１が接続される。演算部１１は、入力部１３を通じて、歪みセンサＳ１によって計測される歪みの計測データを取得する。

　入力部１３には、蓄電素子２の温度を計測する温度センサＳ２、蓄電素子２の環境温度を計測する温度センサＳ３などが接続されてもよい。温度センサＳ２は、蓄電素子２又は蓄電素子２を収容するケース３１の適宜箇所に設けられ、蓄電素子２の温度を計測するセンサである。温度センサＳ３は、蓄電素子２の周囲に設けられ、蓄電素子２の周囲の温度（環境温度）を計測するセンサである。温度センサＳ２，Ｓ３には、熱電対、サーミスタなどの既存のセンサが用いられる。演算部１１は、気象サーバなどの外部サーバから環境温度のデータを取得してもよい。

　入力部１３には、蓄電素子２に流れる電流を計測する電流計Ｓ４、蓄電素子２の電圧を計測する電圧計Ｓ５が接続されてもよい。

　出力部１４は、外部装置を接続する接続インタフェースを備える。出力部１４に接続される外部装置は、液晶ディスプレイなどを備える表示装置１４０である。この場合、演算部１１は、推定した蓄電素子２の内部応力に係る情報を出力部１４から出力することによって、当該情報を表示装置１４０に表示させる。代替的に、推定装置１は、表示装置１４０を備えてもよい。

　更に、出力部１４は、外部装置と通信するための通信インタフェースを備えてもよい。出力部１４に通信可能に接続される外部装置は、蓄電素子２の状態を監視する監視サーバである。代替的に、出力部１４に通信可能に接続される外部装置は、蓄電素子２から供給される電力により動作する携帯端末や電気自動車などの制御装置であってもよい。

　以下、推定装置１が実行する演算処理の内容について説明する。
　推定装置１は、蓄電素子２の内部の力学的状態を表すシミュレーションモデルを用いて、入力部１３を通じて入力される歪みセンサＳ１の計測データに基づき、蓄電素子２の内部応力を推定する。

　蓄電素子２の内部の力学的状態を表すシミュレーションモデルは、力のつり合い式によって表され、例えば、以下の数１が用いられる。

　ここで、Ｆ_ext は拘束部材３によって蓄電素子２に付与される拘束力である。Ｓは蓄電素子２の拘束力に垂直な断面積である。Ｅは蓄電素子２の弾性率である。εは蓄電素子２の弾性歪みである。ε_iso,e は活物質粒子の孤立化による蓄電素子２の固有歪みであり、ε_pre,e は析出物の成長による蓄電素子２の固有歪みである。

　実施の形態１では、蓄電素子２に歪みが発生する要因として、（１）活物質粒子の孤立化、（２）析出物の成長を考慮する。

（１）孤立化による固有歪み
　孤立化とは、充放電の際に活物質粒子に電荷担体（例えばリチウム原子）が挿入離脱することによって、膨張収縮を繰り返し、応力によって活物質粒子が割れる現象を表す。孤立化は、微粉化（Pulverization）とも亀裂進展とも呼ばれる。活物質粒子が割れることによって隙間が生じ、見かけの体積が増すため、蓄電素子２の内部に固有歪みが発生する。

　孤立化による固有歪みの進行速度は、例えば、数２のように記述される。

　ここで、ε_iso,e は孤立化による固有歪みを表す。下付き添え字の「ｉｓｏ」は孤立化（isolation）を表し、下付き添え字の「ｅ」は固有歪みであることを表す。上付き添え字の「ｋ」及び「ｋ＋１」は時間ステップを表す。ｋ_iso,0 ，ｋ_iso,1 は速度係数であり、それぞれ時間と共に孤立化による固有歪みが増加する度合い、通電によって孤立化が進行して固有歪みが増加する度合いを表す。孤立化は、経時変化のみによって進行することは殆どないため、ｋ_iso,0 ＝０．０としても問題はないことが多い。Ｉは蓄電素子２に流れる電流である。α_iso,1 は電流の比例べき乗定数である。ｖ_iso,e は孤立化進行の外乱項を表す。

（２）析出物の成長による固有歪み
　例えば、蓄電素子２の負極にリチウム金属が用いられている場合、長期間にわたる充放電の繰り返しによって、負極表面には析出物が析出することがある。この析出物の成長によって、蓄電素子２の内部に固有歪みが発生する。

　析出物の成長による固有歪みの進行速度は、例えば、数３のように記述される。

　ここで、ε_pre,e は析出物の成長による固有歪みを表す。下付き添え字の「ｐｒｅ」は析出物（Precipitation）を表し、下付き添え字の「ｅ」は固有歪みであることを表す。上付き添え字の「ｋ」及び「ｋ＋１」は時間ステップを表す。ｋ_pre,0 ，ｋ_pre,1 は速度係数であり、それぞれ時間と共に析出物が成長して固有歪みが増加する度合い、通電によって析出物が成長して固有歪みが増加する度合いを表す。σ_inは蓄電素子２の内部応力である。速度係数ｋ_pre,0 ，ｋ_pre,1 は共に内部応力σ_inの関数であることを表す。σ_inは蓄電素子内部の場所の関数であってもよい。Ｉは蓄電素子２に流れる電流である。α_pre,1 は電流の比例べき乗定数である。ｖ_pre,e は析出物の成長による固有歪みの外乱項を表す。

　析出物は疎なリチウム金属や不動体被膜（ＳＥＩ被膜）などであってもよい。どのような析出物が蓄電素子２の内部に析出する場合であっても、析出物の成長に伴う固有歪みの進行速度は、数３と同様の数式によって表される。

　蓄電素子２に付与される拘束力Ｆ_ext は、拘束部材３を構成する側面部３１２の引張力と一致するので、以下の数４によって表される。

　ここで、Ｅ_ref は側面部３１２のヤング率、ε_ref は歪みセンサＳ１によって計測される側面部３１２の歪み、Ｓ_ref は側面部３１２の断面積である。

　蓄電素子２の内部応力σ_inと弾性歪みεとの間には、σ_in＝εＥの関係があるので、数１及び数４の関係式を用いると、内部応力σ_inは、次のように表される。

　数５において、上付きの添え字ｋは時間ステップを表す。ｖ_inは内部応力の外乱項である。数５では、歪みの要因として、孤立化による固有歪みと、析出物の成長による固有歪みとが考慮されている。他の要因が更に存在する場合であっても、固有歪みの和の項に追加することによって、同様の手続きが適用できる。

　２つの固有歪み、歪みセンサＳ１の計測値、及び内部応力の方程式を、状態方程式として表した場合、以下の数６が得られる。

　数６は、孤立化による固有歪みの進行速度を表す式（数２）、析出物の成長による固有歪みの進行速度を表す式（数３）、内部応力に係る方程式（数５）に加え、側面部３１２の歪みの方程式を含む。数６においてε_ref は観測量である。電流Ｉについても電流計Ｓ４によって計測された値が用いられてもよい。数６では、固有歪みを状態量として用いたが、固有歪みにヤング率を乗じて得られる応力を状態量に用いてもよい。

　数６の状態方程式は、数７のようにベクトルを用いた表現に書き換えられる。

　ここで、ｘ^k は状態量を要素に持つベクトル（状態ベクトル）、ｖ^k は外乱量を要素に持つベクトル（外乱ベクトル）である。ｆは数６に示す状態方程式の非線形変換を表す。外乱項は一部又は全部の要素を０にして計算してもよい。

　実施の形態１では、蓄電素子２の歪みは歪みセンサＳ１によって計測されるので、ε_ref が観測量である。観測方程式は次式によって表される。

　ここで、ｙ^k は観測値、Ｃ^T は観測ベクトルである。観測方程式についても外乱ベクトルを加えてもよい。観測量として第３成分のε_ref を取り出す場合、観測ベクトルＣ^T は数９にように表される。

　実施の形態１に係る推定装置１は、非線形フィルタを利用して、数７の状態方程式及び数８の観測方程式によって表されるシミュレーションモデルの時間更新を逐次計算し、内部応力σ_inの時間推移を導出する。

　以下、非線形フィルタの一例として、アンサンブルカルマンフィルタを用いて、時間更新を逐次計算する手法について説明する。

　図５は実施の形態１における内部応力の推定手順を説明するフローチャートである。推定装置１の演算部１１は、ｋ＝１の初期値を与える（ステップＳ１０１）。演算部１１は、歪みセンサＳ１を用いて予め計測した歪みの計測値をε_ref ^kの初期値として与え、孤立化による固有歪みε_iso,e ^k、析出物の成長による固有歪みε_pre,e ^k、内部応力σ_in ^kの初期値として予め設定した仮の値を与えればよい。

　次いで、演算部１１は、各状態変数について、Ｎ個の粒子を発生させる（ステップＳ１０２）。ここで、Ｎは、１０²～１０⁶ 程度の数である。

　次いで、演算部１１は、ｉ＝１，２，…，Ｎについて、ｖ_k に相当する乱数を発生させる（ステップＳ１０３）。ｖ_k は正規分布に従うものとし、分散は既知とする。

　演算部１１は、全てのＮ個の粒子について、数１０に基づく演算を実行し、粒子の状態を次の時間ステップにおける粒子の状態に更新する（ステップＳ１０４）。

　演算部１１は、ｉ＝１，２，…，Ｎの各粒子の状態ベクトルと、全粒子の状態ベクトルの平均値との差ｘ_k ⁽ⁱ⁾_barを算出する（ステップＳ１０５）。ｘ_k ⁽ⁱ⁾_barは、数１１によって表される。

　演算部１１は、全ての粒子に係る状態量予測値の共分散行列Ｐ_k を算出する（ステップＳ１０６）。共分散行列Ｐ_k は、数１２によって表される。

　演算部１１は、入力部１３を通じて歪みセンサＳ１のセンサ出力を取得する（ステップＳ１０７）。取得した歪みセンサＳ１のセンサ出力は、時間ステップｋにおける各粒子の観測値ｙ_k ⁱを与える。

　演算部１１は、ｉ番目の粒子の時間ステップｋにおける観測誤差ｒ_k ⁱを算出する（ステップＳ１０８）。ここで、ｗ_k は観測外乱である。観測誤差ｒ_k ⁱは、数１３によって表される。

　演算部１１は、時間ステップｋにおけるカルマンゲインＫ_k を算出する（ステップＳ１０９）。カルマンゲインＫ_k は、数１４によって表される。

　演算部１１は、ｉ番目の粒子の推定値ｘ_k ⁽ⁱ⁾_hatを算出する（ステップＳ１１０）。推定値ｘ_k ⁽ⁱ⁾_hatは、数１５によって表される。すなわち、演算部１１は、数１０の最初の予測値を、数１３の観測誤差ｒ_k ⁱと数１４のカルマンゲインＫ_k とを用いて修正する。

　演算部１１は、各粒子の平均値ｘ_k _hatを算出する（ステップＳ１１１）。各粒子の平均値ｘ_k _hatは、アンサンブルカルマンフィルタによって得られる状態ベクトル推定値を表し、次式によって算出される。

　数１６によって得られる推定値（各粒子の平均値ｘ_k _hat）には、内部応力σ_inの推定値が含まれる。

　次いで、演算部１１は、演算を終了するか否かを判断する（ステップＳ１１２）。例えば、ユーザから終了指示が与えられた場合、演算部１１は、演算を終了すると判断する。演算を終了しないと判断した場合（Ｓ１１２：ＮＯ）、演算部１１は、処理をステップＳ１０２へ戻し、次の時間ステップにおける演算を実行する。

　演算を終了すると判断した場合（Ｓ１１２：ＹＥＳ）、演算部１１は、推定した内部応力σ_inに係る情報を出力部１４から出力し（ステップＳ１１３）、本フローチャートによる処理を終了する。演算部１１が出力する内部応力σ_inに係る情報は、内部応力の値そのものであってもよく、内部応力を基に導出される物理量（例えば、蓄電素子２の内部抵抗）であってもよい。更に、演算部１１が出力する内部応力σ_inに係る情報は、内部応力σ_inの時間推移を示すグラフであってもよく、応力分布を示す２次元又は３次元のグラフやコンターマップであってもよい。

　以上のように、推定装置１は、アンサンブルカルマンフィルタを利用して蓄電素子２の内部応力σ_inを推定する。アンサンブルカルマンフィルタは、非線形性や非ガウス性を有する状態空間モデルを対象としたフィルタ手法であり、より一般的な状態空間モデルを対象とすることができる。アンサンブルカルマンフィルタは、アルゴリズムが比較的単純であり、推定装置１に容易に実装できる。

　図５のフローチャートは、アンサンブルカルマンフィルタによる演算手法を説明している。代替的に、推定装置１は、粒子フィルタ、拡張カルマンフィルタ、無香料カルマンフィルタなどの非線形フィルタを用いて、蓄電素子２の内部応力σ_inを推定してもよい。

　実施の形態１では、数５を導出する際、固有歪みと内部応力との間の線形性を考慮したが、両者の関係は非線形であってもよい。固有歪みと内部応力との間の関係が非線形であっても、推定装置１は、非線形フィルタを用いた演算を実行することによって、蓄電素子２の内部応力σ_inを推定できる。

（実施の形態２）
　実施の形態２では、温度に起因した固有歪みを更に考慮して内部応力を推定する手法について説明する。
　推定装置１及び蓄電素子２の構成は実施の形態１と同様であるため、その説明を省略する。

　実施の形態２では、蓄電素子２に歪みが発生する要因として、（１）活物質粒子の孤立化、（２）析出物の成長、（３）熱膨張が考慮される。活物質粒子の孤立化による固有歪み、及び、析出物の成長による固有歪みは、実施の形態１と同様であるため、説明を省略する。

（３）熱膨張による固有歪み
　熱膨張は、温度が上昇することによって蓄電素子２の体積が増大する現象である。熱膨張は蓄電素子２の劣化と無関係であり、ある瞬間の温度のみによって定まる。実施の形態２では、熱膨張が温度に比例することを仮定し、基準温度との差異に応じた固有歪みが発生するとしたモデルが記述される。

　温度による固有歪みは、例えば、数１７のように記述される。

　ここで、ε_th,0は熱膨張による固有歪みを表す。α_thは線熱膨張係数である。Ｔはある時刻における温度、Ｔ_ref は基準温度である。ｖ_th,eは熱膨張の外乱項である。上付き添え字のｋは時間ステップを表す。下付きの添え字のｔｈは温度（thermal）を表す。温度Ｔには温度センサＳ２の計測データ、基準温度Ｔ_ref には温度センサＳ３の計測データが用いられる。

　温度を含む状態方程式は、以下の数１８のように記述される。

　実施の形態２では、孤立化による固有歪みの進行速度を表す速度係数ｋ_iso,0 ，ｋ_iso,1 を温度Ｔの関数としている。関数形としては、温度Ｔの単調増加関数が用いられる。単調増加関数として、例えば、高温になる程、孤立化の進行速度が速くなるような特徴を表すアレニウス型の関数が用いられる。実施の形態２では、析出物の成長による固有歪みの進行速度を表す速度係数ｋ_pre,0 ，ｋ_pre,1 は温度Ｔ及び内部応力σ_inの関数としている。関数形は、温度Ｔの単調減少関数が用いられる。単調減少関数として、例えば、低温になる程、析出物の成長の進行が速くなるような特徴を表す関数が用いられる。数１８の第５式には、熱膨張による固有歪みε_th,e ^k に関する項が追加される。

　数１８の状態方程式は、数１９のようにベクトルを用いた表現に書き換えられる。

　観測方程式は、実施の形態１において説明した数８と同様であるが、実施の形態２では、蓄電素子２の歪み及び温度を観測量として取り出すので、観測ベクトルＣ^T は数２０のように表される。

　実施の形態２に係る推定装置１は、非線形フィルタを利用して、数１９の状態方程式及び数８の観測方程式によって表されるシミュレーションモデル（時系列モデル）の時間更新を逐次計算し、温度を考慮した内部応力σ_inの時間推移を導出する。演算手法は実施の形態１と同様であり、推定装置１は、図５に示すフローチャートの手順にて演算を実行することによって、蓄電素子２の内部応力σ_inを推定する。

　実施の形態１及び２に係る推定装置１は、推定プログラムＰＧ１を実行することによって、蓄電素子２の内部応力を推定する構成とした。代替的に、推定装置１は、蓄電素子２の劣化を推定する推定プログラムを併用し、蓄電素子２の使用に伴って電気化学特定が劣化していく様子と、蓄電素子２が膨張していく様子とを同時的にシミュレートしてもよい。蓄電素子２の劣化を推定する推定プログラムとして、例えば、特願２０２０－４８３６９号に記載の手法が使用される。

　実施の形態２では、熱膨張に起因した蓄電素子２の固有歪みが考慮された。代替的に、活物質粒子の挿入脱離に伴う膨張収縮に起因した固有歪みが考慮されてもよい。このような膨張収縮は、蓄電素子２の劣化とは関係なく、通常の充放電時に正極及び負極に活物質粒子が挿入脱離することによって生じる。活物質粒子の挿入脱離による固有歪みは可逆的なものであり、例えばＳＯＣ（State Of Charge）の関数として表される。

　実施の形態２では、熱応力を介した温度の影響を考慮するモデルが記述された。代替的に、熱膨張による固有歪みを考慮せず、孤立化や析出物成長の速度係数の温度依存性のみが考慮されてもよい。その場合、数１８及び数１９の第３式が除外されるとよい。

（実施の形態３）
　実施の形態３では、蓄電素子２の電気化学現象を推定する際、推定装置１によって推定される内部応力σ_inの値を用いる構成について説明する。
　推定装置１及び蓄電素子２の構成は実施の形態１と同様であるため、その説明を省略する。

　蓄電素子２の電気化学現象は、ＮｅｗｍａｎモデルやＲａｎｄｌｅモデルなどの物理モデルによって記述される。観測方程式は、例えば、数２１のように記述される。

　ここで、Ｖは蓄電素子２の端子電圧であり、電圧計Ｓ５により観測される観測値である。ＯＣＰ_p （ｃ_p,1 ）は正極の平衡電位であり、正極活物質粒子の界面での吸蔵リチウムイオン濃度ｃ_p,1 の関数である。ＯＣＰ_n （ｃ_n,1 ）は負極の平衡電位であり、負極活物質粒子の界面での吸蔵リチウムイオン濃度ｃ_n,1 の関数である。Ｒ_ohm （σ_in）は蓄電素子２のオーム抵抗（内部抵抗）を表す。Ｒ_ohm （σ_in）は、オーム抵抗が内部応力σ_inの関数であることを表す。内部応力σ_inの値として推定装置１により推定される値が用いられる。Ｒ_ohm （σ_in）は温度Ｔの関数であってもよい。Ｉは蓄電素子２に流れる電流である。すなわち、Ｒ_ohm （σ_in）Ｉの項は、オーム抵抗による電圧降下を表す。η_act,p （ｃ_p,1 ，Ｉ）は正極活物質粒子の界面での活性化過電圧であり、正極活物質粒子の界面での吸蔵リチウムイオン濃度ｃ_p,1 、電流Ｉ、及び温度Ｔの非線形関数である。η_act,n （ｃ_n,1 ，Ｉ）は負極活物質粒子の界面での活性化過電圧であり、負極活物質粒子の界面での吸蔵リチウムイオン濃度ｃ_n,1 、電流Ｉ、及び温度Ｔの非線形関数である。すなわち、観測値の電圧Ｖは、正極活物質粒子の界面での吸蔵リチウムイオン濃度ｃ_p,1 、負極活物質粒子の界面での吸蔵リチウムイオン濃度ｃ_n,1 、電流Ｉ、温度Ｔの複雑な非線形関数となる。

　図６は蓄電素子２の内部応力σ_inとオーム抵抗Ｒ_ohm との関係を示すグラフである。グラフの横軸は蓄電素子２の内部応力σ_inを表し、縦軸は蓄電素子２のオーム抵抗Ｒ_ohm を表す。図６のグラフに示すように、圧縮応力が高くなる程、オーム抵抗は低くなるという実験事実を考慮して、Ｒ_ohm の関数形は、∂Ｒ_ohm ／∂σ_in ≧０となるように定められる。推定装置１の記憶部１２には、内部応力σ_inをオーム抵抗Ｒ_ohm に変換する関数が記憶されてもよく、内部応力σ_inをオーム抵抗Ｒ_ohm に変換する変換テーブルが記憶されてもよい。

　推定装置１は、上述した手法を用いて推定した内部応力σ_inの値を、予め定めた関数（又はテーブル）に従ってオーム抵抗Ｒ_ohm の値に変換する。推定装置１は、変換後に得られるオーム抵抗Ｒ_ohm の値を用い、数２１に基づき状態推定を実行することによって、正極及び負極の平衡電位や活性化過電圧を含む物理量を推定する。推定手法には、例えば、特願２０２０－１６０９７１号に記載の手法が用いられる。

　全固体電池では、拘束力や内部応力に応じて固体電解質及び活物質粒子の接触面積が変化し、これに伴って電池特性は大きく変化する。実施の形態３では、内部応力σ_inの推定結果を用いて電気化学現象が推定されるので、内部応力に応じて電池特性が大きく変わり得る全固体電池の電気化学現象が精度良く推定される。

　負極に金属リチウムを用いた電池では、内部応力に応じて、オーム抵抗などの内部抵抗や、析出物の成長速度が変化することが知られている。

（実施の形態４）
　実施の形態４では、蓄電素子２の等価回路モデルを用いて電気化学現象を推定する構成について説明する。
　推定装置１及び蓄電素子２の構成は実施の形態１と同様であるため、その説明を省略する。

　図７は等価回路モデルの一例を示す回路図である。蓄電素子２の等価回路モデルは、例えば図７に示すような抵抗器、容量成分、電圧源の組み合わせによって表現されることが多い。

　図７中のＲ₀ はオーム抵抗成分、Ｒ₁ は正極の反応抵抗成分、Ｃ₁ は正極の容量成分、Ｒ₂ は負極の反応抵抗成分、Ｃ₂ は負極の容量成分、Ｅ_eqは開回路電圧（OCV : Open Circuit Voltage）である。ただし、図７は例示であって、直列、並列の組み合わせや電気回路素子の個数や種類に制限はない。

　蓄電素子２の充放電特性は、温度やＳＯＣの影響を受けることが知られている。開回路電圧（ＯＣＶ）はＳＯＣの関数であるとし、Ｒ₀ ～Ｒ₂ 、Ｃ₁ ～Ｃ₂ は温度の関数であるとする。このとき、観測方程式は数２２によって表される。

　ここで、ｙ_U は観測値であり、実施の形態４では蓄電素子２の端子電圧Ｖを表す。上付きの添え字ｋは時間ステップを表す。ＯＣＶ（ＳＯＣ）は開回路電圧であり、ＳＯＣの非線形関数として表される。Ｃ^T は観測ベクトル、ｘ_U は状態ベクトルを表す。Ｒ₀（σ_in）はオーム抵抗であり、内部応力σ_inの関数として表される。Ｒ₀ の関数形は、∂Ｒ₀ ／∂σ_in ≧０となるように定められる。推定装置１の記憶部１２には、内部応力σ_inをオーム抵抗Ｒ₀ に変換する関数が記憶されてもよく、内部応力σ_inをオーム抵抗Ｒ₀ に変換する変換テーブルが記憶されてもよい。ｕは蓄電素子２に流れる電流である。

　推定装置１は、上述した手法を用いて推定した内部応力σ_inの値を、予め定めた関数（又はテーブル）に従ってオーム抵抗Ｒ₀ の値に変換する。推定装置１は、変換後に得られるオーム抵抗Ｒ₀ の値を用い、数２２に基づき状態推定を実行することによって、開回路電圧ＯＣＶを含む物理量を推定する。推定手法には、例えば、特願２０２０－１６０９７１号に記載の手法が用いられる。

　全固体電池では、拘束力や内部応力に応じて固体電解質及び活物質粒子の接触面積が変化し、これに伴って電池特性は大きく変化する。実施の形態４では、内部応力σ_inの推定結果を用いて電気化学現象が推定されるので、内部応力に応じて電池特性が大きく変わり得る全固体電池の電気化学現象が精度良く推定される。

（実施の形態５）
　図８は実施の形態５に係る推定装置１の内部構成を示すブロック図である。推定装置１は、上述した演算部（推定部）１１、記憶部１２、出力部１４に加え、操作部１５を備える。操作部１５は、各種の操作装置を接続するためのインタフェースを備える。操作装置は、キーボード、マウス、タッチパネルなど、ユーザの操作を受付けるための装置である。操作部１５は、操作装置を通じて受付けた操作情報を演算部１１へ出力する。演算部１１、記憶部１２、出力部１４の構成は実施の形態１と同様であるから、その説明は省略される。

　実施の形態５に係る推定装置１は、予め設定された情報又は操作部１５を通じて入力された情報に基づき、蓄電素子２の電気化学現象を推定する。実施の形態５におけるシミュレーション対象の蓄電素子２は、全固体電池である。

　従来の電解液系リチウムイオン電池では、活物質粒子の周囲が電解液によって満たされており、活物質粒子の全表面が電解液と接触した状態にある。この場合、見かけの拡散係数（有効拡散係数）やイオン伝導度（有効イオン伝導度）は、真の拡散係数やイオン伝導度と一致する。

　これに対し、全固体電池では、活物質粒子と固体電解質とが接触する微小面積でしかイオンのやり取りがないため、見かけの拡散係数（有効拡散係数）やイオン伝導度（有効イオン伝導度）は、真の拡散係数やイオン伝導度よりも小さくなると予測される。しかしながら、活物質粒子と固体電解質との間の接触面積と、有効拡散係数や有効イオン伝導度との関係についての実証例は、現時点では存在していない。

　本願発明者らは、活物質粒子と固体電解質との間の接触面積と、有効拡散係数や有効イオン伝導度との関係を数値計算（シミュレーション）により見出した。実施の形態５では、これらの関係を用いて、活物質粒子と固体電解質との間の接触面積（若しくは蓄電素子２の内部応力）を基に、蓄電素子２の電気化学現象を推定するシミュレーション手法を提案する。

　以下では、参考例として、電解液系リチウムイオン電池に関し、電解液の有効拡散係数を求めるための計算手法が説明される。

　図９は有効拡散係数の計算手法を説明する説明図である。図９は電解液系リチウムイオン電池の正極を簡略化したモデルを表す。図９は正極の内部に複数の活物質粒子が存在し、その内部及び周囲が電解液によって満たされている状態を示している。正極の厚みをＬ（ｍ）とし、図中の破線により示される伝導経路の長さをｌ（ｍ）としたとき、間隙・活物質粒子による屈曲度τは、τ＝ｌ／Ｌと記述される。このとき、有効拡散係数Ｄ_l,eff （ｍ² ｓ^-1）と真の拡散係数Ｄ_l （ｍ² ｓ^-1）との関係は、次式により表される。εは電解液の体積占有率を表す。

　同様に、電解液の有効イオン伝導度σ_l,eff （Ｓｍ^-1）と真のイオン伝導度σ_l （Ｓｍ^-1）との関係は、次式により表される。

　図９に示すモデルの上辺を電解液の流入口、下辺を流出口（観測点）とし、適宜の境界条件を与えることにより、観測点での電解液の流速が算出される。境界条件として、例えば、流入口における電解液の濃度（例えば１０００ｍｏｌｍ^-3 ）、及び流出口における電解液の濃度（例えば０ｍｏｌｍ^-3 ）が与えられる。

　観測点における電解液の流速をＪ_l 、電解液の濃度ｃ_l としたとき、流速Ｊ_l と濃度ｃ_l との関係は、次式により表される。Ｄ_l,eff は有効拡散係数を表す。電解液の有効拡散係数Ｄ_l,eff は数２５に基づいて算出される。

　全固体電池では、活物質粒子の一部が固体電解質と接触した状態にあり、その内部及び周囲は電解質によって満たされていない。そのため、電解液系リチウムイオン電池のおける有効拡散係数の計算手法は、そのまま全固体電池に適用することはできない。そこで、本願発明者らは、活物質粒子と固体電解質との間の接触面積という概念を導入し、接触面積に応じて変化する有効拡散係数をシミュレーションにより算出する手法を提案する。

　図１０は全固体電池における有効拡散係数の計算手法を説明する説明図である。実施の形態５において、全固体電池の活物質粒子は球体として仮定される。推定装置１の演算部１１は、この球体において、球体表面を流入口、球体内部に設定した同心球の表面を観測点として設定し、境界状態を与えた場合における、観測点での流速Ｊ_AM（ｍｏｌｍ^-2 ｓ^-1）を算出する。観測点は、例えば、球体の半分の半径を有する同心球の表面として設定される。境界条件として、球体の周囲に例えば１０００ｍｏｌｍ^-3の濃度、球体の内部に例えば０ｍｏｌｍ^-3の濃度が与えられる。流速Ｊ_AMと有効拡散係数Ｄ_AM,effとの関係は、数２６によって表される。ここで、c_AMは活物質の濃度である。演算部１１は、算出した流速Ｊ_AMを数２６に代入し、活物質の有効拡散係数Ｄ_AM,effを算出する。

　図１１Ａ－１１Ｃは全固体電池における有効拡散係数の計算結果を示す図である。図１１Ａは球体として仮定した全固体電池の活物質を表現するために、回転軸を中心に回転させて作成した球体であり、球体表面の全てが電解質と接触している場合の計算結果を示している。球体表面の全てが電解質に接触している場合、観測点に均一な流れ込みが発生するので、活物質の拡散係数とその有効拡散係数とは同一になると予測される。実際に計算を行ったところ、有効拡散係数Ｄ_AM,effは、１．０×１０^-13 （ｍ² ｓ^-1）となり、活物質の拡散係数Ｄ_AMと同一の値となった。これにより、計算手法の妥当性が示された。

　図１１Ｂは球の中心を通る対称軸に対して半頂角１０度で回転させて作成した円錐で切り取られた球冠のうち、球表面に含まれる表面積が電解質と接触している場合の計算結果を示している。半頂角１０度の場合、有効拡散係数Ｄ_AM,effは、３．８×１０^-15 （ｍ² ｓ^-1）となり、活物質の拡散係数Ｄ_AMよりも２桁程度小さな値となった。

　図１１Ｃは半頂角５度の場合の計算結果を示している。半頂角５度の場合、有効拡散係数Ｄ_AM,effは、１．９×１０^-15 （ｍ² ｓ^-1）となり、半頂角１０度の場合と比較して更に小さな値となった。

　図１２は接触面積比率と有効拡散係数との関係を示すグラフである。図１２に示すグラフの横軸及び縦軸は共に対数軸であり、横軸は接触面積比率、縦軸は有効拡散係数（ｍ² ／ｓ）を表す。接触面積比率は、球の表面積に対する球の表面に含まれる球冠の面積の比率であり、活物質と固体電解質との接触の度合いを表すパラメータである。横軸の接触面積比率に活物質を仮定した球の表面積を乗じることによって、接触面積が得られる。図１２に示す対数軸プロットにより、接触面積と有効拡散係数との間の関係は、数２７により定式化されることが分かる。Ｘ線ＣＴ像などを解析して実施に接触面積を見積もることができるほか、接触面積に影響を及ぼすパラメータも使うことができる。例えば、プレス成型して得られる電極の場合には、電極合材の残存空隙率などの値を用いてもよい。

　ここで、Ｄ_AM,effは有効拡散係数、ｘは接触面積を表す。ａ，ｂ，ｃは係数である。係数ａ，ｂ，ｃは、グラフ上の各点を通る近似曲線（図１２に示す対数プロットでは直線）を求めることによって算出される。近似曲線は最小二乗法など公知の手法を用いて求められる。

　演算部１１は、算出した有効拡散係数Ｄ_AM,effの値に基づき、有効イオン伝導度σ_AM,effを算出してもよい。Einsteinの関係式により、イオン伝導度σと拡散係数Ｄとの関係は、数２８により記述される。

　ここで、σはイオン伝導度（Ｓｍ^-1）、ｚはイオンの電荷（無次元）、Ｆはファラデー定数（Ｃｍｏｌ^-1）、Ｄは拡散係数（ｍ² ｓ^-1）、ｃはリチウム濃度（ｍｏｌｍ^-3）、Ｒは気体定数（ｍ² ｋｇｓ^-1 Ｋ^-1 ｍｏｌ^-1）、Ｔは温度（Ｋ）を表す。

　演算部１１は、数２８の拡散係数Ｄに、接触面積から算出した有効拡散係数Ｄ_AM,effの値を代入することにより、有効イオン伝導度σ_AM,effを算出できる。有効イオン伝導度σ_AM,effは、蓄電素子２の電気抵抗に影響を与える物理量である。すなわち、推定装置１は、全固体電池における活物質粒子と固体電解質との接触面積を基に、蓄電素子２における電気化学現象を推定できる。

　以下、推定装置１が実行する演算の実行手順について説明する。
　図１３は実施の形態５に係る推定装置１が実行する演算処理の実行手順を説明するフローチャートである。推定装置１の演算部１１は、全固体電池のモデルとして、球体の活物質粒子と固体電解質とが接触したモデルを設定し、境界条件を与える（ステップＳ５０１）。全固体電池のモデルや境界条件は事前に設定され、記憶部１２にされていてもよい。この場合、演算部１１は、事前に設定されているモデル及び境界条件を記憶部１２から読み出せばよい。代替的に、演算部１１は、操作部１５を通じて、モデル及び境界条件の設定を受け付けてもよい。

　次いで、演算部１１は、接触面積（半頂角）の設定を受付ける（ステップＳ５０２）。接触面積は、蓄電素子２の内部応力を想定して設定されるとよい。蓄電素子２のＳＥＭ画像（ＳＥＭ：Scanning Electron Microscope）を解析することによって得られる現実の値を与えてもよい。シミュレーションに用いる接触面積の値は記憶部１２に事前に記憶されていてもよく、計算時に操作部１５を通じて与えられてもよい。

　次いで、演算部１１は、ステップＳ５０２で設定された接触面積を有する球体の表面部分から観測点に流れ込む流速を算出することによって、有効拡散係数を算出する（ステップＳ５０３）。流速と有効拡散係数との関係は数２６によって表され、有効拡散係数は濃度勾配の係数として算出される。

　次いで、演算部１１は、ステップＳ５０３で算出した有効拡散係数を用いて、有効イオン伝導度を算出する（ステップＳ５０４）。有効拡散係数と有効イオン伝導度との関係は数２８によって表される。

　図１３に示すフローチャートでは、演算部１１は、有効拡散係数及び有効イオン伝導度の双方を算出する構成としたが、何れか一方のみを算出する構成としてもよい。

　演算部１１は、算出した有効拡散係数又は有効イオン伝導度から、蓄電素子２の電気化学現象に関する他の物理量を推定してもよい。例えば、一般的な伝導度σ（Ｓｍ^-1）と内部抵抗Ｒ_ohm （Ω）との間には、Ｌを長さ（ｍ）、Ａを断面積（ｍ² ）として、Ｒ_ohm ＝Ｌ／（σ×Ａ）の関係があるので、演算部１１は、この関係式を用いて、蓄電素子２の内部抵抗を推定してもよい。

　次いで、演算部１１は、算出した有効拡散係数及び有効イオン伝導度の情報を出力部１４より出力する（ステップＳ５０５）。演算部１１は、算出した有効拡散係数及び有効イオン伝導度の値を出力してもよく、数値範囲を出力してもよい。代替的に、演算部１１は、接触面積に対して有効拡散係数及び有効イオン伝導度をプロットしたグラフを出力してもよい。

　以上のように、実施の形態５に係る推定装置１は、活物質粒子と固体電解質との間の接触面積を考慮しながら、有効拡散係数や有効イオン伝導度などの、蓄電素子２の電気化学に関する物理量を推定できる。

（実施の形態６）
　実施の形態６では、蓄電素子２の内部応力を基に、有効拡散係数や有効イオン伝導度を算出する構成について説明する。

　実施の形態５は、活物質粒子と固体電解質との間の接触面積と、有効拡散係数や有効イオン伝導度との関係を説明している。例えば、数２６は、接触面積ｘと有効拡散係数Ｄ_AM,effとの関係を表しているが、接触面積に代えて圧縮応力を用いてもよい。応力は、荷重に対して部材内部に生じる抵抗力であり、荷重を部材面積で除算した単位面積あたりの力である。球状の弾性体同士が接触する部分にかかる応力に関しては、ヘルツ理論が適用される。ヘルツ理論に依れば、接触する部分の圧縮応力は、数２９により記述される。

　ここで、Ｓ_p は圧縮応力（Ｐａ）、Ｆは荷重（Ｎ）、ｒは接触面の半径（ｍ）である。数２９を用いて数２７を書き換えた場合、有効拡散係数Ｄ_AM,effと圧縮応力Ｓ_p との関係が得られる。推定装置１の演算部１１は、接触面積ｘに代えて圧縮応力Ｓ_p を与えることにより、活物質粒子の有効拡散係数Ｄ_AM,effを算出してもよい。

　更に、演算部１１は、算出した有効拡散係数Ｄ_AM,effを数２８に代入し、有効イオン伝導度σ_AM,effを算出してもよい。

　図１４は実施の形態６に係る推定装置１が実行する演算処理の実行手順を説明するフローチャートである。推定装置１の演算部１１は、全固体電池のモデルとして、球体の活物質粒子と固体電解質とが接触したモデルを設定し、境界条件を与える（ステップＳ６０１）。全固体電池のモデルや境界条件は事前に設定され、記憶部１２にされていてもよい。この場合、演算部１１は、事前に設定されているモデル及び境界条件を記憶部１２から読み出せばよい。代替的に、演算部１１は、操作部１５を通じて、モデル及び境界条件の設定を受け付けてもよい。

　次いで、演算部１１は、内部応力の設定を受付ける（ステップＳ６０２）。シミュレーションに用いる内部応力の値は記憶部１２に事前に記憶されていてもよく、計算時に操作部１５を通じて与えられてもよい。

　次いで、演算部１１は、活物質粒子と固体電解質との接触部分から観測点に流れ込む流速を算出することによって、有効拡散係数を算出する（ステップＳ６０３）。流速と有効拡散係数との関係は数２６によって表される。演算部１１は、数２６に基づき、濃度勾配の係数として有効拡散係数を算出できる。

　次いで、演算部１１は、ステップＳ６０３で算出した有効拡散係数を用いて、有効イオン伝導度を算出する（ステップＳ６０４）。有効拡散係数と有効イオン伝導度との関係は数２８によって表される。演算部１１は、数２８に基づき、有効イオン伝導度を算出できる。

　次いで、演算部１１は、算出した有効拡散係数及び有効イオン伝導度の情報を出力部１４より出力する（ステップＳ６０５）。演算部１１は、算出した有効拡散係数及び有効イオン伝導度の値を出力してもよく、数値範囲を出力してもよい。代替的に、演算部１１は、接触面積に対して有効拡散係数及び有効イオン伝導度をプロットしたグラフを出力してもよい。

　図１４のフローチャートにおいて、演算部１１は、有効拡散係数及び有効イオン伝導度の双方を算出する構成としたが、有効拡散係数及び有効イオン伝導度の何れか一方のみを算出する構成としてもよい。

　代替的に、演算部１１は、実施の形態３と同様の手法を用いて、蓄電素子２の内部応力を基に内部抵抗を推定してもよい。

　以上のように、実施の形態６に係る推定装置１は、蓄電素子２の内部応力を考慮しながら、有効拡散係数や有効イオン伝導度などの、蓄電素子２の電気化学に関する物理量を推定できる。

（実施の形態７）
　以下、実施の形態７のシミュレーション対象である蓄電素子２について説明する。
　図１５は蓄電素子２の構成を説明する説明図である。蓄電素子２は、例えば金属リチウム電池である。蓄電素子２は、正極集電体層２１、正極活物質層２２、電解質層２３、及び負極集電体層２４の積層体を含む。

　正極集電体層２１は、金属箔、金属メッシュ等により構成される。正極集電体層２１を構成する金属は、アルミニウム、ニッケル、チタン、ステンレス鋼等の良好な導電性を有する金属である。正極集電体層２１の表面には、接触抵抗を調整するためのコート層が形成されていてもよい。コート層の一例は炭素コートである。

　正極活物質層２２は、リチウムイオンを吸蔵放出させることができるリチウム含有化合物により形成される。リチウム含有化合物として、例えば、Ｌｉx ＣｏＯ₂ 、Ｌｉ_x ＮｉＯ₂ 、Ｌｉ_x Ｍｎ₂ Ｏ₄ 、Ｌｉ_x ＦｅＰＯ₄ などが用いられる。正極活物質層２２は、正極活物質の他に、固体電解質、導電助剤、及びバインダを含んでもよい。

　電解質層２３は、電解液が含浸されたセパレータである。電解液は、例えば非水溶媒と、非水溶媒に溶けたリチウム塩とを含む。非水溶媒の例としては、環状炭酸エステル溶媒、環状エーテル溶媒、鎖状エーテル溶媒、環状エステル溶媒、鎖状エステル溶媒などが挙げられる。リチウム塩の例としては、ＬｉＰＦ₆ 、ＬｉＢＦ₄ 、ＬｉＮ（ＳＯ₂ＣＦ₃ ）₂ 、ＬｉＮ（ＳＯ₂ Ｃ₂ Ｆ₅ ）₂ などが挙げられる。

　負極集電体層２４は、金属箔、金属メッシュ等により構成される。負極集電体層２４を構成する金属は、銅、ニッケル、チタン、ステンレス鋼等の良好な導電性を有する金属である。負極集電体層２４の表面において、例えばリチウム金属又はリチウム合金が溶解又は析出する。ここで、溶解とは、完全に溶解することに限定されず、部分的に溶解することを含む。すなわち、リチウム金属又はリチウム合金は、蓄電素子２の放電状態において、負極集電体層２４の表面に残存してもよい。

　蓄電素子２は、拘束部材３によって拘束されてもよい。拘束部材３は、例えば、蓄電素子２を収容するケース３１と、ケース３１内に圧縮された状態にて配置される弾性部材３２とを備える。ケース３１は、例えば直方体状の容器であり、底面部３１１と側面部３１２とにより構成されるケース本体３１０と、ケース本体３１０の開口を閉塞する蓋体３１３とを備える。ケース本体３１０（底面部３１１及び側面部３１２）並びに蓋体３１３は、例えば、ステンレス鋼、アルミニウム、アルミニウム合金など溶接可能な金属によって形成される。代替的に、ケース本体３１０（底面部３１１及び側面部３１２）並びに蓋体３１３は、樹脂によって形成されてもよい。ケース本体３１０は、蓄電素子２が収容された後、蓋体３１３によって密閉される。

　弾性部材３２は、蓄電素子２の最下層（図１５の例では正極集電体層２１）と底面部３１１との間、及び、蓄電素子２の最上層（図１５の例では負極集電体層２４）と蓋体３１３との間に、圧縮された状態にて配置される。弾性部材３２は例えばゴム状のシートである。弾性部材３２は、その弾性力によって、蓄電素子２に積層方向（図の上下方向）の拘束力を付与する。

　図１５の例では、ケース３１の内部に弾性部材３２が配置されることによって、蓄電素子２に対して拘束力が付与される構成とした。代替的に、ケース３１内に高圧の流体が充填されることによって蓄電素子２に拘束力が付与されてもよい。この場合、流体としては、電池材料に対して不要な反応を生じさせないものが好ましい。例えば、窒素等の不活性ガスや乾燥空気等が用いられる。代替的に、蓄電素子２が積層方向の両側から板部材によって挟み込まれ、蓄電素子２に拘束力が付与された状態にて板部材同士が連結されることによって、蓄電素子２に拘束力が付与される構成としてもよい。

　図１６は析出物の生成量と固有歪みとの関係の一例を示すグラフである。図１６に示すグラフの横軸は析出物の生成量を表し、縦軸は固有歪みを表す。長時間にわたる充放電の繰り返しにより、蓄電素子２の内部には析出物が析出することがある。例えば、蓄電素子２の負極にリチウム金属が用いられている場合、長期間にわたる充放電の繰り返しによって、負極表面にはリチウム金属が疎に析出することがある。析出形態は密である場合もある。

　グラフに示すように、蓄電素子２の固有歪みは、析出物の析出量に比例するが、析出形態によって大きさが異なる。同じ析出量であっても、稠密な金属が析出した場合には固有歪みは比較的小さくなり、疎に析出した場合には固有歪みは比較的大きくなる。苔状の析出物が析出した場合には、固有歪みは両者の中間の大きさとなる。

　図１７は析出物の生成量と固有歪みの関係の他の例を示すグラフである。図１７に示すグラフの横軸は析出物の生成量を表し、縦軸は固有歪みを表す。グラフは、析出の初期段階において稠密な析出物が析出し、次いで疎な析出物が析出している様子を示している。蓄電素子２の固有歪みは、稠密な析出物が析出する段階、及び疎なリチウムが析出する段階の双方において、析出物の生成量に比例する値として与えられる。

　推定装置１の記憶部１２には、析出物の生成量を蓄電素子２の固有歪みに変換する関数若しくはテーブルが記憶される。推定装置１の演算部１１は、記憶部１２に記憶された関数若しくはテーブルを参照し、析出物の生成量を与えたときの蓄電素子２の固有歪みを算出する。

　演算部１１は、算出した固有歪みに基づき、蓄電素子２における応力又は歪みの分布を推定する。

　例えば、変形異方性のない線形弾性体の応力と歪みの関係式は、数３０によって表される。代替的に、対象物の構成材料によって、弾塑性体や脆性材料などの特性を表現する式が用いられてもよい。

　ここで、σ_ii，σ_ijを要素としたテンソルは応力テンソルを表す。σ_iiは法線がｉ方向の面に働く垂直応力、σ_ijは法線がｉ方向の面に働くｊ方向のせん断応力を表す。ラメ定数λ，μによって記述されるテンソルは弾性テンソルを表す。ラメ定数λ，μは、ヤング率及びポアソン比によって表される。ラメ定数λ，μは析出物の生成量の関数として与えてもよい。εは歪みテンソルであり、ε_ii，ε_ijを要素に含む。ε_iiは法線がｉ方向の面に働く垂直歪み、σ_ijは法線がｉ方向の面に働くｊ方向のせん断歪みを表す。ε⁰ は固有歪みテンソルであり、ε_ii ⁰，ε_ij ⁰ を要素に含む。ε_ii ⁰ は法線がｉ方向の面に働く垂直固有歪み、ε_ij ⁰は法線がｉ方向の面に働くｊ方向のせん断固有歪みを表す。実施の形態７において、ε_ii ⁰ ，ε_ij ⁰ は析出物の生成量の関数として与えられる。

　推定装置１の演算部１１は、固有歪みを析出物が析出する部材に与え、任意の拘束条件にて、数３０の応力・歪みの関係式と、数３１の力及びモーメントのつり合い式とを解くことによって、応力又は歪みの分布を算出する。

　以下、推定装置１が実行する処理の手順を説明する。
　図１８は応力・歪み分布の計算手順を示すフローチャートである。演算部１１は、シミュレーション対象の蓄電素子２について析出物の生成量を与え（ステップＳ７０１）、析出物が析出する部材の固有歪みを計算する（ステップＳ７０２）。このとき、演算部１１は、記憶部１２から、析出物の生成量を固有歪みに変換する関数若しくはテーブルを読み出し、読み出した関数若しくはテーブルに従って、生成量を固有歪みに変換すればよい。

　次いで、演算部１１は、蓄電素子２に対して拘束条件を与え（ステップＳ７０３）、数３０及び数３１に基づき、応力・歪みの関係式と、力及びモーメントのつり合い式とに基づき、応力又は歪みの分布を計算する（ステップＳ７０４）。

　以上のように、実施の形態７に係る推定装置１は、蓄電素子２の内部に析出する析出物の生成量を考慮して、蓄電素子２の応力又は歪みの分布を計算できる。

（実施の形態８）
　実施の形態８では、蓄電素子２の内部応力に基づき、析出物の生成速度を算出し、算出した生成速度に基づき、析出物の生成量を算出する構成について説明する。

　蓄電素子２の内部における析出物の生成は、析出物生成反応の場（例えば負極表面）に生じる応力に影響を受ける。析出物の生成速度は、析出物生成反応の場に生じる応力の関数として以下のように表される。

　ここで、Ｒ_p は析出物の生成速度（ｋｇ／（ｓｍ² ））であり、σは応力テンソル（Ｎ／ｍ² ）である。

　析出物の生成は、応力だけでなく、充放電の電流密度や過電圧の経時的な変化の影響を受ける。析出物の生成速度は、数３３のように記述されてもよい。

　ここで、ｉ_react は反応電流密度（Ａ／ｍ² ）であり、ηは過電圧（Ｖ）である。

　析出物の生成速度が数３２（又は数３３）によって与えられる場合、演算部１１は、数３４の演算により、析出物の生成量を算出できる。

　演算部１１は、算出した析出物の生成量に基づき、固有歪みを算出できる。演算部１１は、算出した固有歪みに基づき、応力又は歪みの分布を算出できる。

　演算部１１は、数３５により、析出物の生成量を面密度として分布を算出し、固有歪みも分布として算出してもよい。

　以下、生成速度Ｒ_p を表す式の一例が説明される。
　析出物の生成反応を電極反応とみなした場合、反応電流密度は、バトラーボルマー式を用いて、以下のように表される。

　ここで、ｉ_{react_p} は析出物の生成反応における反応電流密度（Ａ／ｍ² ）であり、ｉ_{0_p} は交換電流密度（Ａ／ｍ² ）である。α_a ，α_c は酸化反応、還元反応それぞれの電荷移動係数、ｎは反応に関与する電子数、Ｆはファラデー定数（Ｃ／ｍｏｌ）、ηは過電圧（Ｖ）、Ｒは気体定数（Ｊ／（ｍｏｌＫ））、Ｔは温度（Ｋ）を表す。φ_s は固相電位（Ｖ）、φ_l は液相電位（Ｖ）、Ｅ_eqは平衡電位（Ｖ）を表す。

　反応電流密度ｉ_{react_p} が数３６によって与えられる場合、析出物の生成速度Ｒ_p は、数３７のように表される。

　ここで、Ｍ_Liは析出物（例えばリチウム金属）のモル質量（ｋｇ／ｍｏｌ）を表す。このとき、電極の膨れ量Δｌ（ｍ）は、数３８によって算出される。ρ_p は析出物の密度（ｋｇ／ｍ³ ）である。

　析出物の析出量がゼロのときの、析出面に垂直な方向に対する電極の寸法をｌ（ｍ）とし、その電極上にΔｌの析出物が生成した場合に、ε_l を以下の数３９のように定義する。数３９は、析出物が析出した場合であっても、析出物が生成する電極部分は全く変形しないと仮定した場合に生じる、析出方向への圧縮歪みを表す。

　析出物の生成による合材膨れ方向（析出面に垂直な方向）の単位ベクトルをｎ＝（ｎ_x ，ｎ_y ，ｎ_z ）としたとき、固有歪みテンソルε⁰ は、数４０のように記述される。

　推定装置１の演算部１１は、例えば、数３６～数４０に基づき、固有歪みテンソルε⁰ を導出し、実施の形態７にて説明した数３０及び数３１を用いることにより、蓄電素子２の応力又は歪みの分布を算出する。

　以上のように、実施の形態８に係る推定装置１は、析出物の生成量と応力・歪み分布とを連成させて計算できる。

（実施の形態９）
　実施の形態９では、析出物の生成量及び応力場を考慮して、蓄電素子２の電気化学現象をシミュレートする構成について説明する。

　蓄電素子２の内部に生じる析出物の生成量及び応力は、蓄電素子２の電池特性に影響を与える。蓄電素子２の固有導電率、液相導電率、交換電流密度は、析出物の生成量、応力、時間の関数として、以下のように表される。

　ここで、ｉ₀ は交換電流密度（Ａ／ｍ² ）、σ_l は液相導電率（Ｓ／ｍ）、σ_s は固相導電率（Ｓ／ｍ）、ｍ_p は析出物の生成量（ｋｇ）又は析出物の面密度（ｋｇ／ｍ² ）、σは応力テンソル、ｔは時間を表す。

　演算部１１は、数４１を用いてＮｅｗｍａｎモデルなどの電気化学モデルを解くことにより、析出物の生成量及び応力の影響を考慮した蓄電素子２の電気化学現象をシミュレートできる。

　Ｎｅｗｍａｎモデルは、以下において説明するＮｅｒｎｓｔ－Ｐｌａｎｃｋ式、電荷保存式、拡散方程式、Ｂｕｔｌｅｒ－Ｖｏｌｍｅｒ式、及びＮｅｒｎｓｔ式により記述される

　Ｎｅｒｎｓｔ－Ｐｌａｎｃｋ式は、電解質や電極におけるイオン拡散を解くための方程式であり、次式により表される。

　ここで、σ_l は液相伝導率（Ｓ／ｍ）、φ_l は液相電位（Ｖ）、Ｒは気体定数（Ｊ／（Ｋ・ｍｏｌ））、Ｔは温度（Ｋ）、Ｆはファラデー定数（Ｃ／ｍｏｌ）、ｆは活量係数、ｃ_l は電解質のイオン濃度（ｍｏｌ／ｍ³ ）、ｔ₊ はカチオン輸率である。数４２のｉ^l _totは液相電流の湧き出し（Ａ／ｍ³ ）を表す。

　電荷保存式は、活物質や集電箔での電子伝導を表す式であり、次式により表される。

　ここで、ｉ_s は固相電流密度（Ａ／ｍ² ）、φ_s は固相電位（ｖ）、σ_s は固相伝導率（Ｓ／ｍ）である。数４３のｉ^s _totは固相電流の湧き出し（Ａ／ｍ³ ）を表す。

　拡散方程式は、活粒子中での活物質の拡散を表す方程式であり、次式により表される。

　ここで、ｃ_s は固相中の活物質濃度（ｍｏｌ／ｍ³ ）、ｔは時間（ｓ）、Ｄ_s は固相中の拡散係数（ｍ² ／ｓ）である。

　Ｂｕｔｌｅｒ－Ｖｏｌｍｅｒ式は、固相と液相との界面で起こる電荷移動反応で生じる反応電流と活性化過電圧との関係を表す式である。Ｎｅｒｎｓｔ式は、活性化過電圧を決定する因子である平衡電位を表す式である。

　ここで、ｉ_react は反応電流密度（Ａ／ｍ² ）、ｉ₀ は交換電流密度（Ａ／ｍ² ）、α_a ，α_c は酸化反応、還元反応それぞれの移行係数、ηは活性化過電圧（Ｖ）、φ_s は固相電位（Ｖ）、φ_l は液相電位（Ｖ）、Ｅ_eqは平衡電位（Ｖ）、Ｅ₀ は標準電極電位（Ｖ）、ｎは酸化還元反応に寄与する電子数、ａ_Ox，ａ_Red は反応前後の化学種の活量である。Ｂｕｔｌｅｒ－Ｖｏｌｍｅｒ式は、代替的に実験値に基づき改変した式が用いられる。例えば、交換電流密度を活物質濃度やイオン濃度の関数にしたり、Ｎｅｒｎｓｔ式を用いて平衡電位Ｅ_eqを算出する代わりに開回路電位としてＳＯＣおよび開回路電位の実測データを用いるなど、適宜に改変され得る。上述した数４２～数４５の式における各パラメータは、他の物理量の関数として記述されてもよい。

　本実施の形態９では、蓄電素子２の物理モデルの一例として、Ｎｅｗｍａｎモデルが示された。充放電特性を表すのであれば、代替的に、等価回路モデルや多項式モデルなどの物理モデル以外のモデルが用いられてもよい。

　図１９は等価回路モデルの一例を示す回路図である。蓄電素子２の等価回路モデルは、抵抗器及び容量成分の組み合わせによって表現される。図１９に示す等価回路モデルにおいて、Ｒ₀ はオーム抵抗成分、Ｒ_pos は正極の反応抵抗成分、Ｃ_pos は正極の容量成分、Ｒ_neg は負極の反応抵抗成分、Ｃ_neg は負極の容量成分を表す。ただし、図１９は例示であって、直列、並列の組み合わせや電気回路素子の個数や種類に制限はない。

　推定装置１は、蓄電素子２の電気化学現象を等価回路モデルを用いて１Ｄで評価する場合、以下の数４６に基づき、各抵抗・電気容量を推定すればよい。

　ここで、ｍ_p は析出物の生成量（ｋｇ）又は析出物の面密度（ｋｇ／ｍ² ）、σは応力、ｔは時間を表す。

　以上のように、実施の形態９に係る推定装置１は、析出物の生成量及び応力場に基づき、蓄電素子２の電気化学現象をシミュレートできる。

（実施の形態１０）
　実施の形態１０では、析出物の生成量及び析出形態に基づき、蓄電素子２の熱現象をシミュレートする構成について説明する。

　蓄電素子２の内部に生じる析出物の生成量及び析出形態は、蓄電素子２の熱現象に影響を与える。蓄電素子２の発熱反応速度、熱伝導率、比熱、密度は、析出物の生成量及び析出形態の関数として、以下のように表される。

　ここで、Ｑ_TRは発熱（Ｗ／ｍ³ ）、Ｃ_p は比熱（Ｊ／（ｋｇＫ））、ｋは熱伝導率（Ｗ／（ｍＫ））、ρは密度（ｋｇ／ｍ³ ）を表す。ｍ_p は析出物の生成量（ｋｇ）であり、αは析出形態に応じて設定される係数である。係数αは、例えば、析出物の比表面積に比例する値として与えられる。代替的に、係数αは（固有歪みの絶対値／析出物の生成量）に比例する値として与えられてもよい。

　演算部１１は、数４７を用いて熱伝導方程式を解くことにより、析出物の生成量及び析出形態の影響を考慮した蓄電素子２の熱現象をシミュレートできる。熱伝導方程式は数４８により表される。

　ここで、Ｑは発熱（Ｗ／ｍ³ ）であり、数４７におけるＱ_TRに相当する。

　演算部１１は、実施の形態９において説明した電気化学モデルと組み合わせ、以下の式に基づきジュール発熱を算出してもよい。

　ここで、ｉ_l は液相電流密度（Ａ／ｍ² ）、φ_l は液相導電率（Ａ／ｍ² ）、ｉ_s は固相電流密度（Ａ／ｍ² ）、φ_s は固相導電率（Ａ／ｍ² ）である。

　以上のように、実施の形態１０に係る推定装置１は、析出物の生成量及び析出形態に基づき、蓄電素子２の熱現象をシミュレートできる。

（実施の形態１１）
　実施の形態１１では、蓄電素子２の内部に発生するガスの量を考慮して、蓄電素子２の応力又は歪みの分布を算出する構成について説明する。

　充放電の繰り返しにより、蓄電素子２の内部にガスが発生する場合がある。蓄電素子２の内部にガスが発生した場合、蓄電素子２は膨張し、膨張することによって歪みが生じる。実施の形態１１では、ガスの発生量を考慮して、蓄電素子２の応力又は歪みの分布を計算する手法について説明する。

　蓄電素子２のガスの発生量は様々な式で与えられる。例えば、反応むらがあり、電流の集中があるとガス発生が進行しやすい場合、上述した反応電流密度ｉ_react に基づき、ガス発生量ｎ_gas （ｍｏｌ）が算出される。

　ここで、Ｊは局所的な反応電流密度が平均よりも高い場合（電流集中がある場合）に１、そうでない場合に０となる変数である。ａｖｅ＿ｉ_react は平均反応電流密度（Ａ／ｍ² ）である。

　このとき、ガス発生による蓄電素子２の内圧ｐ（Ｎ／ｍ² ）は、以下の数５１により算出される。

　ここで、ｎ₀ は初期にあったガス量（ｍｏｌ）であり、ｖは蓄電素子内部の隙間の体積（ｍ³ ）である。厳密には、ギブズ＝デュエムの関係（ギブズ＝デュエムの式）に従うが、計算単純化のために数５１を用いてもよい。理想気体と仮定して数５１を用いてもよいが、分子間力を考慮した気体の状態方程式を用いてもよい。

　このとき応力テンソルσは、蓄電素子２の内圧ｐを用いて、次式のように記述される。

　演算部１１は、数５２に示す応力を加味し、数３０の応力・歪みの関係式と、数３１の力及びモーメントのつり合い式とを解くことによって、応力又は歪みの分布を算出する。

　以上のように、実施の形態１１に係る推定装置１は、蓄電素子２の内部に発生するガスの量を考慮して、蓄電素子２の応力又は歪みの分布を算出できる。

　今回開示された実施形態は、全ての点において例示であって、制限的なものではないと考えられるべきである。本発明の範囲は、上述した意味ではなく、請求の範囲によって示され、請求の範囲と均等の意味及び範囲内での全ての変更が含まれることが意図される。

　例えば、蓄電素子２は、複数のセルを直列に接続したモジュール、複数のモジュールを直列に接続したバンク、複数のバンクを並列に接続したドメイン等であってもよい。

　１　推定装置
　２　蓄電素子
　３　拘束部材
　１１　演算部
　１２　記憶部
　１３　入力部
　１４　出力部
　２１　正極集電体層
　２２　正極活物質層
　２３　固体電解質層
　２４　負極活物質層
　２５　負極集電体層
　３１　ケース
　３２　弾性部材
　３１０　ケース本体
　３１１　底面部
　３１２　側面部
　３１３　蓋体
　ＭＤ１　シミュレーションモデル
　ＰＧ１　推定プログラム
　Ｓ１　歪みセンサ
　Ｓ２，Ｓ３　温度センサ
　Ｓ４　電流計
　Ｓ５　電圧計

Claims

　蓄電素子に発生する歪みに係るデータを取得する取得部と、
　前記蓄電素子の内部の力学的状態を表すシミュレーションモデルを用いて、前記取得部が取得したデータに基づき、前記蓄電素子の内部応力を推定する推定部と
　を備える推定装置。
　前記シミュレーションモデルは、前記蓄電素子の固有歪みと、前記蓄電素子に対する拘束力とをパラメータに含み、前記歪みに係るデータの入力に応じて、前記蓄電素子の内部応力に係るデータを出力するよう構成される
　請求項１に記載の推定装置。
　前記固有歪みは、活物質粒子の孤立化、析出物の成長、及び熱膨張の少なくとも１つに起因して生じる前記蓄電素子の歪みである
　請求項２に記載の推定装置。
　前記推定部は、非線形フィルタを用いた状態推定器を備える
　請求項１から請求項３の何れか１つに記載の推定装置。
　前記推定部は、前記蓄電素子の内部抵抗を前記内部応力の関数として推定する
　請求項１から請求項４の何れか１つに記載の推定装置。
　前記蓄電素子は、電解質が固体である全固体電池である
　請求項１から請求項５の何れか１つに記載の推定装置。
　前記蓄電素子は、負極に金属リチウムを用いた電池である
　請求項１から請求項６の何れか１つに記載の推定装置。
　活物質粒子と固体電解質との接触面積をパラメータに含むシミュレーションモデルを用いて、前記固体電解質を含む全固体電池の電気化学現象をシミュレートする演算部
　を備える推定装置。
　充放電に応じて析出物が生じる蓄電素子に関して、前記析出物の生成量及び析出形態に基づき、前記蓄電素子の固有歪みを算出し、
　算出した固有歪みに基づき、前記蓄電素子に生じる応力又は歪みの分布を算出する
　演算部を備える推定装置。
　蓄電素子に発生する歪みに係るデータを取得し、
　前記蓄電素子の内部の力学的状態を表すシミュレーションモデルを用いて、取得したデータに基づき、前記蓄電素子の内部応力を推定する
　処理をコンピュータによって実行する推定方法。
　活物質粒子と固体電解質との接触面積をパラメータに含むシミュレーションモデルを用いて、前記固体電解質を含む全固体電池の電気化学現象をシミュレートする
　処理をコンピュータによって実行する推定方法。
　充放電に応じて析出物が生じる蓄電素子に関して、前記析出物の生成量及び析出形態に基づき、前記蓄電素子の固有歪みを算出し、
　算出した固有歪みに基づき、前記蓄電素子に生じる応力又は歪みの分布を算出する
　処理をコンピュータによって実行する推定方法。
　コンピュータに、
　活物質粒子と固体電解質との接触面積をパラメータに含むシミュレーションモデルを用いて、前記固体電解質を含む全固体電池の電気化学現象をシミュレートする
　処理を実行させるためのコンピュータプログラム。
　前記シミュレーションモデルは、前記接触面積と前記活物質粒子の有効拡散係数との関係を規定してあり、
　前記接触面積の値に基づき、前記活物質粒子の有効拡散係数を推定する処理を前記コンピュータに実行させるための請求項１３に記載のコンピュータプログラム。
　前記シミュレーションモデルは、前記接触面積と前記活物質粒子の有効イオン伝導度との関係を規定してあり、
　前記接触面積の値に基づき、前記活物質粒子の有効イオン伝導度を推定する処理を前記コンピュータに実行させるための請求項１３に記載のコンピュータプログラム。
　前記接触面積は、前記全固体電池の内部応力の関数であり、
　前記内部応力の値に基づき、前記全固体電池の電気化学現象をシミュレートする処理を前記コンピュータに実行させるための請求項１３から請求項１５の何れか１つに記載のコンピュータプログラム。
　前記全固体電池に発生する歪みを計測する歪みセンサから、前記歪みに係る計測データを取得し、
　前記全固体電池の内部の力学的状態を表すモデルを用いて、取得した計測データに基づき、前記全固体電池の内部応力を推定する
　処理を前記コンピュータに実行させるための請求項１６に記載のコンピュータプログラム。
　前記全固体電池の内部抵抗は、前記内部応力の関数であり、
　前記内部応力の値に基づき、前記内部抵抗の値を推定する処理を前記コンピュータに実行させるための請求項１６又は請求項１７に記載のコンピュータプログラム。
　コンピュータに、
　充放電に応じて析出物が生じる蓄電素子に関して、前記析出物の生成量及び析出形態に基づき、前記蓄電素子の固有歪みを算出し、
　算出した固有歪みに基づき、前記蓄電素子に生じる応力又は歪みの分布を算出する
　処理を実行させるためのコンピュータプログラム。
　前記析出物の生成速度を、生成反応の場に生じる応力の関数として記述してあり、
　前記関数により算出される前記析出物の生成速度に基づき、前記生成量を算出する
　処理を前記コンピュータに実行させるための請求項１９に記載のコンピュータプログラム。
　前記生成量及び応力場に基づき、前記蓄電素子の電気化学現象をシミュレートする
　処理を前記コンピュータに実行させるための請求項１９又は請求項２０に記載のコンピュータプログラム。
　前記生成量及び析出形態に基づき、前記蓄電素子における熱現象をシミュレートする
　処理を前記コンピュータに実行させるための請求項１９から請求項２１の何れか１つに記載のコンピュータプログラム。
　前記蓄電素子内部のガス発生量を算出し、
　算出したガス発生量に基づき、前記蓄電素子に生じる応力又は歪みの分布を算出する
　処理を実行させるための請求項１９から請求項２２の何れか１つに記載のコンピュータプログラム。