Solving the Savage–Hutter equations for granular avalanche flows with a second‐order Godunov type method on GPU

In this article, we apply Davis's second‐order predictor‐corrector Godunov type method to numerical solution of the Savage–Hutter equations for modeling granular avalanche flows. The method uses monotone upstream‐centered schemes for conservation laws (MUSCL) reconstruction for conservative variables and Harten–Lax–van Leer contact (HLLC) scheme for numerical fluxes. Static resistance conditions and stopping criteria are incorporated into the algorithm. The computation is implemented on graphics processing unit (GPU) by using compute unified device architecture programming model. A practice of allocating memory for two‐dimensional array in GPU is given and computational efficiency of two‐dimensional memory allocation is compared with one‐dimensional memory allocation. The effectiveness of the present simulation model is verified through several typical numerical examples. Numerical tests show that significant speedups of the GPU program over the CPU serial version can be obtained, and Davis's method in conjunction with MUSCL and HLLC schemes is accurate and robust for simulating granular avalanche flows with shock waves. As an application example, a case with a teardrop‐shaped hydraulic jump in Johnson and Gray's granular jet experiment is reproduced by using specific friction coefficients given in the literature. Copyright © 2014 John Wiley & Sons, Ltd.     

关于采用粗粒化提高颗粒材料多尺度模拟守恒特性的研究

连续体-颗粒耦合方法常用来描述连续-非连续颗粒行为或解决颗粒材料与其他可变形构件间相互作用问题。粗粒化coarse-graining (CG)是基于统计力学的均匀化方法,由离散的颗粒运动定义连续的宏观物理场。本文利用粗粒化(CG)推导有限元-离散元(FEM-DEM)表面和体积耦合的一般性表达式。对于表面耦合,CG可以将耦合力分布到颗粒-单元接触点以外的位置,如相邻的积分点;对于体积耦合,CG可以将颗粒尺度的运动均匀化到耦合单元上。由粗粒化推导出的耦合项仅包含一个参数,即粗粒化宽度,为均匀化后的宏观场定义了一个可调整的空间尺度。当粗粒化宽度为零时,表面和体积耦合表达式简化为常规局部耦合。本文通过弹性立方体冲击颗粒床和离散-连续介质间波传播两个数值算例,展示使用粗粒化方法提高耦合系统能量守恒的优势,并结合其他耦合参数(如体积耦合深度)讨论了粗粒化参数对数值稳定性和计算效率的影响。     

Performance of a three‐dimensional unstructured mesh compressible flow solver on NVIDIA Fermi‐class graphics processing unit hardware

We describe the performance of Chicoma , a 3D unstructured mesh compressible flow solver, on graphics processing unit (GPU) hardware. The approach used to deploy the solver on GPU architectures derives from the threaded multicore execution model used in Chicoma , and attempts to improve memory performance via the application of graph theory techniques. The result is a scheme that can be deployed on the GPU with high‐level programming constructs, for example, compiler directives, rather than low‐level programming extensions. With an NVIDIA Fermi‐class GPU (NVIDIA Corp., Sta. Clara, CA, USA) and double precision floating point arithmetic, we observe performance gains of 4–5 × on problem sizes of 10⁶– 10⁷ tetrahedra. We also compare GPU performance to threaded multicore performance with OpenMP and demonstrate hybrid multicore‐GPU calculations with adaptive mesh refinement. Published 2012. This article is a US Government work and is in the public domain in the USA.     

The locally conservative Galerkin (LCG) method for solving the incompressible Navier–Stokes equations

In this paper, the locally conservative Galerkin (LCG) method (Numer. Heat Transfer B Fundam. 2004; 46 :357–370; Int. J. Numer. Methods Eng. 2007) has been extended to solve the incompressible Navier–Stokes equations. A new correction term is also incorporated to make the formulation to give identical results to that of the continuous Galerkin (CG) method. In addition to ensuring element‐by‐element conservation, the method also allows solution of the governing equations over individual elements, independent of the neighbouring elements. This is achieved within the CG framework by breaking the domain into elemental sub‐domains. Although this allows discontinuous trial function field, we have carried out the formulation using the continuous trial function space as the basis. Thus, the changes in the existing CFD codes are kept to a minimum. The edge fluxes, establishing the continuity between neighbouring elements, are calculated via a post‐processing step during the time‐stepping operation. Therefore, the employed formulation needs to be carried out using either a time‐stepping or an equivalent iterative scheme that allows post‐processing of fluxes. The time‐stepping algorithm employed in this paper is based on the characteristic‐based split (CBS) scheme. Both steady‐ and unsteady‐state examples presented show that the element‐by‐element formulation employed is accurate and robust. Copyright © 2007 John Wiley & Sons, Ltd.     

Acceleration of iterative Navier-Stokes solvers on graphics processing units

