pp. X–XX SIMULATING BINARY FLUID-SURFACTANT DYNAMICS BY A PHASE FIELD MODEL

In this paper, the dynamics of a binary fluid-surfactant system 
described by a phenomenological phase field model is investigated 
through analytical and numerical computations. We first consider 
the case of one-dimensional planar interface and prove the 
existence of the equilibrium solution. Then we derive the 
analytical equilibrium solution for the order parameter and the 
surfactant concentration in a particular case. The results show 
that the present phase field formulation qualitatively mimics the 
surfactant adsorption on the binary fluid interfaces. We further 
study the time-dependent solutions of the system by numerical 
computations based on the pseudospectral Fourier computational 
framework. The present numerical results are in a good agreement 
with the previous theoretical study in the way that the surfactant 
favors the creation of interfaces and also stabilizes the 
formation of phase regions.

1. Introduction.Due to the amphiphilic nature, surfactant molecules in binary fluid tend to migrate to fluid interfaces, and consequently alter the interfacial tension causing significant changes on the interfacial properties of the fluid mixture.Such a delicate property has many applications, for instance, facilitating the breakup of large droplets into smaller ones [6,10], preventing the coalescence of smaller droplets [10], and reducing the risk from bubbles formed in blood due to rapid decompression [1].
Since 1990's, a number of research studies [11,12,13,14,18] related to the surfactant in mixture have been conducted.Investigating such a complex system requires modelling the evolution behavior of the binary fluid, the mass transportation of surfactant in fluids, and the interactions between surfactant and binary fluid interfaces.Generally, the fluid mixture is described by the Cahn-Hilliard type of energy [2,17] while the behavior of surfactant in fluids can be modelled either by a microscopic or macroscopic approach.In [11,13], cell dynamical systems are used to account for the positions and the amphiphilic nature of surfactant molecules in fluids.The macroscopic approach of modelling the role of surfactant in fluids is by specifying a mixing energy law to provide a mechanism of the adsorption of surfactant on binary fluid interfaces such as in [12,19].
In this paper, we investigate the equilibrium and dynamical behavior of the binary fluid-surfactant system through a phase field model.The present formulation is adopted from the model introduced by Perkins et.al. [7], with a slight modification.The resulting energy law consists of three parts: (1) the Cahn-Hilliard free energy [2,17] describing the behavior of binary fluid, (2) a surfactant-interface coupling energy adopted from the energy law in [7], and (3) an entropy term specifying the ideal mixing manner of surfactant in fluids.Both analytical and numerical approaches are employed to study the solutions of the system.
In the theoretical part, we first consider the case of one-dimensional planar interface and prove the existence of the equilibrium solution.Then we derive the analytical equilibrium solution for the order parameter and the surfactant concentration in a particular case.The effect on the interfacial tension due to the presence of surfactant is also discussed for this particular solution.The results show that the present phase field formulation qualitatively mimics the surfactant adsorption on the binary fluid interfaces observed in experiments.
Numerical methods are also employed to investigate the time evolution of the system formulated by a system of partial differential equations consisting of a Cahn-Hilliard type of equation for the order parameter and a Fickian-diffusion type of equation for the surfactant.These time-dependent equations are discretized by the pseudospectral Fourier method in space and the backward Euler method in time.It is shown that the scheme preserves the mass conservation of the discrete solutions, and the discrete energy of the scheme decreases as time evolves.The present numerical results are in a good agreement with the previous theoretical study [7] in a way that the surfactant favors the creation of interfaces and also stabilizes the formation of phase regions.
The rest of the paper is organized as follows.In Section 2, we present the formulation of the proposed phase field model followed by the analysis of the equilibrium solution for the one-dimensional planar interface with surfactant.We then describe the time evolutional equations to study the dynamics of the binary fluid-surfactant system.Section 3 describes the numerical method for the time-dependent equations in detail.Numerical results and discussions are provided in Section 4 followed by some concluding remarks in Section 5.

