Simulation of Dual Mixed Refrigerant Natural Gas Liquefaction Processes Using a Nonsmooth Framework

Abstract

Natural gas liquefaction is an energy intensive process where the feed is cooled from ambient temperature down to cryogenic temperatures. Different liquefaction cycles exist depending on the application, with dual mixed refrigerant processes normally considered for the large-scale production of Liquefied Natural Gas (LNG). Large temperature spans and small temperature differences in the heat exchangers make the liquefaction processes difficult to analyze. Exergetic losses from irreversible heat transfer increase exponentially with a decreasing temperature at subambient conditions. Consequently, an accurate and robust simulation tool is paramount to allow designers to make correct design decisions. However, conventional process simulators, such as Aspen Plus, suffer from significant drawbacks when modeling multistream heat exchangers. In particular, no rigorous checks exist to prevent temperature crossovers. Limited degrees of freedom and the inability to solve for stream variables other than outlet temperatures also makes such tools inflexible to use, often requiring the user to resort to a manual iterative procedure to obtain a feasible solution. In this article, a nonsmooth, multistream heat exchanger model is used to develop a simulation tool for two different dual mixed refrigerant processes. Case studies are presented for which Aspen Plus fails to obtain thermodynamically feasible solutions.

1. Introduction

Natural gas plays a major role in the global shift towards new environmentally friendly energy sources. Low CO₂ emissions and no particle emissions upon combustion means that natural gas provides a cleaner alternative to oil and coal. However, a significant challenge with using natural gas is related to its transportation, especially over long distances. Alternative technologies for the transportation of natural gas exist, where the conventional approach is to use pipelines. However, pipeline transportation requires a large initial investment in infrastructure. Moreover, it ties the seller of the gas to a small set of customers at the receiving terminals. Excessive infrastructure and transportation costs for long distances also make pipeline gas difficult to export to the global energy market. Because of this, the recent trend has been towards natural gas liquefaction. The liquefaction of natural gas is a very energy intensive process that requires cooling to about −162 C. Investments in expensive, customized, and proprietary technology are necessary. Together with high operating costs, liquefaction accounts for about 30–40% of the total cost of the Liquefied Natural Gas (LNG) chain [1]. LNG production plants are frequently divided into three main categories: base-load, peak-shaving, and small-scale plants. Single mixed refrigerant (SMR) liquefaction processes are normally considered for small-scale and peak-shaving LNG production where capital costs rather than operational costs are the main concern. In the case of base-load plants, high production volumes with accompanying high operating costs advocate more efficient designs. A popular alternative is the Dual Mixed Refrigerant (DMR) processes due to their high efficiency and flexible design. The added flexibility comes at a cost computationally, however, as DMR processes are significantly more complex to model, simulate, and optimize than SMR processes.

A large temperature span and small driving forces in the heat exchangers at cryogenic temperatures make LNG processes particularly difficult to analyze. The small temperature driving forces are a direct result of the exergetic losses from irreversible heat transfer increasing exponentially with a decreasing temperature at subambient conditions. As a result, even small inaccuracies in the process model at cryogenic temperatures may propagate to a significant amount of additionally required actual compression power. Furthermore, the high operating and capital costs in LNG processes favor an accurate and robust simulation tool for engineers to make the best design decisions. Nevertheless, state-of-the-art process simulation tools, such as Aspen Plus [2] and Aspen HYSYS [3], suffer from limitations in the modeling of multistream heat exchangers (MHEXs), in particular, the lack of any rigorous checks to prevent temperature crossovers and thus, the possibility that the process simulator could converge to an infeasible design [4,5].

Different modeling approaches for MHEXs that perform these checks have been proposed in the literature. A MHEX model for optimization of LNG processes was developed by Kamath et al. [4]. The model uses concepts from Pinch Analysis by treating the MHEX as a heat integration problem with no external utilities. Phase changes in the MHEX are handled using complementarity constraints and do not need to be determined a priori. Phase transitions are a known issue in modeling MHEXs as phase boundaries vary dynamically during simulation and optimization. Consequently, the location of the phase boundary as well as the phase path of a process stream cannot readily be determined a priori. Models that are capable of handling unknown phase states are therefore necessary to fully analyze LNG processes. The result is a fully equation-oriented (EO) model that was solved to local optimality for the single mixed refrigerant PRICO process. Another optimization model for MHEXs was developed by Hasan et al. [6,7] using a superstructure approach, where the MHEX is represented by a network of two-stream heat exchangers. Optimal operation is then determined by solving a heat integration (HI) problem with no external utilities. However, the model can handle phase changes in the MHEX only as long as the inlet and outlet states are known a priori. The result is a mixed-integer nonlinear program (MINLP) that is computationally expensive and requires global optimization tools to solve. Later, Rao and Karimi [8] proposed an alternative superstructure that can handle unknown phase states without including binary variables. In the superstructure, the MHEX is represented as a series of stream bundles, where each bundle is modeled as a set of two-stream heat exchangers between the hot and cold streams. Nonlinear constraints are included for the phase changes, which are modeled to occur at the endpoints of the two-stream heat exchangers. In this way, each heat exchanger operates in a single phase-regime and can be solved using a corresponding property model. In the article, property evaluations were done using a conventional process simulator (here Aspen Plus). However, tailored property models can also be added with accompanying binary variables for the phase transitions. The final model is a nonconvex NLP (or MINLP with custom property models) that is computationally expensive to solve, mostly because of the repeated property evaluations done by the process simulator, and again, requires global optimization methods. Pattison and Baldea [9] developed a MHEX model using their own pseudo-transient EO approach for modeling and simulation. The model requires no prior knowledge of the stream paths and phase boundaries. However, a relative sequence of stream temperatures must be known and fixed prior to optimization to define the enthalpy intervals for the composite curves. Temperatures are then calculated at the endpoint of each interval by introducing a nonphysical time-dependent temperature variable and then solving a system of differential-algebraic equations (DAEs). No Boolean variables or disjunctions are used for handling phase transitions. Instead, the dummy transient variable is perturbed across the phase boundary while keeping the temperature constant, and the property models are resolved using the solution from the previous time step as a starting point.