While new power-efficient computer architectures exhibit spectacular theoretical peak performance, they require specific conditions to operate efficiently, which makes porting complex algorithms a challenge. Here, we report results of the semi-implicit method for pressure linked equations (SIMPLE) and the pressure implicit with operator splitting (PISO) methods implemented on the graphics processing unit (GPU). We examine the advantages and disadvantages of the full porting over a partial acceleration of these algorithms run on unstructured meshes. We found that the full-port strategy requires adjusting the internal data structures to the new hardware and proposed a convenient format for storing internal data structures on GPUs. Our implementation is validated on standard steady and unsteady problems and its computational efficiency is checked by comparing its results and run times with those of some standard software (OpenFOAM) run on central processing unit (CPU). The results show that a server-class GPU outperforms a server-class dual-socket multi-core CPU system running essentially the same algorithm by up to a factor of 4.     

Contactless optical extensometer for textile materials

In this paper we present a contactless extensometer. For some flexible materials, with great displacements and deformations, contact during measurement is not acceptable. In fact, contact measurement can modify the tensile behavior, as is the case for fibrous materials. Contactless extensometers usually have to print or glue some marks on the sample, which may cause problems during measurement. These extensometers typically use digital image processingto obtain deformation data. The principle used in this study uses the natural periodicity or surface patterns inherent in most textile materials without any image processing. During deformation the distance between two periods or pattern elements changes, allowing this method to measure the real-time modification of this in-plane distance. The extensometer consists of two parts: an optical device and a signal processing unit performing a Fourier analysis. Some results obtained during a tensile test on woven fabrics and non-wovens are presented here.     

A parallel image encryption algorithm based on the piecewise linear chaotic map and hyper-chaotic map

This paper proposes a parallel digital image encryption algorithm based on a piecewise linear chaotic map (PWLCM) and a four-dimensional hyper-chaotic map (FDHCM). Firstly, two decimals are obtained based on the plain-image and external keys, using a novel parallel quantification method. They are used as the initial value and control parameter for the PWLCM. Then, an encryption matrix and four chaotic sequences are constructed using the PWLCM and FDHCM, which control the permutation and diffusion processes. The proposed algorithm is implemented and tested in parallel based on a graphics processing unit device. Numerical analysis and experimental results show that the proposed algorithm achieves a high encryption speed and a good security performance, which provides a potential solution for real-time image encryption applications.     

A massively parallel GPU‐accelerated model for analysis of fully nonlinear free surface waves

We implement and evaluate a massively parallel and scalable algorithm based on a multigrid preconditioned Defect Correction method for the simulation of fully nonlinear free surface flows. The simulations are based on a potential model that describes wave propagation over uneven bottoms in three space dimensions and is useful for fast analysis and prediction purposes in coastal and offshore engineering. A dedicated numerical model based on the proposed algorithm is executed in parallel by utilizing affordable modern special purpose graphics processing unit (GPU). The model is based on a low‐storage flexible‐order accurate finite difference method that is known to be efficient and scalable on a CPU core (single thread). To achieve parallel performance of the relatively complex numerical model, we investigate a new trend in high‐performance computing where many‐core GPUs are utilized as high‐throughput co‐processors to the CPU. We describe and demonstrate how this approach makes it possible to do fast desktop computations for large nonlinear wave problems in numerical wave tanks (NWTs) with close to 50/100 million total grid points in double/single precision with 4 GB global device memory available. A new code base has been developed in C++ and compute unified device architecture C and is found to improve the runtime more than an order in magnitude in double precision arithmetic for the same accuracy over an existing CPU (single thread) Fortran 90 code when executed on a single modern GPU. These significant improvements are achieved by carefully implementing the algorithm to minimize data‐transfer and take advantage of the massive multi‐threading capability of the GPU device. Copyright © 2011 John Wiley & Sons, Ltd.     

Topology of the Integral Manifolds of the N-Body Problem with Positive Energy

The integral manifolds of the N-body problem are the level sets of energy and angular momentum. For positive energy and non-zero angular momentum, all level sets are diffeomorphic to a non-zero level set of angular momentum on the unit tangent bundle of the configuration space. The one complication that arises in attempting to describe this level set explicitly is the degeneracy at the syzygies of the equations that define angular momentum. In this work, we analyze the behavior of the angular momentum near syzygies, and show how to construct local coordinates near the syzygies. In particular, we show that the projection of the integral manifold onto the configuration space _c is a homotopy equivalence, and use this to compute the homology of the integral manifolds.     

A geometric mass‐preserving redistancing scheme for the level set function