2.1.
Free energy of the system.Consider a binary fluid-surfactant system in a domain Ω ⊆ R 2 or R 3 formulated by a phase field model.Let u : Ω → R be an order parameter where u = 1 and u = 0 correspond to different fluid phases, and ρ : Ω → R be the surfactant concentration.The total free energy of the phase field system is defined as where ε, α, and β are all small positive parameters.To have the equally competing effect, here we choose both α and β to be O(ε).Notice that, the present energy model is motivated by [7] but differs from that by adding the last energy term in Eq. (1).We now briefly describe how those energy terms make the contributions to the whole system.
Cahn-Hilliard free energy [2]: f (u)/ε + ε|∇u| 2 /2.The double-well potential f (u)/ε, termed as the bulk energy density, takes on the minimum value (f = 0) at u = 0 or u = 1, which drives the system to two phases.This effect leads to the separation of fluid domains into pure phases, and thus creates diffusive interfaces of width characterized by ε, connecting different phase domains.The interfacial energy is characterized by the gradient term, ε|∇u| 2 .Minimizing the gradient energy density, thus, minimizes the total perimeter of interfaces.These two combined energetic mechanism leads to the coarsening of separated phase domains.
Entropy for ideal mixing of surfactant: The energy density h(ρ) is an entropy term and the purpose of this term is twofold.The first one is to model the ideal mixing of surfactant in fluids which induces a Fickian type of equation for ρ by choosing the mobility M ρ ∝ ρ(1 − ρ) [16], see the time evolutional equation (22b) in detail.The second purpose is to restrict the value of ρ to be in the range (0, 1); ρ ln ρ ensures the value of ρ to be positive and (1 − ρ) ln(1 − ρ) enforces ρ < 1.Notice that the upper bound for ρ also models the saturation of surfactant at interfaces.This entropic part of free energy is also adopted in literature such as in [19,4] to study the kinetics of surfactant adsorption.
Surfactant-interface coupling energy: (ρ − |∇u|) 2 .In order to describe the adsorption behavior of the surfactant near the binary fluid interface, we must add a coupling energy in the whole energy system.This surfactant-interface coupling term is often constructed as −ρ|∇u| 2 in previous literature such as [12,18,19].At an interface where |∇u| takes on a larger value, the more surfactant is adsorbed on the interface, so the more the energy decreases.However, on a pure phase region where |∇u| vanishes, this coupling term does not imply whether the surfactant should move or not since no matter what value of ρ takes this surfactant-interface coupling term vanished.Thus, an additional energy term needs to be introduced to the total energy to prevent surfactant clustering in pure phase regions [12,18].In [7], Fonseca, Morini, and Slastikov added the present surfactant energy into the Cahn-Hilliard free energy and studied the role of surfactant in the formation of bubbles in foam.Through a sophisticated mathematical analysis, it was qualitatively shown that the surfactant segregates to the interfaces and enhances the creation of fluid interfaces.Here, we shall quantitatively demonstrate the above surfactant segregation behavior by numerical simulations.
2.2.Equilibrium solution.When the system is in equilibrium, the chemical potentials defined as the functional derivatives of energy G with respect to u and ρ become some constants.From the calculus of variations, those chemical potentials denoted by δG/δu and δG/δρ can be written as where denotes the differentiation with respect to the argument.
Let us assume the interface to be one dimensional planar interface (in yz plane) with normal along the x direction so that u = 0 and u = 1 while x tends to negative and positive infinity, respectively.We seek the solution profiles of u and ρ, with u x ≥ 0 so the equations ( 2)-(3) can be simplified as where the subscript x denotes the derivative with respect to x, and µ u = δG δu and µ ρ = δG δρ are the constant chemical potentials.The far field conditions of u and ρ are where 0 < ρ b < 1 denotes the surfactant concentration in the bulk regions.

