Scenario-Free Optimal Design under Uncertainty of the PRICO Natural

Subscriber access provided by READING UNIV

Article

Scenario-Free Optimal Design under Uncertainty of the PRICO® Natural Gas Liquefaction Process Calvin Tsay, and Michael Baldea Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.7b03634 • Publication Date (Web): 16 Jan 2018 Downloaded from http://pubs.acs.org on January 20, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Industrial & Engineering Chemistry Research is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Scenario-Free Optimal Design under Uncertainty of the PRICO Natural Gas Liquefaction Process R

Calvin Tsay and Michael Baldea∗ McKetta Department of Chemical Engineering, The University of Texas at Austin 200 East Dean Keeton St., Stop C0400, Austin, TX 78712 November 20, 2017

Abstract Liquefaction accounts for a large portion of the cost of liquefied natural gas (LNG), motivating significant efforts in the simulation and optimization of liquefaction processes. While most optimization studies typically assume a fixed feed stream concentration, the properties of natural gas feeds can vary in time or between different sources. This constitutes a significant exogenous uncertainty, particularly for field or offshore liquefaction facilities with limited feed-blending capabilities. This paper examines the R plant as a representative mixed-refrigerant liquefaction process and presents PRICO

a design optimization formulation for natural gas liquefaction processes with varying feed composition. We introduce models for the multistream heat exchanger and for the compressor that allow for the size of key equipment to be optimized along with the process operating conditions. Using a recently-developed scenario-free framework for optimization under parametric uncertainty, our study provides new insights (including some counterintuitive results) concerning optimal liquefaction equipment sizing for dealing with uncertainty in feed composition. ∗

corresponding author: [email protected]

1 Environment ACS Paragon Plus

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Keywords Optimization under uncertainty, process design, pseudo-transient modeling, natural gas liquefaction 1

Introduction

Increased interest in natural gas as a safe and clean fuel source (U.S. liquefaction capacity is expected to grow by over 500% from 2016-2019 1) has led to considerable efforts in improving the efficiency of natural gas liquefaction processes. Liquefied natural gas (LNG) has approximately 600 times the density of gaseous natural gas 2 , allowing for economic transportation of this resource over long distances. Nevertheless, cooling and liquefying natural gas require significant refrigeration and represent around 52% of the cost of LNG 3 . Accordingly, many efforts have been dedicated to improving the energy efficiency of liquefaction processes, including LNG flowsheet simulation and optimization advances 4;5;6;7;8 by the process systems engineering community. Modeling, simulation, and optimization of natural gas liquefaction processes are complicated by the need to account for phase transitions and recycle streams present in refrigeration cycles 4 . Moreover, such processes typically incorporate multistream heat exchangers (MHEXs), which allow for thermal contact between multiple streams in a single process unit, but further complicate numerical simulation 7;8 . Due to the discontinuities inherently present in MHEX models 8 (e.g. capturing phase transitions, transport and transfer correlations dependent on flow regime), flowsheet convergence remains an obstacle in the optimization of liquefaction processes, and many proposed approaches employ simplified models in an attempt to balance simulation accuracy and robustness 6 . Many such simplified approaches to modeling MHEXs involve the introduction of integer variables, resulting in mixed-integer, nonlinear models that may be challenging to use in optimization. R R Single mixed refrigerant (SMR) liquefaction processes, such as the PRICO (PRICO

is a registered trademark of Black & Veatch Holding Company) process, have relatively simple designs and small footprints 6 , making them cost-effective for offshore and stranded LNG plants. Many previous works have studied the optimization of SMR processes: Lee et al. 9 found that the composition of the refrigerant is the most significant variable for SMR


Page 2 of 35

Page 3 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


process design; Aspelund et al. 10 used the Nelder-Mead algorithm for local searches in a simulation-based optimization scheme; Tak et al. 11 compared the optimized designs for various compression configurations; Khan and Lee 6 implemented a particle-swarm optimization paradigm; and Jensen 12 considered a nonlinear objective function to balance capital and operating cost tradeoffs. In most of these studies, the specific work of the liquefaction cycle (i.e., the amount of energy required per unit mass of LNG produced) plays a large role in the objective function. These studies generally consider steady-state optimization for a known natural gas feed 5;7 or a few scenarios involving different natural gas feed compositions 10;11 . While assuming fixed natural gas composition(s) facilitates the comparison among specific works associated with optimal SMR process designs reported previously 5;7;8;13 , in reality natural gas feed compositions can be uncertain and/or change in time, depending on the reservoir, and liquefaction systems may not always operate at nominal design conditions. Consideration of variations in feed composition is especially important in the natural gas supply chain, where products are typically sent to customers with little or no processing 14 . Recognizing these challenges, Li et al. 14 considered the optimization of natural gas production systems at the infrastructure level, using a scenario-based stochastic programming approach and a nonconvex decomposition method to maximize the expected process profit. The benefits that were identified for designing the natural gas infrastructure for variable feed compositions suggest that the design of LNG liquefaction processes should also consider this uncertainty. Optimization under uncertainty often results in process designs that can differ considerably from the nominal case 15;16 . Nevertheless, feed concentration uncertainty has rarely been considered in the optimization of SMR liquefaction processes and remains an open research question. R Motivated by the above, in this work we investigate the design of the PRICO natural

gas liquefaction process under feed composition uncertainty, which we model by treating the concentration of the natural gas to be liquefied as a parametric uncertainty. We use R a recently proposed scenario-free dynamic optimization-based scheme 15 and the PRICO

process as a representative, single mixed-refrigerant liquefaction system. 3 Environment ACS Paragon Plus


The novelty of this contribution consists of: • A pseudo-transient compressor model that can be seamlessly integrated with other pseudo-transient process unit models 17 . The model is based on compressor dynamics/compressor curves to account for variable efficiency and instability (i.e. operation past the surge point) at different operating points • A flowsheet representation that allows for equipment (compressor and MHEX) sizing and operation to be considered simultaneously during optimization, which is not typically available in conventional steady-state flowsheet models. • An investigation, using a scenario-free dynamic optimization approach for parametric uncertainty 15 , of the effects of feed concentration deviations from nominal values on R the optimal design of the PRICO process.

We consider the optimization of the process flowsheet with both single-stage and twostage compression configurations. In addition, we consider optimization of i) operating cost only (for comparison with other studies 5;7 ) and ii) a nonlinear objective function considering both operating and capital costs, which has been shown to produce more realistic solutions 18 . 2

Liquefaction Process Model

NG liquefaction processes convert natural gas into LNG using a combination of heat exchangers, compressors, flash tanks, and valves 4;10 . They are generally categorized as (i) mixed-refrigerant processes, based on adjusting coolant composition to match the LNG cooling curve, (ii) cascade processes, employing multiple refrigeration cycles to limit the mean temperature difference, and (iii) mixed-fluid cascade processes, combining both techniques to further improve energy efficiency. In the first two categories, Figure 1 shows a cascaded LEC (Liquefied Energy Chain) al. 13;19 and Figure 5 shows two flowsheets adapted from the R industrial liquefaction process. As illustrated by these exammixed-refrigerant PRICO

ples, compressors and MHEXs are typically the central units for LNG processing plants, and their careful, optimization-oriented design is crucial to maximize process energy efficiency (minimize specific work). 4 Environment ACS Paragon Plus

Page 4 of 35

Page 5 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