In this paper we describe and evaluate a geometric mass‐preserving redistancing procedure for the level set function on general structured grids. The proposed algorithm is adapted from a recent finite element‐based method and preserves the mass by means of a localized mass correction. A salient feature of the scheme is the absence of adjustable parameters. The algorithm is tested in two and three spatial dimensions and compared with the widely used partial differential equation (PDE)‐based redistancing method using structured Cartesian grids. Through the use of quantitative error measures of interest in level set methods, we show that the overall performance of the proposed geometric procedure is better than PDE‐based reinitialization schemes, since it is more robust with comparable accuracy. We also show that the algorithm is well‐suited for the highly stretched curvilinear grids used in CFD simulations. Copyright © 2010 John Wiley & Sons, Ltd.     

Finite element methods for a class of continuum models for immiscible flows with moving contact lines

In this paper, we present a finite element method for two‐phase incompressible flows with moving contact lines. We use a sharp interface Navier–Stokes model for the bulk phase fluid dynamics. Surface tension forces, including Marangoni forces and viscous interfacial effects, are modeled. For describing the moving contact lines, we consider a class of continuum models that contains several special cases known from the literature. For the whole model, describing bulk fluid dynamics, surface tension forces, and contact line forces, we derive a variational formulation and a corresponding energy estimate. For handling the evolving interface numerically, the level‐set technique is applied. The discontinuous pressure is accurately approximated by using a stabilized extended finite element space. We apply a Nitsche technique to weakly impose the Navier slip conditions on the solid wall. A unified approach for discretization of the (different types of) surface tension forces and contact line forces is introduced. Results of numerical experiments are presented, which illustrate the performance of the solver. Copyright © 2016 John Wiley & Sons, Ltd.     

Reordering and incomplete preconditioning in serial and parallel adaptive mesh refinement and coarsening flow solutions

The effects of reordering the unknowns on the convergence of incomplete factorization preconditioned Krylov subspace methods are investigated. Of particular interest is the resulting preconditioned iterative solver behavior when adaptive mesh refinement and coarsening (AMR/C) are utilized for serial or distributed parallel simulations. As representative schemes, we consider the familiar reverse Cuthill–McKee and quotient minimum degree algorithms applied with incomplete factorization preconditioners to CG and GMRES solvers. In the parallel distributed case, reordering is applied to local subdomains for block ILU preconditioning, and subdomains are repartitioned dynamically as mesh adaptation proceeds. Numerical studies for representative applications are conducted using the object‐oriented AMR/C software system libMesh linked to the PETSc solver library. Serial tests demonstrate that global unknown reordering and incomplete factorization preconditioning can reduce the number of iterations and improve serial CPU time in AMR/C computations. Parallel experiments indicate that local reordering for subdomain block preconditioning associated with dynamic repartitioning because of AMR/C leads to an overall reduction in processing time. Copyright © 2011 John Wiley & Sons, Ltd.     

Running unstructured grid‐based CFD solvers on modern graphics hardware

Techniques used to implement an unstructured grid solver on modern graphics hardware are described. The three‐dimensional Euler equations for inviscid, compressible flow are considered. Effective memory bandwidth is improved by reducing total global memory access and overlapping redundant computation, as well as using an appropriate numbering scheme and data layout. The applicability of per‐block shared memory is also considered. The performance of the solver is demonstrated on two benchmark cases: a NACA0012 wing and a missile. For a variety of mesh sizes, an average speed‐up factor of roughly 9.5 × is observed over the equivalent parallelized OpenMP code running on a quad‐core CPU, and roughly 33 × over the equivalent code running in serial. Copyright © 2010 John Wiley & Sons, Ltd.     

A new volume of fluid method in three dimensions—Part II: Piecewise‐planar interface reconstruction with cubic‐Bézier fit

A new interface reconstruction method in 3D is presented. The method involves a conservative level‐contour reconstruction coupled to a cubic‐Bézier interpolation. The use of the proposed piecewise linear interface calculation (PLIC) reconstruction scheme coupled to a multidimensional time integration provides solutions of second‐order spatial and temporal accuracy. The accuracy and efficiency of the proposed reconstruction algorithm are demonstrated through several tests, whose results are compared with those obtained with other recently proposed methods. An overall improvement in accuracy with respect to other recent methods has been achieved, along with a substantial reduction in the central processing unit time required. Copyright © 2008 John Wiley & Sons, Ltd.     

Aerodynamic design using the truncated Newton algorithm and the continuous adjoint approach

In this paper, the so‐called ‘continuous adjoint‐direct approach’ is used within the truncated Newton algorithm for the optimization of aerodynamic shapes, using the Euler equations. It is known that the direct differentiation (DD) of the flow equations with respect to the design variables, followed by the adjoint approach, is the best way to compute the exact matrix, for use along with the Newton optimization method. In contrast to this, in this paper, the adjoint approach followed by the DD of both the flow and adjoint equations (i.e. the other way round) is proved to be the most efficient way to compute the product of the Hessian matrix with any vector required by the truncated Newton algorithm, in which the Newton equations are solved iteratively by means of the conjugate gradient (CG) method. Using numerical experiments, it is demonstrated that just a few CG steps per Newton iteration are enough. Considering that the cost of solving either the adjoint or the DD equations is approximately equal to that of solving the flow equations, the cost per Newton iteration scales linearly with the (small) number of CG steps, rather than the (much higher, in large‐scale problems) number of design variables. By doing so, the curse of dimensionality is alleviated, as shown in a number of applications related to the inverse design of ducts or cascade airfoils for inviscid flows. Copyright © 2011 John Wiley & Sons, Ltd.     

