Unconditionally convergent nonlinear solver for hyperbolic conservation laws with S-shaped flux functions Patrick Jenny a, * , Hamdi A. Tchelepi b , Seong H. Lee c a Institute of Fluid Dynamics, ETH Zurich, Sonneggstrasse 3, CH-8092 Zurich, Switzerland b Dept. of Energy Resources Engineering, Stanford University, USA c Chevron Company, 6001 Bollinger Canyon Rd., San Ramon, 94583-2324 CA, USA article info Article history: Received 5 January 2009 Received in revised form 25 May 2009 Accepted 10 June 2009 Available online 3 July 2009 Keywords: Implicit discretization Non-convex flux Conservation laws Nonlinear analysis Upwinding Multiphase flow abstract This paper addresses the convergence properties of implicit numerical solution algorithms for nonlinear hyperbolic transport problems. It is shown that the Newton–Raphson (NR) method converges for any time step size, if the flux function is convex, concave, or linear, which is, in general, the case for CFD problems. In some problems, e.g., multiphase flow in porous media, the nonlinear flux function is S-shaped (not uniformly convex or concave); as a result, a standard NR iteration can diverge for large time steps, even if an implicit dis- cretization scheme is used to solve the nonlinear system of equations. In practice, when such convergence difficulties are encountered, the current time step is cut, previous itera- tions are discarded, a smaller time step size is tried, and the NR process is repeated. The criteria for time step cutting and selection are usually based on heuristics that limit the allowable change in the solution over a time step and/or NR iteration. Here, we propose a simple modification to the NR iteration scheme for conservation laws with S-shaped flux functions that converges for any time step size. The new scheme allows one to choose the time step size based on accuracy consideration only without worrying about the conver- gence behavior of the nonlinear solver. The proposed method can be implemented in an existing simulator, e.g., for CO 2 sequestration or reservoir flow modeling, quite easily. The numerical analysis is confirmed with simulation studies using various test cases of nonlinear multiphase transport in porous media. The analysis and numerical experiments demonstrate that the modified scheme allows for the use of arbitrarily large time steps for this class of problems. Ó 2009 Elsevier Inc. All rights reserved. 1. Introduction A wide range of CFD problems, including multiphase flow dynamics in porous media, are described by nonlinear hyper- bolic conservation laws. Various methods are used to solve these conservation equations numerically [7]. Explicit time inte- gration schemes offer accuracy and computational efficiency as long as the limit on the stable time step size is not a major concern [5]. In some problems, however, the time step restriction associated with an explicit scheme is quite severe, and the use of implicit schemes is necessary. In reservoir simulation problems [1], where the evolution of the saturation field in the geologic porous formation, as a function of space in time is sought, it is often the case that for a given global time step size, the CFL numbers in the computational domain can vary by orders of magnitude [6]. In such cases, the use of explicit time integration schemes is simply not feasible, and implicit time integration is required [1,3]. 0021-9991/$ - see front matter Ó 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.jcp.2009.06.032 * Corresponding author. Tel.: +41 44 632 6987. E-mail address: jenny@ifd.mavt.ethz.ch (P. Jenny). Journal of Computational Physics 228 (2009) 7497–7512 Contents lists available at ScienceDirect Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp