Uncertainty and Global Sensitivity Analysis of a Membrane Biogas Upgrading Process Using COCO Simulator

1. Introduction

1.1. Uncertainty and Global Sensitivity Analysis

In chemical engineering, modelling and simulation allow for studying process performance. Process simulators are helpful for process design because they can rapidly assess different process configurations and operating parameters.

It should be noted that process simulators provide a deterministic solution. However, in real processes, some parameters of the process units are not known with precision. Similarly, the feed characteristics and the influencing environment are subjected to variations. The uncertainty of these factors introduces, in turn, uncertainty in the process response, leading to a mismatch with the simulator results. Because numerical predictions are often the basis of engineering decisions, uncertainty quantification is a subject of concern [1]. Uncertainty analysis allows for quantifying the impact of the variability of factors on process outcomes, enabling the determination of potential limits for the process response. This allows for the robustness and confidence of the simulation results to be assessed. Many process simulators allow parametric sensitivity studies where one or several parameters can be modified to study their effect on critical process outcomes. However, a further step can be taken by applying global sensitivity analysis to know how the uncertainty in the output of a model can be apportioned to different sources of uncertainty in the model input [2]. Thus, this technique can answer which of the uncertain input factors are more critical in determining the uncertainty in the outputs of interest.

1.2. CAPE-OPEN Process Simulators

Process simulators combine physical properties databases, thermodynamic models, unit operation models, and numerical solvers with the capability to develop flowsheets in which model units are connected to form a process. As the best modelling options may come from different sources, many researchers point out the advantages of using CAPE-OPEN standards for their integration. These standards facilitate the implementation of new unit operations in different coding software. This can be crucial for particular unit operations like membrane stages. The simulator chosen for this study, the COCO simulator, offers excellent flexibility for this purpose. This simulator presents the advantages of being a free-of-charge CAPE-OPEN compliant steady-state simulation environment that demonstrates strong reliability in performing material and energy balances, with results closely aligned with those obtained from commercial process simulators [3]. The software allows interaction with user unit operations through add-ons developed by the COCO developer [4]. These units can be implemented in Excel, MATLAB, or free numerical software like Scilab and Python [5].

1.3. Scope of This Work

This work aims to demonstrate the applicability of global uncertainty and sensitivity analysis using the COCO simulator. For this purpose, a biogas upgrading process was chosen. The increasing interest in biomethane as a substitute for natural gas has led to the development of biogas upgrading techniques. Among these technologies, gas permeation based on membranes stands out. In recent years, the development efforts for these processes have focused on the design and analysis of stage configurations with the aim of energetic and economic optimization [6,7,8,9]. We have not found studies addressing the global sensitivity and uncertainty of this specific case. However, some studies deal with local sensitivity through parameter variation [10,11,12]. In these cases, the procedures used were aimed at optimization. In contrast, in our case, global sensitivity was applied to study the uncertainty and global sensitivity of the responses of an optimized system.

This work is structured in the following steps: i) A case definition based on an optimized CO₂/CH₄ separation process existing in the literature. ii) The development of the deterministic model of the process in the COCO simulator, adding a user unit to model gas membrane stages. iii) Identifying and characterizing the uncertainty factors through a bibliographic study. iv) Preparing the process model to accept uncertainty inputs and generating model results for specified input combinations. v) The subsequent uncertainty and global sensitivity analysis and their posterior interpretation.

2. Materials and Methods

2.1. Case Study

To demonstrate the application of uncertainty calculations using the COCO simulator to a biogas-upgrading process, we defined a case study based on the process configuration developed by Evonik [13]. This process was intended to obtain high gas purities using a three-membrane stage configuration with recirculation using a single compressor. Figure 1 shows the process flowsheet implemented in the simulator. In the process, two membrane stages are placed in series; the permeate of the first stage is sent to the third stage, and the permeate stream of the third stage, after decompression, is mixed with the permeate of the second stage to form a stream that is recirculated to be mixed with the process feed. Structural optimization for this process was performed in a thesis [14] for a biogas feed composed of 40% CO₂ and 60% CH₄ upgraded using a typical polyimide membrane with a CO₂ permeance of 60 GPU, obtaining optimal stage areas and pressure for the permeate stream from stage 1 feeding stage 3 (Table 1).