Cahn–Hilliard modeling of particles suspended in two‐phase flows

In this paper, we present a model for the dynamics of particles suspended in two‐phase flows by coupling the Cahn–Hilliard theory with the extended finite element method (XFEM). In the Cahn–Hilliard model the interface is considered to have a small but finite thickness, which circumvents explicit tracking of the interface. For the direct numerical simulation of particle‐suspended flows, we incorporate an XFEM, in which the particle domain is decoupled from the fluid domain. To cope with the movement of the particles, a temporary ALE scheme is used for the mapping of field variables at the previous time levels onto the computational mesh at the current time level. By combining the Cahn–Hilliard model with the XFEM, the particle motion at an interface can be simulated on a fixed Eulerian mesh without any need of re‐meshing. The model is general, but to demonstrate and validate the technique, here the dynamics of a single particle at a fluid–fluid interface is studied. First, we apply a small disturbance on a particle resting at an interface between two fluids, and investigate the particle movement towards its equilibrium position. In particular, we are interested in the effect of interfacial thickness, surface tension, particle size and viscosity ratio of two fluids on the particle movement towards its equilibrium position. Finally, we show the movement of a particle passing through multiple layers of fluids. Copyright © 2011 John Wiley & Sons, Ltd.     

Discontinuous Galerkin methods on graphics processing units for nonlinear hyperbolic conservation laws

We present a novel implementation of the modal DG method for hyperbolic conservation laws in two dimensions on graphics processing units (GPUs) using NVIDIA's Compute Unified Device Architecture. Both flexible and highly accurate, DG methods accommodate parallel architectures well as their discontinuous nature produces element‐local approximations. High‐performance scientific computing suits GPUs well, as these powerful, massively parallel, cost‐effective devices have recently included support for double‐precision floating‐point numbers. Computed examples for Euler equations over unstructured triangle meshes demonstrate the effectiveness of our implementation on an NVIDIA GTX 580 device. Profiling of our method reveals performance comparable with an existing nodal DG‐GPU implementation for linear problems. Copyright © 2014 John Wiley & Sons, Ltd.     

A confrontation of 1D and 2D RKDG numerical simulation of transitional flow at open‐channel junction

In this study, a comparison between the 1D and 2D numerical simulation of transitional flow in open‐channel networks is presented and completely described allowing for a full comprehension of the modeling water flow. For flow in an open‐channel network, mutual effects exist among the channel branches joining at a junction. Therefore, for the 1D study, the whole system (branches and junction) cannot be treated individually. The 1D Saint Venant equations calculating the flow in the branches are then supplemented by various equations used at the junction: a discharge flow conservation equation between the branches arriving and leaving the junction, and a momentum or energy conservation equation. The disadvantages of the 1D study are that the equations used at the junction are of empirical nature due to certain parameters given by experimental results and moreover they often present a reduced field of validity. On the contrary, for the 2D study, the entire network is considered as a single unit and the flow in all the branches and junctions is solved simultaneously. Therefore, we simply apply the 2D Saint Venant equations, which are solved by a second‐order Runge–Kutta discontinuous Galerkin method. Finally, the experimental results obtained by Hager are used to validate and to compare the two approaches 1D and 2D. Copyright © 2009 John Wiley & Sons, Ltd.     

DEM prediction of particle flows in grinding processes

Particle size reduction is a critical unit process in many industries including mineral processing, cement, food processing, pigments and industrial minerals and pharmaceuticals. The aim is to take large feed material and as efficiently as possible reduce the size of the particles to a target size range. Over time, a very large range of equipment has been developed to perform this for many materials and in many different conditions. Discrete element method (DEM) modelling is a computational tool that can allow detailed exploration of the particle flow and breakage processes within comminution equipment and can assist in developing a clearer and more comprehensive understanding of the detailed processes occurring within. In this paper, we examine the particle and energy flows in representative examples of the equipment used in many grinding processes. We study a 36′ semi‐autogenous mill used in primary grinding for mineral processing, a ball mill used for cement clinker grinding, a grinding table also used for cement grinding, a tower mill used for fine grinding both in mineral processing and for industrial minerals and finally for an Isamill, which is used for ultra‐fine grinding in mineral processing. Copyright © 2008 John Wiley & Sons, Ltd.