MHEX modeling and optimization have received considerable attention in the literature. However, optimization with MHEX models is challenging because the stream phases and phase boundaries can change between flowsheet optimization iterations, especially for mixed-refrigerant processes. Previous approaches include capturing phase change effects and variable heat capacities through disjunctive programming and smoothing functions 5 or through superstructures of two-stream heat exchangers 13 . For accurate simulation outside of an optimization framework, a nonsmooth formulation approach 20 was proposed, which allows phase changes and nonideal thermodynamics to be modeled and solved robustly. In our recent contributions 7;8 , we showed that MHEX mathematical models (of various levels of complexity) can be seamlessly incorporated into equation-oriented flowsheet simulation and optimization using a pseudo-transient continuation approach.

LNG -164.1 °C

LN2

MHEX 1

CO2 MHEX 2 CO2

LN2

HX

NG 70 bar 15 °C

Figure 1: Process flow diagrams for a LEC liquefaction process 19 . The LEC process employs CO2 from industrial sources and N2 from a connected air separation plant as refrigerants. The outlet CO2 stream is typically used for enhanced oil recovery; N2 is released to the atmosphere.

2.1

Pseudo-Transient Process Modeling

In this work, process models are solved using a pseudo-transient continuation approach 17 , which has shown promise in improving convergence properties of flowsheets incorporating detailed models or complex process units such as multistream heat exchangers (MHEXs) 7;8 . This section contains a brief overview of pseudo-transient process modeling concepts as applicable to the modeling of LNG liquefaction process units and the optimization of the 5 Environment ACS Paragon Plus


Page 6 of 35

R PRICO liquefaction process. A more thorough discussion, including stability analyses 7;17

and the reformulation of distributed process models 8;21;22 , is available elsewhere. Originally proposed for deterministic, steady-state process simulation and optimization 7;17;21;23 , pseudo-transient process models improve the convergence properties of process flowsheets by converting a subset of the model algebraic variables x into differential variables xdif f , with the remaining variables xalg remaining algebraic. Providing initial conditions for the differential variables thus replaces providing initial guesses for all variables to an algebraic solver. Mathematically, the steady-state model of a process flowsheet is a set of (nonlinear) algebraic equations usually consisting of mass balances, energy balances, equations of state, reaction/transport kinetics, and other constitutive relations:

f(d, z, x, θ) = 0

(1)

where d is the vector of design variables, z is the vector of recourse variables (that can be changed in response to the realization of uncertainty), x is the vector of process variables, and θ is the vector of process parameters. In a pseudo-transient model, the algebraic equations (1) are systematically 17;21 converted to a system of differential-algebraic equations (DAEs) of the form: ˆf(x˙ dif f , xdif f , xalg , z, θ, t) = 0

(2)

g(xdif f , xalg , z, θ) = 0

(3)

where ˆf comprise the differential equations and g the algebraic equations of the DAE system, defined such that, given initial conditions for the differential variables xdif f (t = 0) = x0dif f , the consistent initialization of the DAE system only requires the solution of a full-rank, linear system of equations 17 . The pseudo-transient dynamics described by ˆf have the same steady-state solution as the original process flowsheet (1), i.e., the models are statically equivalent 17 . This is achieved by modulating the dynamic mass/energy balances with artificial, user-defined time constants 21


Page 7 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


or by introducing dynamics proportional to the residual of equality constraints in (1), resulting in solution of the original equality constraint at steady state 8;17 . Based on the above, the process steady state is obtained by integrating the pseudo-transient model to steady state; the pseudo-transient dynamics are not representative of physical phenomena, but rather are solely a mathematical tool for solving the original system (1). With these general reformulation concepts in mind, mathematical models for MHEXs and compressors based on pseudo-transient continuation will be presented below. In particular, we discuss a pseudo-transient MHEX modeling approach 7 that allows for approximation of the heat-exchange area, and we then present a novel pseudo-transient compressor model based on dynamics that allows for compressor size to be considered during flowsheet optimization. 2.2

Pseudo-Transient Multistream Heat Exchanger Model

Multistream heat exchangers (MHEXs) are central to natural gas liquefaction processes, as they improve process efficiency by allowing multiple hot and cold streams to exchange heat simultaneously in a single physical device. MHEXs introduce several challenges into the modeling and simulation of LNG processes, including: i) the minimum temperature approach constraints between any two streams cannot be written explicitly in terms of the inlet and outlet stream temperatures, ii) the temperature driving force along MHEXs is often very small, and iii) the physical properties of the streams (notably, heat capacity) cannot be assumed to be constant over the typically wide temperature range of the streams present in the exchanger. The above challenges can be overcome by discretizing the exchanger into a finite number of enthalpy segments and reformulating the resulting mathematical model using pseudotransient techniques 7;20 . As is often true in practical application, the ordering of MHEX inlet and outlet streams in terms of temperature magnitude is assumed to be invariant and known prior to simulation and optimization. For example, if the inlet temperature of stream A is initially lower than the inlet temperature of stream B, then it will always (during a simulation or a full optimization procedure) be lower. The cold stream inlet and outlet temperatures comprise the cold temperature sequence, and those of the hot streams comprise 7 Environment ACS Paragon Plus


Page 8 of 35

the hot temperature sequence. Once the number of points (where each point corresponds to an individual inlet or outlet stream temperature) in the cold temperature sequence NC and the number of points in the hot temperature sequence NH are known, the number of heat exchange intervals NHX is given by

NHX = NC + NH − 3

(4)

Each heat exchange interval contains a fixed set of cold and hot streams, separated from its surrounding intervals by a stream inlet or outlet. The continuity of stream properties between consecutive intervals is ensured by additional equations (6)–(7), written for each heat exchange interval (enthalpy interval) comprising inlet and outlet conditions for the set of streams (Si ) present in each enthalpy interval (i ∈ INT ). INT denotes the set of all enthalpy intervals.

Qi =

X

c∈Si

out (Hc,i − Hcini ) =

X

in (Hh,i − Hhout ) i

(5)

in Hcout = Hc,i+1 , ∀ c ∈ (Si ∩ Si+1 ) i

(6)

out Hhini = Hh,i+1 , ∀ h ∈ (Si ∩ Si+1 )

(7)

h∈Si

where index c ∈ Ci represents the set of cold streams in enthalpy interval i, index h ∈ Hi represents the set of hot streams (Hi = Si \ Ci ), and Qi represents the total enthalpy exchanged in enthalpy interval i. Each heat exchange interval (enthalpy interval) i ∈ INT is then discretized into a finite number of enthalpy segments Ni , as shown in Figure 2. Previous works 7;20 suggest that the model accuracy is dependent on the number of enthalpy segments selected for each interval, and more enthalpy segments should be used when the dependence of heat capacity on temperature and/or the likelihood of phase changes are significant. The enthalpies of the cold and hot composite curves, denoted respectively as


Page 9 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


II

I

C1in

H2out

C1out

H1out

H1in H2in Segment boundary Interval boundary

Figure 2: Segmentation of an example MHEX into enthalpy intervals and enthalpy segments. The MHEX has a single cold stream and two hot streams. The hot streams share a common inlet temperature, and the resultant MHEX has two associated enthalpy intervals (I & II).

Hci and Hhi , can then be calculated at each segment zi = [0, 1, ..., Ni ] by X Qi in (zi ) + Hc,i Ni c∈Si X Qi out Hhi (zi ) = (zi ) + Hh,i Ni h∈S Hci (zi ) =

(8) (9)

i