For this study, we considered two operating scenarios for the same process configuration: (i) with active control of the compressor pressure to achieve a constant CH₄ mole fraction of 0.96 in the biomethane product, and (ii) without pressure control. The mole fraction set point of the first scenario was chosen to be very close to the output obtained for the second scenario at the reference values. These scenarios were compared regarding their sensitivity to input uncertainties and process variability. The setpoint of composition for the scenario with regulation was chosen to be close to that obtained for the reference parameters in the unregulated system.

2.2. Deterministic Simulation of the Process in COCO

The deterministic simulation of the case study was performed before the uncertainty analysis, considering the parameters of Table 1 as deterministic. The process model was developed using the COCO simulator (version 3.9.0.1). The thermodynamic model for the whole process was based on the Peng-Robinson equation of state, which is commonly used in modelling membrane gas separation, and which handles non-ideal behavior of CO₂ reasonably well [15].

Figure 1 shows the flow diagram of the case study implemented in the simulator’s flowsheet environment (COFE). The feed process (BG) is separated into an upgraded biomethane stream (BM) and an off-gas stream (OG), rich in carbon dioxide. The stage membranes (units MB1, MB2, and MB3) used the reference parameter values from Table 1. The compressor CP, placed after mixing the feed and recirculation, was assigned an isentropic efficiency of 72%. An aftercooler HX was set to reduce the temperature of the feed of the membrane stage MB1 to 40 ℃. The output pressure of the valve after the retentate of MB1 was set at the same pressure as the permeate of MB2 (1 bar).

The process simulator included some predefined units needed: compressor, heat exchanger, valve, and mixer units. Still, it does not include a membrane gas separation stage. The code for this unit operation was implemented in MATLAB and connected to the add-on “MATLAB unit operation” developed by Amsterchem [4] to define a user unit of a membrane gas separation stage able to communicate with the COCO environment. The equations used to create this user unit can be found in Appendix A. This subroutine was included as an additional user MATLAB unit operation file. The MATLAB unit operation also required an interface to receive the default stage parameters (stage area, component permeances, and permeate pressure). Ports were also added for the input and output streams of the stage. Additional ports for information were also included. These ports were a key component for the uncertainty analysis, as they allowed modifying the default parameters of the unit according to the information sent from an uncertainty generation unit.

Note that Figure 1 includes the control of the compressor output pressure to fix the methane product composition using the information streams xCH4 y Pout (not active in the unregulated scenario).

2.3. Application of Uncertainty and Sensitivity Methodology

Uncertainty and sensitivity analyses require characterizing the uncertainty of the most essential uncontrollable factors in the process. Appropriate ranges of variation must be determined, and probability density distribution functions must be proposed, with statistical parameters to be fitted. The uncertainty distributions of the factors were used to perform a global sensitivity analysis with the Morris method to rank the importance of the factors. The uncertainty analysis was also applied to determine the uncertainty expected for the product flow and composition and the compression energy requirements.

2.3.1. Characterization of the Uncertainty Sources

Feed variability can be significant in terms of flow and composition. When periodical measurements are available, an empirical distribution can be obtained for the flow rate, but typically, global sensitivity methods use log-normal, normal (truncated to zero), or uniform distributions. The choice between these options depends on the approximate knowledge of the empirical distribution. However, for an exploratory sensitivity study, a uniform distribution can be appropriate when a range can be defined, if we want to avoid introducing unsupported assumptions. As we are interested in studying the effect of variability, we will use uncertainty ranges of flow at two levels of variation, ±5% and ±10%, with respect to the reference value of the study. In the real system, the magnitude of the fluctuations will depend on the size of the gas holder tank placed before the membrane upgrading system.

