Search

We present an integral equation approach to solving the Cahn-Hilliard equation equipped with boundary conditions that model solid surfaces with prescribed Young's angles. The discretization of the system in time using convex splitting leads to a modified biharmonic equation at each time step. To solve it, we split the solution into a volume potential computed with free space kernels, plus the solution to a second kind integral equation (SKIE). The volume potential is evaluated with the help of a box-based volume-FMM method. For non-box domains, the source density is extended by solving a biharmonic Dirichlet problem. The near-singular boundary integrals are computed using quadrature by expansion (QBX) with FMM acceleration. Our method has linear complexity in the number of surface/volume degrees of freedom and can achieve high order convergence in space with adaptive refinement to manage error from function extension.

In this paper, we consider numerical approximations for a three-phase phase-field model, where three fourth-order Cahn-Hilliard equations are nonlinearly coupled together through a Lagrange multiplier term and a sixth-order polynomial bulk potential. By combining the recently developed SAV approach with the linear stabilization technique, we arrive at a novel stabilized-SAV scheme. At each time step, the scheme requires solving only four linear biharmonic equations with constant coefficients, making it the first, to the best of the author's knowledge, totally decoupled, second-order accurate, linear, and unconditionally energy stable scheme for the model. We further prove the unconditional energy stability rigorously and demonstrate the stability and the accuracy of the scheme numerically through the comparisons with the non-stabilized SAV scheme for simulating numerous benchmark numerical examples in 2D and 3D.

In this paper we propose and analyze a weakly nonlinear, energy stable numerical scheme for the strongly anisotropic Cahn-Hilliard model. In particular, a highly nonlinear and singular anisotropic surface energy makes the PDE system very challenging at both the analytical and numerical levels. To overcome this well-known difficulty, we perform a convexity analysis on the anisotropic interfacial energy, and a careful estimate reveals that all its second order functional derivatives stay uniformly bounded by a global constant. This subtle fact enables one to derive an energy stable numerical scheme. Moreover, a linear approximation becomes available for the surface energy part, and a detailed estimate demonstrates the corresponding energy stability. Its combination with an appropriate treatment for the nonlinear double well potential terms leads to a weakly nonlinear, energy stable scheme for the whole system. In particular, such an energy stability is in terms of the interfacial energy with respect to the original phase variable, and no auxiliary variable needs to be introduced. This has important implications, for example, in the case that the method needs to satisfy a maximum principle. More importantly, with a careful application of the global bound for the second order functional derivatives, an optimal rate convergence analysis becomes available for the proposed numerical scheme, which is the first such result in this area. Meanwhile, for a Cahn-Hilliard system with a sufficiently large degree of anisotropy, a Willmore or biharmonic regularization has to be introduced to make the equation well-posed. For such a physical model, all the presented analyses are still available; the unique solvability, energy stability and convergence estimate can be derived in an appropriate manner. In addition, the Fourier pseudo-spectral spatial approximation is applied, and all the theoretical results could be extended for the fully discrete scheme. Finally, a few numerical results are presented, which confirm the robustness and accuracy of the proposed scheme.

Recent applications (e.g. active gels and self-assembly of elastic sheets) motivate the need to efficiently simulate the dynamics of thin elastic sheets. We present semi-implicit time stepping algorithms to improve the time step constraints that arise in explicit methods while avoiding much of the complexity of fully-implicit approaches. For a triangular lattice discretization with stretching and bending springs, our semi-implicit approach involves discrete Laplacian and biharmonic operators, and is stable for all time steps in the case of overdamped dynamics. For a more general finite-difference formulation that can allow for general elastic constants, we use the analogous approach on a square grid, and find that the largest stable time step is two to three orders of magnitude greater than for an explicit scheme. For a model problem with a radial traveling wave form of the reference metric, we find transitions from quasi-periodic to chaotic dynamics as the sheet thickness is reduced, wave amplitude is increased, and damping constant is reduced., 22 pages, 10 figures

