Page 129 - Advances in Biomechanics and Tissue Regeneration

7.2 EQUATIONS THAT GOVERN THE ELECTRICAL ACTIVITY OF THE HEART 125

k
V k +1 ¼ V + ΔtV k + θ , (7.38)
k + θ k k +1
V ¼ð1 θÞV + θV , (7.39)
where θ 2 [0, 1] is a scalar parameter. Eqs. (7.37)–(7.39) can be combined to obtain an algebraic system of equations to
determine V k+1 .
When using the operator splitting algorithm for solving the monodomain model, the equations are solved in two
steps. First, the electrophysiological cellular model

k k (7.40)
∗
V ¼ V ΔtJ ion ðV ,uÞ + J stm ðtÞ
is solved at each mesh point to obtain an intermediate transmembrane potential vector V* (Step I). Even though a
forward Euler scheme has been used in Eq. (7.40), any other ODE solver can be used to calculate V*. With this inter-
k k
mediate solution at hand, along with Eqs. (7.38), (7.39) and MV ¼ KV from the previous converged time increment,
Eq. (7.37) becomes
V k +1 V ∗ k +1 k
M ¼ K θ V + ð1 θÞ V , (7.41)
Δt
or alternatively
^
^
KV k +1 ¼ b, (7.42)
^
where ^ K is everything that multiplies onto V k+1 , and b contains the other terms in Eq. (7.41). Eq. (7.42) is solved for the
entire domain to the obtained V k+1 (Step II).
Hence, the basic algorithm at time t k+1 can be summarized, as:
k
• Step I:Use V as the initial condition to integrate Eq. (7.40) to obtain V*.
• Step II: Use the result obtained in Step I to solve Eq. (7.42) for V k+1 .
For different values of the parameter θ, different time integration schemes are obtained for integrating the discre-
tized homogeneous parabolic equation system (7.41):
θ ¼ 0 Forward Euler (conditionally stable)
θ ¼ 0.5 Crank-Nicolson scheme (unconditionally stable)
2 Galerkin scheme (unconditionally stable)
θ ¼ 3
θ ¼ 1 Backward Euler (unconditionally stable)
The Crank-Nicolson scheme is second-order accurate in time, whereas the others are first-order accurate in time.
However, for θ 0:5, integration schemes are unconditionally stable.
As mentioned before, Step I can be performed using a backward difference approximation in time (implicit inte-
gration) or a forward difference approximation in time (explicit integration). Implicit integration requires the solution
of a nonlinear system of equations at each point of the mesh, making it computationally costly. However, it ensures the
stability of the numerical solution. On the contrary, explicit integration is computationally cheaper but imposes more
stringent conditions on the size of the time step in order to avoid numerical instabilities.

7.2.4.2 Integration of the Mass Matrix
e
For the standard finite element formulation, the elemental mass matrix, M of Eq. (7.41), is given by [51]
Z
e
ij
M ¼ N i N j dx,
Ω e
where N j is the shape function of node j of the element e and C m ¼ 1 has been assumed without loss of generality.
Remember that C m is the membrane capacitance, as defined in Eq. (7.29).
e
When the shape functions N, used to compute M , are the same used to approximate the potential V, the resulting
nondiagonal matrix is known as the consistent mass matrix. Using a consistent mass matrix implies that a linear system
of equations has to be solved for V k+1 when an explicit scheme is used in Eq. (7.41). In order to improve numerical
e
efficiency, the proposed algorithm evaluates M using a mass preserving nodal quadrature [52]. Nodal quadrature
is based on the use of different base functions to those used to approximate the transmembrane potential, V.

I. BIOMECHANICS