The pressure drop for each stream ∆Ps is assumed to be known at the design stage and to vary linearly with heat duty (the pressure drop for a segment i is computed as ∆Ps,i /∆Ps = P Qi / Qi ). The model as presented above only provides the enthalpies of the hot and cold streams at each enthalpy segment. To calculate the composite curve temperatures at each

segment, we assume that a physical property package is available to evaluate enthalpy as a function of stream temperature, pressure, and composition (noting that this functionality is available in practically all commerical flowsheet simulators and physical property packages):    H L (T, P, F, x) T ≤ Tbub    H pp (T, P, F, x) = H 2p (T, P, F, x) Tbub < T < Tdew      H V (T, P, F, x) T ≥ Tdew

(10)

where T , P , F , and x are respectively the temperature, pressure, flowrate, and composition for which enthalpy is evaluated. The superscripts L, 2p, and V respectively denote the liquid, two-phase, and vapor regimes. 9 Environment ACS Paragon Plus


Page 10 of 35

To evaluate composite curve temperature for each segment as a function of enthalpy, a pseudo-transient equation is introduced. Specifically, the composite curve temperatures, denoted respectively as T c(zi ) and T h(zi ), are modeled using a set of first-order differential equations X dT c(zi ) Href τ = Hc (z ) − Hcpp (T c(zi ), Pc (z, i), Fc , xc ) T i i T c0 dt c∈Si X dT h(zi ) Href τT = Hhi (zi ) − Hhpp (T h(zi ), Ph (z, i), Fh , xh ) 0 Th dt h∈S

(11) (12)

i

where

Href T c0

and

Href T h0

are scaling factors that ensure consistent units (T c0 and T h0 denote

the stream inlet temperatures). The dynamics of this system have no physical significance and are goverened by the selection of τ , an arbitrary time constant. The selection of time constants for pseudo-transient models is driven by solution robustness and stability considerations and has been examined in previous works 7;17 . Finally, the model is initialized by setting the initial temperature profiles to be constant along the MHEX, with temperature values equal to the inlet conditions T c(zi ) = T c0 , ∀ zi = [0, 1, ..., Ni ]

(13)

T h(zi ) = T h0 , ∀ zi = [0, 1, ..., Ni ]

(14)

While most MHEX models only account for one side of the heat exchanger 5;24 , this pseudo-transient representation provides both the hot and cold composite curves. In turn, this allows for approximating the heat exchanger area by using the required heat duty and the temperature driving force between the composite curves. Assuming a heat transfer coefficient Ui for each enthalpy interval, the area of the heat exchanger is computed as 7

Ai =

Qi X 1 Ui z T h(zi ) − T c(zi ) i X AHX = Ai i


(15) (16)

Page 11 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where Ai is the required area for enthalpy interval i, and AHX is the total area of the heat exchanger. In turn, AHX can be used to approximate the capital investment for the designed MHEX. 2.3

Pseudo-Transient Compressor Model

Compressors account for the majority of energy used in natural gas liquefaction, and therefore predictive compressor models are crucial for optimization of energy consumption. Specifically in the context of optimization under uncertainty, compressor models should be able to describe changes in efficiency and stability (e.g. operation past the surge point) for a range of operating conditions. Compressors are pressure change units that can be described by the following energy balance equation at steady state: ˙ − Qext 0 = Hin − Hout + W

(17)

˙ is the work added to the system and Qext includes losses due to inefficiency. We where W assume that the pressure at the inlet and outlet of the compressor (Pin and Pout ) are design parameters, and the energy balance equation above is used to compute the outlet enthalpy (and temperature) given the inlet enthalpy. The efficiency of a compressor is defined as ise Hout − Hin ηc = Hout − Hin

(18)

ise where Hout is the enthalpy of the outlet stream for an isentropic process:

ise ise Sin = S(Tin , Pin , F, x) = S(Tout , Pout , F, x) = Sout

(19)

ise ise Hout = H pp (Tout , Pout , F, x)

(20)

H pp denotes enthalpy evaluated by a physical properties package (10). A general concept in pseudo-transient process modeling involves introducing a continuation parameter 7;17;21 α ˆ to modulate sink/source terms such that they are not present at initialization, yet are fully



Page 12 of 35

enforced at steady state. dα ˆ = 1−α ˆ dt

(21)

α(t ˆ = 0) = 0

(22)

τα

Treating the nonideality (ηc < 1) of the compressor as a sink in the energy balance (17), we multiply the term containing ηc by the continuation parameter α ˆ and formulate a pseudotransient compressor model, whereby first-order differential equations are used to calculate ise the unknown temperatures. The isentropic outlet temperature Tout is calculated by matching

the inlet and (isentropic) outlet entropies (19), and the actual outlet temperature Tout is calculated using the energy balance (17) and the definition of efficiency (18): ise dTout Sin ise τT = Sin − Sout 0 T dt α ˆ ise = Hin − Hout + (Hout − Hin ) ηc

Hin dTout τT 0 T dt

(23) (24)

where T 0 is the inlet temperature to the compressor. Most previous works regarding steadystate LNG process optimization 5;7;8;12 assumed the compressor effiency to be known and constant. While this assumption may be appropriate for the nominal-case operation, compressor efficiency is in effect dependent on the operating regime, and varying the natural gas feed composition under uncertainty may result in sub-optimal compressor performance. In order to capture variations in efficiency, we use the cubic equation given by Jensen 12 to model the characteristic curve of a compressor. Given a desired pressure ratio P r = Pout /Pin , the required rotational speed N of the compressor can be computed from 3 mr ˙ mr ˙ − 1 − 0.5 −1 P r = P r0 + H 1 + 1.5 W W

(25)

H = H0 − 1.2(H0 + 0.5P r0 − 1)(1 − Nr)

(26)

W = W0 (Nr)1/3 p Nr RT 0 /MW N= D

(27)


(28)

Page 13 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where P r0, H0 and W0 are parameters obtained from fitting compressor characteristic curves, D is the diameter of the compressor wheel, and mr ˙ is the non-dimensionalized flow, defined p −1 . H and W are the semi-height and semi-width of the as mr ˙ = m ˙ RT 0 /MW D −2 Pin

compressor characteristic curve, MW is the molar weight of the compressor stream, and Nr

is a non-dimensional speed ratio given by Saravanamuttoo et al. 25 . It is noteworthy that in many empirical expressions 4;12 , compressor rotational speed is given as a percentage (of maximum speed), rather than as an absolute velocity. The compressor efficiency is then computed 12

ηc = η0

H − H0 2 1− − 1000(mr ˙ − 2W )2 H0

(29)

where η0 , the nominal efficiency, is a parameter of the compressor. 8 7 6 N=100%

5 4 3 2 N=10% 1 0 0

0.002 0.004 0.006 0.008

0.01

0.012 0.014 0.016 0.018

0.02

Figure 3: Compressor curves over a range non-dimensionalized flow for various compressor speeds. Black dots are plotted at the peak the curve for each speed, with operation on a curve to the left of the peak representing compressor surge.

A set of compressor curves predicted using (25)–(28) for various rotational speeds is shown in Figure 3. The surge point, where dynamic instability occurs, is located near the peak value of the pressure ratio for a given compressor rotational speed 12 , plotted as black



Page 14 of 35

dots in Figure 3. Operation to the left of this point on a curve may be possible with active (feedback) control, but generally it is preferable to operate to the right of the surge line 26 . As proposed by Jensen 12 , we define mr ˙ peak as the value of mr ˙ for which the pressure ratio of the compressor characteristic curve reaches its maximum, and we can then express a dimensionless surge margin:

Surge Margin = mr ˙ − mr ˙ peak

3

(30)

Steady-State Design of LNG Plants for Uncertainty in Feed Composition

Dynamic market conditions and diverse portfolios of available feedstocks have led to considerable uncertainty in the design and operation of petrochemical processes 27;15 , such as natural gas liquefaction. The problem of process flowsheet optimization under uncertainty is generally formulated with some uncertain parameter(s) assumed to be drawn from known continuous probability distributions 28 , with probabilistic approaches considering the problem of optimizing the expected value of the objective function. Scenario-based formulations are a widely used approach for solving the infinite-dimensional problem of optimization under uncertainty, relying on sampling the uncertain parameters at a finite number of points (scenarios) to estimate expected values. In this work, a specific case of parametric uncertainty is employed to capture the situation of a varying-composition feed stream to a liquefaction process. For a nominal natural gas feed stream 5;7 (Table 1), we consider an uncertain parameter δ drawn from a probability distribution J(δ), which represents the additional amount of methane present in the feed (in kmol/s). The flow rates of the other components are scaled accordingly, such that ethane, propane, n-butane, and nitrogen remain in the same molar proportions to each other and the total molar mass of the natural gas stream remains constant for all values of δ. We assume that δ is bounded between -0.1 kmol/s and 0.1 kmol/s. The feed gas compositions at these bounds are shown in Table 1, and we note that they are very similar to the lean and rich gas compositions found in other works 10;11 . A multiple-scenario optimization formulation with an equal-width discretization (δ is 14 Environment ACS Paragon Plus

Page 15 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1: Natural gas feed composition (in mol %) under uncertainty

Lower Bound Base Case 5;7 Upper Bound

δ Methane -0.1 84.0 0 89.7 0.1 95.0

Ethane 8.6 5.5 2.7

Propane 2.8 1.8 0.9

n-Butane Nitrogen 0.2 4.4 0.1 2.8 0.0 1.4

Molar Mass 17.6 g/mol 17.6 g/mol 17.6 g/mol

sampled at n equally spaced intervals between δ L and δ U , including the bounds) results in n realizations of δ, where the value of n is defined prior to the optimization procedure. Scenarios (δ1 , δ2 ..., δn ) are equally spaced, such that δi = δ L + (i − 1)∆δ, where

∆δ = δi+1 − δi =

δU − δL n−1

(31)

The multiple-scenario optimization problem associated with the above discretization can be expressed as 15 : n