As regards the variability of composition, some general ranges can be found in the literature. Guerrero et al. found that in biogas obtained from anaerobic digestion of biomass, the methane range can be 35-85%, and the carbon dioxide range 10-65%, the rest of the components being much less significant [16]. The type and mixture of the feedstocks, operating conditions, and temporal fluctuations cause the variability observed [17,18]. Our reference feed values fall inside the mentioned ranges. However, to assume the wide ranges shown as uncertainty ranges in the plant design may not be realistic, as a plant is fed with a restricted class of organic sources. In our case, only the uncertainty range for CO₂ molar fraction was used as a factor, as the analysis did not consider minor biogas components. Consequently, the CH₄ molar fraction is entirely dependent. Like in the flow analysis, for the CO₂ molar fraction, a uniform distribution will be used at two levels of variation, ±5% and ±10%, with respect to the reference value of the study, which falls inside the bounds previously commented.

Membrane performance is another great source of uncertainty in many processes. Permeability and selectivity of membranes for gas separation can be affected by compaction of the polymeric membrane due to pressure [19,20]. For example, CO₂ permeability can be reduced by 10% to 20% with respect to its nominal value. Besides, variations in the performance of different batches of membrane fibers have been observed [21]. The same membrane design can vary significantly depending on the manufacturing scale or production batch. However, it must also be considered that in a membrane stage, there is more than one membrane module; this means that the uncertainty in the performance of the stage is smaller than that of a single module. Permeability of CO₂ and CH₄ generally increases with temperature, showing an Arrhenius behavior [22]. Permeability increases because the diffusion through the polymer increases as the temperature rises, and it dominates over the possible solubility decrease. However, for CH_4, this effect is more critical than for CO₂, resulting in a decline in CO₂/CH₄ selectivity. For example, for polyimide membrane samples, an increase in the operating temperature from 35 ℃ to 45 ℃ increased CO₂ permeabilities between 7-12% but CO₂/CH₄ selectivity was reduced by 9-21% [23]. The temperature after the compressor plus cooling subsystem will have low variability, as a well-instrumented system can control the temperature of the stream between ±1 ℃. Therefore, the temperature variability of the membranes will be caused by the environment surrounding the membrane stages, and the magnitude of the variation will depend on the system’s insulation. Considering the ranges exposed for variability of CO₂ permeability, either due to manufacturing or affected by temperature, we will use for this permeability a uniform distribution with a variation of ±10% with respect to the reference value of the study. Membrane permeability and selectivity are not independent variables, as there is a trade-off between them. Robeson found an upper limit with a potential relationship between the permeability of CO₂ and the CO₂/CH₄ selectivity [24]. The fitting of this limit (Eq. 6) is a superior bound that relates both properties to the same material. We will use this expression to relate selectivity (S_CO2/CH4) and permeability to CO₂ (pmb_CO2) using the study’s reference values, assuming a variability of ±5% with respect to the calculated selectivity value.

$$S_{C{O}_2/C{H}_4}=S_{C{O}_2/C{H}_{4,ref}}·{\left({\displaystyle \frac{{pmb}_{CO2}}{{pmb}_{CO2,ref}}}\right)}^{-0.379}$$

(6)

Based on this discussion, the uncertainty factors and their associated uncertainties used in the uncertainty and sensitivity study are summarized in Table 2. For the feed, the variability of the flow and composition will be studied jointly through two levels of variation (FV=±5% and FV= ±10%).

2.3.2. Response Variables Studied

The effects of the uncertainty were studied on the flow and composition of the biomethane product stream (BM) and off-gas (OG), as well as on the utilities of energy and heat duty of the compressor. Table 3 shows the response variables considered.

2.3.3. User Module for Uncertainties

Using the COCO simulator, it is possible to perform local sensitivity analysis in which one or more variables of the process (stream variables or unit parameters) can be modified to study the effect on product streams or the energy requirements of the process units. However, global sensitivity analysis requires performing calculations for combinations of variables following probability distributions. To address this issue, we created a MATLAB unit that transmits the combinations of factors with uncertainty to the simulator for calculation (unit UG of Figure 2). This unit uses the deterministic values ​​of the problem as parameters when a case variable has a flag value equal to 0. A parametric study modifies the flag value from 1 to the number of executions required. When the flag value is 1, the values of all the factor combinations are loaded into the computer’s random-access memory from a comma-delimited file previously prepared in MATLAB using its statistical functions. The values ​​for the subsequent factor combinations are sent through information streams to the process units. A MATLAB user unit was incorporated for the inlet stream to modify its flow and composition around the inlet stream reference values. Another MATLAB user unit was created to distribute the uncertain permselective parameters to the three process stages (information dividers ID1 and ID2).

