Turn Excel into a Physics-Informed Simulator — Free Algebraic Reduction Method for Electrochemical Interfaces (VBA Code Included)

【無料公開】ExcelがマルチフィジックスFEMに化ける——「電気化学×流体渦×交流応答」を一瞬で解く代数縮約ソルバー（VBAコード付）

1. なぜ界面シミュレーションは「現場で回らない」のか？

めっき、半導体配線（ダマシン）、リチウムイオン電池、電解銅箔。

ものづくりの現場で「固液界面」を扱ったことのある理系技術者や研究者なら、誰もが一度はこの壁にぶつかります。

> 「市販のマルチフィジックスFEM（有限要素法）は、重すぎて現場のライン制御や日常の試行錯誤では回らない」

> 

界面近傍の現象を真面目に追及しようとすると、複数の物理が複雑に連立します。

 * 流体力学（Navier-Stokes）：攪拌や噴流ノズルが作る対流と渦

 * 物質移動（Nernst拡散層）：界面近傍の数十μmに生じる急峻な濃度勾配

 * 電荷移動反応（Butler-Volmer則）：過電圧と活性化障壁に支配される非線形反応

 * 界面化学（吸着平衡）：添加剤分子（にかわ、PEG、SPS等）による活性サイトの被覆・阻害

これらをメッシュで細かく切り、非線形方程式をニュートン法で反復収束させると、1条件の計算に数十分から数時間かかります。論文発表なら良くても、「現場でリアルタイムにインライン制御したい」「数千通りの条件を一瞬でスクリーニングしたい」という要求には到底耐えられません。

筆者はこれまで、電気化学の世界的泰斗である逢坂哲彌 先生（元ECSプレジデント）、半導体・無電解めっきの基盤を築かれた沖中裕 先生（元AT&Tベル研究所フェロー）、ならびに機能薄膜工学をリードされる小岩一郎 教授らと学術研究・共著論文の執筆に携わってきました。

学術の最前線から現場実装までを経験して導き出された結論は一つです。

「現場が真に求めているのは、超巨大な微分方程式を重く解くことではなく、物理の本質を代数的に縮約した『1ステップで動く陽解法』である」ということ。

本記事では、この連立非線形系を1本の閉じた代数解へと縮約し、Excel VBAだけで瞬時に界面状態を予測するオープンな計算モデルを全コード付きで無料公開します。

2. 数理モデルの核心：ニュートン法を消し去る「特異境界代数縮約」

① Butler-Volmer式を1本の陽解へ

カソード析出電流密度 i_c は、界面濃度 C_s に依存する物質移動限界と、電荷移動反応が連立した非線形方程式です。

定常拡散層厚み \delta を横切る流束の釣り合い：

から、界面濃度比 C_s / C_b = 1 - i_c / i_L をButler-Volmer式へ代入して整理すると、反復収束（ニュートン法）なしで求まる陽解が得られます。

 * K_c = \exp\left( \frac{\alpha_c z F}{R T} \eta_c \right)、K_a = \exp\left( - \frac{\alpha_a z F}{R T} \eta_c \right)

 * i_L = \frac{z F D C_b}{\delta}（物質移動限界電流密度）

 * \theta：添加剤による被覆率（0.0 \sim 0.99）

流体渦の角速度 \omega_{\text{vort}} からLevich型の境界層厚み \delta を動的に更新することで、流動攪拌 \rightarrow 拡散層収縮 \rightarrow 限界電流増大 \rightarrow 界面電流の追従が一瞬で定まります。

② 交流インピーダンス（EIS）の「半球の右肩」とミクロ移動の直結

早稲田大学の門間聰之教授（逢坂・門間研）らが切り拓かれた導電性高分子電極（ポリピロール膜等への過塩素酸基・硫酸基導入）の解析知見では、交流インピーダンス（EIS）のナイキスト線図における「高周波円弧（半球）の右肩」が、「電子移動から膜内イオンの吸脱着・鎖間ホッピング律速へ移り変わる境界」を表します。

本モデルでは、分子内移動距離や吸脱着の緩和時間 \tau_s に基づく複素緩和弾性 G_{\text{re}} を導入。