X 1 J(δ L + (i − 1)∆δ)P (d, zi , xi , δ L + (i − 1)∆δ) L + (j − 1)∆δ) J(δ j=1 i=1

max Pn

d,z1 ,...,zn

(32)

s.t. f(d, zi , xi , δ L + (i − 1)∆δ) = 0, ∀ i = 1, ..., n d ∈ D, z1 , ..., zn ∈ Z, x1 , ...xn ∈ X

In our recent work 15 , we showed that the computational expense of solving a scenariobased formulation such as (32) can be reduced by forgoing the discretization of continuous probability density functions, instead representing uncertain parameters as a set of (continuous) disturbance variables in a fictitious (pseudo) time domain. This scenario-free approach, illustrated in Figure 4, maintains the steady-state flowsheet equations as path constraints throughout the pseudo-time integration. A sequential dynamic optimization routine is then carried out, wherein the objective function, inequality constraints, and their respective sensitivities are integrated in each iteration and used to update the decision variables. The process design decision variables d are held constant for each iteration, while the control/recourse decision variables z are treated using a control vector parameterization 29 . We demonstrated that order-of-magnitude reductions in computational time and memory usage can be achieved using the proposed scenario-free approach, mainly owing to the fact 15 Environment ACS Paragon Plus


Page 16 of 35

P(d,z(t),θ1(t)) Objective Function θ1U

θ1(t)

Uncertain Parameter

θ1L z(t2)

z(t3)

z(t4)

z(t1)

...

z(tm) Recourse Variable

d Design Variable

Pseudo-time

Figure 4: A single iteration of optimization under uncertainty (single uncertain parameter) using sequential dynamic optimization 15 . The design variables d are held constant, and the control variables z(t) are handled using a piecewise-constant control vector parametrization with m pieces.

that only a single instance of the process flowsheet must be initialized and stored during an optimization iteration. For a single uncertain parameter δ, as in (32), the trajectory is defined to be linear in pseudo-time, efficiently and completely searching the parameter uncertainty space between the lower bound δ L and the upper bound δ U . The corresponding optimization problem is written as follows 15 : Z tf 1 δU − δL J(δ(t))P (d, z(t), x(t), δ(t))dt max R δU d,z(t) tf 0 J(δ)dδ L δ s.t.f(d, z(t), x(t), δ(t)) = 0, ∀t ∈ [0, tf ]

(33)

d ∈ D, z(t) ∈ Z, x(t) ∈ X, ∀t ∈ [0, tf ] δ(t) =

3.1

δU − δL t + δL tf

Practical Considerations for the Dynamic Optimization Approach

The resulting scenario-free optimization problem is equivalent to the scenario-based formulation (32) as n → ∞ if an infinite-dimensional control vector parametrization is used for the recourse variables z 15 . In practice, an infinite-dimensional control vector cannot be implemented, and we instead solve the optimization with increasingly fine finite-dimensional control vector parameterizations until the difference in objective function between two suc-


Page 17 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


cessive refinements drops below a predefined threshold 15 . While this dynamic optimization approach for dealing with uncertainty can be generally applied, we note that the DAE solvers used for sequential dynamic optimization are implicit and rely on algebraic solvers. Thus, the application of the proposed formulation to chemical processes may lead to DAE solver difficulties in handling the constraint system of the proposed formulation (i.e. the algebraic, steady-state process models) for some values of the uncertain parameter and/or control variables. We previously showed 15 that this issue can be largely avoided by using a pseudo-transient representations of the process models 17;7 to assist the time integration. Since the solution of the original algebraic process flowsheet model (1) is recovered at the steady state of a pseudo-transient model (per the static equivalance condition), we incorporate the pseudo-transient models into the proposed optimization approach under uncertainty by selecting time constants such that the dynamics of the pseudo-transient process model evolve in a much faster time scale than the dynamics of the uncertain parameters. With this judicious choice of the time constants τ , the imposed dynamics of the resulting system are thus in a standard singularly perturbed form, and the (fast) pseudo-transient process model can be assumed to be at a (pseudo)steady-state at all times 30 . Further, we showed that the error (deviation from steady-state solution) can be tuned by increasing the difference between the time constants of the pseudo-transient model and dynamics of the uncertain parameter (33) 15 . 4

R PRICO Process Description

R process 31 is a representative single-stage, mixed refrigerant natural gas liqThe PRICO

uefaction process that has been studied extensively in the recent past 5;7;8;12 . In this work, we ignore the removal of heavy components from and pretreatment of the natural gas feed stream, including removal of water, carbon dioxide, etc. Further, we assume that the natural gas stream passes through the full length of the MHEX, differing slightly from our previous work 7 . This is motivated by the fact that the optimal solution found there involved both hot streams exiting the MHEX at nearly identical temperatures. The process flow diagram



R for the PRICO process is shown in Figure 5.

A 1 kmol/s natural gas stream (S1), nominally composed of 89.7% methane, 5.5% ethane, 1.8% propane, 0.1% butane, and 2.8% nitrogen, enters the MHEX at 25◦ C and 55 bar and exits the MHEX liquefied and subcooled at -155◦C. The mixed-refrigerant stream is composed of the same five components and is employed in a single-stage refrigeration cycle. Hot mixed refrigerant is cooled in the MHEX, expanded through a throttle valve, repassed through the MHEX to liquefy the natural gas, and then compressed and chilled to 25◦ C in a salt water (SW) cooler. The pseudo-transient models described in the previous section are used to model all process units, and the Soave Redlich Kwong cubic equation of state is used to evaluate thermodynamic properties, with fluid properties calculated using mixing rules. Several numerical challenges arise during optimization of this process flowsheet 7 : • the phase of the mixed-refrigerant stream is unknown throughout the flowsheet • the pressures and temperatures of the refrigerant stream can vary throughout the process • the composition and flow rate of the refrigerant stream are allowed to vary during optimization, resulting in changing phase boundaries in the MHEX • a small minimum approach temperature in the MHEX (1.2◦ C) is imposed to ensure maximum energy recovery In addition, natural gas feedstocks can be of various compositions, and variations in natural gas composition affect the specific work, as the specific enthalpy of each component is different from the others. As a consequence, the gas compositions considered in past studies that focused on simulation and optimization of LNG plants vary considerably 11;10 4.1

Design Optimization Problem Formulation

R liquefaction process (FigThe design decision variables considered here for the PRICO

ure 1) include compressor diameter (as described in Section 2.3), refrigerant flowrate, refrigerant composition, and the high and low pressures of the refrigeration cycle. The single degree of freedom of the MHEX is the exit temperature of the cold refrigerant stream (S6 18 Environment ACS Paragon Plus

Page 18 of 35

Page 19 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


in Figure 5), which must be completely vaporized to avoid damaging the compressor. With these considerations, the following constraints are included in the optimization problem: • The temperature of the cold refrigerant exiting the MHEX (S6) must be at least 10◦ C greater than the dew point • The minimum approach temperature in the MHEX must be at least 1.2◦ C • The compressor does not operate past 100% of its maximum rotational speed • The compressor surge margin, as given by (30), is positive The compressor surge margin is merely constrained to be positive because operation beyond this point is still possible with feedback control 12;26 ; however, the bound can easily be increased to provide a more conservative design. The pseudo-transient models from Section 2 are used to simulate and optimize the flowR sheet. For optimization of the PRICO process, we considered two objective functions: (i)

minimization of the compressor work 7 and (ii) minimization of a nonlinear objective function incorporating the size of the MHEX 12;18 . The objective functions are respectively denoted as OBJ1 or OBJ2. As will be described in the following sections, optimization solutions were computed by minimizing either OBJ1 or OBJ2

OBJ1 = Wcomp

(34)

OBJ2 = Wcomp + C0 (AM HEX )0.65

(35)

The size of the MHEX for OBJ2 is calculated assuming an average heat transfer coefficient of U = 500 W m−2 K −1 , and the cost factor (C0 ) was selected to be 2135 W m−1.3 , as given by Jensen and Skogestad 18 . 4.2

Compressor and MHEX Considerations

Flowsheet optimization was carried out for design configurations with both single- and two-staged compression (the latter with intercooling, similar to the configuration proposed by Tak et al. 11 ). The process flowsheets for both configurations are shown in Figure 5. 19 Environment ACS Paragon Plus


Page 20 of 35

Throttle Valve Throttle Valve

S6

S7

S6

S7 (CMR)

(CMR) LNG -155 °C MHEX

LNG -155 °C MHEX

S2

S2 SW Cooler 2

SW Cooler

S5

S5

25 °C (HMR)

S4

25 °C (HMR)

S1

NG 55 bar, 25 °C

Comp 2

NG 55 bar, 25 °C

S4

S1

S3

S3

Comp 1

Comp

SW Cooler 1

R Figure 5: Process flow diagram for the simplified PRICO natural gas liquefaction process. The configuration with a single compressor is shown on the left, and the configuration with two compression stages and intercooling is shown on the right.

The nominal efficiency of all compressors was assumed to be 82.2%, and compressor curve parameters were regressed from empirical data 4 (Table 2). We note that although the fitted pressure ratio at zero flow rate (P r0 ) is negative, the non-surge regime is to the right side of the peak, and the parameter does not have physical meaning as modeled 12;4 . Compressor curves with these parameters are plotted for increasing speeds in Figure 3. The first value of W0 is used for the first compressor, and the second value is used for the second compressor in the two-stage, intercooled compressor configuration. Table 2: Compressor Model Parameters

Parameter P r0 H0 W0

Physical Significance Pressure ratio at zero flow Reference semiheight of compressor curve Reference semiwidth of compressor curve

Value -22.0 14.0 0.007785, 0.00213991

Using the modeling framework presented in Section 2.2, the MHEX is discretized into 50 segments, for which enthalpies are established and the pseudo-transient energy balances, given by (11)–(12), are used to compute temperature. The pressure drops across the process MHEX are assumed to be 5 bar for the natural gas stream, 4 bar for the hot refrigerant stream, and 1 bar for the cold refrigerant stream. Additionally, a 0.1 bar pressure drop is 20 Environment ACS Paragon Plus

Page 21 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


assumed for each SW cooler. To facilitate satisfying the constraint of maintaining constant heat exchanger area for various concentrations of natural gas feed, we fix the heat exchanger area A to a constant Asetpoint using an integral-only controller, as proposed in our previous work 21 . The controller eliminates a degree of freedom from the flowsheet model, which we select as the outlet pressure of the compressor, and replaces it with a new decision variable, Asetpoint .

τI

Asetpoint dPout = A − Asepoint 0 Pout dt

(36)

0 Pout (t = 0) = Pout

(37)

Here τI is the integral time constant of the controller, which should be selected such that the dynamics of the controller are at least as slow as the dominant dynamic mode of the pseudo-transient MHEX model 32 . The sign of the gain between the controlled variable (A) and the selected manipulated variable must also be known (at least locally) in order to maintain closed-loop stability. We note that a PI control law could potentially be implemented; while this may provide some convergence benefits, stability of the resulting closed-loop dynamical system is more difficult to ensure. 5

Base Case Optimization Results

The process flowsheet was optimized seperately for both objective functions (34)–(35) and both compression sequences using a time-relaxation-based algorithm 17 . The problems were all solved with the NLPSQP (sequential quadratic programming) solver in gPROMS 33 5.0.1 using a 64-bit Windows 7 desktop system with a 3.40GHz Intel Core i7 processor and 16GB of RAM. Several randomly generated initial guesses were provided for each optimization to find the best locally optimal solution for each case. The problem formulations are described in Table 3 and results of all optimization calculations are presented in Table 4. Although the same mathematical models were used in every optimization, it was found that minimization of OBJ2 instead of OBJ1 always results in a smaller MHEX, reduced mixed refrigerant flow, higher total compression ratio, and increased compressor work (Ta-



Page 22 of 35

Table 3: Optimization Problem Formulations

P1a # of Compressors Compressor Model Objective Function CPU Time (s)

P1b

P2a

P2b

P3a

P3b

P4a

P4b

1

2 Fixed ηc Equal Pr OBJ1 OBJ2 OBJ1 OBJ2 OBJ1 OBJ2 OBJ1 OBJ2 719 204 413 89 460 347 393 232 -

ble 4). Optimization of OBJ2 also required less CPU time in every case, likely because the steeper objective function allowed local optima to be reached in fewer iterations. A comparison of optimization problems P1 and P2 reveals that the assumption of constant compressor isentropic efficiency is adequate for nominal-case flowsheet optimization (a constant ηc = 80% was assumed, as was used in previous works 7;8 ), and optimization with the new compressor model results in nearly the same operating pressures, MHEX size, and mixed-refrigerant (MR) stream. Table 4: Nominal-Case Optimization Results. C1 and C2 refer to compressors 1 and 2.

Pressure S3 (bar) Pressure S4 (bar) MR Flow (kmol/s) N2 (mol%) CH4 (mol%) C2 H6 (mol%) C3 H8 (mol%) n-C4 H6 (mol%) C1 D (m) C2 D (m) C1 Pr C2 Pr C1 Work (MW) C2 Work (MW) Total Pr AM HEX (m2 ) OBJ1 (MW) OBJ2 (MW)

P1a P1b 9.97 5.69 35.72 33.85 3.889 2.838 16.93 13.11 39.41 38.64 27.14 26.76 0.00 0.00 16.52 21.48 1.425 1.654 3.58 5.95 15.54 16.06 3.58 5.95 79871 48137 15.54 16.05 18.87 18.44

P2a 9.99 35.76 3.889 16.95 39.40 27.14 0.00 16.51 3.58 15.89 3.58 79894 15.89 19.22

P2b 5.90 33.89 2.889 13.37 38.63 26.83 0.00 21.15 5.74 16.44 5.87 49567 16.44 18.88

P3a 10.95 39.95 3.823 18.72 38.57 27.09 0.00 15.59 1.318 1.667 2.06 1.78 8.35 5.92 3.67 78123 14.27 17.55

P3b 3.47 42.63 2.010 10.44 39.76 25.40 0.00 24.40 1.803 1.481 4.63 2.67 9.76 5.06 12.36 34058 14.82 16.73

P4a 11.41 42.36 3.777 19.62 38.11 27.08 0.00 15.18 1.286 1.685 1.93 1.93 7.44 6.85 3.72 76479 14.29 17.52

P4b 5.18 44.47 2.376 13.97 38.99 25.30 0.00 21.75 1.543 1.652 2.93 2.93 7.78 7.09 8.58 38163 14.87 16.93

Aside from the solutions for problems P2a and P2b, every compressor operates near its nominal efficiency of 82.2% (the lowest value is 81.5%). As expected, when the compressor size is simultaneously optimized with the nominal flowsheet design, the compressor is sized 22 Environment ACS Paragon Plus

Page 23 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


such that it operates near its peak efficiency. The configuration with two compressors in sequence (with intercooling) reduces the total compression work required by the system, and it also allows for the refrigerant stream flowrate and heat exchanger size to be further reduced when heat exchanger size is penalized (P3b and P4b). Although the results found through elimination of the intermediate pressure in the two-stage compression sequence (P4) and the results with the intermediate pressure as an additional decision variable (P3) have similar objective function values, the designs are noticeably different. This is most apparent from the solution to problem P3b, where the first compressor accounts for much more of the compression work than the second, and the overall pressure ratio is also much higher. The total refrigerant flow rate differs greatly between designs, but the optimal refrigerant stream compositions for all design cases are relatively similar, comprising mostly methane and ethane, and with no propane. The MHEX composite curves for the optimal designs P1a and P1b are shown in Figure 6. The composite curves for P2a, P3a, and P4a closesly resemble the ones for P1a, and those found in the solutions for P2b, P3b, and P4b closesly resemble the ones P1b. The designs are tightly heat-integrated, with the MHEXs found in the minimization of OBJ1 reaching the minimum allowable temperature driving force of 1.2 K at multiple points. We also note here the importance of using a distributedparameter mathematical model in decribing the MHEX: all MHEX designs involve reaching the minimum allowable temperature driving force at one or more interior points. 6

Optimization Under Uncertainty Results

Upon completing the nominal-case optimization calculations presented above, the process flowsheet was optimized for an uncertain feedstock, with δ drawn from a uniform distribution δ ∼ U(−0.1, 0.1) as in Table 1. Along with the MHEX area and the compressor size(s), the R refrigerant flow and composition were selected as design decision variables for the PRICO

process, such that the pressurized refrigeration cycle would not have to be altered during process operation. The rest of the decision variables were treated as recourse variables, as described in Section 3. The partitioning of decision variables into design and recourse/control variables is shown in Table 5.



P1a

300

P1b

300

280

280

260

260

240

240

Temperature (K)

Temperature (K)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 35

220 200 180

220 200 180

160

160

140

140

120

120

100

Hot Composite Cold Composite

100 0

10

20

30

40

50

60

70

80

0

Enthalpy (MW)

10

20

30

40

50

60

Enthalpy (MW)

Figure 6: Optimal temperature-enthalpy diagram for the MHEX designs found in P1a and P2a. The vertical gray lines denote enthalpy segments where the temperature driving force exactly meets the constraint of 1.2 K. Table 5: Optimization Under Uncertainty Decision Variables

Design MR Flow MR Composition MHEX Area C1 Diameter C2 Diameter

Recourse C1 Pr C2 Pr

The problems were again solved with the NLPSQP (sequential quadratic programming) solver in gPROMS 33 5.0.1 using a 64-bit Windows 7 desktop system with a 3.40GHz Intel Core i7 processor and 16GB of RAM. Pseudo-transient models were used to initialize the process model for each optimization iteration and to smooth the algebraic process constraints 15 , and the maximum residual in any steady-state algebraic equation or deviation from contoller setpoint during integration through pseudo-time was less than 0.004%. The values of the time constants were selected to be τT = 10 and τα = 100, and the integration horizon (tf in (33)) was selected to be 1000000. The optimal result for the nominal case (Table 4) for each problem formulation was used as the initial guess for optimization. The problem was not solved with fixed compressor efficiency (P2a and P2b), nor was the problem solved with equal compression ratios (P4a and P4b), as smaller objective function values were found by leaving the intermediate pressure as a decision variable. The results of the optimization calculations for problems P1a∗ , P1b∗ , 24 Environment ACS Paragon Plus

Page 25 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


P3a∗ , and P3b∗ , where the asterisk denotes optimization under uncertainty, are presented in Table 6, along with deviations of the respective values from the base case results. The solution of the optimization problems under uncertainty required slightly more memory, with the peak memory usage in the optimizations reaching 121.23 MB, compared to 116.78 MB in the nominal case optimizations. The solutions for P1a∗ , P1b∗ , P3a∗ , and P3b∗ were found in 36743, 26930, 36030, and 58638 seconds respectively. As found in previous work 15 , the bulk of the time per iteration was spent initializing the model at the initial condition of the uncertain parameter, or solving the equality constraints of the algebraic process model (1) using a pseudo-transient time integration to steady-state. Comparably, less time was spent on the integration of the model (with sensitivities) through the trajectory of the uncertain parameter. While the nominalcase optimization problems were solved using a time relaxation-based algorithm, wherein each optimization iteration was initialized from the steady-state solution of the previous iteration, dynamic optimization in gPROMS does not allow for this time-saving initialization strategy. We estimate the solution times if this capability were available by calculating the average percentage of time within each optimization iteration spent on initialization and replacing it with the amount of time required to initialize each iteration of the same nominal-case problem (since the models are approximately the same size). We found that initialization took approximately 81% and 91% of the total time of an iteration for the single compressor and two-stage flowsheets, respectively. Using this finding, we estimate that the solution times for P1a∗ , P1b∗ , P3a∗ , and P3b∗ could be respectively reduced to 8429, 5610, 3889, and 5635 seconds. The number of degrees of freedom available for MHEX operation is very limited, as is characteristic of many integrated processes 34 . In a related area, the literature surrounding the design of air separation processes demonstrates that MHEXs represent a significant design challenge when variable operation is considered, from both steady-state design 16 and dynamic operation 35 perspectives. In the results presented in Table 6, the optimal MHEX areas found for P1a∗ and P3a∗ are significantly reduced from those found in the nominal case optimization calculations, and the optimal refrigerant stream flowrate found for P1bb is 25 Environment ACS Paragon Plus


Page 26 of 35

Table 6: Optimization Under Uncertainty Results

MR Flow (kmol/s) N2 (mol%) CH4 (mol%) C2 H6 (mol%) C3 H8 (mol%) n-C4 H6 (mol%) C1 D (m) C2 D (m) E[C1 Work] (MW) E[C2 Work] (MW) AM HEX (m2 ) E[OBJ1] (MW) E[OBJ2] (MW)

P1a∗ 3.963 17.11 39.41 27.06 0.00 16.42 1.434 15.72 75186 15.72 18.92

∆ P1b∗ 1.9% 2.976 1.1% 13.83 0.0% 38.77 -0.3% 26.95 0.0% 0.00 -0.6% 20.45 0.6% 1.618 1.2% 16.05 -5.9% 48970 1.2% 16.05 0.3% 18.44

∆ P3a∗ 4.6% 3.534 5.5% 18.50 0.3% 38.80 0.7% 26.70 0.0% 0.00 -4.8% 15.99 -2.2% 1.333 1.650 0.0% 8.38 6.02 1.7% 66173 0.0% 14.40 0.0% 17.30

∆ P3b∗ -7.6% 2.005 -1.2% 10.40 0.6% 39.78 -1.4% 25.44 0.0% 0.00 2.6% 24.37 1.1% 1.841 -1.0% 1.485 0.4% 9.78 1.7% 5.04 -15.3% 34035 0.9% 14.83 -1.5% 16.71

∆ -0.3% -0.4% 0.0% 0.2% 0.0% -0.1% 2.1% 0.3% 0.2% -0.4% 0.0% 0.0% -0.1%

increased. As suggested by Jensen and Skogestad 18 , the exclusion of MHEX area from the objective function results in a very tightly integrated process with a large heat exchanger and refrigerant stream. While intuition might suggest that equipment should be oversized for more flexible operation, we find that in reality the MHEX should be undersized to reduce the degree of heat integration. The temperature driving force profiles along the MHEXs in the optimal designs are plotted in Figure 7, where it can be seen that the designs found in P1a∗ and P3a∗ deviate slightly more from the minimum driving temperature force than the nominal case optimization results. To compare the level of heat integration between designs, we also compute the mean temperature difference ∆T of each MHEX 18 1 ∆T = A

Z

∆T dA

(38)

The expected value of ∆T found in P1a∗ is 2.133 K, throughout the range of operation for all possible feed concentrations, while in the nominal case, the mean temperature difference in the MHEX is 1.974 K. This confirms that the larger MHEX design is more tightly heat integrated and less flexible, as the ∆T is smaller than the minimum found in the optimization under uncertainty case (2.131 K). The same increase in ∆T is seen in the solution to P3a∗ , where the expected and minimum values are respectively 2.141 K and 2.136 K, compared to 26 Environment ACS Paragon Plus

Page 27 of 35

P1a *

10

15

6

T (K)

T (K)

P1b *

20

8

10

4 5

2 0

0 0

0.2

0.4

0.6

0.8

1

0

0.2

Q/Q tot

0.4

0.6

0.8

1

0.8

1

Q/Q tot

P3a *

10

P3b *

30

Range of Operation

25

8

Expected Value Nominal Case Result

20 6

T (K)

T (K)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


15

4 10 2

5

0

0 0

0.2

0.4

0.6

0.8

1

0

Q/Q tot

0.2

0.4

0.6

Q/Q tot

Figure 7: Temperature driving force along the MHEXs found by optimization under uncertainty. The shaded areas show the range of values found in the entire range of uncertain feed concentrations, and the red dashed lines mark the 1.2 K lower bound for ∆Tmin .

the ∆T in the optimal nominal-case design of 1.950 K. Similarly, in the solution to P1b∗ , the expected and minimum values of ∆T are 2.594 K and 2.590 K, compared to the nominal result of 2.541 K. We observe that the solution found for P3b∗ does not deviate much from the nominal case optimization (P3b), and the MHEX has an identical expected ∆T as the nominal case optimization, indicating that the design found in the solution to P3b is already flexible enough to accomodate the range of methane feed compositions considered. The required work profiles, divided between compressors, found in P3a∗ and P3b∗ are shown in Figure 8. As control of tightly integrated processes require quick computation of the optimal pro27 Environment ACS Paragon Plus


13 *

P3a C1

12

*

P3a C2 *

11

P3b C1 *

P3b C2

10

Work (MW)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 35

9 8 7 6 5 4 84

85

86

87

88

89

90

91

92

93

94

95

Feed Methane Concentration (mol%)

Figure 8: Compressor work profiles found in P3a∗ and P3b∗ . The dotted lines represent the values found in the nominal-case optimization.

cess inputs 34 , identification of the operating regimes for different feedstock concentrations along with the optimal design specifications is crucial. Therefore, the recourse variable profiles are a second important result from the solution of the optimization problem under uncertainty using the proposed sequential dynamic optimization approach. The trajectories of the recourse variables were parameterized 29 using piecewise-linear profiles, giving the optimal operating conditions for the explored range of feed conditions 15 . The profiles of the high and low pressures, along with their values found in the nominal case optimization, are shown in Figure 9, where it can be seen that the dependence of the optimal high and low pressures on feed composition for various designs is highly nonlinear. 7

Conclusions

In this paper, we present a novel approach for optimization of natural gas liquefaction processes under feedstock uncertainty. Our work provides a general, scenario-free framework for considering varying steady-state operating regimes at the design stage. We present mathematical models based on novel pseudo-transient continuation concepts that allow for key equipment (MHEXs and compressors) to be sized simultaneously with the flowsheet optimization. We then introduce variations in feedstock concentration as a specific parametric R uncertainty and optimized the representative PRICO mixed-refrigerant NG liquefaction

process gas over a range of feed methane concentrations. 28 Environment ACS Paragon Plus

Page 29 of 35

48 P1a

46

*

*

P3b

44

Pressure S4 (bar)

*

P1b P3a

*

42 40 38 36 34 32 84

85

86

87

88

89

90

91

92

93

94

95

Feed Methane Concentration (mol%) 12 11 10 9

Pressure S3 (bar)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


P1a * P1b *

8

P3a *

7

P3b *

6 5 4 3 2 84

85

86

87

88

89

90

91

92

93

94

95


Figure 9: High and low pressure profiles found in the optimization under uncertainty results. The dotted lines represent the values found in the nominal-case optimization.

In the nominal case, we find that the use of two-stage compression allows for a significant reduction of the total required compressor power and that the inclusion of MHEX area in the objective function results in smaller, less tightly heat-integrated MHEX designs. The same flowsheet is then optimized while accounting for uncertainty in feed methane concentration, showing notable differences in the MHEX sizing when MHEX area is not included in the objective function. We find that, counterintuitively, the optimal MHEXs designed for the nominal-case operation (when size is not penalized) are too large – and thus too tightly heat integrated – to operate flexibly, and that they should actually be designed to be “undersized” for an uncertain feed stream.



Acknowledgements The authors gratefully acknowledge funding from the National Science Foundation (NSF) through the CAREER Award 1454433 and Award CBET-1512379.


Page 30 of 35

Page 31 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Literature Cited 1. Zaretskaya, V. and Dyl, K. United States expected to become a net exporter of natural gas this year. Technical report, U.S. Department of Energy: Energy Information Administration, August 2017. 2. Jensen, J. B. and Skogestad, S. Optimal operation of a mixed fluid cascade LNG plant. Computer Aided Chemical Engineering, 21:1569–1574, 2006. 3. Initiative, M. E. The future of natural gas: An interdisciplinary MIT study. Cambridge, MA: Massachusetts Institute of Technology, 2011. 4. Michelsen, F. A., Halvorsen, I. J., Lund, B. F., and Wahl, P. E. Modeling and simulation for control of the TEALARC liquified natural gas process. Industrial & Engineering Chemistry Research, 49(16):7389–7397, 2010. 5. Kamath, R. S., Biegler, L. T., and Grossmann, I. E. Modeling multistream heat exchangers with and without phase changes for simultaneous optimization and heat integration. AIChE Journal, 58(1):190–204, 2012. 6. Khan, M. S. and Lee, M. Design optimization of single mixed refrigerant natural gas liquefaction process using the particle swarm paradigm with nonlinear constraints. Energy, 49:146–155, 2013. 7. Pattison, R. C. and Baldea, M. Multistream heat exchangers: Equation-oriented modeling and flowsheet optimization. AIChE Journal, 61(6):1856–1866, 2015. 8. Tsay, C., Pattison, R. C., and Baldea, M. Equation-oriented simulation and optimization of process flowsheets incorporating detailed spiral-wound multistream heat exchanger models. AIChE Journal, 63(9):3778–3789, 2017. 9. Lee, G. C., Smith, R., and Zhu, X. X. Optimal synthesis of mixed-refrigerant systems for low-temperature processes. Industrial & Engineering Chemistry Research, 41(20):5016–5028, 2002.



10. Aspelund, A., Gundersen, T., Myklebust, J., Nowak, M. P., and Tomasgard, A.

Page 32 of 35

An

optimization-simulation model for a simple LNG process. Computers & Chemical Engineering, 34(10):1606–1617, 2010. 11. Tak, K., Lee, I., Kwon, H., Kim, J., Ko, D., and Moon, I. Comparison of multistage compression configurations for single mixed refrigerant processes. Industrial & Engineering Chemistry Research, 54(41):9992–10000, 2015. 12. Jensen, J. B. Optimal operation of refrigeration cycles. PhD thesis, Norwegian University of Science and Technology, 2008. 13. Rao, H. N. and Karimi, I. A. A superstructure-based model for multistream heat exchanger design within flow sheet optimization. AIChE Journal, 63(9):3764–3777, 2017. 14. Li, X., Tomasgard, A., and Barton, P. I. Natural gas production network infrastructure development under uncertainty. Optimization and Engineering, 18(1):35–62, 2017. 15. Tsay, C., Pattison, R. C., and Baldea, M. A dynamic optimization approach to probabilistic process design under uncertainty. Industrial & Engineering Chemistry Research, 56(30):8606– 8621, 2017. 16. Pattison, R. C. and Baldea, M. Optimal design of air separation plants with variable electricity pricing. Foundations of Computer Aided Process Design (FOCAPD), pages 393–398, 2014. Cle Elum, WA. 17. Pattison, R. C. and Baldea, M. Equation-oriented flowsheet simulation and optimization using pseudo-transient models. AIChE Journal, 60(12):4104–4123, 2014. 18. Jensen, J. B. and Skogestad, S. Problems with specifying ∆Tmin in the design of processes with heat exchangers. Industrial & Engineering Chemistry Research, 47(9):3071–3075, 2008. 19. Wechsung, A., Aspelund, A., Gundersen, T., and Barton, P. I. Synthesis of heat exchanger networks at subambient conditions with compression and expansion of process streams. AIChE Journal, 57(8):2090–2108, 2011. 20. Watson, H. A. J. and Barton, P. I. Modeling phase changes in multistream heat exchangers. International Journal of Heat and Mass Transfer, 105:207–219, 2017.


Page 33 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


21. Pattison, R. C., Tsay, C., and Baldea, M. Pseudo-transient models for multiscale, multiresolution simulation and optimization of intensified reaction/separation/recycle processes: Framework and a dimethyl ether production case study. Computers & Chemical Engineering, 105: 161–472, 2017. 22. Tsay, C., Pattison, R. C., and Baldea, M. A pseudo-transient optimization framework for periodic processes: Pressure swing adsorption and simulated moving bed chromatography. AIChE Journal, 2017. doi:10.1002/aic.15987. 23. Ma, Y., Luo, Y., and Yuan, X. Simultaneous optimization of complex distillation systems with a new PTC model. Industrial & Engineering Chemistry Research, 56(21):6266–6274, 2017. 24. Dowling, A. W. and Biegler, L. T. A framework for efficient large scale equation-oriented flowsheet optimization. Computers & Chemical Engineering, 72:3–20, 2015. 25. Saravanamuttoo, H. I. H., Rogers, G. F. C., Cohen, H., and Straznicky, P. V. Gas turbine theory. Pearson Education: New York, NY, 2001. 26. Gravdahl, J. T. and Egeland, O. Centrifugal compressor surge and speed control. IEEE Transactions on Control Systems Technology, 7(5):567–579, 1999. 27. Baldea, M. and Daoutidis, P. Dynamics and Nonlinear Control of Integrated Process Systems. Cambridge University Press: Cambridge, United Kingdom, 2012. 28. Wang, S. and Baldea, M. Identification-based optimization of dynamical systems under uncertainty. Computers & Chemical Engineering, 64:138–152, 2014. 29. Vassiliadis, V. S., Sargent, R. W. H., and Pantelides, C. C. Solution of a class of multistage dynamic optimization problems. 1. Problems without path constraints. Industrial & Engineering Chemistry Research, 33(9):2111–2122, 1994. ˇ 30. Ladde, G. S. and Siljak, D. D. Multiparameter singular perturbations of linear systems with multiple time scales. Automatica, 19(4):385–394, 1983. 31. Price, B. C. and Mortko, R. A. PRICO: A simple, flexible proven approach to natural gas liquefaction. In GASTECH International LNG/LPG/Natural Gas Conference, pages 1–14, 1996. Vienna, Austria.



Page 34 of 35

32. Seborg, D. E., Mellichamp, D. A., Edgar, T. F., and Doyle III, F. J. Process dynamics and control, 3rd ed. John Wiley & Sons: Hoboken, NJ, 2011. 33. Process

Systems

Enterprise.

general

PROcess

Modeling

System

(gPROMS).

www.psenterprise.com/gproms, 1997-2017. 34. Baldea, M. From process integration to process intensification. Computers & Chemical Engineering, 81:104–114, 2015. 35. Cao, Y., Swartz, C. L. E., and Baldea, M. Design for dynamic performance: Application to an air separation unit. In American Control Conference (ACC), 2011, pages 2683–2688. IEEE, 2011. San Francisco, CA.


Page 35 of 35

13 P3a * C1

12 11

P3a * C2 P3b * C1 P3b * C2

10

Work (MW)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


9 8 7 6 5 4 84

85

86

87

88

89

90

91

92


For Table of Contents Only.


93

94

95

Scenario-Free Optimal Design under Uncertainty of the PRICO Natural

Recommend Documents