A system of Fredholm second kind integral equations (SKIEs) is constructed for the modified biharmonic equation in two dimensions with gradient boundary conditions. Such boundary value problem arises naturally when the incompressible Navier-Stokes equations are solved via a semi-implicit discretization scheme and the resulting boundary value problem at each time step is then solved using the pure stream-function formulation. The advantages of such an approach (Greengard and Kropinski, 1998) [14] are two fold: first, the velocity is automatically divergence free, and second, complicated (nonlocal) boundary conditions for the vorticity are avoided. Our construction, though similar to that of Farkas (1989) [10] for the biharmonic equation, is modified such that the SKIE formulation has low condition numbers for large values of the parameter. We illustrate the performance of our numerical scheme with several numerical examples. Finally, the scheme can be easily coupled with standard fast algorithms such as FFT, fast multipole methods (Greengard and Rokhlin, 1987) [15], or fast direct solvers (Ho and Greengard, 2012; Martinsson and Rokhlin, 2005) [17,25] to achieve optimal complexity, bringing accurate large-scale long-time fluid simulations within practical reach.

In this paper, we describe a new class of fast solvers for separable elliptic partial differential equations in cylindrical coordinates $(r,\theta,z)$ with free-space radiation conditions. By combining integral equation methods in the radial variable $r$ with Fourier methods in $\theta$ and $z$, we show that high-order accuracy can be achieved in both the governing potential and its derivatives. A weak singularity arises in the Fourier transform with respect to $z$ that is handled with special purpose quadratures. We show how these solvers can be applied to the evaluation of the Coulomb collision operator in kinetic models of ionized gases., Comment: 20 pages, 5 figures

A Fredholm second kind integral equation (SKIE) formulation is constructed for the Dirichlet problem of the biharmonic equation in three dimensions. A fast numerical algorithm is developed based on the constructed SKIE. Its performance is illustrated via several numerical examples.