2.2.2.
Quantitative behavior of the solution.Although Eq. ( 10) relates the values of u and ρ in equilibrium implicitly, it does not show the profiles of u and ρ quantitatively.In the following, we first derive the exact solution for the case of α > 0, β = 0, and then we compute the numerical profiles of u(x) and ρ(x) for the general case α > 0, β > 0. Before we proceed, we need the following equality which is derived in detail in Appendix A.
Case(1), α > 0, β = 0. Substituting the far field condition (7) into Eqs.( 11) and ( 5) with β = 0, we obtain One can further simplify Eq. ( 12) by substituting the equation of ρ into the equation of u, so we have This gives the exact form of u as which immediately leads the solution of ρ as It is interesting to note that under this case, despite the presence of surfactant, the order parameter u still has the same profile as the case of without surfactant (α = 0).Nevertheless, the absorption of the surfactant does appear in the diffusive interfacial region (near x = 0) where the concentration ρ has significantly larger values than the buck value ρ b .
To see how the surfactant at the interface alters the interfacial tension, we adopt the approaches from [9,15] and define the interfacial tension as where g(ρ b , u = 0) and g ρ (ρ b , u = 0) denote the value of g and g ρ as x → −∞, respectively.Substituting the solutions of u and ρ into Eq.( 18), we have Therefore, the interfacial tension can be calculated as One can immediately see that the presence of surfactant in the diffusive interfacial region reduces the interfacial tension by the magnitude of αρ b .Case(2), α > 0, β > 0. Substituting the far field condition (7) into Eqs.( 11) and ( 5), we obtain Then, the equations of u and ρ become For this case, it is unlikely to find analytical solutions for u and ρ since the above equations are fully nonlinear and coupled.Thus, we solve the equations numerically.We simply discretize Eqs.(20-21) by finite difference method so that the resultant nonlinear equations are solved by Newton's method.The computational domain is chosen as Ω = [−5, 5], and the mesh size is chosen fine enough to resolve the interfacial region.The parameters ε = α = β = 0.02, and ρ b = 0.1.The tolerance of stopping criteria for Newton's iteration is 10 −8 .Figure 1 shows the initial guess (dashed line) and the final computed equilibrium solution (solid line) for the order parameter u (left column) and the surfactant concentration ρ (right column).We perform three different runs (shown by different rows in the figure) by using three different initial guesses for ρ with the same initial guess for u.One can see that the sharp interface eventually forms near the center of the domain x = 0 and the surfactant absorption appears almost in the same interfacial region despite the different initial surfactant distributions.
We also compare the above equilibrium solutions (as shown in the first row of Figure 1) with the one obtained by solving the algebraic equation (10) with the same parameters ε = α = β = 0.02, and ρ b = 0.1.Figure 2 shows the u − ρ plot of those solutions in which we can see that both solutions coincide with each other very well.

