Numerical implementation
With the view of solving the two partial differential equations (1.16)
and (1.17) numerically, the 1-D domain [0, Ω] (see Figure 3) is
discretized into J+1 nodes numbered with superscript j . Every
node ω j is assigned a length interval
[ω j -½Δω , ωj +½Δω ], with Δω = Ω /J. Half of
the length interval associated with each of the two boundary nodes is
inside the domain [0, Ω]. To allow for imposing the above boundary
conditions, one ghost node is added at either side of the domain, with
indices j =1 and j =J+3, such that the plane ω =0 is
in ω 2 and the plane ω = Ω is inω J+2.
The equations are solved with an Euler-forward finite-difference scheme
implemented with MATLAB R2014b. To impose numerical stability in our
explicit scheme, the time step Δt should obey the criterion
where Fo Δ is a local Fourier number andCe 0 is a constant much larger
than the maximum value of the consolidation coefficientCe of Eq. (1.18), i.e.Ce 0 should be much larger than
the constant factor in Eq. (118).
The discretisation of Eqs. (1.16) and (1.17) is pretty straightforward.
The same applies to the boundary conditions, except that for,
calculating a value of e 2 at nodeω 2 (i.e., at ω =0) from Eqs. (1.17) and
(1.19), a value of e 1 is needed at ghost nodeω 1. It is found by extrapolating from thee 1 values at nodes ω 2,ω 3and ω 4 by using
equal ratios of differences between these nodes, given that thee 1 profile is found to be square root shaped.
Care must be taken that the e 1 values at nodesω 1 should not become negative. To realize the
Neumann boundary condition at ω = Ω by applying a central
differencing scheme to Eq. (1.20), the value ofe 1 at node J +3 is taken equal to that at
node J +1. More details can be found in Hazelhoff Heeres’ MSc
thesis [29].
The time dependent filter cake thickness is calculated with
in which the superscript i denotes the time step of a variable
and =1 unless j =2 or j =J +2: then it is ½.
The outflow velocity through the filter cloth then follows from