Fourth-order analysis of a diffusive lattice Boltzmann method for barrier coatings

TL;DR

This paper analyzes a diffusive lattice Boltzmann method for barrier coatings, identifying higher-order corrections and boundary condition issues, ultimately demonstrating high accuracy in simulations for most parameter ranges.

Contribution

It provides a fourth-order analysis of the lattice Boltzmann method, improves boundary conditions, and assesses the method's accuracy for simulating fluid transport in barrier coatings.

Findings

Higher-order corrections improve accuracy in certain parameter ranges.

Standard boundary conditions cause significant errors, which can be mitigated.

The method achieves below 0.1% error in most practical scenarios.

Abstract

We examine the applicability of diffusive lattice Boltzmann methods to simulate the fluid transport through barrier coatings, finding excellent agreement between simulations and analytical predictions for standard parameter choices. To examine more interesting non-Fickian behavior and multiple layers of different coatings, it becomes necessary to explore a wider range of parameters. However, such a range of parameters exposes deficiencies in such an implementation. To investigate these discrepancies, we examine the form of higher-order terms in the hydrodynamic limit of our lattice Boltzmann method. We identify these corrections to fourth order and validate these predictions with high accuracy. However, it is observed that the validated correction terms do not fully explain the bulk of observed error. This error was instead caused by the standard finite boundary conditions for the contact of the coating with the imposed environment. We identify a self-consistent form of these boundary conditions for which these errors are dramatically reduced. The instantaneous switching used as a boundary condition for the barrier problem proves demanding enough that any higher-order corrections meaningfully contribute for a small range of parameters. There is a large parameter space where the agreement between simulations and analytical predictions even in the second-order form are below 0.1%, making further improvements to the algorithm unnecessary for such an application.

Tables

Table 1. Table 1: Values of τ 𝜏 \tau and θ 𝜃 \theta used in simulations, with corresponding diffusion constant D 𝐷 D and time scale T 𝑇 T corresponding to four hours of macroscopic equivalent exposure with F = 5.76 × 10 − 2 𝐹 5.76 superscript 10 2 F=5.76\times 10^{-2} (all in lattice units).

$\tau$	$\theta$	$D$	$T$
0.55	0.5	0.025	23040
0.70	0.5	0.10	5760
1.0	0.5	0.25	2304
1.5	0.5	0.50	1152
2.0	0.5	0.75	768
10.0	0.5	4.75	121

Equations

$${\partial }_t\rho +\nabla \mathbf{j}=0,$$

$$f_i(x+v_i,t+1)=f_i(x,t)+\sum _j{\Lambda }_{ij}\left[f_j^0(\rho (x,t))-f_j(x,t)\right].$$

$$\rho (x,t)=\sum _i{f}_i(x,t).$$

$${\Lambda }_{ij}=\frac{1}{\tau }{\delta }_{ij}$$

$$\sum _i{f}_i^0$$

$$\sum _i{f}_i^0{v}_{i\alpha }$$

$$\sum _i{f}_i^0{v}_{i\alpha }{v}_{i\beta }$$

$${\partial }_t\rho ={\nabla }_{\alpha }\left(\tau -\frac{1}{2}\right){\nabla }_{\alpha }(\rho \theta )$$

$$D=\left(\tau -\frac{1}{2}\right)\theta .$$

$$f_i^0=\rho w_i,$$

$$w_0$$

$$w_1$$

$$w_2$$

$$f_i(0,t)=f_i^0({\rho }^b(t))$$

$${\rho }^{th,1}(x,t)={\rho }_0\left(1-\mathrm{erf}\left(\frac{x}{\sqrt{4Dt}}\right)\right).$$

$$\rho^{th}(x,t)=\rho_{0}\sum_{i=0}^{\infty}(-1)^{i}\biggl{[}2+\operatorname{erf}\left(\frac{x-2(i+1)L_{x}}{\sqrt{4Dt}}\right)-\\ \operatorname{erf}\left(\frac{x+2iL_{x}}{\sqrt{4Dt}}\right)\biggr{]}$$

$$t$$

$$x$$

$$\rho$$

$$F\equiv \frac{TD}{X^2}$$

$$ϵ(x)\equiv \left|\rho (x,t)-{\rho }^{th}(x,t)\right|$$