We analyze the source of the self-force errors in the node-centered adaptive-mesh-refinement particle-in-cell (AMR-PIC) algorithm and propose a method for reducing those self-forces. Our approach is based on a method of charge deposition due to Mayo [A. Mayo, The fast solution of Poisson's and the biharmonic equations on irregular regions, SIAM Journal of Numerical Analysis 21(2) (1984) 285-299] that can reduce the self-force error to any specified degree of accuracy.

We present a nonstiff, fully adaptive mesh refinement-based method for the Cahn-Hilliard equation. The method is based on a semi-implicit splitting, in which linear leading order terms are extracted and discretized implicitly, combined with a robust adaptive spatial discretization. The fully discretized equation is written as a system which is efficiently solved on composite adaptive grids using the linear multigrid method without any constraint on the time step size. We demonstrate the efficacy of the method with numerical examples. Both the transient stage and the steady state solutions of spinodal decompositions are captured accurately with the proposed adaptive strategy. Employing this approach, we also identify several stationary solutions of that decomposition on the 2D torus.

The evaluation of sums (matrix-vector products) of the solutions of the three-dimensional biharmonic equation can be accelerated using the fast multipole method, while memory requirements can also be significantly reduced. We develop a complete translation theory for these equations. It is shown that translations of elementary solutions of the biharmonic equation can be achieved by considering the translation of a pair of elementary solutions of the Laplace equations. The extension of the theory to the case of polyharmonic equations in R^3 is also discussed. An efficient way of performing the FMM for biharmonic equations using the solution of a complex valued FMM for the Laplace equation is presented. Compared to previous methods presented for the biharmonic equation our method appears more efficient. The theory is implemented and numerical tests presented that demonstrate the performance of the method for varying problem sizes and accuracy requirements. In our implementation, the FMM for the biharmonic equation is faster than direct matrix-vector product for a matrix size of 550 for a relative L"2 accuracy @e"2=10^-^4, and N=3550 for @e"2=10^-^1^2.

It is shown that a discrete delta function can be constructed using a technique developed by Anita Mayo [The fast solution of Poisson's and the biharmonic equations on irregular regions, SIAM J. Sci. Comput. 21 (1984) 285-299] for the numerical solution of elliptic equations with discontinuous source terms. This delta function is concentrated on the zero level set of a continuous function. In two space dimensions, this corresponds to a line and a surface in three space dimensions. Delta functions that are first and second order accurate are formulated in both two and three dimensions in terms of a level set function. The numerical implementation of these delta functions achieves the expected order of accuracy.

A high-order harmonic spectral filter (HSF) is further studied using double Fourier series (DFS), which performs filtering in terms of successive inversion of tridiagonal matrices with complex-valued elements. The high-order harmonics filter equation is split into multiple Helmholtz equations. It is found that the filter provides the same order of accuracy as the spectral filter in [J. Comput. Phys. 177 (2002) 313] that consists of the pentadiagonal matrices with real-valued elements. The advantage of the filter over the previous one lies on the simplicity and easiness of numerical implementation or computer coding, just requiring the same complexity as Poisson's equation solver. However, the operation count associated with the filter increases by a factor of about 2. To circumvent the inefficiency while preserving the simplicity, an easy way to construct pentadiagonal matrices associated with the biharmonic equation is presented in which the tridiagonal matrices related with Poisson's equation are manipulated. Computational efficiency of the spectral filter is discussed in terms of the relative computing time to the spectral transform. It is revealed that the computing cost (requiring O(N2) operations with N being the truncation) for the spectral filtering, even with the complex-valued matrices, is not significant in the DFS spectral model that is characterized by O(N2 log2 N) operations. Filtering with different DFS expansions is discussed with a focus on the accuracy and pole condition. It is shown that the DFS violating the pole conditions produces a discontinuity at poles in case of wave truncation, and its influence spreads over the globe. The spectral filter is applied to two kinds of uniform-grid data, the regular and the shifted grids, and the results are compared with each other. The operator splitting (or spherical harmonics factorization) makes it feasible to apply the finite difference method to the high-order harmonics filter with ease because only the five-point stencil computations are required. The application could also be extended to other numerical methods only if the Helmholtz equation solver is available.

We present a new method for the solution of the Stokes equations. The main features of our method are: (1) it can be applied to arbitrary geometries in a black-box fashion; (2) it is second-order accurate; and (3) it has optimal algorithmic complexity. Our approach, to which we refer as the embedded boundary integral method (EBI), is based on Anita Mayo's work for the Poisson's equation: "The Fast Solution of Poisson's and the Biharmonic Equations on Irregular Regions", SIAM Journal on Numerical Analysis, 21 (1984) 285-299. We embed the domain in a rectangular domain, for which fast solvers are available, and we impose the boundary conditions as interface (jump) conditions on the velocities and tractions. We use an indirect boundary integral formulation for the homogeneous Stokes equations to compute the jumps. The resulting equations are discretized by Nystrom's method. The rectangular domain problem is discretized by finite elements for a velocity-pressure formulation with equal order interpolation bilinear elements (Q1-Q1). Stabilization is used to circumvent the inf-sup condition for the pressure space. For the integral equations, fast matrix-vector multiplications are achieved via an N logN algorithm based on a block representation of the discrete integral operator, combined with (kernel independent) singular value decomposition to sparsify low-rank blocks. The regular grid solver is a Krylov method (conjugate residuals) combined with an optimal two-level Schwartz-preconditioner. For the integral equation we use GMRES. We have tested our algorithm on several numerical examples and we have observed optimal convergence rates.

A fast Schur complement algorithm is presented for computing the piecewise Hermite bicubic orthogonal spline collocation solution of the biharmonic Dirichlet problem on a rectangular region. On an N × N uniform partition, the algorithm, which involves the preconditioned conjugate gradient method and fast Fourier transforms, requires O(N2 log2 N) arithmetic operations.

In this paper, we study the convergence property of PHYSALIS when it is applied to incompressible particle flows in two-dimensional space. PHYSALIS is a recently proposed iterative method which computes the solution without imposing the boundary conditions on the particle surfaces directly. Instead, a consistency equation based on the local (near particle) representation of the solution is used as the boundary conditions. One of the important issues needs to be addressed is the convergence properties of the iterative procedure. In this paper, we present the convergence analysis using Laplace and biharmonic equations as two model problems. It is shown that convergence of the method can be achieved but the rate of convergence depends on the relative locations of the cages. The results are directly related to potential and Stokes flows. However, they are also relevant to Navier-Stokes flows, heat conduction in composite media, and other problems.

This article provides accurate spectral solutions of the driven cavity problem, calculated in the vorticity-stream function representation without smoothing the corner singularities--a prima facie impossible task. As in a recent benchmark spectral calculation by primitive variables of Botella and Peyret, closed-form contributions of the singular solution for both zero and finite Reynolds numbers are subtracted from the unknown of the problem tackled here numerically in biharmonic form. The method employed is based on a split approach to the vorticity and stream function equations, a Galerkin-Legendre approximation of the problem for the perturbation, and an evaluation of the nonlinear terms by Gauss-Legendre numerical integration. Results computed for Re = 0, 100, and 1000 compare well with the benchmark steady solutions provided by the aforementioned collocation-Chebyshev projection method. The validity of the proposed singularity subtraction scheme for computing time-dependent solutions is also established.

We present new, highly accurate, and efficient methods for computing the motion of a large number of two-dimensional closed interfaces in a slow viscous flow. The interfacial velocity is found through the solution to an integral equation whose analytic formulation is based on complex-variable theory for the biharmonic equation. The numerical methods for solving the integral equations are spectrally accurate and employ a fast multipole-based iterative solution procedure, which requires only O(N) operations where N is the number of nodes in the discretization of the interface. The interface is described spectrally, and we use evolution equations that preserve equal arclength spacing of the marker points. We assume that the fluid on one side of the interface is inviscid and we discuss two different physical phenomena: bubble dynamics and interfacial motion driven by surface tension (viscous sintering). Applications from buoyancy-driven bubble interactions, the motion of polydispersed bubbles in an extensional flow, and the removal of void spaces through viscous sintering are considered and we present large-scale, fully resolved simulations involving O(100) closed interfaces.

A high-order harmonic spectral filter (HSF) is implemented to smooth out a field variable defined on a spherical surface using the double Fourier series (DFS) as orthogonal basis functions. The filter consists of the solution of the high-order harmonic diffusion equation with the implicit method, where the high-order harmonic operator is split into second- or lower-order harmonic operators. The second-order harmonic operator is replaced by a pentadiagonal matrix whose elements are the spectral coefficients. First, a biharmonic spectral filter (BiHSF), the prototype of the high-order HSF, is developed where only the second-order harmonic operator is included. It is found that the computational error for the inversion of a pentadiagonal matrix remains in the order of machine rounding. Compared to the mixed-order HSF with DFS used in the previous study, which contains the second- and third-order harmonic operators, the BiHSF can provide a sharper cutoff of high wavenumbers as well as improved computational efficiency. These advantages come from the fact that for each set of spectral coefficients the BiHSF needs only a single inversion of the pentadiagonal matrix whereas the mixed-order HSF requires triple inversions and an auxiliary operation of the tridiagonal matrix. Based on the BiHSF, the high-order HSF up to the sixth order, which is composed of a multiple inversion of tri- or pentadiagonal matrices, is designed. Tests with the rotated Gaussian fields revealed that the HSF with DFS is isotropic. Application to various problems, including a time-dependent strongly nonlinear case and the observed flow, indicates that the high-order HSF is well capable of smoothing out selectively high-wavenumber components without significantly affecting the conserved quantity, such as total (kinetic) energy.

This paper describes a variational formulation for solving the time-dependent Navier–Stokes equations expressed in terms of the stream function and vorticity around multiple airfoils. This approach is an extension to the case of multiply connected domains of the weak formulation based on explicit viscous diffusion recently proposed by Guermond and Quartapelle. In their method the momentum equation was interpreted as a dynamical equation governing the evolution of the (weak) Laplacian of the stream function, while the Poisson equation for the latter was used as an expression to evaluate the distribution of the vorticity. Time discretizations with the viscous term made explicit were used, leading to the viscosity being split from the incompressibility, similarly to the primitive variable fractional-step method. In the present work the multiconnectedness is addressed by introducing an influence matrix to determine the constant values of the stream function on the airfoils in a noniterative fashion. The explicit treatment of the viscous term leads to an influence matrix rooted in the harmonic problem instead of in the biharmonic problem occurring in methods enforcing integral conditions on the vorticity, such as the Glowinski–Pironneau method. The influence matrix changes at each time step or is constant depending on whether a semi-implicit or fully explicit treatment is adopted for the nonlinear term. The resulting split method is implemented using a first-order Euler backward difference or a second-order BDF scheme and linear finite elements. Numerical results are given and compared with the solutions obtained by means of the biharmonic formulation for multiply connected domains.

We present new methods for computing the motion of two-dimensional closed interfaces in a slow viscous flow. The interfacial velocity is found through the solution to an integral equation whose analytic formulation is based on complex-variable theory for the biharmonic equation. The numerical methods for solving the integral equations are spectrally accurate and employ a fast multipole-based iterative solution procedure, which requires only O(N) operations where N is the number of nodes in the discretization of the interface. The interface is described spectrally, and we use evolution equations that preserve equal spacing in arclength of the marker points. A small-scale decomposition is performed to extract the dominant term in the evolution of the interface, and we show that this dominant term leads to a CFL-type stability constraint. When in an equal arclength frame, this term is linear and we show that implicit time-integration schemes that are explicit in Fourier space can be formulated. We verify this analysis through several numerical examples.

It is proved that a discrete delta function introduced by Smereka [P. Smereka, The numerical approximation of a delta function with application to level set methods, J. Comput. Phys. 211 (2006) 77-90] gives a second-order accurate quadrature rule for surface integrals using values on a regular background grid. The delta function is found using a technique of Mayo [A. Mayo, The fast solution of Poisson's and the biharmonic equations on irregular regions, SIAM J. Numer. Anal. 21 (1984) 285-299]. It can be expressed naturally using a level set function.

We present a class of integral equation methods for the solution of biharmonic boundary value problems, with applications to two-dimensional Stokes flow and isotropic elasticity. The domains may be multiply-connected and finite, infinite or semi-infinite in extent. Our analytic formulation is based on complex variables, and our fast multipole-based iterative solution procedure requiresO(N) operations, whereNis the number of nodes in the discretization of the boundary. The performance of the methods is illustrated with several large-scale numerical examples.

Compact computational schemes for the first biharmonic problem in a rectangle a × b with fourth- and second-order truncation errors and expressed in matrix form are presented. The matrix formulation of the fourth-order schemes is based on a kind of (p, q, r, s) non-coupled approach which may be viewed as an extension of the known non-coupled (p, q) method--see 6]--that uses 9-point stencils to approximate the Laplacian at the l × m interior nodes of a grid, with l and m standing for the number of equidistant subdivision points of the edges a and b, respectively. The final matrix equation which serves as a fourth-order discrete equivalent of the problem may be expressed in conventional formulation as a symmetric system with l × m unknowns. For second-order schemes, the matrix equation in its initial form is based on the non-coupled (p, q) approach employing 5-point stencils to approximate the Laplacian at the l × m interior grid nodes, while the final system in its conventional formulation is again symmetric. For the specific values of p, q, r, and s used in this paper, namely p = 1, q = 2, r = 3, and s = 4, it is possible to solve the problem by means of a quasi direct method after reducing the solution of both fourth- and second-order schemes to the solution of two symmetric and positive definite linear systems each of order equal to min(l, m). In addition, for the same values of p, q, r, and s, the employment of the SOR iterarive method leads, after a reasonable number of iterations, to results which agree with the ones obtained by the already mentioned quasi-direct method of solution. The experimentally observed accuracy for schemes with fourth-order truncation error (at least for the problems considered in this paper) was also of fourth order. For schemes with second-order truncation errors the observed accuracy was also of second order, as it should be, since the schemes in question express in effect the (p, q) approach for which--see 5, 7]--a formal proof of accuracy exists.

The method of fundamental solutions (MFS) is a relatively new technique for the numerical solution of certain elliptic boundary value problems. It falls in the class of methods generally called boundary methods, and, like the well-known boundary integral equation method, is applicable when a fundamental solution of the differential equation is known. In the MFS, the approximate solution is a linear combination of fundamental solutions with singularities placed outside the domain of the problem. The locations of the singularities are either preassigned or determined along with the coefficients of the fundamental solutions so that the approximate solution satisfies the boundary conditions as well as possible. In many applications, these quantities are determined by a least squares fit of the boundary conditions, a nonlinear problem, which is solved using standard software. In this paper, the MFS is formulated for biharmonic problems and is applied to a variety of standard test problems as well as to problems arising in elasticity and fluid flow. © 1987. 69 2 434 459 Cited By :138

The biharmonic boundary integral solution of two-dimensional creeping flow of an incompressible Newtonian fluid is extended to problems with free surfaces. Iterative construction of the free surface with the normal stress condition is achieved with a shooting method that is most appropriate for thin films and menisci of widely varying curvature. A continuation scheme is also introduced to systematically construct the solution as a function of a given parameter. These new techniques, combined with the known storage and computational advantages of the biharmonic boundary integral method, allow construction of difficult free surfaces. We demonstrate this by solving the Bretherton problem of a two-dimensional air bubble traveling in a channel. Agreement with asymptotic theory and the finite-difference method is shown. Other problems on formation and transport of bubbles in channels are also tackled.

The repair of small blood vessels and the pathological growth of internal blood clots involve the formation of platelet aggregates adhering to portions of the vessel wall. Our microscopic model represents blood by a suspension of discrete massless platelets in a viscous incompressible fluid. Platelets are initially noncohesive; however, if stimulated by an above-threshold concentration of the chemical ADP or by contact with the adhesive injured region of the vessel wall, they become cohesive and secrete more ADP into the fluid. Cohesion between platelets and adhesion of a platelet to the injured wall are modeled by creating elastic links. Repulsive forces prevent a platelet from coming too close to another platelet or to the wall. The forces affect the fluid motion in the neighborhood of an aggregate. The platelets and secreted ADP both move by fluid advection and diffusion. The equations of the model are studied numerically in two dimensions. The platelet forces are calculated implicitly by minimizing a nonlinear energy function. Our minimization scheme merges Gill and Murray's ( Math. Programming 7 (1974) , 311) modified Newton's method with elements of the Yale sparse matrix package. The stream-function' formulation' of the Stokes' equations for the fluid motion under the influence of platelet forces is solved using Bjorstad's biharmonic solver (“Numerical Solution of the Biharmonic Equation,” Ph. D. Thesis, Stanford University, 1980) . The ADP transport equation is solved with an alternating-direction implicit scheme. A linked-list data structure is introduced to keep track of changing platelet states and changing configurations of interplatelet links. Results of calculations with healthy platelets and with diseased platelets are presented.

The biharmonic or Navier-Stokes problems in the form of a coupled pair of Dirichlet problems (J. Smith, SIAM J. Numer. Anal . 5 , 323 (1968) are numerically solved by using a two parameter point SOR method. We emphasize the dependance of the convergence domain on the discrete boundary formulae. Optimization of this SOR method is heuristic but can be foreseen with a satisfactory precision. The optimal region is rather large and although using imprecisely optimal parameters, we can greatly improve the classical block SOR method (L. W. Ehrlich, SIAM J. Numer. Anal . 8 , 278 (1971); L. W. Ehrlich and M. M. Gupta, SIAM J. Numer. Anal . 12 , 773 (1975); M. M. Gupta and R. P. Manohar, J. Comput. Phys . 31 , 265 (1979); M. Khalil, Thesis, Universite Paul Sabatier, Toulouse 1983 (unpublished)).

A biharmonic boundary integral equation (BBIE) method is used to solve a two dimensional contained viscous flow problem. In order to achieve a greater accuracy than is usually possible in this type of method analytic expressions are used for the piecewise integration of all the kernel functions rather than the more time-consuming method of Gaussian quadrature.

New difference formulas are derived for solving the biharmonic problem in two dimensions over a rectangular domain. These methods use only the nine grid points of a single mesh cell and do not require fictitious points in order to approximate the boundary conditions. Derivatives of the solution are obtained as a by-product of the methods. Second order formulas are derived for both the first and second biharmonic problems. In numerical experiments, the new second order formulas compare favourably with the standard second order methods. Extensions to fourth order formulas are given. The method of deriving these formulas can be used to derive similar formulas for arbitrarily shaped regions.

The biharmonic equation transformed in the computational domain is solved for the generation of boundary-orthogonal curvilinear coordinate systems. The method permits direct and complete control of the mesh point location on the boundary as well as the angle of intersection of the coordinate lines with the boundary. The method may also be used for the generation of meshes in segmented fields. Finally, the method can be easily extended in three dimensions.

The Dirichlet problem for the biharmonic equation is solved using the thirteen-point formula. The prescribed normal derivative on the boundary is replaced by two classes of boundary approximations in order to define the solution at certain fictitious node points. A direct method is used to solve the resulting system of algebraic equations. It is found that the accuracy of the numerical solution strongly depends upon the boundary approximation used, as in the coupled-equation approach. However, the cost of obtaining the solution is almost independent of the boundary approximation, unlike the coupled-equation approach.

This paper describes a fourth-order finite difference model of the equatorial ocean that is designed to study dynamic and thermodynamic processes on time scales of a decade or less. It is a primitive equation model employing the reduced gravity assumption so that the deep ocean is at rest below the active upper ocean. The model consists of a surface mixed layer and an active layer below, which is divided into an arbitrary number of numerical layers by means of a sigma coordinate. The model can be used in an unstratified version, when temperature acts as a passive tracer, as well as in the full stratified version. The numerical formulation of the model is described in detail. Experiments comparing three different horizontal smoothers: Shapiro filter, Laplacian friction, and biharmonic friction are presented. It is concluded that, at the level needed to maintain computational stability, the Shapiro filter damps the fields least; in addition, it is faster and easier to implement when the horizontal finite difference grid is stretched.

A strongly implicit procedure is described which solves the system of 13 point finite difference equations associated with the biharmonic and similar fourth order elliptic equations. No factorization of the equation is required, and for the majority of problems, a universal set of iteration parameters provide rapid rates of convergence. In a comparison with another solution procedure for the biharmonic equation, the new method appears to reduce the computation required to about one-third.

With the aim of efficiently simulating three-dimensional multiphase turbulent flows with a phase-field method, we propose a new discretization scheme for the biharmonic term (the 4th-order derivative term) of the Cahn-Hilliard equation. This novel scheme can significantly reduce the computational cost while retaining the same accuracy as the original procedure. Our phase-field method is built on top of a direct numerical simulation solver, named AFiD (www.afid.eu) and open-sourced by our research group. It relies on a pencil distributed parallel strategy and a FFT-based Poisson solver. To deal with large density ratios between the two phases, a pressure split method [1] has been applied to the Poisson solver. To further reduce computational costs, we implement a multiple-resolution algorithm which decouples the discretizations for the Navier-Stokes equations and the scalar equation: while a stretched wall-resolving grid is used for the Navier-Stokes equations, for the Cahn-Hilliard equation we use a fine uniform mesh. The present method shows excellent computational performance for large-scale computation: on meshes up to 8 billion nodes and 3072 CPU cores, a multiphase flow needs only slightly less than 1.5 times the CPU time of the single-phase flow solver on the same grid. The present method is validated by comparing the results to previous studies for the cases of drop deformation in shear flow, including the convergence test with mesh refinement, and breakup of a rising buoyant bubble with density ratio up to 1000. Finally, we simulate the breakup of a big drop and the coalescence of O(10^3) drops in turbulent Rayleigh-B\'enard convection at a Rayleigh number of $10^8$, observing good agreement with theoretical results., Comment: 32 pages, 14 figures