Recent literature on MHEX modeling for natural gas liquefaction has primarily focused on flowsheet optimization. The models presented in the paragraph above, require solving nonconvex optimization problems, sometimes to guaranteed global optimality, where feasible heat transfer is guaranteed at the solution through minimum approach temperature inequality constraints. Furthermore, the models only consider optimization of the simple PRICO process, which consists of a single MHEX for natural gas liquefaction. Single mixed refrigerant processes normally feature relatively simple designs with minimal auxiliary equipment, and are favored for small-scale applications or offshore production where size and capital cost are of primary concern. The large-scale production of LNG, however, demands the use of highly complex and optimized processes with additional MHEXs and refrigerants to reduce the driving forces between the hot and cold composite curves. Dual mixed refrigerant processes are examples of such designs. However, the added complexity makes DMR processes difficult to simulate and optimize with the custom MHEX models presented in the literature. Instead, these processes are normally studied using conventional process simulation tools, such as Aspen HYSYS or Aspen Plus, for simulation and a stochastic search algorithm for optimization [10].

Although the custom MHEX models focus on flowsheet optimization, process simulation remains an important step in the design and analysis of liquefaction processes. Whereas optimization ostensibly provides the designer a best known design or operation, simulation models can yield useful information on existing designs or operating points that are not necessarily optimal. Furthermore, process simulation allows for probing the behaviour of a system in the neighborhood of the current operating point, to investigate possible improvements in design that can be achieved with little effort. In addition, it can be used as a first step towards flowsheet optimization by using the simulation results as starting points for the optimizer, usually improving reliability of the convergence [11]. However, flowsheet convergence to a thermodynamically feasible operation is more challenging in simulation than optimization due to the absence of minimum temperature difference inequality constraints in the problem formulation. This is a particular issue when using conventional process simulators such as Aspen HYSYS or Aspen Plus, which frequently converge to a solution with temperature crossovers [4,5]. Furthermore, the MHEX block included in Aspen solves the overall energy balance for a single outlet temperature from the heat exchanger. Other process variables (i.e., refrigerant pressure levels and compositions) must be held fixed during simulation, making Aspen inflexible to use for LNG simulation, and often requiring the user to resort to a manual iterative procedure for locating feasible solutions.

A MHEX model that has been developed for both the simulation and optimization of LNG processes was presented by Watson et al. [5,12]. The model uses new developments in nonsmooth analysis [13] and pinch analysis concepts to solve the MHEX as a heat integration problem with no external utilities. A reformulation of the Duran and Grossmann [14] pinch location algorithm was developed to calculate the minimum temperature difference in the MHEXs and to prevent temperature crossovers during simulation. The model size is independent of the number of streams exchanging heat, which is particularly advantageous when many substreams are needed to model accurately nonlinear cooling curves [15]. Furthermore, the model includes area calculations for economic analysis. No Boolean variables or disjunctive representations are used for handling phase transitions. Instead, the model uses the nonsmooth operators max, min and mid to correctly determine the phase state, enthalpy, and temperature of the process streams. Sensitivities are computed using lexicographic directional (LD-)derivatives, which are extensions of the classical directional derivative to certain classes of nonsmooth functions that provide useful sensitivity information at nonsmooth points. This nonsmooth MHEX model has previously been used for the simulation and optimization of different single mixed refrigerant processes [11,16]. A hybrid modeling approach was used to develop the flowsheet models, where auxiliary equipment and the two-phase stream variables are solved sequentially in nested subroutines. The result is a reduced model size and increased robustness, making it possible to simulate and optimize larger and more complex SMR models. It is, therefore, interesting to investigate whether the novel modeling approach is also capable of simulating DMR processes. In this article, the nonsmooth modeling approach is used to simulate two dual mixed refrigerant processes. The first process is a relatively simple DMR process with cascading PRICO processes, whereas the second is a version of the commercial AP-DMR process [17] with single stage compression. Case studies are performed, each starting from initial points for which the process simulator Aspen Plus fails to obtain thermodynamically feasible solutions.

2. The Nonsmooth Multistream Heat Exchanger Model

The standard two-stream countercurrent heat exchanger is completely described by Equations (1)–(3), which represent the energy balance, minimum approach temperature in the heat exchanger, and total heat exchanger conductance, respectively.

$$m{C}_{\mathrm{p},\mathrm{H}}\left(T_{\mathrm{H}}^{\mathrm{IN}}-T_{\mathrm{H}}^{\mathrm{OUT}}\right)=m{C}_{\mathrm{p},\mathrm{C}}\left(T_{\mathrm{C}}^{\mathrm{OUT}}-T_{\mathrm{C}}^{\mathrm{IN}}\right),$$

