Physics and Mathematical Formulation

This document describes the physical models and mathematical formulations implemented in the fuel_cell_3D code.

Reproducing Kernel Particle Method (RKPM)

RKPM is a meshfree method that constructs shape functions using a reproducing kernel approximation. The displacement field is approximated as:

\[u^h(\mathbf{x}) = \sum_{I=1}^{n_p} \Psi_I(\mathbf{x}) \, d_I\]

where \(\Psi_I(\mathbf{x})\) is the shape function for node \(I\) and \(d_I\) are the nodal coefficients.

Shape Function Construction

The RKPM shape function is defined as:

\[\Psi_I(\mathbf{x}) = \mathbf{H}^T(\mathbf{0}) \, \mathbf{M}^{-1}(\mathbf{x}) \, \mathbf{H}(\mathbf{x} - \mathbf{x}_I) \, \phi_a\left(\frac{\|\mathbf{x} - \mathbf{x}_I\|}{a_I}\right)\]

where:

• \(\mathbf{H}(\mathbf{x})\) is the polynomial basis vector

• \(\mathbf{M}(\mathbf{x})\) is the moment matrix

• \(\phi_a(z)\) is the kernel function

• \(a_I\) is the support size for node \(I\)

For 3D problems, the linear polynomial basis is:

\[\begin{split}\mathbf{H}(\mathbf{x}) = \begin{bmatrix} 1 \\ x/h \\ y/h \\ z/h \end{bmatrix}\end{split}\]

where \(h = 10^{-6}\) m is a scaling factor for numerical stability.

The moment matrix is computed as:

\[\mathbf{M}(\mathbf{x}) = \sum_{I=1}^{n_p} \mathbf{H}(\mathbf{x} - \mathbf{x}_I) \, \mathbf{H}^T(\mathbf{x} - \mathbf{x}_I) \, \phi_a(z_{xI})\]

Kernel Function

The cubic B-spline kernel is used:

\[\begin{split}\phi(z) = \begin{cases} \frac{2}{3} - 4z^2 + 4z^3 & 0 \leq z < 0.5 \\ \frac{4}{3} - 4z + 4z^2 - \frac{4}{3}z^3 & 0.5 \leq z \leq 1 \\ 0 & z > 1 \end{cases}\end{split}\]

where \(z = \|\mathbf{x} - \mathbf{x}_I\| / a_I\) is the normalized distance.

Diffusion Equation

The steady-state diffusion equation governs species transport:

\[\nabla \cdot (D \nabla c) + f = 0 \quad \text{in } \Omega\]

where:

• \(c\) is the concentration (or potential)

• \(D\) is the diffusion coefficient

• \(f\) is the source term

Weak Form

The weak form is:

\[\int_\Omega D \nabla w \cdot \nabla c \, d\Omega - \int_{\Gamma_N} w q \, d\Gamma - \int_\Omega w f \, d\Omega = 0\]

where \(w\) is the test function and \(q\) is the prescribed flux.

Nitsche’s Method for Boundary Conditions

Dirichlet boundary conditions are enforced weakly using Nitsche’s method:

\[\int_{\Gamma_D} \beta (c - g) w \, d\Gamma - \int_{\Gamma_D} D \nabla c \cdot \mathbf{n} \, w \, d\Gamma = 0\]

where:

• \(g\) is the prescribed boundary value

• \(\beta\) is the Nitsche penalty parameter

• \(\mathbf{n}\) is the outward normal vector

The discretized stiffness matrix has three components:

\[\mathbf{K} = \mathbf{K}_1 + \mathbf{K}_2 + \mathbf{K}_3\]

where:

• \(\mathbf{K}_1\): Domain diffusion integral

• \(\mathbf{K}_2\): Consistency term (flux boundary)

• \(\mathbf{K}_3\): Nitsche penalty term

Butler-Volmer Electrochemistry

The electrochemical reaction rate at the electrode-electrolyte interface follows the Butler-Volmer equation:

\[i = i_0 \left[ \exp\left(\frac{\alpha_a F \eta}{RT}\right) - \exp\left(-\frac{\alpha_c F \eta}{RT}\right) \right]\]

where:

• \(i\) is the current density

• \(i_0\) is the exchange current density

• \(\eta\) is the overpotential

• \(F\) is Faraday’s constant (96,485 C/mol)

• \(R\) is the gas constant (8.314 J/mol·K)

• \(T\) is the temperature

• \(\alpha_a, \alpha_c\) are the anodic and cathodic transfer coefficients

For small overpotentials, the linearized form is:

\[i \approx i_0 \frac{F \eta}{RT}\]

Linear Elasticity

The mechanical equilibrium equation is:

\[\nabla \cdot \boldsymbol{\sigma} + \mathbf{b} = 0 \quad \text{in } \Omega\]

where \(\boldsymbol{\sigma}\) is the Cauchy stress tensor and \(\mathbf{b}\) is the body force.

Constitutive Relation

The stress-strain relationship for linear elasticity:

\[\boldsymbol{\sigma} = \mathbf{C} : (\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^*)\]

where:

• \(\mathbf{C}\) is the fourth-order elasticity tensor

• \(\boldsymbol{\varepsilon}\) is the total strain

• \(\boldsymbol{\varepsilon}^*\) is the chemical expansion strain

Elasticity Tensor

For isotropic materials, the elasticity tensor in Voigt notation is:

\[\begin{split}\mathbf{C} = \begin{bmatrix} \lambda + 2\mu & \lambda & \lambda & 0 & 0 & 0 \\ \lambda & \lambda + 2\mu & \lambda & 0 & 0 & 0 \\ \lambda & \lambda & \lambda + 2\mu & 0 & 0 & 0 \\ 0 & 0 & 0 & \mu & 0 & 0 \\ 0 & 0 & 0 & 0 & \mu & 0 \\ 0 & 0 & 0 & 0 & 0 & \mu \end{bmatrix}\end{split}\]

where the Lamé constants are:

\[\lambda = \frac{E \nu}{(1+\nu)(1-2\nu)}, \quad \mu = \frac{E}{2(1+\nu)}\]

with \(E\) being Young’s modulus and \(\nu\) Poisson’s ratio.

Chemical Expansion Strain

The chemical expansion strain due to concentration changes:

\[\varepsilon^*_{ij} = \frac{\beta}{3} \Delta c \, \delta_{ij}\]

where \(\beta\) is the chemical expansion coefficient (m³/mol) and \(\Delta c\) is the concentration change.

Gauss Quadrature

Domain integrals are evaluated using Gauss quadrature:

\[\int_\Omega f(\mathbf{x}) \, d\Omega \approx \sum_{g=1}^{n_g} f(\mathbf{x}_g) \, |J_g| \, w_g\]

where:

• \(\mathbf{x}_g\) are the Gauss point coordinates

• \(|J_g|\) is the Jacobian determinant

• \(w_g\) are the Gauss weights

For 3D hexahedral cells, 2×2×2 Gauss integration is used with points at:

\[\xi_i = \pm \frac{1}{\sqrt{3}}, \quad w_i = 1\]

Damage Model

When enabled, a damage variable \(D\) modifies the stiffness:

\[\mathbf{C}_{eff} = (1 - D) \mathbf{C}\]

The damage variable evolves based on equivalent strain:

\[\begin{split}D = \begin{cases} 0 & \varepsilon_{eq} < \kappa_i \\ 1 - \frac{\kappa_i}{\varepsilon_{eq}} \exp\left(-\frac{\varepsilon_{eq} - \kappa_i}{\kappa_f - \kappa_i}\right) & \varepsilon_{eq} \geq \kappa_i \end{cases}\end{split}\]

where \(\kappa_i\) and \(\kappa_f\) are the initial and final damage thresholds.