Time evolutional equations.
As in [3], we assume a generalized Fick's law so that the mass fluxes of u and ρ are proportional to the gradient of the corresponding chemical potentials.Thus, the Cahn-Hilliard type of equations in the domain Ω can be written as where u 0 and ρ 0 are the initial conditions for u and ρ, respectively.The mobility M ρ is a function of ρ and is chosen to obtain a Fickian equation for ρ [16].Under suitable boundary conditions for the u, ρ and their corresponding chemical potentials δG δu and δG δρ , one can easily derive that both u and ρ satisfy the mass conservation.Throughout the rest of this paper, we assume the periodic boundary conditions are used.
Moreover, the above Cahn-Hilliard formulation leads to the total free energy G defined in (1) is decreasing as time evolves.This can be shown by taking the time derivative of G as Applying integration by parts on the second integral and using the periodic boundary conditions, we have Substituting Eqs.(22) into the above equation and using the divergence theorem, we obtain provided that M ρ ≥ 0. Before closing this section, we shall address an issue related to the range of the surfactant concentration ρ as time evolves.From the energy law, the range of ρ in equilibrium is restricted in (0, 1) which we have verified this numerically in Figure 1.Whether the value of ρ can still be in the range (0, 1) for the time evolution of Eq. ( 22b) is beyond the scope of present study and needs further theoretical investigation.Nevertheless, as shown later in Section 4, our numerical results in Figure 4-6 indicate that the range of ρ is indeed in (0, 1) which ensures the mobility function M ρ is positive.
3.1.1.Basic Concepts.We start by reviewing some basic concepts of the pseudospectral Fourier method [8].Let I = [0, 2π], Z be the set of integers, and N > 0 be an even integer.Denote the 2π periodic function space by H = span{e inx |n ∈ Z}.
Introduce the grid points To approximate u ∈ H, we seek a function I N u ∈ B N of the form where l j (x) is the Lagrange interpolating function satisfying l j (x i ) = δ ij with δ ij being the Kronecker delta function.The symbol and I N is interpreted as an interpolation operator since To approximate the first and high-order derivatives of u, we take the following approach The pseudospectral differentiation formulas ensure that The order reduction is due to the fact that sin(N x/2) resulting from differentiating cos(N x/2), the basis function of B N , is projected out of B N −1 , after applying I N on d dx I N u.
In addition to approximating the derivatives of u, we also need a quadrature rule to compute the integral of u: where 3.1.2.Two-dimensional framework.The above one-dimensional pseudospectral Fourier method can be extended directly into a two-dimensional framework for functions defined on Ω = [0, 2π] × [0, 2π].Let N x and N y be positive even integers.Define the two-dimensional grid points x ij = (x i , y j ) by Let u(x ij ) = u ij .To approximate u, we seek an approximation of the form .
where I N is the interpolation operator on the grid points x ij .One can also define the discrete version of the gradient, Laplace and the bi-harmonic operators as where î and ĵ are unit vectors along x and y directions, respectively.These discrete vector differential operators satisfy the following properties.Consider a vectorvalued function F (x) and scalar functions u(x) and v(x) defined on Ω.Then, from Eq. ( 25) and Eq. ( 24) one can show that which is the summation by parts for 2π-periodic grid functions.Note that, the above formula immediately leads to the following identities The first one is obtained by taking u = 1 in Eq. ( 26), while the second one is obtained by substituting F = ∇ N v into Eq.(26).
3.2.The semi-discrete scheme.To solve the phase field model formulated by Eqs.(22), we seek numerical solutions u N and ρ N of the form satisfying the collocation scheme where Multiplying Eq. (28a) by ∆V = 4π 2 /(N x N y ) and performing a summation over the whole grid points, then applying the first identity of Eq. ( 27), we have Thus, we derive the semi-discrete analogue of the mass conservation for the order parameter u N .The mass conservation for the surfactant concentration ρ N can also be derived similarly.
The collocation scheme of (28) leads to the semi-discrete analogue of energy decreasing as in Eq. ( 23).The discrete energy functional is defined as where Using the second identity of Eq. ( 27) and then employing Eqs.(28a-b), S 1 and S 2 become Then, applying second identity of Eq. ( 27) again, we obtain indicating that G D is decreasing in time provided that 0 ≤ ρ N ≤ 1.
3.3.Fully-discrete scheme and computational techniques.To advance the solution in time, we adopt the semi-implicit time discretization as follows.For clarity, we omit the subscript N and denote u n ij and ρ n ij the numerical solutions u N and ρ N evaluated at the collocation point x ij and at time t n .The fully-discrete scheme for Eq. ( 22) is where Generally speaking, to march an implicit scheme in time often requires the matrix inversion, which can be expensive.However, for the present semi-implicit scheme, there is an efficient way to advance the solution in time.Recall that the collocation points are chosen so that the discrete Fourier coefficients of an approximate function can be efficiently computed through the Fast Fourier Transform (FFT).Taking the discrete Fourier transform of the scheme Eq. ( 29), we have where |ω| 2 = p 2 + q 2 and the symbol (  Notice that, the backward Euler method is only applied to the fourth-order derivative and other lower order derivatives of nonlinear terms are still treated explicitly.Thus, a CF L number constraint for numerical stability is required.In general, since the diffusion operator applied on the nonlinear term f (u) is treated explicitly, we have to choose the time step size like N −2 , i.e., where CF L = O(1).However, for the present study where the parameter ε ranges from 0.02 to 0.05, we found that (after a series numerical experiments) the time step size chosen as still works for stable computations.
4.1.Convergence test.We start our numerical simulations by performing the convergence test for the present scheme.The convergent rate between two consecutive L 2 errors is computed as log 2 , where u N is the numerical solution obtained by the grid number N = N x = N y while u ref is the reference solution obtained by a larger grid number N = 4096.The initial conditions are given as u 0 (x, y) = 0.5 + 0.1 cos(6x) cos(6y), ρ 0 (x, y) = 0.1, and the solutions are computed up to time T = 0.01.Table 1 shows the discrete L 2 -errors and the convergent rates for different values of N .One can see that the  rate of convergence approaches to one, which indicates the first-order of accuracy of the temporal discretization.
In the previous section, we have shown that the present semi-discrete scheme does preserve the mass conservations of u N and ρ N , and the discrete total energy of the system decreases as time evolves.To demonstrate those properties, we define the discrete mass of u and ρ as and the discrete free energy as The time evolution of |U (t) − U (0)|, |Θ(t) − Θ(0)| and G(t), are shown in Fig. 3.One can see that the discrete mass of u and ρ are conserved and the errors are within the machine accuracy.It is also shown that the energy decreases during the time evolution and eventually reaches an equilibrium state.
4.2.One-dimensional phase separation and surfactant absorption.In this subsection, we simulate the one-dimensional phase separation and surfactant absorption for the binary fluid-surfactant system.The initial profiles for u and ρ are chosen as u 0 (x) = 0.3 + 0.01 cos(6x), ρ 0 (x) = 0.2.
For comparison, we also simulate the case of without surfactant in which only the equation of the order parameter u needs to be solved with the parameter α = β = 0. Figure 4 shows different time plots for the order parameter u and the surfactant concentration ρ.
For the case of without surfactant (left column of Fig. 4), we observe that the order parameter u forms multiple separated phase regions of u = 0 and u = 1 immediately.For convenience, in what follows, we refer the regions of u = 1 and u = 0 to the 1-phase and 0-phase regions, respectively.As time evolves, the coarsening behavior continues and eventually (at t = 100) decomposes the domain into the one with only two 1-phase regions left.
For the case of with surfactant (right column of Fig. 4), we also observe that multiple separated 1-and 0-phase regions are formed within a very short time interval.Meanwhile, the ρ profile starts to form peaks indicating that the surfactant is absorbed into the binary fluid interfaces.Around t = 3.2, u starts to change its profile again.As the result, the surfactant located at those gradually disappeared fluid interfaces diffuse and merge into the neighboring two lumps in the regions of x ∈ [1.6, 2.6] and x ∈ [4.8, 5.8] (see the second and third panels of t = 11.8 and t = 20.0).During the time period 50 ≤ t ≤ 100, we have observed no significant changes on the solution profiles, indicating that the system arrives at a low energy state.
From the simulation results, it is shown that the surfactant-interface coupling term, indeed, models the surfactant absorption near the interfacial region, and the entropy term restricts the value of ρ to be in the range of (0, 1).In addition, we also observe that the surfactant delays the coarsening process eventually so that the number of interfaces for the system with surfactant may be larger than the one without surfactant as shown in last panel of Fig. 4. The above numerical results are in a good agreement with the theoretical study [7] regarding that the surfactant favors the creation of interfaces and also stabilizes the formation of phase regions.4.3.Two-dimensional phase separation and surfactant absorption.4.3.1.Surfactant uniformly distributed initially.Figure 5 illustrates the evolution of u and ρ fields with the surfactant uniformly distributed over a two-dimensional domain initially.The initial profiles for u and ρ are chosen as u 0 (x, y) = 0.3 + 0.01 cos(6x) cos(6y), ρ 0 (x, y) = 0.1 + 0.01 cos(6x) cos(6y).
As in the 1D result, we can see that the u field forms multiple separated regions in the beginning (see Fig. 5(a)).Then the coarsening dynamics inherent in the Cahn-Hilliard equation soon drives the apparent phase separation which results in a smaller 1-phase region generally shrinks and disappears, while a larger 1-phase domain expands its area, due to the mass conservation of u.As time evolves, the u field eventually becomes a circular 1-phase region surrounded by 0-phase region, as shown in Fig. 5(j).Driven by the interfacial energy term, the surfactant is absorbed into the binary fluid interfaces so that higher concentration of ρ appears near the interfaces.When a binary fluid interface disappears due to further coarsening, the surfactant adherent on that interface diffuses to other nearby interfaces which causes the maximum value of the surfactant concentration inside the remaining area to be temporally higher due to the mass conservation of the surfactant.Nevertheless, the surfactant will be diffused away and eventually absorbed in the interfaces.

4.3.2.
Surfactant locally distributed initially.The evolution of the system with the surfactant distributed locally at the center of the domain initially is shown in Fig. 6.The initial profiles for u and ρ are chosen as u 0 (x, y) = 0.3 + 0.01 cos(6x) cos(6y), Similar to the previous case, the u field is gradually coarsen into a circular 1-phase region surrounded by the 0-phase region.However, the absorption behavior of the surfactant is completely different from the previous one.Since the surfactant is initially concentrated at the center region of the domain, it takes time for the surfactant to diffuse away from this center region.Consequently, during the early stage of the evolution, the higher concentration of surfactant only appears around the center area of the domain.When there is enough amount of surfactant diffused over the domain, the surfactant concentration field starts to form multiple rings at binary fluid interfaces, as shown in Fig. 6(g-j).This behavior of surfactant diffusing into binary fluid interfacial regions is qualitatively in agreement with the result shown in [12].In this simulation, we have seen many interesting phenomena caused by the highly concentrated surfactant, and it is worth to discuss in more details.As pointed out in [7], the surfactant favors the creation of interfaces.Besides, the surfactant also stabilizes the formation of a phase region, in a way of resisting phase separation.This creation of interfaces can be seen in Fig. 6(a).Notice that, the initial profile for the order parameter u is nearly a flat surface (0.29 ≤ u ≤ 0.31).At t = 0.01 which is just in the beginning of the evolution, we see a fast creation of binary fluid interfaces resulting from the rapid formation of small multiple 1-phase regions around the center area, where the surfactant is highly concentrated but not yet diffused away.The stabilizing phenomenon is visualized from the evolution of the center most located at 1-phase region shown in Fig. 6(d-j).In these u field plots, we see that the center 1-phase region residing in the area containing a large amount of surfactant, maintains its phase state and gradually expands its area, during the entire evolution.
The stabilizing argument has the following consequence.For two 1-phase regions of nearly equal size, the one in the domain with more surfactant is more stable.In other words, the one with less surfactant is likely to shrink and disappear first.The u field illustrated in Fig. 6(d) shows nine 1-phase regions (a 1-phase region surrounded by eight 1-phase regions) of nearly equal size in the center area of the domain, which the central 1-phase region contains more surfactant.The subsequent u field plots (see Fig. 6(e-f)) show that the surrounding 1-phase regions shrink gradually and disappear one by one but the central one remains still.This 1-phase region finally resides at the center of the domain due to the chain reaction of the creation mechanism followed by the stabilizing one, but not due to the initial u profile.Hence, for a random initial profile for u, distributing the surfactant in a localized manner can be used to control the final stable phase regions.

Conclusion.
In this paper, the dynamics of a binary fluid-surfactant system formulated by a phenomenological phase field model is investigated through analytical and numerical computations.We first consider the case of one-dimensional planar interface and prove the existence of the equilibrium solution.Then we derive the analytical equilibrium solution for the order parameter and the surfactant concentration in a particular case.The results show that the present phase field formulation qualitatively mimics the surfactant adsorption on the binary fluid interfaces to that the interfacial tension is reduced.We further study the time-dependent solutions of the system by numerical computations based on pseudospectral Fourier method.The present numerical results are in a good agreement with the previous theoretical study in a way that the surfactant favors the creation of interfaces and also stabilizes the formation of phase regions.

Figure 1 .
Figure 1.The equilibrium solutions of u and ρ.Initial guess: dashed line.Final solution: solid line.
Eqs. (30), then the grid functions u n+1 ij and ρ n+1 ij can be obtained by taking the inverse FFT.

Figure 4 .
Figure 4. Evolution of u and ρ fields.Right and left columns correspond to the evolutions with and without surfactant, respectively.The plots on the left column are taken at t=0.5, 11.8, 20.0, 100.0 (top to bottom) and the plots on the right column are taken at t=0.5, 3.2, 5.0, 100.0 (top to bottom).N = 512, ε = 0.02, α = 0.02, β = 0.02.

Figure 5 .
Figure 5. Field evolution plots with uniformly distributed surfactant initially.The u field plots are placed in the left two columns, and the associated ρ field plots are placed in the right two columns.ε = 0.02, α = 0.02, β = 0.02, N = 256.

Figure 6 .
Figure 6.Field evolution plots with surfactant given initially at the center of the domain.The u field plots are placed in the left two columns, and the associated ρ field plots are placed in the right two columns.ε = 0.02, α = 0.02, β = 0.02, N = 256 Calculate the Fourier coefficients ûn pq by the FFT from the values of u n ij .2. Multiply ip and iq to ûn pq and perform the inverse FFT to obtain ∇ N u. 3. Use the grid functions for ∇ N u, and ρ to construct the grid function F (u, ρ). 4. Compute the discrete divergence ∇ N • F , where the derivatives are computed as in Step 2. 5. Perform FFT on ∇ N • F and multiply the coefficients by −α|ω| 2 to get F n pq .Once the Fourier coefficients ûn+1

Table 1 .
Grid convergence test for u N and ρ N .