(1)

$$\Delta T_{\min }=\min \left\{T_{\mathrm{H}}^{\mathrm{IN}}-T_{\mathrm{C}}^{\mathrm{OUT}},T_{\mathrm{H}}^{\mathrm{OUT}}-T_{\mathrm{C}}^{\mathrm{IN}}\right\},$$

(2)

$$UA=\frac{Q}{\Delta T_{\mathrm{LM}}},$$

(3)

where $m{C}_{\mathrm{p},\mathrm{H}/\mathrm{C}}$ is the heat capacity flowrate for the hot (H) and cold (C) streams, respectively, $UA$ is the heat transfer conductance, $Q\equiv m{C}_{\mathrm{p},\mathrm{H}}\left(T_{\mathrm{H}}^{\mathrm{IN}}-T_{\mathrm{H}}^{\mathrm{OUT}}\right)$ is the total heat exchange duty, and $\Delta T_{\mathrm{LM}}$ is the log-mean temperature difference.

The energy balance can readily be extended to the case of $n_{\mathrm{H}}$ hot and $n_{\mathrm{C}}$ cold streams as follows:

$$\sum _{i=1}^{n_{\mathrm{H}}}mC{p}_{\mathrm{H},i}\left(T_{\mathrm{H},i}^{\mathrm{IN}}-T_{\mathrm{H},i}^{\mathrm{OUT}}\right)=\sum _{j=1}^{n_{\mathrm{C}}}mC{p}_{\mathrm{C},j}\left(T_{\mathrm{C},j}^{\mathrm{OUT}}-T_{\mathrm{C},j}^{\mathrm{IN}}\right).$$

(4)

Also, the equation for the total heat exchanger conductance can be applied to the case of multistream heat exchangers by assuming vertical heat exchange between the hot and cold composite curves

$$UA=\sum _{k=1}^{K-1}\frac{\Delta Q^k}{\Delta T_{\mathrm{LM}}^k},$$

(5)

where K is the total number of enthalpy intervals and $\Delta Q^k$ is the enthalpy change of interval k.

As for the minimum approach temperature constraint, however, it cannot readily be extended to that of multistream heat exchangers. The minimum temperature difference in two-stream heat exchangers occurs at the physical endpoints of the heat exchangers, assuming constant $mCp$ values for the streams. For multistream heat exchangers, on the other hand, the pinch point can occur at any of the stream inlet temperatures and are not necessarily situated at the physical endpoints. Instead, the approach temperature can be calculated using concepts from pinch analysis and heat integration through a pinch location algorithm. Several pinch location algorithms exist in the literature, most of which require a disjunctive optimization problem to be solved to global optimality. However, in order to avoid solving a separate optimization problem for the minimum temperature calculations, Watson et al. [5] developed a reformulation of the Duran and Grossmann model for simultaneous optimization and heat integration [14]. The reformulation solves the pinch problem through the nonsmooth equation

$$\underset{p\in P}{\min \{}EB{P}_{\mathrm{C}}^p-EB{P}_{\mathrm{H}}^p\}=0,$$

(6)

where P is the (finite) set of candidate pinch points and $EB{P}_{\mathrm{H}/\mathrm{C}}^p$ is the enthalpy of the extended hot/cold composite curves for pinch candidate p, as defined in Watson et al. [5]. This equation can be solved using nonsmooth numerical equation solvers, and therefore, no optimization problem needs to be solved.

Multistream heat exchangers are particularly useful in natural gas liquefaction processes where the streams normally undergo phase changes. Phase changes in the heat exchanger present a commonly known modeling issue, as phase boundaries and the phases traversed in the heat exchanger are not known a priori and may change during the simulation. Instead of using Boolean variables or disjunctions to detect the correct phase behavior, the model uses the nonsmooth operators max, min and mid where the mid operator is a function mapping to its median argument [12]. This is done by partitioning each process stream into superheated (sup), two-phase (2p) and subcooled (sub) substreams whose inlet and outlet temperatures are determined by the following equations:

$$T_{\sup }^{\mathrm{in}/\mathrm{out}}=\max \left(T_{\mathrm{DP}},T^{\mathrm{IN}/\mathrm{OUT}}\right),$$

(7)

$$T_{2\mathrm{p}}^{\mathrm{in}/\mathrm{out}}=\mathrm{mid}\left(T_{\mathrm{DP}},T_{\mathrm{BP}},T^{\mathrm{IN}/\mathrm{OUT}}\right),$$

(8)

$$T_{\mathrm{sub}}^{\mathrm{in}/\mathrm{out}}=\min \left(T_{\mathrm{BP}},T^{\mathrm{IN}/\mathrm{OUT}}\right),$$

(9)

where $T^{\mathrm{IN}/\mathrm{OUT}}$ is the inlet or outlet temperature of the process stream, $T_{\mathrm{sub}/2\mathrm{p}/\sup }^{\mathrm{in}/\mathrm{out}}$ is the corresponding inlet or outlet temperature of the substreams, and $T_{\mathrm{DP}}$ and $T_{\mathrm{BP}}$ are the dew and bubble point temperatures of the process stream. Additional stream segments may be used to improve the accuracy of the calculations, which is particularly important for the two-phase region where enthalpy and temperature vary highly nonlinearly due to phase changes. Watson et al. showed, by using the PRICO process [18] as an illustrative example, that 20 segments provided a sufficient level of accuracy for representing the two-phase region [15].

Stream temperatures in the two-phase region are calculated using successive pressure-enthalpy (PQ)-flash operations for the stream segments. As stream propertiesm and thus phase boundaries change during the simulation, the PQ-flash algorithm must be capable of handling instances of single phase flow. This is also an issue for models of auxiliary equipment in LNG processes, such as compressors and valves, which may experience instances of single phase flow during the iterations of the nonsmooth solver. For this application, a nonsmooth extension of the well-known Boston and Britt [19] flash algorithm was developed that handles instances of single phase flow without relying on post-processing methods. The flash algorithm was summarized in a three-paper series by Watson et al. [15,20] and Watson and Barton [21], and employs a mid-function for detecting the correct phase state at the outlet. The algorithm was shown to handle instances near the phase-boundaries robustly, where the conventional inside-out algorithm is prone to failure [20]. A methodology for calculating correct sensitivity information from the flash equations was presented, which allows the nonsmooth inside-out algorithm to be integrated into the flowsheet models. Rather than using a fully equation-oriented framework where the MHEX model and flash calculations are solved simultaneously, the flash calculations are embedded in the model as subroutine calls and solved with the specialized algorithm. The result is a reduced model size and increased robustness, particularly for larger and more complex models.

Solving the MHEX Model

The MHEX model requires a nonsmooth algebraic equation system to be solved. The inclusion of the nonsmooth operators max, min and mid, means that points of nondifferentiability exist where the Jacobian is undefined. A common approach for handling nonsmoothness in the literature is to formulate a smooth approximation for the nonsmooth function around the kinks. However, formulating these approximating functions is nontrivial, and sometimes, poor selection may lead to either an ill-conditioned approximation or the loss of accuracy [22].

Alternatively, extensions of the conventional derivative to certain classes of nonsmooth functions exist. There are several of these generalized derivatives described in the literature, where the Clarke Generalized Jacobian [23] for locally Lipschitz continuous functions is particularly well-known. However, the Clarke Jacobian is challenging to compute, especially for composite functions, as it follows calculus rules (e.g., the chain rule) as inclusions rather than as equations. Therefore, new interest in nonsmooth analysis has centered around a different generalized derivative known as the lexicographic (L-)derivative [24]. The L-derivative is a generalized derivative for functions satisfying the conditions for lexicographic (L)-smoothness as formulated by Nesterov [24]. It was shown later by Khan and Barton [25] that L-derivatives are as useful as elements of the Clarke Jacobian itself in many nonsmooth numerical methods. The renewed attention on the L-derivative was fueled by the development of an automatic procedure for calculating its elements for composite functions. Rather than computing elements of the L-derivative directly, Khan and Barton [13] introduced a nonsmooth extension to the directional derivative, known as the lexicographic directional (LD)-derivative. Unlike the Clarke Jacobian and L-derivatives, the LD-derivative follows calculus rules sharply, and can be computed for composite functions using a nonsmooth analog of the vector forward mode of automatic differentiation [13]. Furthermore, the L-derivative can be obtained directly from the LD-derivative. An extensive review of computing LD-derivatives and their applications is provided by Barton et al. [26].

Nonsmooth Newton-type solvers, where generalized derivative information is used in place of the function’s Jacobian, exist for solving nonsmooth equation systems:

$$\mathbf{G}\left({\mathbf{x}}^k\right)({\mathbf{x}}^{k+1}-{\mathbf{x}}^k)=-\mathbf{f}\left({\mathbf{x}}^k\right).$$

(10)

Here, $\mathbf{G}\left({\mathbf{x}}^k\right)$ is an element of a generalized derivative of $\mathbf{f}$ at ${\mathbf{x}}^k$. Equation (10) solves for the next iterate ${\mathbf{x}}^{k+1}$ provided $\mathbf{G}\left({\mathbf{x}}^k\right)$ is nonsingular at ${\mathbf{x}}^k$. However, singular generalized derivative elements may occur in the MHEX model from residuals of the form $\min \left\{0,y\right\}$ used to solve Equation (6). A Newton-type solver that is also applicable to singular generalized derivative elements is the linear programming (LP) Newton method by Facchinei et al. [27]:

$$\begin{array}{cc}\hfill \underset{\gamma ,\mathbf{x}}{\min }& \gamma \hfill \\ \hfill \mathrm{s}.\mathrm{t}.& {∥\mathbf{f}\left({\mathbf{x}}^k\right)+\mathbf{G}\left({\mathbf{x}}^k\right)(\mathbf{x}-{\mathbf{x}}^k)∥}_{\infty }\le \gamma min\left({∥\mathbf{f}\left({\mathbf{x}}^k\right)∥}_{\infty },{∥\mathbf{f}\left({\mathbf{x}}^k\right)∥}_{\infty }^2\right),\hfill \\ & {∥(\mathbf{x}-{\mathbf{x}}^k)∥}_{\infty }\le \gamma {∥\mathbf{f}\left({\mathbf{x}}^k\right)∥}_{\infty },\hfill \\ & \mathbf{x}\in X,\hfill \end{array}$$

(11)

where X is a polyhedral set (e.g., bounds) on the problem and $\gamma$ is a supplementary variable that drives convergence towards the solution. The LP-Newton solves a linear program upon each iteration, where the next iterate ${\mathbf{x}}^{k+1}$ is given by the $\mathbf{x}$ part of the solution. Although the program in Equation (11) is also applicable for ill-conditioned generalized derivatives, previous experience with simulating single mixed refrigerant processes has shown that including a backtracking line search [28] significantly improves the step quality of the LP-Newton method at singular points [16]. Nevertheless, solving a linear program on every iteration is comparatively more expensive than computing a step with the method in Equation (10). Therefore, the model includes a hybrid solution where the LP-Newton is invoked only when the generalized derivative is either singular or ill-conditioned.

3. Simulation Cases and Results

Previously, the developed nonsmooth simulation tool was used to analyse different SMR processes under conditions for which the commercial simulator Aspen Plus failed to obtain results [16]. In particular, the additional two unknowns computed by Equations (5) and (6) add versatility, making it possible to obtain feasible operating points in cases where more than one operating parameter is unknown to the designer. This article develops simulation models for two different DMR processes. The first design constitutes a simple design where the PRICO process is used for both the warm mixed refrigerant (WMR) and the cold mixed refrigerant (CMR) cycles. Furthermore, a natural gas liquid (NGL) separator is added for the extraction of heavier hydrocarbons. The second DMR process is a version of the AP-DMR process. The main focus here is on the liquefaction part of the process, and thus, not the compression scheme. Therefore, only a single-stage compression with aftercooling is included in the model. Different case studies are presented, each solving for different sets of unknown variables. The variables considered in the analysis are the high and low pressure levels, refrigerant flowrates, inlet and outlet temperatures from the MHEXs, and refrigerant compositions for the warm and cold mixed refrigerant, as well as the MHEX specifications, i.e., the minimum temperature difference and the heat exchanger conductance. The following nomenclature is used for the parameters and unknown variables in the models:

• Pressure level of the (warm/cold) high pressure refrigerant: $P_{\mathrm{HP}(\mathrm{W}/\mathrm{C})}$.
• Pressure level of the (warm/cold) low pressure refrigerant: $P_{\mathrm{LP}(\mathrm{W}/\mathrm{C})}$.
• Inlet/outlet temperatures of the high pressure refrigerants (equal to the natural gas stream): $T_{\mathrm{HP}}^{\mathrm{IN}/\mathrm{OUT}}$.
• Inlet/outlet temperatures of the low pressure refrigerant streams: $T_{\mathrm{LP}}^{\mathrm{IN}/\mathrm{OUT}}$.
• Molar flowrate of the (warm/cold) refrigerants: $F_{(\mathrm{W}/\mathrm{C})}$.
• Molar flowrate of component i in (warm/cold) refrigerants: $f_{(\mathrm{W}/\mathrm{C}),i}$.

The models were written using Julia v0.6.0 (Julia Lab – Massachusetts Institute of Technology, Cambridge, MA, USA) and run on a Dell Latitude E5470 laptop in the Ubuntu v16.10 environment with an Intel Core i7-6820HQ CPU at 2.7 GHz and 8.2 GB RAM. Simulations were done for the Peng–Robinson Equation of State where the property parameters were taken from Aspen Plus [2]. To ensure an accurate representation of the process streams, the single-phase substreams were partitioned into five stream segments for each MHEX. Similarly, two-phase substreams were represented by 20 stream segments, which was shown to be sufficiently accurate for capturing the nonlinearities in this region by Watson et al. [15] The overall flowsheet convergence tolerance was set to ${∥\mathbf{y}∥}_{\infty }<{10}^{-5}$, where $\mathbf{y}$ represents the equation residuals, whereas the tolerance for the individual flash calculations was selected to be ${∥\mathbf{y}∥}_{\infty }<{10}^{-8}$.

4. Example 1

The first DMR process studied in this paper is a simple configuration with cascading PRICO cycles for the warm and cold mixed refrigerants. This cycle consists of two MHEXs as well as an NGL separator for the extraction of heavier hydrocarbons. Heavier hydrocarbons freeze out at cryogenic temperatures and can cause plugging in process equipment. Furthermore, they represent a valuable commodity if sold separately. The feed gas is sent to the process at 295.15 K where it is precooled by a warm mixed refrigerant cycle consisting of ethane, propane, and n-butane. The feed gas and the refrigerants exit the heat exchanger at a temperature of $T_{\mathrm{HP},1}^{\mathrm{OUT}}$. The feed gas is then sent to the NGL separator, where heavier hydrocarbons are extracted for further fractionation and/or export, before the gas enters MHEX 2 for liquefaction. The cold mixed refrigerant consists of a lighter refrigerant mixture consisting of nitrogen, methane, ethane, and propane for the liquefaction of the natural gas. Along with the feed gas, the CMR is precooled in the WMR PRICO cycle before it enters the cold heat exchanger. A process flowsheet of the DMR process is presented in Figure 1.

The parameter values and initial guess values for the unknown variables are provided in

Table 1. Multistream heat exchanger (MHEX) and refrigerant stream data for Example 1. For unknown variables, the value listed is an initial guess.

Property	Value	Property	Value
$\eta$	0.8	$U{A}_2$ (MW/K)	3.0
$\Delta T_{\min ,1}$ (K)	3.0	$\Delta T_{\min ,2}$ (K)	3.0
$F_{\mathrm{W}}$ (kmol/s)	1.65	$F_{\mathrm{C}}$ (kmol/s)	1.55
$P_{\mathrm{HPW}}$ (MPa)	1.67	$P_{\mathrm{HPC}}$ (MPa)	4.30
$P_{\mathrm{LPW}}$ (MPa)	0.42	$P_{\mathrm{LPC}}$ (MPa)	0.25
$T_{\mathrm{HP},1}^{\mathrm{OUT}}$ (K)	245.15	$T_{\mathrm{HP},2}^{\mathrm{OUT}}$ (K)	120.15
$T_{\mathrm{LP},1}^{\mathrm{OUT}}$ (K)	290.15	$T_{\mathrm{LP},2}^{\mathrm{OUT}}$ (K)	240.15
Composition (mol %):
Ethane	47.83	Nitrogen	10.00
Propane	34.17	Methane	43.80
n-Butane	18.00	Ethane	35.20
		Propane	11.00

. The parameter values were selected such that Aspen Plus failed to converge to a feasible solution using its standard MHEX model with one equation. Essentially, the Aspen Plus model only solves the overall energy balance in Equation (4) for a single unknown temperature (chosen here as the inlet inlet temperatures to the compressors). Refrigerant compositions and pressures cannot be handled as unknown variables in these models. This results in less versatility in cases where pressure, compositions, and multiple temperatures must be adjusted to find a feasible design.

The model contains 61 unknowns, five of which are provided by the solution of the MHEX equations. The remaining 56 variables are the temperatures of the intermediate stream segments for the single phase regions. As flash calculations are decoupled and solved separately, the model size is independent of the number of stream segments in the two-phase region. A two-equation model is used for MHEX 1 as the $UA$-value depends on the stream results from MHEX 2. Fixing the area prior to simulation can therefore be challenging. Therefore, the heat exchanger conductance value is instead calculated through post-processing. The simulations were carried out for a rich feed gas composition with 1.00 mol % nitrogen, 85.60 mol % methane, 4.93 mol % ethane, 3.71 mol % propane, 2.90 mol % n-butane, 1.30 mol % i-butane, and 0.56 mol % n-pentane at a pressure of 4 MPa and flowrate of 1.0 kmol/s where both the refrigerant mixtures and the feed gas enter the precooling MHEX at a temperature of 295.15 K, as indicated in Figure 1. Two simulation cases were constructed, solving for different sets of unknown variables:

• Case I:$P_{\mathrm{LPW}}$, $P_{\mathrm{HPC}}$, $f_{\mathrm{W},\mathrm{propane}}$, $T_{\mathrm{LP},2}^{\mathrm{OUT}}$, $U{A}_2$.
• Case II:$P_{\mathrm{HPW}}$, $P_{\mathrm{LPC}}$, $F_{\mathrm{C}}$, $\Delta T_{\min ,1}$, $\Delta T_{\min ,2}$.

Case I. Solved for a variable WMR composition, an unknown inlet temperature to the CMR compressor, the heat exchanger conductance in MHEX 2, as well as the low pressure level $P_{\mathrm{LPW}}$ and high pressure level $P_{\mathrm{HPC}}$ of the warm and cold mixed refrigerants, respectively. The refrigerant composition was changed in the model by varying the component molar flowrate of propane $f_{\mathrm{W},\mathrm{propane}}$. A solution was obtained after four iterations and a total simulation time of 62.7 s, including initialization. The model converged to a solution with $P_{\mathrm{LPW}}=0.57$ MPa, $P_{\mathrm{HPC}}=6.53$ MPa and $T_{\mathrm{LP},2}^{\mathrm{OUT}}=242.15$ K. The design resulted in a total compression work of 21.33 MW, with heat exchanger conductance values of $U{A}_1=2.69$ MW/K and $U{A}_2=1.89$ MW/K. The work distribution of the two compressors was 16.79 MW for compressing the CMR and 4.54 MW for compressing the WMR. A new WMR composition was obtained consisting of 53.30 mol % ethane, 26.64 mol % propane and 20.06 mol % n-butane with a corresponding molar flowrate of 1.48 kmol/s. Figure 2 presents the composite curves and driving force plot for the solution.

Case II. Solved for the flowrate of the CMR, the minimum approach temperatures in both MHEXs, and the high pressure and low pressure levels of the warm and cold refrigerant mixtures, respectively. A solution was obtained after 13 iterations with $P_{\mathrm{HPW}}=1.58$ MPa, $P_{\mathrm{LPC}}=0.388$ MPa, $F_{\mathrm{C}}=1.952$ kmol/s, and minimum approach temperatures of 5.00 K and 3.15 K for MHEX 1 and 2, respectively. The design resulted in a total compressor work of 20.56 MW, and a heat exchanger conductance value $U{A}_1$ of 2.04 MW/K. Compressing the CMR required a total of 14.10 MW, whereas compressing the WMR required only 6.46 MW. The total simulation time, including initialization of the model, was 102.5 s. The composite curves and driving force plots for the process are presented in Figure 3. As can be seen, both solutions have a similar trend, although the cold low pressure refrigerant superheating is noticeably larger in Case I. This results in lower driving forces in MHEX 1, but also introduces a larger temperature difference at the cold end of the process, leading to a higher compression power. A discussion on superheating and its effect on design of DMR processes was made by Kim and Gundersen [29].

5. Example 2

A more complex DMR design is presented in Figure 4. Instead of having two cascading PRICO cycles, this design features a spiral-wound heat exchanger (SWHX) for the CMR where the refrigerant is separated after precooling to provide cooling at different temperature levels. The vapor product, which consists mainly of the light components nitrogen and methane, is liquefied and subcooled to provide cooling in the cold MHEX 3. At the same time, the liquid product is subcooled and mixed with the low pressure refrigerant from MHEX 3, where it is used to cool the feed gas in the intermediate MHEX 2. Since part of the refrigerant mixture only circulates in the warm end of the SWHX, the overall molar flowrate and thus the required heat transfer area decrease.

The refrigerant streams and MHEX data for the DMR model are given in

Table 2. MHEX and refrigerant stream data for Example 2. For unknown variables, the value listed is an initial guess.

Property	Value	Property	Value
$\eta$	1.0	$\Delta T_{\min ,1}$ (K)	4.0
$U{A}_3$ (MW/K)	0.3	$\Delta T_{\min ,2}$ (K)	11.0
		$\Delta T_{\min ,3}$ (K)	4.0
$F_{\mathrm{W}}$ (kmol/s)	1.55	$F_{\mathrm{C}}$ (kmol/s)	1.45
$P_{\mathrm{HPW}}$ (MPa)	1.67	$P_{\mathrm{HPC}}$ (MPa)	4.85
$P_{\mathrm{LPW}}$ (MPa)	0.50	$P_{\mathrm{LPC}}$ (MPa)	0.25
$T_{\mathrm{HP},1}^{\mathrm{OUT}}$ (K)	240.15	$T_{\mathrm{HP},2}^{\mathrm{OUT}}$ (K)	170.15
$T_{\mathrm{HP},3}^{\mathrm{OUT}}$ (K)	120.15	$T_{\mathrm{LP},1}^{\mathrm{OUT}}$ (K)	280.15
$T_{\mathrm{LP},2}^{\mathrm{OUT}}$ (K)	230.15	$T_{\mathrm{LP},3}^{\mathrm{OUT}}$ (K)	145.15
Composition (mol %):
Ethane	47.83	Nitrogen	7.00
Propane	34.17	Methane	41.80
n-Butane	18.00	Ethane	33.20
		Propane	18.00

. Again, the parameter values were selected such that a solution could not be obtained using the commercial simulation tool Aspen Plus. A feed gas composition with 1.00 mol % nitrogen, 91.60 mol % methane, 4.93 mol % ethane, 1.71 mol % propane, 0.35 mol % n-butane, 0.40 mol % i-butane, and 0.01 mol % i-pentane at a pressure of 5.5 MPa and a flowrate of 1.0 kmol/s was used in the simulations.

The simulation model with three MHEXs consists of 96 variables and exhibits seven unknowns. Again these unknowns may be used to solve for any process stream variable such as the pressure, composition, flowrate, or temperature, as well as important MHEX data such as the minimum temperature difference and the heat exchanger conductance. Specifying the $UA$-value in MHEXs 1 and 2 is challenging as it depends on the solution of MHEX 3. Therefore, as for Example 1, the $UA$-values are calculated during post-processing to make problem specification easier. Two simulation cases were constructed solving for the following sets of unknown variables:

• Case I:$P_{\mathrm{LPW}}$, $P_{\mathrm{HPW}}$, $P_{\mathrm{LPC}}$, $T_{\mathrm{HP},2}^{\mathrm{OUT}}$, $T_{\mathrm{LP},3}^{\mathrm{OUT}}$, $F_{\mathrm{C}}$, $\Delta T_{\min ,2}$.
• Case II:$P_{\mathrm{LPW}}$, $P_{\mathrm{HPC}}$, $T_{\mathrm{HP},2}^{\mathrm{OUT}}$, $T_{\mathrm{LP},3}^{\mathrm{OUT}}$, $f_{\mathrm{W},\mathrm{ethane}}$, $\Delta T_{\min ,2}$, $U{A}_3$.

Case I: Solved for both pressure levels of the WMR, the low pressure level and refrigerant flowrate of the CMR, the feed gas and high pressure refrigerant temperatures out of MHEX 2, the low pressure refrigerant temperature out of MHEX 3, as well as the minimum approach temperature in MHEX 2. A solution was obtained after six iterations and a total simulation time of 83.0 s with $P_{\mathrm{LPW}}=0.43$ MPa, $P_{\mathrm{HPW}}=1.62$ MPa, $P_{\mathrm{LPC}}=0.27$ MPa, $T_{\mathrm{HP},2}^{\mathrm{OUT}}=155.34$ K, $T_{\mathrm{LP},3}^{\mathrm{OUT}}=151.34$ K, $F_{\mathrm{C}}=1.42$ kmol/s, and $\Delta T_{\min ,2}=6.68$ K. The $UA$-values were calculated to be 1.99 MW/K and 2.12 MW/K for MHEXs 1 and 2, respectively. The obtained feasible design resulted in a combined compression work of 14.40 MW, where 9.76 MW was needed to compress the CMR, and 4.64 MW was used to compress the WMR. Figure 5 presents the composite curves and driving force distribution in the MHEXs at the solution.