2.3.4. Morris’ Analysis of Global Sensitivity

Morris’ global sensitivity analysis was applied to the two scenarios of the case study. This method can be used as a screening method to determine the relative importance of the factors by computing elementary effects from distinct perturbations.

Figure 3 shows the workflow diagram used to perform the analysis. The information on reference values and uncertainty ranges was processed by the subroutine “morris” of the package sensitivity (version 1.30.1) [25] of R (version 4.5.1) to generate a design matrix with the combinations of factors to be calculated. The COFE flowsheet of the process, including the MATLAB unit operations for uncertainty generation and distribution, permits COCO to perform a parametric study. The results of this parametric study were manually exported to a comma-delimited file. The file was processed with the design matrix, again calling the subroutine “morris” in R software to obtain the sensitivity analysis results.

The design was of type OAT (one-at-a-time) with r=30 trajectories to have a relatively dense sampling. The number of evaluations for the k=4 uncertainty factors was N=(r+1)×k=150, which was considered an acceptable computing cost. Using a grid jump of less than or equal to half of the number of levels is recommended, so we started with a typical number of 4 levels with a grid jump of 2. Still, after analysis of the design matrix, these values were increased to avoid repeated combinations.

It is important to note that the method was applied to normalized responses obtained by dividing each response by its reference value. The reference values of the responses were the deterministic values obtained for the reference values of the factors. This allows for better comparison of relative sensitivity across responses as sensitivity results are interpreted in terms of percent change of the response.

2.3.5. Uncertainty Analysis for Determination of Tolerance Limits

In addition to determining the importance of the uncertainty factors, assessing their effects on the probability distribution of the response variables is interesting. This information can be used to obtain tolerance values that capture a certain probability content. Specifically, we are interested in symmetric tolerance limits (two-sided limits) to cover a probability content of 95%, which is the usual standard when no safety issues are involved.

The Monte Carlo propagation method was used to estimate the probability distribution of the outputs. A design matrix was generated in MATLAB by combining factors with random values that follow a uniform distribution within the ranges of Table 2. These data were arranged in a comma-delimited file and processed similarly to the workflow shown in Figure 3, but using MATLAB as the analysis tool.

The number of computer evaluations to achieve a specific confidence level can be obtained using Wilks’ formula. Using rank 1 (tolerance limits placed in the first and last values of the ordered vector) is usually acceptable when the tolerance limits are not close to safety limits. However, according to Porter [26], a higher rank statistic gives a more accurate and precise answer, and the improvement obtained by increasing the rank from 1 to 2 is greater than that obtained by increasing from rank 2 to 3. Therefore, we used the formula considering rank 2, which means that the lower tolerance limit is defined by the value in the second position of the ordered vector of results, and the upper tolerance limit is the value in the last but one position. A confidence level of 99% requires performing 198 simulator evaluations (68 more than for rank 1).

3. Results

3.1. Deterministic Simulation Results

Table 4 shows the results obtained using the COCO simulator for the process configuration and the reference values in Table 1. Comparing the values of the flows of biomethane and off-gas for the non-controlled scenario with those obtained by Scholz [14], we observed differences of 0.12% and −0.16% for the same configuration and reference parameters. These minor differences are probably a consequence of the most elaborate model of the mentioned work in terms of the internal structure of the stage. However, compared with the variability posteriorly observed in the uncertainty study, they are small enough to validate the process model for capturing the uncertainty characteristics of the process.

3.2. Sensitivity Analysis

3.2.1. Refined Matrix Design

Using four levels with a grid jump of 2 yielded a design matrix with repeated combinations of factors. This is not desirable when evaluating a computer code, as redundant effects are obtained. Therefore, we increased to 8 levels and a grid jump of 4 for an experimental plan with no repetitions (Figure 4).

Using an Intel Core i9 at 3.50GHz, the computing time in the simulator for the 150 combinations to calculate was 5.91 min for the non-controlled scenario and 75.3 min for the controlled case.

3.2.2. Morris Sensitivity Analysis of the Flow and Composition of the Product Streams