交流サイクル時間と分子緩和時間が一致する「半球の右肩」の変曲点を代数的に捉え、膜質変化や劣化の兆候を捉えます。

③ 断面組織・SEM観察と直結する平滑性指数 \Pi_{\text{cat}}

Mullins-Sekerka形態不安定性と添加剤効果に基づき、皮膜の健全性を無次元指数で判定します。

 * 界面濃度比（C_s/C_b）：イオン枯渇が進むと急落（デンドライト・ヤケ発生の危険域）。

 * 添加剤効果（1 + 2.5\theta）：吸着分子が突起をマスキングして平滑化を促す効果を反映。

3. コピペで動くExcel VBAコード

Excelを開き、Alt + F11 ➔「挿入」➔「標準モジュール」を作成し、以下をそのまま貼り付けてください。

Option Explicit

'========================================================================================

' 界面ダイナミクス・電気化学マスター方程式 多重スケール代数縮約ソルバー

'========================================================================================

Sub SolveMultiScaleInterfaceSolver()

    Dim ws As Worksheet

    Set ws = ActiveSheet

    ' --- 1. 基本物理定数 ---

    Const F As Double = 96485.332 ' ファラデー定数 [C/mol]

    Const R As Double = 8.314462 ' 気体定数 [J/(mol*K)]

    Const PI As Double = 3.14159265358979

    Const c_light As Double = 299792458# ' 真空中光速 [m/s]

    Const h_bar As Double = 1.0545718e-34 ' ディラック定数 [J*s]

    Const G_grav As Double = 6.6743e-11 ' 万有引力定数 [m^3/(kg*s^2)]

    

    Dim l_Planck As Double

    l_Planck = Sqr((h_bar * G_grav) / (c_light ^ 3))

    ' --- 2. パラメータ読み込み (G列) ---

    Dim T As Double, z As Double, Cb As Double, D_diff As Double

    Dim i0 As Double, alpha_c As Double, alpha_a As Double

    Dim M_mol As Double, rho_metal As Double, theta As Double, eta_eff As Double

    Dim nu_fluid As Double, gamma_0 As Double, tau_s As Double, lambda_0 As Double

    Dim scale_mode As Integer

    T = ws.Range("G2").Value ' 液温 [K] (例: 298.15)

    z = ws.Range("G3").Value ' イオン価数 (例: Cu = 2)

    Cb = ws.Range("G4").Value ' バルク濃度 [mol/m^3] (例: 1000 = 1.0 M)

    D_diff = ws.Range("G5").Value ' 拡散係数 [m^2/s] (例: 1.0e-9)

    i0 = ws.Range("G6").Value ' 交換電流密度 [A/m^2] (例: 20)

    alpha_c = ws.Range("G7").Value ' カソード移動係数 (例: 0.5)

    alpha_a = 1# - alpha_c

    M_mol = ws.Range("G8").Value ' 原子量 [g/mol] (例: 63.55)

    rho_metal = ws.Range("G9").Value ' 密度 [g/cm^3] (例: 8.96)

    theta = ws.Range("G10").Value ' 添加剤被覆率 [0.0 - 0.99]

    eta_eff = ws.Range("G11").Value ' 電流効率 [0.0 - 1.0] (例: 0.98)

    

    nu_fluid = IIf(ws.Range("G12").Value > 0, ws.Range("G12").Value, 1e-6)

    gamma_0 = IIf(ws.Range("G13").Value > 0, ws.Range("G13").Value, 1.3)

    tau_s = IIf(ws.Range("G14").Value > 0, ws.Range("G14").Value, 1e-4)

    lambda_0 = IIf(ws.Range("G15").Value > 0, ws.Range("G15").Value, 1e-6)

    scale_mode = CInt(ws.Range("G16").Value) ' 1:Nano, 2:Meso[推奨], 3:Macro, 0:Limit

    ' --- 3. 現象スケール別スケーリング ---

    Dim c_eff As Double, l_min As Double

    Select Case scale_mode

        Case 1 ' Nano (核生成・極薄膜)

            c_eff = 50# ' 表面弾性波 [m/s]

            l_min = 0.25e-9 ' 原子半径スケール [m]

        Case 2 ' Meso (電解銅箔・ビア配線・添加剤吸着) ※標準

            c_eff = 1500# ' 音響波速度 [m/s]

            l_min = 5e-9 ' 高分子吸着鎖長 [m]

        Case 3 ' Macro (大型めっき槽・主流流体)

            c_eff = 2# ' 攪拌主流流速 [m/s]

            l_min = 50e-6 ' コルモゴロフ渦長 [m]

        Case Else

            c_eff = c_light: l_min = l_Planck

    End Select

    Dim Gamma_min As Double, Omega_c As Double

    Gamma_min = 2# * PI * l_min * c_eff

    Omega_c = lambda_0 / (2# * PI * c_eff * tau_s)

    ' --- 4. 時系列入力の読み込み ---

    Dim lastRow As Long

    lastRow = ws.Cells(ws.Rows.Count, "A").End(xlUp).Row

    Dim N As Long: N = lastRow - 1

    If N < 2 Then

        MsgBox "A列(時間), B列(過電圧), C列(渦角速度 rad/s) を入力してください。", vbExclamation

        Exit Sub

    End If

    Dim t() As Double, eta_c() As Double, w_vort() As Double

    ReDim t(1 To N), eta_c(1 To N), w_vort(1 To N)

    Dim i As Long

    For i = 1 To N

        t(i) = ws.Cells(i + 1, 1).Value

        eta_c(i) = ws.Cells(i + 1, 2).Value

        w_vort(i) = ws.Cells(i + 1, 3).Value

    Next i

    ' --- 5. 出力ヘッダー生成 ---

    ws.Cells(1, 5).Value = "限界電流 iL (A/m^2)"

    ws.Cells(1, 6).Value = "界面電流 ic (A/m^2)"

    ws.Cells(1, 7).Value = "表面濃度比 Cs/Cb"

    ws.Cells(1, 8).Value = "積算膜厚 H (nm)"

    ws.Cells(1, 9).Value = "複素緩和弾性 G_re (-)"

    ws.Cells(1, 10).Value = "膜応力補正因子 m_eff/m0"

    ws.Cells(1, 11).Value = "循環比 Gamma/Gamma_min"

    ws.Cells(1, 12).Value = "平滑性判定指数 Pi_cat (-)"

    ' --- 6. 一括代数計算ループ ---

    Dim outData() As Variant

    ReDim outData(1 To N, 1 To 8)

    Dim f_factor As Double, v_factor As Double

    f_factor = (z * F) / (R * T)

    v_factor = ((M_mol / 1000#) / ((rho_metal * 1000#) * z * F)) * eta_eff

    Dim H_acc As Double: H_acc = 0#

    Dim dt As Double, delta As Double, iL As Double, ic As Double, Cs_ratio As Double

    Dim Kc As Double, Ka As Double, v_nms As Double

    Dim omega_hat As Double, G_re As Double, warp_factor As Double

    Dim Gamma_vort As Double, tau_shear As Double, Pi_cat As Double

    For i = 1 To N

        dt = IIf(i = 1, 0#, t(i) - t(i - 1))

        ' ① 流体渦による拡散層厚み delta の動的評価

        If Abs(w_vort(i)) > 0.01 Then

            delta = 1.61 * (D_diff ^ (1# / 3#)) * (nu_fluid ^ (1# / 6#)) * (Abs(w_vort(i)) ^ (-0.5))

        Else

            delta = 150e-6

        End If

        iL = (z * F * D_diff * Cb) / delta

        ' ② 複素緩和弾性

        omega_hat = Abs(w_vort(i)) * (lambda_0 / c_eff)

        G_re = (Omega_c ^ 2) / (Omega_c ^ 2 + omega_hat ^ 2)

        ' ③ 特異境界直接陽解法

        If eta_c(i) > 0.00001 Then

            Kc = Exp(alpha_c * f_factor * eta_c(i))

            Ka = Exp(-alpha_a * f_factor * eta_c(i))

            ic = ((1# - theta) * i0 * (Kc - Ka)) / (1# + (1# - theta) * (i0 / iL) * Kc)

            If ic > iL Then ic = iL

            If ic < 0# Then ic = 0#

            Cs_ratio = 1# - (ic / iL)

            If Cs_ratio < 0# Then Cs_ratio = 0#

            v_nms = (ic * v_factor) * 1000000000#

        Else

            ic = 0#: Cs_ratio = 1#: v_nms = 0#

        End If

        H_acc = H_acc + v_nms * dt

        ' ④ せん断応力と張力バランス

        tau_shear = (1000# * nu_fluid * (Abs(w_vort(i)) * lambda_0)) / delta + 1e-6

        warp_factor = 1# / (1# + (tau_shear / (gamma_0 / lambda_0)))

        ' ⑤ 循環評価

        Gamma_vort = 2# * PI * (lambda_0 ^ 2) * Abs(w_vort(i))

        ' ⑥ 平滑性・相転移判定指数 Pi_cat

        Pi_cat = Cs_ratio * G_re * (1# + 2.5 * theta) * warp_factor

        outData(i, 1) = iL

        outData(i, 2) = ic

        outData(i, 3) = Cs_ratio

        outData(i, 4) = H_acc

        outData(i, 5) = G_re

        outData(i, 6) = warp_factor

        outData(i, 7) = Gamma_vort / Gamma_min

        outData(i, 8) = Pi_cat

    Next i

    ws.Range(ws.Cells(2, 5), ws.Cells(N + 1, 12)).Value = outData

    MsgBox "解析完了（スケールモード: " & scale_mode & "）", vbInformation

End Sub

4. ワークシートの入力設定

G列：パラメータ設定（硫酸銅めっき標準例）

 * G2（液温 T）：298.15 [K]

 * G3（価数 z）：2 (Cu²⁺)

 * G4（濃度 C_b）：1000 [mol/m³] (1.0 M)

 * G5（拡散係数 D）：1.0E-9 [m²/s]

 * G6（交換電流 i_0）：20 [A/m²]

 * G7（移動係数 \alpha_c）：0.5

 * G8（原子量 M）：63.55 [g/mol]

 * G9（密度 \rho）：8.96 [g/cm³]

 * G10（被覆率 \theta）：0.3

 * G11（電流効率）：0.98

 * G12（動粘性係数）：1.0E-6 [m²/s]

 * G13（界面張力）：1.3 [N/m]

 * G14（緩和時間 \tau_s）：1.0E-4 [s]

 * G15（代表長さ \lambda_0）：1.0E-6 [m]

 * G16（スケール選択）：2（1:Nano, 2:Meso[推奨], 3:Macro）

A〜C列：時系列プロファイル

 * A列（時間 t）：0, 1, 2, 3... [秒]

 * B列（過電圧 \eta_c）：0.15（定電位）やパルス波形

 * C列（渦角速度 \omega_{\text{vort}}）：10, 50, 150... [rad/s]（攪拌流速）

5. 現場・投資家・研究者が見るべきポイント

 * 理系技術者・研究者へ：

   何千行もの時系列プロファイルを入力しても、数秒で計算が完了します。条件スクリーニングや、実験で測定した電位波形を食わせて「どの瞬間に限界電流に達していたか」の逆算検証にそのままご活用ください。

 * 製造業DX・投資家へ：

   現在の先端製造業（半導体PKG・電池）の課題は「ブラックボックスな物理現象をどうソフトウェアで制御するか」です。単なる統計的AI（深層学習）では解けない界面物理を、こうした数理モデルを組み込んだ「Physics-Informed」な制御プログラムにすることで、極少データでのリアルタイム自動補正が可能になります。

> 【知的財産・産業実装に関するご案内】

> ※本記事で公開しているVBAコードは、基礎物理および電気化学の公知理論を代数縮約したオープンな学術研究・事前検討用ツールです。

> 当研究所（Sakaqi Interface）が出願済みの特許ポートフォリオ（特願2026-120730、特願2026-206793等）に含まれる「等価回路を用いずに実測EISデータから直接状態を逆算し、インラインで自動補正を行うリアルタイム制御プログラムおよびハードウェア統合システム」のコアアルゴリズムは本コードには含まれておりません。

> 工場ラインへのインライン自動制御の実装、添加剤濃度のリアルタイム特定、有償PoC（共同評価）に関するご相談は、X（@interfacesakaqi）のDM等よりお気軽にお問い合わせください。

>