Case II: Solved for the low pressure level of the WMR, high pressure level of the CMR, the natural gas and high pressure refrigerant temperatures out of MHEX 2, the low pressure refrigerant temperature out of MHEX 3, the composition of the WMR, the minimum temperature difference in MHEX 2, and the heat exchanger conductance value for MHEX 3. The model converged after three iterations and a total simulation time of 64.2 s to a solution with $P_{\mathrm{LPW}}=0.44$ MPa, $P_{\mathrm{HPC}}=4.59$ MPa, $T_{\mathrm{HP},2}^{\mathrm{OUT}}=161.03$ K, $T_{\mathrm{LP},3}^{\mathrm{OUT}}=157.03$ K, $\Delta T_{\min ,2}=7.88$ K, and $U{A}_3=0.33$ MW/K. A new WMR composition was obtained with 49.24 mol % ethane, 33.24 mol % propane, and 17.51 mol % n-butane and a total molar flowrate of 1.59 kmol/s. The feasible design required a total compression power of 14.85 MW, where 10.08 MW was spent compressing the CMR and 4.76 MW was used to compress the WMR. The heat exchanger conductance values were calculated during post-processing to be $U{A}_1=2.02$ MW/K and $U{A}_2=1.84$ MW/K, respectively. The composite curves and driving force plots for the process are presented in Figure 6.

Case III: Included an NGL separator for the extraction of heavier hydrocarbons (see Figure 7). The case solved for the same set of variables as in Case I and with the same initial guesses and parameter values as given in Table 2. As for the previous examples, the model was solved from an initial guess at which Aspen Plus obtained no feasible solutions with the built-in MHEX module. Rich feed gas compositions at a reduced pressure of 4.0 MPa were used to ensure adequate separation. Simulations were carried out at three different compositions with varying methane contents (

Table 3. Natural gas compositions for Case III.

	Compostion I:	Composition II:	Composition III:
Nitrogen (mol %)	2.00	2.00	2.00
Methane (mol %)	85.60	87.60	89.60
Ethane (mol %)	6.93	5.93	4.93
Propane (mol %)	3.71	2.71	1.71
n-Butane (mol %)	1.35	1.35	1.35
i-Butane (mol %)	0.40	0.40	0.40
i-Pentane (mol %)	0.01	0.01	0.01

).

Driving force distributions for the three solutions are provided in Figure 8. Solutions were obtained for all three cases within a few iterations. The first two feed gas compositions converged after seven iterations and total simulation times of 85.6 s and 86.4 s for compositions I and II, respectively. The third case converged after six iterations and a total simulation time of 82.1 s. All three solutions exhibited similar driving force profiles, with temperature differences varying mainly in the intermediate MHEX. The same trend can also be seen from the simulation results in

Table 4. Simulation results for the DMR process with NGL extraction.

	Composition I:	Composition II:	Composition III:
Total work (MW)	14.99	15.11	15.13
$W_{\mathrm{C}}$ (MW)	10.26	10.40	10.43
$W_{\mathrm{H}}$ (MW)	4.73	4.71	4.70
$U{A}_1$ (MW/K)	2.07	2.05	2.03
$U{A}_2$ (MW/K)	2.63	2.87	3.41
$P_{\mathrm{LPW}}$ (MPa)	0.43	0.43	0.43
$P_{\mathrm{HPW}}$ (MPa)	1.66	1.66	1.65
$P_{\mathrm{LPC}}$ (MPa)	0.26	0.25	0.25
$T_{\mathrm{HP},2}^{\mathrm{OUT}}$ (K)	158.22	157.95	157.50
$T_{\mathrm{LP},3}^{\mathrm{OUT}}$ (K)	154.22	153.95	153.50
$F_{\mathrm{C}}$ (kmol/s)	1.46	1.47	1.48
$\Delta T_{\min ,2}$ (K)	3.69	3.14	2.20
LNG composition (mol %):
Nitrogen	2.08	2.05	2.03
Methane	88.04	89.24	90.61
Ethane	6.50	5.68	4.81
Propane	2.67	2.15	1.48
n-Butane	0.52	0.65	0.81
i-Butane	0.19	0.23	0.27
i-Pentane	0.00	0.00	0.00
LNG flowrate (kmol/s)	0.96	0.97	0.98

, where the main differences between the three solutions are the $UA$ and $\Delta T_{\min }$ values for the intermediate MHEX.

6. Conclusions

This paper successfully simulated two different dual mixed refrigerant processes using recent modeling advances in the field of nonsmooth analysis. Natural gas liquefaction processes are particularly difficult to simulate with conventional process simulators, such as Aspen Plus. Lacking rigorous checks for preventing temperature crossovers means that the user must often resort to an iterative approach to obtain feasible designs. A multistream heat exchanger model that features rigorous checks for preventing temperature crossovers during simulation was included in the models. It includes nonsmooth formulations for pinch point location and phase regime detection in the heat exchangers. As a result, it offers a compact formulation where the flowsheet models can be simulated by solving a set of nonsmooth equations, rather than solving a local or global optimization problem. All the cases successfully converged within a total simulation time of 100 s, including initialization of the models. Furthermore, each case was constructed such that Aspen Plus failed to obtain a feasible solution using the same starting point and initial guess values as the nonsmooth model. The use of additional unknowns and the possibility of varying ostensibly difficult process stream variables, such as compositions makes the tool more versatile and capable of handling even complex LNG liquefaction processes.