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