The Morris sensitivity analysis provides three main indices: the mean of elementary effects (μ), the mean of the absolute values of the elementary effects (μ*), and the standard deviation of the elementary effects (σ). The μ* index is usually used over μ as it avoids the cancellation between positive and negative effects, and indicates, on average, how much the input changes the output. The σ index measures the nonlinearity, i.e., the interaction between effects; low σ values indicate that the effect is linear or additive, high σ values indicate interaction and nonlinearity. Parameters with low values of μ* and σ are unimportant. A diagonal (σ=μ*) is also included in the figures to identify if the effects tend to be predominantly linear (σ<<μ*) or nonlinear (σ>>μ*).

Figure 5 shows the analysis for the response of the biomethane product flow (FBM). Let us consider the non-controlled scenario at a feed variability FV = ±5. The effects of the membrane parameter uncertainty (given by factors PCO2 and S) are negligible because of their low μ* and σ values. The flow of biogas FBG is the dominant factor, strong and with a linear effect; the biogas composition xCO2BG is lower, but not negligible. Comparing this situation with that obtained at FV=±10%, we can see that μ* has almost doubled, and that σ has increased notably, which suggests non-linear effects and interaction with other variables. The situation is similar to that in the controlled scenario.

Figure 6 shows the analysis for the response of the off-gas product flow (FOG). Thanks to the normalization of this variable, we observed similar index values. The result of the analysis is very similar to that performed for the FBM variable, except that now the effect of the CO2 fraction in the feed becomes somewhat more important than the biogas flow.

Figure 7 shows the analysis for the response of the methane fraction in the biomethane product (xCH4BM) for a non-controlled scenario. This response was not analyzed in the controlled scenario because it is set to a fixed value. We can see that for low feed variability (FV=±5), the membrane parameter uncertainty predominates over the feed variability as the permeance factor PCO2 has the higher μ* value; however, the variability of the selectivity with respect to the value expected for each permeance is not significant. Nevertheless, for FV=±10, the effect of the flow of biogas becomes practically as important as that of the membrane permeance. The smaller μ* values currently observed in the analysis of the composition response compared to those observed for the flow responses do not imply less significant effects. Note now that they represent relative composition variations with respect to a composition fraction value close to unity.

Figure 8 shows the analysis for the response of the CO₂ fraction in the off-gas (xCO2OG). The results for the non-controlled and controlled scenarios are very similar. Like in the case of the xCH4BM response, at FV= ±5% the effect of the CO₂ membrane permeance is the most important. Similarly, at the higher feed variability of FV= ±10%, feed variability becomes as important as membrane permeance uncertainty.

3.2.3. Morris’ Analysis of the Process Utilities

Figure 9 shows the analysis for the EDC response. In all situations, the uncertainty of biogas flow is the most critical factor (higher μ*), doubling when the feed variability is doubled. For the controlled system, the non-linear effects are higher (higher σ).

An additional study was performed on the mechanical power ratio to the flow of processed biogas (EDC/FBG). In Figure 10, we can see that, compared with the EDC analysis, a large part of the effect of the variability due to FBG has been cancelled, and the μ* values are smaller. Now, we appreciate different behaviors for the non-controlled and the controlled scenarios. In the non-controlled scenario, the effects of all factors are similar; however, for the controlled scenario, the impact of the permeance of CO₂ is higher than the effect of FBG at FV= ±5% but it takes on a comparable value at FV= ±10%.

Morris’ analysis results of the thermal load to be removed in the compressor (HDX and HDX/FBG) were virtually identical to those of the mechanical power (Supplementary Materials figures S1 and S2). This is logical since both variables depend almost exclusively on the flow through the compression system.

3.3. Determination of Tolerance Limits of the Responses

Table 5 shows the uncertainty analysis results to obtain the lower and upper tolerance limit values for a 95% probability content and a 99% confidence level (LTL and UTL). The approximate form of each distribution can be represented, but for the sake of concision, two coefficients, RLG and RUG, were created as the relative gaps between lower and upper tolerance limits and the mean value. The coefficient of variation (CV), defined as the ratio of the standard deviation to the mean, was also included.