$$\tau=\tau(x)\equiv\left\{\begin{array}[]{llr}\tau_{\operatorname{out}}&,&0\leq x\leq L_{x}/2\\ \tau_{\operatorname{in}}&,&L_{x}/2\leq x\leq L_{x}\end{array}\right.$$

$$\begin{gathered} ({\partial }_t+v_{i\alpha }{\nabla }_{\alpha })f_i^0-\left(\tau -\frac{1}{2}\right)({\partial }_t+v_{i\alpha }{\nabla }_{\alpha }{)}^2{f}_i^0+\left({\tau }^2-\tau +\frac{1}{6}\right)({\partial }_t+v_{i\alpha }{\nabla }_{\alpha }{)}^3{f}_i^0 \\ -\left({\tau }^3-\frac{3}{2}{\tau }^2+\frac{7}{12}\tau -\frac{1}{24}\right)({\partial }_t+v_{i\alpha }{\nabla }_{\alpha }{)}^4{f}_i^0+O({\partial }^5)=\frac{1}{\tau }(f_i^0-f_i). \end{gathered}$$

$$\sum _i{f}_i^0$$

$$\sum _i{v}_{i\alpha }{f}_i^0$$

$$\sum _i{v}_{i\alpha }{v}_{i\beta }{f}_i^0$$

$$\sum _i{v}_{i\alpha }{v}_{i\beta }{v}_{i\gamma }{f}_i^0$$

$$\sum _i{v}_{i\alpha }{v}_{i\beta }{v}_{i\gamma }{v}_{i\delta }{f}_i^0$$

$$\partial_{t}\rho-A(\tau)(\partial_{t}^{2}\rho+\nabla_{\alpha}\nabla_{\beta}\rho\theta\delta_{\alpha\beta})+\\ B(\tau)(\partial_{t}^{3}\rho+\partial_{t}\nabla_{\alpha}\nabla_{\beta}\rho\theta\delta_{\alpha\beta}+\partial_{t}\nabla_{\alpha}\nabla_{\gamma}\rho\theta\delta_{\alpha\gamma}+\partial_{t}\nabla_{\beta}\nabla_{\gamma}\rho\theta\delta_{\beta\gamma})-\\ C(\tau)\bigg{(}\partial_{t}^{4}\rho+\partial_{t}^{2}\nabla_{\alpha}\nabla_{\beta}\rho\theta\delta_{\alpha\beta}+\partial_{t}^{2}\nabla_{\alpha}\nabla_{\gamma}\rho\theta\delta_{\alpha\gamma}+\partial_{t}^{2}\nabla_{\alpha}\nabla_{\delta}\rho\theta\delta_{\alpha\delta}+\partial_{t}^{2}\nabla_{\beta}\nabla_{\gamma}\rho\theta\delta_{\beta\gamma}\\ +\partial_{t}^{2}\nabla_{\beta}\nabla_{\delta}\rho\theta\delta_{\beta\delta}+\partial_{t}^{2}\nabla_{\gamma}\nabla_{\delta}\rho\theta\delta_{\gamma\delta}+\nabla_{\alpha}\nabla_{\beta}\nabla_{\gamma}\nabla_{\delta}\left[\frac{\rho\theta}{3}(\delta_{\alpha\beta}\delta_{\gamma\delta}+\delta_{\alpha\gamma}\delta_{\beta\delta}+\delta_{\alpha\delta}\delta_{\beta\gamma})\right]\bigg{)}+O(\partial^{5})=0$$

$$A(\tau )$$

$$B(\tau )$$

$$C(\tau )$$

$${\partial }_t\rho =\left(\tau -\frac{1}{2}\right){\nabla }_{\alpha }^2\rho \theta +O({\nabla }^3).$$

$${\partial }_t^2\rho ={\left(\tau -\frac{1}{2}\right)}^2{\nabla }_{\alpha }^2{\nabla }_{\beta }^2\rho {\theta }^2+O({\nabla }^5).$$

$$\begin{gathered} {\partial }_t\rho -\left(\tau -\frac{1}{2}\right){\nabla }^2\rho \theta -{\left(\tau -\frac{1}{2}\right)}^3{\nabla }^4\rho {\theta }^2+\left({\tau }^2-\tau +\frac{1}{6}\right)\left(\tau -\frac{1}{2}\right)3{\nabla }^4\rho {\theta }^2 \\ -\left({\tau }^3-\frac{3}{2}{\tau }^2+\frac{7}{12}\tau -\frac{1}{24}\right){\nabla }^4\rho \theta =0. \end{gathered}$$

$${\partial }_t\rho =D{\nabla }^2\rho +\alpha {\nabla }^4\rho$$

Full text

I Introduction

Coating systems are used heavily in industry for the protection of materials and infrastructure. Common examples include the paints on cars, bicycles, and houses; the layered coating systems used on boats and airframes; and the coatings used to protect bridges. In all cases, the goal of the coating system is to protect the underlying substrate from ingress by aggressive particulate, gaseous, or fluid materials while remaining aesthetically intact.

Crosslinked polymer networks, also called thermosets, are typically chosen in protective applications due to their net-like structure. In such a structure, precursor materials are chemically bonded through a crosslinking and curing process to form a three-dimensional structural network. This network acts as a physical and chemical barrier that attempts to prevent permeation by water, salt, particulate matter, and other environmental contaminants. Since substrates are often materials susceptible to corrosion, like aluminum or steel, it is essential that moisture not be permitted to reach the substrate in appreciable quantities.

As is known, most coatings permit, to some degree, moisture ingress Taylor2012169 . This can be due to imperfections in the preparation process Kroll201582 , the formation of void space during curing or cooling Zee201555 , or because of damage in service. To help detect coating formulations that may be unsuitable for use in the field, accelerated weathering testing is used to determine failure rates and modes in the lab. Much research has been devoted to the relationship between accelerated testing, comparable real-world testing, and service life, but no complete and predictive model exists that accurately correlates a coating’s performance in lab testing, performance in field testing, and failure modes or lifespan that is likely to occur in service doi:10.1021/bk-2002-0805.ch001 .

Since moisture entering the coating is conserved, the dynamics of fluid density $\rho$ obey the continuity equation

[TABLE]

where $\mathbf{j}$ is the mass current. Assuming an isotropic coating, mass current will be in the direction of negative density current. We denote the proportionality between the current and negative gradient by $D$, which in the simplest case is a constant. Later, we consider a more general $D(\rho)$. We therefore have $\mathbf{j}=-D\nabla\rho$. With this constitutive relation for the mass current, we recover the well-known diffusion equation.

Several methods exist to model idealized diffusion. Early work focused primarily on precise mathematical modeling and numerical solutions to boundary-matched differential equations governing diffusion crank1979mathematics ; de1998monitoring . Modern approaches include network connectivity models AIC:AIC690480104 , Monte Carlo simulations SAHIMI19912225 , and finite-element analysis Mu2007 for more complex structures like porous media where an effective diffusivity is desired. However, approaches dealing with pore structures may depend on the structure and porosity of the material in question, quantities that most often unknown a priori. Additionally, finite-element models tend to be computationally complex and often rely on commercial closed codes. As a whole, there is comparatively little known about the precise dynamics of diffusion through polymeric coatings Pathania2017149 .

Additionally, different approaches exist when multiple layers are considered. In the case of multiple hydrophobic barrier coatings, boundary-matching Fickian solutions can be used de1998monitoring and matched to experiment via electrochemical methods. When a base coating is hydrophilic, as is the case with some primers, an alternative approach couples Fickian diffusion for any overlying barrier coatings with the assumption of an instantaneous reservoir for the base layer Baukh20123304 . In either case, different coatings in a multi-layer stackup differ in their effective diffusivity.

Any numerical technique used to model the progression of moisture in such a stackup must stably account for a wide range of diffusion constants. Since laboratory testing of candidate barrier coating systems typically includes cyclic exposure to moisture and dry ambient air over long periods of time, simulations of cyclic processes must maintain numerical stability over correspondingly longer time scales.

In this paper, we use lattice Boltzmann numerical techniques to determine the accuracy of modeling moisture ingress through a finite coating system exposed to a reservoir and adhered to an ideal substrate. Because of the necessity of modeling a wide range of saturation levels and diffusivity in the case of multi-layer systems, we analyze the error introduced in the traditional second-order approximation to the diffusion equation used in lattice Boltzmann approaches. To investigate the nature of this error, we introduce a fourth-order correction and perform a Fourier component analysis to confirm the correctness of our results. We show that the bulk of the second-order error in such a system arises from the boundary conditions used, and comment on the proper use of periodic systems to remove this error. Applications to multi-layer systems with variable diffusivity are discussed in the context of our analysis.

II Lattice Boltzmann methods

The lattice Boltzmann approach models densities $f_{i}$ defined on a discrete lattice space associated with discrete lattice velocities $v_{i}$. After being displaced to a new lattice position $x+v_{i}$, the densities at each lattice point are redistributed in a collision step. This method has been used extensively to model hydrodynamic behavior qian1992lattice ; PhysRevLett.75.830 ; PhysRevLett.56.1505 , diffusion wolfdiffusion ; shan , electrostatics electrostatics , and similar systems with high accuracy and computational efficiency. Notably, the hydrodynamic partial differential equations underlying such systems are not the starting point for the method, but rather emerge from it. Choices like the number of quantities conserved in the collision allow for the freedom of recovering the governing equations for a variety of different systems, as mentioned above.

A popular collision term defines a local equilibrium $f_{i}^{0}$ that only depends on the conserved quantities and then relaxes the actual density towards the local equilibrium. In this form the lattice Boltzmann equation can be written as

[TABLE]

Here $f_{j}^{0}$ is the local equilibrium density, $\Lambda_{ij}$ is a collision matrix, and $\rho(x,t)$ is the local density of the system, given by

[TABLE]

The form of the collision matrix allows for further control of the algorithm, but this freedom is not explored in this paper. Most examples where this freedom has shown to be useful relate to simulations of hydrodynamic systems with very low viscosity. Such low viscosities may give rise to instabilities that can be controlled by a careful choice of the collision matrix. For diffusive systems like the one considered here, the advantages of utilizing multiple relaxation times are less well established (See Ginzburg PREGinzburg ), so we will employ the particularly simple collision matrix

[TABLE]

that was originally proposed by Qian qian1992lattice , using a single relaxation time $\tau$.

It is necessary to impose moments on the equilibrium distribution, following the method of wolfdiffusion . While not considered here, formulations of this method in the case of multiple components shan and multiple relaxation times LiPRE2014 ; LePRE2015 . In particular, we impose the following (non-unique) moments on the distribution:

[TABLE]

where the Greek indices are spatial dimensions and follow the Einstein notation.

Local density conservation is assured by Eqn. (5), while Eqn. (7) introduces a spatially uniform imposed temperature $\theta$. Following PhysRevE.94.033302 , a second-order Taylor approximation using this choice of moments leads to the lattice diffusion equation

[TABLE]

and, if the temperature is constant, this recovers a diffusion equation with the diffusion constant

[TABLE]

Testing for coating applications usually applies moisture somewhat homogeneously on the sample, either in soak testing or weathering chambers; drying also proceeds evenly. This reduces the problem of interest to an effectively one-dimensional case. For simulations, we use the simplest one-dimensional lattice Boltzmann model with the velocities $\{v_{i}\}=\left\{0,1,-1\right\}$. This one-dimensional lattice with the given velocities is known as a D1Q3 scheme. For this implementation of a diffusive system, the local equilibrium distribution can be written as

[TABLE]

where $w_{i}$ are the weights related to the magnitude of the velocities $\{v_{i}\}$. To recover the necessary moments, the weights are

[TABLE]

The D1Q3 implementation then allows for a full and self-contained simulation method for a diffusive system.

III Application to water content of coatings

We wish to model the wetting of a single-layer coating via Fickian diffusion. Since coatings are frequently examined in the laboratory on test panels using weathering chambers that subject the coating to moisture, we will consider the case where the coating, represented by a lattice from $0\leq x\leq L_{x}$, is exposed to a reservoir of varying concentration $\rho^{b}(t)$ at $x=0$ and an impermeable substrate at the right end of the simulation lattice. The meaning of $\rho^{b}$ is the amount of water that will be absorbed just inside the coating as it is exposed to the environment. For an immersion in water, this corresponds to the maximal water content the coating can absorb, and we scale the density so that this value corresponds to $\rho=1$.

We must account for these two boundary conditions in our numerical simulation. We implement the source term by setting

[TABLE]

and by replacing the streaming step at the right end by a bounceback algorithm, where the right-moving $f_{1}(L_{x})$ is reinserted as an $f_{2}$ in the streaming step. The result for a step function $\rho^{b}(t)=\Theta(t)$ in the exposure is shown in Fig. 1. We used a system with $L_{x}=100$ lattice points, $\tau=1$, $\theta=0.5$, $\rho_{0}=1$, and ran the simulation for a variable number of iterations $T$. As expected, moisture is at first located closely to the surface and then penetrates the sample.

To verify the correctness of the simulation results, we construct an analytical solution for the concentration over time, using linear combinations of the well-known error function solution crank1979mathematics . These are solutions of the diffusion equation for the initial condition of a step function in an infinite system. If the initial step goes from $2\rho_{0}$ to zero, then the solution is

[TABLE]

This solution has a fixed point $\rho_{\operatorname{bath}}(0,t)=\rho_{0}$ at $x=0$, which corresponds to our boundary condition. So $\rho^{th,1}(x,t)$ for $x\geq 0$ and $t>0$ is the analytical solution for an infinite dry coating exposed to a reservoir starting at time $t=0$. Note that the long-time behavior gives $\rho^{th,1}(x,t\to\infty)=\rho_{0}$ as expected.

Suppose now that we have a finite one-dimensional coating extending from $0\leq x\leq L_{x}$. At $x=0$, the coating is exposed to a reservoir with fixed concentration $\rho(x=0,t)=\rho_{0}$. At $x=L_{x}$ is an impermeable substrate where $\nabla\rho(x=L_{x},t)=0$.

To account for the vanishing current at the substrate, we use an image source reservoir at $x=2L_{x}$. This will ensure a vanishing gradient at $x=L_{x}$ and, by symmetry, a vanishing current. However, when the reflected concentration becomes nonzero at the reservoir again, we must subtract another image source reservoir at $x=-2L_{x}$ to maintain the correct boundary condition. Repeating this process infinitely, we arrive at the final solution that includes both reservoir and substrate:

[TABLE]

For practical purposes, we find ten terms of the infinite sum in Eqn. (14) are entirely sufficient for most cases.

It is instructive to determine the correspondence between these numerical parameters and a laboratory case. A typical barrier coating might have thickness $X=50\ \mu$m, diffusion constant in water $D\sim 10^{-14}\ \text{m}^{2}/\text{s}$, and be exposed to moisture in a weathering chamber for $T=4$ hours at a time for testing. We can introduce reduced time, length, and density scales $t^{\prime},x^{\prime},\rho^{\prime}$ such that

[TABLE]

and $0\leq\{t^{\prime},x^{\prime},\rho^{\prime}\}\leq 1$. Since the unit relationship $T=X^{2}/D$ holds by dimensional analysis, for any given experimental setup the quantity

[TABLE]

is dimensionless and we have the scaled diffusion equation $\partial_{t^{\prime}}\rho^{\prime}=-\nabla_{x^{\prime}}F\nabla\rho^{\prime}$. Using the experimental parameters suggested above gives $F=5.76\times 10^{-2}$. In our simulations, we use total length $X=L_{x}=100$ lattice sites, reservoir concentration $\rho_{0}=1$, $\theta=0.5$, and $\tau=1$. Since this gives a time scale $T\approx 2300$ iterations, this means one hour of equivalent macroscopic exposure corresponds to approximately $575$ simulation iterations. Further, the choice of $\tau=1$ yields immediate relaxation of local distributions, so we would expect excellent agreement to theory.

We are now in a position to comment on the accuracy of this simulation method in comparison to the analytical solution $\rho^{th}$ in Eqn. (14). For each of the exposure times in Fig. 1, we compute the absolute error

[TABLE]

across the lattice space profile. The result is plotted logarithmically in Fig. 2, showing excellent agreement. It is interesting to observe how the error changes over time; initially, the error drops substantially since moisture has not yet permeated through the entire coating lattice. This tail increases as the entire lattice becomes wet, but then uniformly decays as the numerical solution approaches saturation and agrees with the corresponding analytical solution.

While this method provides efficient and stable numerical modeling of a single coating, a given coating system might consist of two or more barrier layers in a stackup, each with a different diffusion constant that permits moisture ingress and egress at different rates from its neighbors. To extend this method to the simplest multi-layer case, we might wish to model a two-layer stackup consisting of idealized barrier coatings with different physical properties. To do so, our reservoir model is modified slightly, with the outer barrier coating represented at lattice sites $0\leq x\leq L_{x}/2$ and the inner barrier coating at $L_{x}/2\leq x\leq L_{x}$. Since the diffusion constant is controlled by the parameter $\tau$ in Eqn. (8), the presence of two diffusion constants requires that $\tau$ be position-dependent:

[TABLE]

Incidentally, changing the value of $\theta$ between the two regions will lead to different maximum water uptake in the layers, an important relationship that will be explored elsewhere.

Although such a two-layer system is not investigated in this paper, it is essential to determine the range of $\tau$ values for which numerical and analytical solutions agree sufficiently over time. For efficient simulations, it is advantageous to choose $\tau$ as large as feasible, since this corresponds to a large diffusion constant and hence a shorter simulation time. For a quick initial evaluation, we run a series of lattice Boltzmann simulations with varying values of $\tau$ to the same macroscopic equivalent time of four hours of exposure. After that time, we compute the absolute error $\epsilon$ between numerical and analytical solutions across the entire lattice profile. Results are shown in Fig. 3. The choices of $\tau$, along with the corresponding time scale $T$, are shown in Table 1.

As shown earlier, the solutions agree very well for $\tau=1$. However, the error may be orders of magnitude larger for $\tau\neq 1$. Depending on the particular application, we may require ratios of diffusion constants that vary significantly (such as in multi-layer systems); however, the errors indicated here may cause the numerical method to appear less than ideal. We discuss the lower asymptotic limit $\tau\to 0.5$ later.

We therefore wish to examine the origin and nature of the $\tau$-dependent error. Of note is that the derivation of the lattice diffusion equation given above (and used heavily in the literature) is done with only a second-order Taylor approximation. To determine the degree to which this approximation leads to the errors shown, we next perform a fourth-order correction to this diffusion equation.

IV Fourth-order limit of diffusion equation

In order to introduce a correction to the diffusion equation, we perform a Taylor expansion of the lattice Eqn. (2) to account for higher orders. As shown by Wagner PhysRevE.74.056703 , this equation expanded to the fourth order takes the form

[TABLE]

Since we have now introduced higher-order powers into this expansion, we must utilize moments up to the fourth-order. Using the form of the equilibrium distribution in Eqn. (10), we calculate the higher-order moments:

[TABLE]

Summing over all indices of Eqn. (21) using these revised moments, we are left with

[TABLE]

where we have defined the $\tau$-dependent prefactors

[TABLE]

for brevity.

This form is not particularly useful since there are mixed spatial and temporal derivatives in the higher-order powers. In our one-dimensional implementation, we can drop our indices. We use the diffusion equation to write the temporal derivatives in terms of the spatial derivatives as

[TABLE]

It immediately follows that

[TABLE]

We can then introduce these two substitutions into Eqn. (27) and we have

[TABLE]

We then obtain the form of a corrected diffusion equation

[TABLE]

with corrections up to the fourth power in spatial derivatives, where we define

[TABLE]

This definition of $\alpha$ represents the expected error between the second-order diffusion equation and the corrected fourth-order equation. For certain parameter values, such as $\tau=1$ and $\theta=1/3$, we have $\alpha=0$, which accounts for higher accuracy observed for such parameters. We plot a density field representation of the relative error quantity $\alpha(\tau,\theta)/D(\tau\theta)$ in Fig. 4. We indicate a contour where this quantity vanishes, as well as additional contours whose numerical importance will be explained in later sections.

The correction term has a similar form to a surface tension term in a Cahn-Hilliard equation. In this case, positive values of $\alpha$ would correspond to a negative surface free energy. This implies that simulations with positive $\alpha$ should be unstable for high frequency perturbations. This equation can be solved in Fourier space, allowing us to verify our analytical predictions with lattice Boltzmann simulations. In the subsequent section, we perform this analysis.

V Fourier analysis of correction term

A Fourier transform of Eqn. (31) yields

[TABLE]

Here $k$ is any specific Fourier mode and $\widehat{\rho}(k,t)$ is the $k$-space density represented by

[TABLE]

where $L_{x}$ is the system size in the $x$-direction. Even though $x$ is continuous, the finite periodicity of $2\pi$ causes $k$ to be discrete. This allows for our system to contain a finite number of $k$ modes which can be now examined independently. The form of Eqn. (33) is simple since different $k$ modes do not couple. In $k$-space, the initial profile at $t=0$ is chosen by defining $\rho(x,0)$, which for $\widehat{\rho}(k,0)$ gives Eqn. (34) and

[TABLE]

We reproduce the uncorrected diffusion equation by setting $\alpha=0$, obtaining

[TABLE]

These predictions are implemented on a discrete lattice which implies that there will be a finite number of $k$ modes. From Eqn. (34), we have

[TABLE]

which implies a maximum allowed $k$ mode when $k=\pi$ and a minimum lattice dimension of $L_{x}=2$. In this finite system, we have the back transform

[TABLE]

It is now possible to verify this theoretical prediction by examining the decay of specific Fourier modes by imposing an initial profile

[TABLE]

Using this profile, the uncorrected and corrected $k$-space densities become, respectively,

[TABLE]

In practice, we change $k$ by varying the system size $L_{x}$. An interesting point to note is that when $\alpha<-\frac{D}{\pi^{2}}$, it is predicted that the numerical simulations would be unstable. This is predicted due to the fact that in Eqn. (40), the negative $\alpha$ term leads to a positive exponent and causes $\widehat{\rho}(k,t,\alpha)$ not to decay.

VI Numerical verification of correction term

To determine the validity of the prediction for the correction term shown in Eqn. (32), we define a ratio between the two forms of k-space density in Eqn. (40) in a simple form such that

[TABLE]

We can use this relation to measure $\alpha$ from numerical simulations. We do this by initializing our probability distributions by $f_{i}(x,0)=f_{i}^{0}(\sin(2\pi x/L_{x}))$ and then varying $L_{x}$. Our first prediction is that $\ln(R(k,t))$ is a linear function of $t$. We can find $\alpha$ from the time evolution of the density through

[TABLE]

where we numerically calculate the temporal derivative using a finite difference method.

The numerical evaluation of Eqn. (41) using the numerical results is shown in Fig. 5. At $t=0$ we have $R=1$ by construction, but for all $\tau\neq 1$ we observe a rapid transient change which manifests itself as a near instantaneous jump in Fig. 5. After this transient period, the behavior of $\ln(R)$ is indeed linear, as expected. We then calculate

[TABLE]

where we take $t_{1}$ when $\widehat{\rho}(k,t,\alpha)=0.5$ and $t_{2}$ when $\widehat{\rho}(k,t,\alpha)=0.01$ to avoid any difficulties with the offset. Eqn. (42) gives our correction polynomial as a function of any Fourier mode $k$. Using this form, we can compare our predicted correction term in Eqn. (32) to a numerical representation. Fig. 6 shows simulation data for $\frac{d}{dt}\left[\ln R(k,t)\right]$ for $\tau=1$ and $\theta=0.1$. We see a good fit for all $k$ modes between simulation and the prediction in Eqn. (42).

We first test the prediction comparing $\alpha_{\operatorname{exp}}$ in Eqn. (42) to our theoretical prediction for $\alpha$ from Eqn. (32). Fig. 7 shows a comparison between $\alpha_{exp}$ and our theoretical prediction for $\alpha$ for various values of $\tau$ and $\theta$ as a function of $L_{x}$. For this analysis, we chose a known stable value for either $\tau$ or $\theta$ and set the other parameter as a more extreme value. For a choice of $\theta=1/3$, we set $\tau=0.51$ as the extreme value. In these cases, we see very good agreement between $\alpha_{\operatorname{exp}}$ and our prediction. In the cases of $\theta=1/3$ with $\tau=1.5$ and $\theta=0.9$ and $\tau=1$ we observe good agreement for $L_{x}>40$, but as $L_{x}$ becomes smaller, deviations begin to increase. This suggests that there is a discrepancy in $\alpha_{\operatorname{exp}}$ for large $k$ modes.

In the case where $\alpha=0$, it is interesting to note that the results match a $\frac{1}{k^{6}}$ rather than the predicted $\frac{1}{k^{4}}$ fit. This implies that there are additional correction terms which may be relevant at specific values of $\tau$ and $\theta$. These higher-order corrections are not considered in the present analysis.

As discussed previously, Eqn. (40) predicts numerical instability when $\alpha<-\frac{D}{\pi^{2}}$. The density representation shown in Fig. 4 implies that this will happen as we increase $\tau$ and decrease $\theta$ to extreme values ($\tau\gtrsim 4$ and $\theta\lesssim 0.3$ simultaneously). A contour showing $\alpha(\tau,\theta)/D(\tau,\theta)=-1/\pi^{2}$, the start of the region of instability, is shown in that figure.

It is instructive to examine $\alpha$ while holding either $\tau$ or $\theta$ fixed. Setting $\theta=\frac{1}{3}$, we examine $\alpha$ as a function of $\tau$ alone in Fig. 8, which shows excellent agreement to theory over 100 independent $k$ modes. We set $\tau=1$ and examine $\alpha$ as a function of $\theta$ alone in Fig. 9, with similarly excellent agreement.

VII Application of correction to reservoir diffusion

With the fourth-order correction term in hand and its correctness assured, we next determine its applicability to our reservoir coating system. Fig. 10 shows the absolute error profile between lattice Boltzmann simulation results and a fourth-order corrected analytical solution. This solution is produced by first setting up an appropriate initial step function

[TABLE]

in a periodic lattice. This is entirely equivalent to the boundary conditions implied by the derivation of the second-order error function solution in Eqn. (14). We transform this step function into $k$ space via a discrete Fourier transform, use the fourth-order correction to perform a time evolution, and then transform the result back into real space. Strictly speaking, the method of Eqn. (14) generates a continuous solution, while the Fourier transform approach yields a discrete solution. We discuss the ramifications of this difference in the Appendix and conclude that the difference in solution discretization is very small and of the same order as the error produced in our best numerical results.

We note no consistent improvement over the second-order shown in Fig. 3 from the introduction of the fourth-order correction. However, the magnitude of the error, especially for high values of $\tau$, remains troubling. The nature of the finite simulation lattice is such that the boundaries are treated independently of other lattice sites. In particular, the reservoir density $\rho$ is set manually and not strictly determined by local distributions. Since we have seen that $\tau$-dependent errors tend to accumulate near the reservoir boundary over an order of magnitude higher than at the substrate boundary, the nature of using such a finite lattice is suspect. The case when $\tau=1$ yielded excellent agreement throughout the finite lattice, but this is consistent with the immediate relaxation of local equilibrium distributions and does not apply to other values of $\tau$.

This $\tau$-dependent error is consistent with the jump observed in Fig. 5, where setting $f_{i}(x,0)=f_{i}^{0}(\rho(x))$ led to deviations. Indeed, Eqn. (21) implies that

[TABLE]

which suggests an approach that would allow us to increase the accuracy of our boundary conditions.

In our current case, however, we can avoid the cumbersome issue of the boundary condition altogether by simply embedding the system into the periodic lattice used for establishing the initial step function condition of the analytical Fourier solution. This permits a more standard lattice Boltzmann approach that does not rely on manual density adjustment at the reservoir (here at $x=3L_{x}$) and uses symmetry to establish the substrate (at $x=4L_{x}$) with no bounceback. With this setup, we have removed the need for boundary conditions altogether; we therefore expect that the $\tau$-dependent error should be substantially reduced, especially at the reservoir boundary. Fig. 11 shows a diagram of the periodic step function from Eqn. 44 used for this analysis.

We again run two sets of simulations for our range of $\tau$ values to the same scaled time, both using the periodically-embedded lattice simulation. The first set of simulations uses only the traditional second-order approximation and is shown in Fig. 12. The second set applies our fourth-order correction and is shown in Fig. 13.

As was hoped, the error at the reservoir is reduced by orders of magnitude when compared to the finite system with imposed boundaries. This confirms that the accumulated error from the finite system is due to the presence of boundary conditions that are only guaranteed to match at $\tau=1$ when relaxation is immediate during collisions. However, contrary to expectation, there is almost no benefit from the fourth-order $\alpha(\tau,\theta)$ correction, even though its validity was verified via Fourier analysis.

It is of value to mention here the effects of using a progressively lower of $\tau$, since our results suggest that errors increase with $\tau$, regardless of the choice of boundary conditions considered. Of course, a lower choice of this parameter leads to increasingly long simulations, which must be balanced with the desired numerical accuracy. We examined a range of $\tau$ values as low as $\tau=0.5001$ and found essentially no change in error from the $\tau=0.55$ lower limit presented throughout this paper. This suggests that (to within machine accuracy), there is likely no theoretical limit to the ratio of diffusion constants possible. This implies that a study of a multi-layer coating stack, which fixes the diffusitivity ratio through a selection of $\tau$ values, is possible for a situation where one coating’s diffusion constant is orders of magnitude higher than the other.

It is natural at this point to wonder if there are any choices of parameters $\tau$ and $\theta$ for which the fourth-order correction provides substantial benefit in our reservoir problem, especially since its use in simulations incurs additional computational burden. Naturally, any such error analysis depends heavily on the particular problem of interest, and therefore on the initial profile and desired time evolution. For our system, we examine the parameter space $0.5<\tau\leq 2.5$ and $0.1\leq\theta\leq 1.0$. For each point in this space, we run a lattice Boltzmann periodic reservoir system simulation to the same scaled time. After this time, we compute the ratio

[TABLE]

where

[TABLE]

is the root mean square error between numerical concentration $\rho$ and second- or fourth-order Fourier analytical concentration $\rho_{2,4}$. If $\overline{\epsilon}\approx 1$, there is no appreciable correction from using the fourth-order solution; as $\overline{\epsilon}\to 0$, the correction becomes more substantial. From a computational perspective, there is a trade-off between the computational burden of the correction and the benefit (if any) from using it. We do not comment on the appropriate balance for any particular situation.

Contrary to expectation, there are no regions of the given parameter space where $\overline{\epsilon}<0.9$ during long times, indicating no appreciable benefit to the correction. Further, the fourth-order analysis predicts numerical instability in the bulk region of parameter space where $\alpha<-1/\pi^{2}$, although the numerical simulations and second-order analysis remain stable. This is a surprising result overall: a fourth-order correction is not only unhelpful in increasing the accuracy of solutions at long times, it is often worse than the second-order approximation and predicts numerical problems incorrectly.

If we instead run the same analysis for a much shorter time (in the equivalent macroscopic system, just $3.5$ seconds), the results are more promising and shown in Fig. 14. For small values of both $\tau$ and $\theta$, the fourth-order correction increases accuracy by an order of magnitude. This is largely due to the fact that the fourth-order theory accurately predicts some early time oscillations at the sharp reservoir interface, as shown by the error reduction in Fig. 15. This discussion of higher-order effects gives the rather surprising result that for our barrier coating application, there is no noticeable improvement. This may also arise because even the second-order results are accurate enough that any resulting errors are of the same order of magnitude as the difference between continuous and discrete analytical solutions, as shown in Fig. 17 in the Appendix.

There is another prediction we can obtain from our derivation of the correction term. In Figure 14, we see bold black lines indicating the domain for $\alpha(\tau,\theta)=0$. For these values, our theory predicts that the original method (derived only to second order) is actually a fourth-order method. To test this prediction, we choose $\tau=1$ and $\theta=1/3$ for which we have $\alpha=0$, and examine the convergence of the method to the Fourier term analytical solution $\epsilon_{2}$ of Eqn. (47). In particular, we examine the same periodic system as before by choosing $F=0.0576$ for different lattice sizes. Keeping $F$ constant implies that the number of iterations scales as the square of the lattice size, which is sometimes known as diffusive scaling. This is shown in Figure 16, where we see that we indeed find fourth-order convergence when $\alpha=0$. For pairs of $\tau$ and $\theta$ for which we have nonzero correction terms (e.g. we show the case $\tau=0.55$ and $\theta=1/3$), we see that we have a second-order convergence instead.

VIII Conclusions

In this paper we have examined whether a diffusive lattice Boltzmann method is an effective tool for examining problems related to Fickian water diffusion in barrier coatings. This validation was assisted by our ability to derive an analytical solution for a simple, but not trivial, coatings problem. In Sec. III we presented a real-space solution for the water content of a dry coating that is initially exposed to a constant moisture reservoir on the surface. A second analytical solution in terms of Fourier components was presented in Sec. V that can be used both for the standard Fickian diffusion case already examined in Sec. III, as well as the more complex fourth-order diffusion equation we derived as part of a higher-order hydrodynamic limit of the lattice Boltzmann equation. The two equivalent analytical solutions differ slightly because our Fourier series corresponds to a discrete system with only a finite number of Fourier terms.

For a simple initial implementation of the inlet boundary, we found excellent agreement only for a relaxation time $\tau=1$. Our analysis revealed that the disagreement for $\tau\neq 1$ was caused by assuming an equilibrium distribution as the reservoir boundary condition. Eventually we were able to define a “perfect” boundary condition by doing away with the boundary altogether through an embedding of the system in a large periodic system that only requires periodic boundary conditions.

Along the way of our examination, we discovered that we can indeed identify a fourth-order accurate hydrodynamic limit of the diffusion equation. However, this higher-order correction was found to be irrelevant for the coatings problem considered here, as we could only identify a small region in parameter space where the fourth-order predictions were significantly more accurate. This may act as a cautionary tale that validating a higher-order correction does not guarantee that such predictions will always be more accurate for specific applications.

However, for the best cases, the numerical solutions agree with our analytical solutions almost as well as the two analytical solutions agree with each other, suggesting that the proposed method is indeed an excellent candidate to be applied to coatings problems.

In the future we expect to extend this from one-dimensional coatings problems, corresponding to full immersion of the coating, to the more complicated problem of droplets sitting on a coating. This case will require a full three-dimensional simulation, and wetting and drying problems then occur in one and the same simulation, spatially separated. Furthermore a droplet sitting on a coating would add a pressure gradient caused by the Laplace pressure in the drop. This may lead to an additional transport mode of advection, driven by the pressure gradient. This will require an extension of the current model to allow for some amount of advection as well. This advection would be expected to be highly overdamped, so that the local advection velocity would be simply proportional to the local pressure gradient. Technically doing this will require a replacement of the equilibrium distribution to one which allows for a non-zero first moment. Such simulations will significantly extend the current state of the art for coatings research which remains firmly focused on one-dimensional problems.

Acknowledgments

The authors thank Prof. Stuart Croll, Kent Ridl, and Reza Parsa of North Dakota State University for ongoing discussions and helpful insights. The second author was supported by the Strategic Environmental Research and Development Program under contract W912HQ-15-C-0012. Views, opinions, and/or findings contained in this report are those of the authors and should not be construed as an official Department of Defense position or decision unless so designated by other official documentation.

Appendix A Continuous and discrete solutions

The error function solution in Eqn. (14) solves the second-order diffusion equation with the given boundary conditions in continuous real space (see crank1979mathematics for a basic form of the derivation). However, we later compute a solution by transforming the appropriate initial condition into Fourier space, performing a second-order time evolution, and then transforming back into real space. This process uses a finite number of $k$ modes in each transform, and necessarily implies a discrete lattice sampling of both the initial condition in real space and the time-evolved form in Fourier space. We therefore expect a discrepancy when directly comparing the two solutions: the first is a solution to the continuous diffusion equation that is examined at discrete lattice points for comparison to the simulation, while the second is a sampled solution to the discrete lattice diffusion equation, the continuous form of which would require (in theory) an infinite number of $k$ modes to match the continuous case.

To examine the extent to which these solution forms differ from each other, we compute both at the same scaled four-hour time at each lattice site, and plot the absolute value of the difference, $\epsilon$, in Fig. 17. The two solutions agree to within $10^{-5}$ of each other. Since this error is on the order of the remaining error for the periodically-embedded simulation, we conclude that any further correction of simulation results renders any error obscured by differences between the discrete and continuous solutions to the diffusion equation, and is of no practical consequence.