Calculating reaction rate

[marbris]

Similarly to issue #121, the particle volume fraction seems to be interpreted as a fraction of the electrode volume, rather than the particles, when calculating R_Vp (mod_cell.py Lines 181-196):

$$R_{Vp}(v) = \beta (1-\epsilon) P_L \sum_{p} \frac{\mathrm{d}\bar{c}_{v,p}}{\mathrm{d}t} V_{v,p}$$

In addition, including $\beta (1-\epsilon) P_L$ seems to imply that $V_{v,p}$ is the volume fraction of the electrode, rather than of the particles. Like in issue #121, the total particle volume seems to be assumed to be equal in all electrode volumes. If R_Vp is the reaction rate in each volume (in mol/s) the code should be replaced by

# Define dimensionless R_Vp for each electrode volume
for trode in trodes:
    for vInd in range(Nvol[trode]):
        eq = self.CreateEquation(
            "R_Vp_trode{trode}vol{vInd}".format(vInd=vInd, trode=trode))
        RHS = 0
        for pInd in range(Npart[trode]):
            Vj = config["psd_vol"][trode][vInd,pInd]
            RHS += - Vj * self.particles[trode][vInd,pInd].dcbardt())
        eq.Residual = self.R_Vp[trode](vInd) - RHS

However, I'm not fully understanding how to think about the particles. Currently there seems to be two independent ways the total capacity could be defined. On the one hand it could be defined by adding up the particle volumes (rho_s, mean_c, shape, Npart_c, Nvol_c), on the other it could be defined with the geometry (L_c, P_L_c, poros_c, rho_s). The latter is what is currently used in derived_values.py lines 173-177. But the particle concentration is what actually determines muR, which in turn determines phi_m, which determines when the simulation ends.

Given this, I don't understand what P_L_c and is used for? It only shows up in scaling and when calculating R_Vp, but then it's divided from R_Vp when the total current is calculated in mod_cell.py Lines 347-357.