The mean values obtained for all variables did not differ significantly from those obtained for the reference factors (Table 4). The coefficients RLG and RUG had similar values for most responses, implying high symmetry for most distributions. For these responses, we can see that CV is approximately half of the relative gap values because, for a normal distribution, twice the CV is approximately 95% of the probability content. The highest difference in RLG and RUG was found for the response xCH4BM under the non-controlled scenario, with RLG representing 59% of the total tolerance gap. A reason can be found in the fact that the values of this variable are near an upper physical bound, and the left tail of the distribution predominates. The values of relative tolerance gaps were much higher for the flows (FBM and FOG), and utility variables related to them through the recirculated flow (EDC and HDX). For these responses, doubling the feed variability approximately doubles the gap value. The tolerance gap and CV were smaller for the composition-related variables (xCH4BM and xCO2OG), and the impact on the tolerance gap of doubling the feed variability was smaller.

The controlled and non-controlled scenarios yield similar results for many variables, but differences were observed in some specific responses. The tolerance interval was null for the controlled variable xCH4BM, as it is a system property, but this fact did not substantially increase the tolerance gap for the response xCO2OG. However, the tolerance gaps for the ratios EDC/FBG and HDX/FBG are much higher for the controlled scenario because of the necessary regulation effort of the pressure to achieve the set point of methane composition in the biomethane product.

4. Discussion

This paper shows the feasibility of conducting sensitivity and uncertainty analysis of processes combining free tools such as the COCO simulator and R. However, it is still necessary to model some specific units using programming software like MATLAB or other free software. This analysis can be largely automated using parameterization and reading from bridge files.

Studies on uncertainty and sensitivity in areas related to our process, such as anaerobic digestion [27], or gas membrane plants [28,29], indicate the usefulness of Morris’ analysis as a screening method to rank the importance of the factors. In this way, factors with a low influence on the uncertainty of a response can be potentially fixed without affecting the prediction.

However, Morris’ analysis has limitations as it does not provide variance-based quantitative indices. Other global sensitivity methods, like Sobol or eFAST, are more robust and can be used in a second phase on a subset of the influential variables to get quantitative sensitivity indices and evaluate interactions. Nevertheless, it should be considered that these methods involve a high number of computing evaluations, and difficulties may arise when an evaluation fails in the calculation process.

In our case, we have seen that the effects of feed variability are most influential in the uncertainty of product flows and power and heat demand responses. However, for product composition, the parameter uncertainty of membrane permeance can be most influential when feed variability is under ±10%. The important effect observed of the flow and composition variability of the feed agrees with what has been determined using parameter variation in other studies in a biogas processing plant [11].

The uncertainty analysis yielded lower and upper tolerance limits to cover a probability content of 95%. A study about the life cycle assessment of a biomethane production recommends its use to test the robustness of the results [30]. According to [31], simplified models are likely to overpredict performance compared to experimental results; therefore, the uncertainty quantification can help quantify the estimates of uncertain designs. In conclusion, knowing the tolerance limits is beneficial information for process design. Let us suppose the obtained tolerance limits are bound to be unacceptable. In that case, we can modify the process (unit dimensions or structure) or the operating parameters and include additional control to ensure the correct functioning of the process. In our case, stream flow showed high variability, so the information of their probability distributions can be used to size the gas-holding tanks before the upgrading process. Seemingly, the lower tolerance limit of the power demand of the compressor could be used to specify its nominal power. For the composition of the product stream, smaller tolerance limits were obtained. However, modifications of the process, like increasing the membrane stage areas or system pressures, may be necessary if a tolerance limit is close to an undesirable value. One crucial observation must be considered when applying the uncertainty analysis procedure to a process. For simplicity, we have considered in the case study that the variability of feed flow rates is independent of the composition variability, which may not be the case in some real-life situations. As Porter [26] discussed, when there is a correlation between variables, the joint distribution should be propagated, since otherwise, the uncertainty is overpredicted.

For the case study, different uncertainty and sensitivity characteristics were obtained when comparing the controlled and non-controlled scenarios of the same system. Similar situations can happen for changes in the configuration or the relative size of the different process units. This encourages us to focus future studies on understanding how the process configuration influences the uncertainty of the outputs.