Development of engine activity cycles for the prime movers of unconventional natural gas well development

ABSTRACT

With the advent of unconventional natural gas resources, new research focuses on the efficiency and emissions of the prime movers powering these fleets. These prime movers also play important roles in emissions inventories for this sector. Industry seeks to reduce operating costs by decreasing the required fuel demands of these high horsepower engines but conducting in-field or full-scale research on new technologies is cost prohibitive. As such, this research completed extensive in-use data collection efforts for the engines powering over-the-road trucks, drilling engines, and hydraulic stimulation pump engines. These engine activity data were processed in order to make representative test cycles using a Markov Chain, Monte Carlo (MCMC) simulation method. Such cycles can be applied under controlled environments on scaled engines for future research. In addition to MCMC, genetic algorithms were used to improve the overall performance values for the test cycles and smoothing was applied to ensure regression criteria were met during implementation on a test engine and dynamometer. The variations in cycle and in-use statistics are presented along with comparisons to conventional test cycles used for emissions compliance.

Implications: Development of representative, engine dynamometer test cycles, from in-use activity data, is crucial in understanding fuel efficiency and emissions for engine operating modes that are different from cycles mandated by the Code of Federal Regulations. Representative cycles were created for the prime movers of unconventional well development—over-the-road (OTR) trucks and drilling and hydraulic fracturing engines. The representative cycles are implemented on scaled engines to reduce fuel consumption during research and development of new technologies in controlled laboratory environments.

Introduction

The United States has experienced growth in the natural gas industry over the past decade. Annual gross withdrawals of natural gas in the United States have increased 25% since 2005 from 23 million to over 31 million cubic feet (U.S. Energy Information Administration [EIA], 2016). The U.S. Energy Information Administration (EIA) also estimates that there are about 2.276 trillion cubic feet (Tcf) of recoverable resources of dry natural gas in the United States, which at the current rate of consumption is enough to last about 84 years (EIA, 2015). These advances are a product of the new technologies of horizontal drilling and hydraulic fracturing. In order to utilize these resources, strategies are being implemented to reduce costs and emissions. One strategy involves displacing diesel fuel with natural gas in the prime movers associated with well development via either dual-fuel or dedicated natural gas applications. These prime movers, or major diesel consumers, include heavy-duty over-the-road (OTR) trucks and high-horsepower drilling and fracturing engines.

The process of unconventional well development consists of four major steps: planning and preparation of the well pad site, drilling, casing and cementing, and completion and stimulation (Canadian Society for Unconventional Resources, 2012). The primary consumer of fuel during preparation of the site is OTR trucks. Heavy-duty OTR trucks transport gravel, dirt, and other materials to and from the site. Vertical and horizontal drilling applications use high-horsepower engines to power rigs. Heavy-duty OTR trucks also haul water, drilling fluids, and other supplies during drilling. The drilling company often performs casing and cementing that involve the same engines as drilling, with an increase in truck traffic. Hydraulic fracturing is the most energy-intensive step in the process and involves the use of high-horsepower fracturing engines and large amounts of truck traffic for the delivery of water, sand, and fracturing chemicals. It is estimated that a Marcellus Shale well pad with four wells requires 20,000–30,000 truck trips during the completion of those wells—yielding up to 4 million truck trips per year for the entire Marcellus region (Carr et al., 2012). Drilling engines operate continuously from rig setup until drilling, casing, and cementing are completed. The Texas Commission on Environmental Quality, in a survey of 1261 wells, found that the average drilling duration for a well was between 8 and 27 days depending on the depth of the vertical section and length of the horizontal sections. On average, the drilling rigs consisted of 2.15 engines, with an average size of 1381 horsepower (hp) each. These engines operate on average 62.6 hr per 1000 feet drilled at an average load of 48.5% (Baker and Pring, 2009). A report by Global Hunter Securities in collaboration with Prometheus Energy states that the average drilling rig uses between 62.5 and 83.3 gallons of diesel fuel per hour, resulting in 1500–2200 gallons of diesel per day (Kelly, 2012). A study by Argonne National Laboratory (ANL) assumed that a rig ran on one 2000-hp engine at 45% load and consumed 0.06 gal/hp/hr in its analysis. The study focused on four major shale plays and the fuel consumed in these plays depended mostly on the depth of the well and the amount of time required for drilling. Diesel fuel consumption by shale play ranged from 43,041 gallons for the Fayetteville Shale to 85,845 gallons for the Haynesville, with the Barnett and Marcellus within this range (Clark et al., 2012).

High-horsepower diesel engines also power hydraulic fracturing pumps. A typical fracturing fleet features engine capacities over 20,000 hp. An average site requires 8–12 fracturing pumps, but some require up to two dozen, each rated between 1500 and 2500 hp. Fracturing a well can consume between 2000 and 8000 gallons of diesel fuel per day. The University of Michigan performed site surveys in the Marcellus and Eagle Ford shale plays and showed average fleet sizes of 14 and 16 engines, respectively. The average fuel consumption during stimulation of a well in these shales was 20,800 and 22,100 gallons, respectively. The engines at both sites were Caterpillar 3512Cs (Peoria, IL) rated at 2250 hp.

Prime mover engines are subject to emissions compliance regulations as tested over specific engine cycles. The primary engine test cycle over which on-road heavy-duty diesel engines are certified is the transient Federal Test Procedure (FTP) engine dynamometer test cycle defined in Code of Federal Regulation (CFR) Title 40 Subchapter C Part 86 Appendix I. The new greenhouse gas (GHG) standards also use the FTP cycle for certification of vocational engines. This cycle is a series of speed and torque points normalized to a percentage for applicability to all engines (DieselNet, 1999). As a result of the consent decrees, engine certification also involves the Supplemental Emission Test (SET) and not-to-exceed (NTE) testing, in addition to the FTP cycle (DieselNet, 2015). Certification of nonroad engines is also based on the emissions standards outlined by the CFR. Engine tests occur on dynamometers based on the procedures of CFR Title 40 Chapter I Subchapter U Part 1039. Certification of Tier 1–3 engines requires operation over steady-state tests. The International Organization for Standardization ISO 8178 is the standard for nonroad engine emissions testing. This test, also known as the Non-Road Steady Cycle (NRSC), consists of five modes of steady-state operation and is applicable to both drilling and hydraulic fracturing engines (DieselNet, 2001).

Although these cycles are well defined, they may not be representative of in-use operation. Additionally, research and development on full-scale engines is energy intensive and cost prohibitive from a fueling standpoint. Therefore, our goal was to collect representative in-use activity data and develop a method to create representative cycles for application on scaled engines within our test cell.

Data collection

OTR trucks

Collection of data from OTR trucks servicing the unconventional well development industry was performed over a 6-month period, beginning June 6, 2014, and concluding January 9, 2015. The vehicles targeted for data collection were heavy-duty diesel vehicles that were travelling to and from unconventional natural gas wells in the Marcellus Shale region. These heavy-duty OTR trucks consisted of those hauling water, sand, and gravel. Over 100 companies in the Marcellus and other shale plays were contacted. Seven different companies located throughout Ohio, Pennsylvania, and West Virginia provided access to 25 vehicles for data collection. Vehicle instrumentation was random within each fleet, and no specifications on vehicle age, mileage, or other variables were required. Of the 25 vehicles, 18 were water haulers, 6 were sand haulers, and 1 carried gravel for site preparation. The engines populating the vehicles consisted of 18 Cummins ISX-15 (Columbus, IN) engines, 4 Caterpillar C-15 (Peoria, IL) engines, 2 Volvo D13 (Greensboro, NC) engines, and 1 Mack MP8 (Greensboro, NC). Table 1 provides a summary of the vehicles instrumented. It is noted that the vehicles all operated within the greater Marcellus region and their activity may not be representative of all unconventional natural gas shale plays. The number of truck trips varies by well pad, region, infrastructure, and other logistical variables. Three studies have estimated that a majority of truck trips per well development are attributed to hauling water to and from sites (67–78% of all truck trips) (Maryland Department of the Environment [MDE], 2014; Felsburg Holt and Ullevig, 2013; Massachusetts Institute of Technology [MIT], 2011). Of the 25 random vehicles in this program, 68% were water-hauling trucks.

Table 1. Well service vehicle data collection information.

Truck no.	Comp.	Hauling service	Make	Model	Size (hp)	Chassis	Vehicle year	Start date	End date
1	A	Water	Cummins	ISX-15	500	Peterbilt	2012	5-Jun	30-Jun
2	A	Water	Cummins	ISX-15	500	Peterbilt	2012	7-Jul	16-Jul
3	A	Water	Mack	MP8	505	Mack	2011	7-Jul	16-Jul
4	A	Water	Cummins	ISX-15	500	Freightliner	2012	7-Jul	16-Jul
5	A	Water	Cummins	ISX-15	500	International	2005	7-Jul	16-Jul
6	B	Gravel	Volvo	D13	435	Volvo	2002	29-Jul	31-Jul
7	C	Sand	Cummins	ISX-15	500	Kenworth	2013	1-Aug	13-Aug
8	C	Sand	Cummins	ISX-15	500	Kenworth	2013	1-Aug	13-Aug
9	C	Sand	Cummins	ISX-15	500	Kenworth	2009	1-Aug	13-Aug
10	C	Sand	Cummins	ISX-15	500	Kenworth	2008	1-Aug	13-Aug
11	D	Water	Caterpillar	C-15	550	Peterbilt	2012	15-Aug	28-Aug
12	D	Water	Cummins	ISX-15	500	Kenworth	2013	15-Aug	26-Aug
13	D	Water	Caterpillar	C-15	550	International	2009	15-Aug	28-Aug
14	E	Sand	Caterpillar	C-15	550	Peterbilt	2006	18-Sep	24-Sep
15	E	Sand	Cummins	ISX-15	500	Peterbilt	2009	26-Sep	22-Oct
16	F	Water	Cummins	ISX-15	500	Peterbilt	2013	19-Nov	2-Dec
17	F	Water	Cummins	ISX-15	500	Peterbilt	2013	19-Nov	2-Dec
18	F	Water	Cummins	ISX-15	500	Peterbilt	2013	19-Nov	2-Dec
19	F	Water	Cummins	ISX-15	500	Peterbilt	2013	19-Nov	2-Dec
20	G	Water	Cummins	ISX-15	500	Freightliner	2012	17-Dec	5-Jan
21	G	Water	Caterpillar	C-15	550	Peterbilt	2010	17-Dec	5-Jan
22	G	Water	Cummins	ISX-15	500	Kenworth	2011	17-Dec	5-Jan
23	G	Water	Volvo	D13	500	Kenworth	2012	17-Dec	5-Jan
24	G	Water	Cummins	ISX-15	500	Peterbilt	2010	5-Jan	9-Jan
25	G	Water	Cummins	ISX-15	500	Peterbilt	2013	5-Jan	9-Jan

Data collection occurred at a rate of 1 Hz using DAWN J1939 Mini Loggers (HEM Data, 2015). Although dozens of engine control unit (ECU) parameters were available and ultimately collected, the only parameters relevant to the creation of engine test cycles under this work were engine speed and percent load at current speed. The resulting data set included 4,724,800 data points, or over 54 days of data. The International Society of Automotive Engineers (SAE) under J1939 protocol defines the parameters as follows:

• Engine percent load at current speed—SPN 92: The ratio of actual engine percent torque (indicated) to maximum available at the current engine speed, clipped to zero torque during engine braking.

• Engine speed—SPN 190: Actual engine speed calculated over a minimum crankshaft angle of 720 degrees divided by the number of cylinders (“Vehicle Application Layer”; SAE, 2015).

Drilling

Drilling activity data collection occurred during the drilling of two horizontal wells in Westover, West Virginia. The collection of data occurred over a period of 11 days from September 9, 2015, to September 21, 2015. The rig used three Caterpillar 3512C generator units, rated at 1101 kW. These engines were outfitted with dynamic gas blending (DGB) kits allowing the engines to run in either dual fuel or diesel only mode. Data collected occurred on a single engine that operated continuously during rig setup, preparation, pipe tripping (PT), and steady-state drilling (SSD).

The engine control unit broadcasted data, which were recorded using a nine-pin Deutsch connector and a VIA HDV100A1 from B&B Electronics (Ottawa, IL). An in-house software recorded SAE J1939 parameters. After conversion of data to CSV format, it was determined that only 233 of the 311 hr were valid.

Use of activity logs from the drilling company assisted in sorting data. Table 2 presents a breakdown of the engine activity. This breakdown shows the amount of time spent performing each type of activity necessary for the completion of drilling. Note that data were obtained during the horizontal portions of the wells and may not be representative of the vertical or conventional drilling portion. In addition, drilling activity varies by well, depth, operator, and shale play; but in discussions with the site operators, the activities were typical for the Marcellus region.

Table 2. Drilling activity breakdown.

	ECU collection		Overall
Activity type	Hours	Percent	Hours	Percent
Rig up	4	1%	101	12%
Drilling	160	51%	216.5	25%
Reaming	3	1%	8	1%
Coring	0	0%	19	2%
Condition and circulate	16.5	5%	56.5	6%
Trips	21	7%	126.5	14%
Lubricate rig	10.5	3%	14.5	2%
Repair rig	9	3%	17.5	2%
Cut off drilling line	1	0%	3	0%
Wireline logs	0	0%	105.5	12%
Run casing and cement	27	9%	77	9%
Wait on cement	8.5	3%	8.5	1%
Nipple up B.O.P.	5	2%	16.5	2%
Test B.O.P.	12	4%	25.5	3%
Directional work	13	4%	21	2%
Other	20.5	7%	58.5	7%
Total	311	100%	875	100%

Note. B.O.P. = blowout preventer.

Hydraulic fracturing

Collection of hydraulic fracturing activity occurred during two separate campaigns with two different fracturing fleets. Collection of data was similar to drilling: with a nine-pin Deutsch connector, serial cable, VIA model HDV100A1 from B&B Electronics (Ottawa, IL), and a laptop computer. The first data collection effort occurred over 21 days between July 19, 2015, and August 8, 2015, from a well pad in central West Virginia. Data collection occurred on four different hydraulic fracturing pumps during the stimulation of six wells. Cummins QSK50 engines, rated at 1678 kW, powered these pumps.

The second data collection effort occurred at a well pad near Morgantown, West Virginia. Data collection occurred on a Caterpillar 3512B-HD engine rated at 1678 kW over a period of 7 days from November 5, 2015, to November 12, 2015. The engine included a dynamic gas blending kit allowing it to operate in either diesel only or dual fuel mode. We collected data from 12 fracturing stages. We eliminated data from three stages due to erratic values. The hydraulic fracturing engines operated at rated speed except for when the engines were idling between stages, which was limited. Table 3 presents a summary of the data. Since the engines typically operated at rated speed, the most important engine parameter was percent load at current speed. Note that all OTR and drilling data were collected at 1 Hz. In some cases of hydraulic fracturing, data were collected at 10 Hz (due to hardware), but a 10-point average was used to down-sample the data for use as 1 Hz. Both data sets were collected in the Marcellus Shale region. In addition, hydraulic stimulation activity varies by well, depth, operator, and shale play; but in discussions with the site operators, the activities were typical for the Marcellus region.

Table 3. Hydraulic fracturing data collection information.

Stage	Site	Engine make	Model	Data rate (Hz)	Length of stage (sec)	Average speed (rpm)	Average load (%)	Average power (kW)
1	A	Caterpillar	3512B	1	5209	1842	66.37	1114
2	A	Caterpillar	3512B	1	4672	1873	52.32	878
3	A	Caterpillar	3512B	1	4001	1937	79.57	1335
4	A	Caterpillar	3512B	1	4914	1866	60.45	1014
5	A	Caterpillar	3512B	1	5226	1936	79.69	1337
6	A	Caterpillar	3512B	1	4638	1936	84.81	1423
7	A	Caterpillar	3512B	1	6271	1842	66.57	1117
8	A	Caterpillar	3512B	1	5134	1823	73.76	1238
9	A	Caterpillar	3512B	1	5234	1936	77.28	1297
10	A	Caterpillar	3512B	1	4879	1602	73	1223
11	B	Cummins	QSK50	10	4860	1937	56.31	945
12	B	Cummins	QSK50	10	5810	1948	50.71	851
13	B	Cummins	QSK50	10	10770	1466	49.5	831
14	B	Cummins	QSK50	10	6130	1947	50.47	847
15	B	Cummins	QSK50	10	5130	1950	49.27	827
16	B	Cummins	QSK50	10	2380	1948	52.82	886
17	B	Cummins	QSK50	10	5910	1949	52.5	881
18	B	Cummins	QSK50	10	5890	1946	48	805
19	B	Cummins	QSK50	1	2507	1949	54.48	914
20	B	Cummins	QSK50	1	3093	1947	53.01	890
21	B	Cummins	QSK50	1	1555	1950	50.09	841
22	B	Cummins	QSK50	10	6850	1944	50.26	843
23	B	Cummins	QSK50	10	5080	1950	54.73	918
24	B	Cummins	QSK50	1	3657	1942	51.31	861

Cycle development

Numerous approaches are available to create engine or chassis dynamometer test cycles, including Markov chain Monte Carlo simulations, micro-trips, or day-in-the-life cycles (Clark et al., 2007; Steven, 2001; Ullman, 1999). The technique selected was a stochastic second-by-second concatenation, which utilized Markov chain theory. We applied the same method to all three types of cycle development. We selected a starting period of engine idle and from there the Markov chain technique was allowed to progress freely based on the transition probability matrix of each of the respective data sets. We imposed a time limit of 30 min (1800 sec) on each of the cycles. The FTP cycle is 20 min long (DieselNet, 1999); however, research suggests that longer Markov chains were more representative when attempting to match a distribution (Geyer, 1992; Smith, 2003). Previous studies on cycle development used 30 min as the time of the cycle (Steven, 2001). We selected 30 min to align with Steven (2001) and to provide a longer than conventional engine cycle length. Longer cycles should be examined under future work, but due to fuel restrictions, we retained 30 min as the upper limit such that a single cold start and three repeated hot starts could be conducted using the same compressed natural gas cylinder bank.

Review

Wide varieties of sciences use Markov chain Monte Carlo (MCMC) simulations as a way to estimate probabilities and perform numerical integration. MCMC is “a general method for the simulation of stochastic processes having probability densities known up to a constant proportionality” (Geyer, 1992).

The theory behind MCMC is subdivided into Markov chain techniques and Monte Carlo simulations. A Markov chain is a stochastic process that has stationary transition probabilities. The Markov chain is a string of concatenated values from a defined “state space,” or data being sampled (Smith, 2003). A Markov chain is often described by saying that there exists some state, or event, and that state contains a transitional probability matrix, which gives the likelihood of the next state.

The next state in the chain depends only on the current state and its transition probability matrix. The Markov chain is “stationary” if any state in the chain has a transitional probability distribution that is independent of the location of the state. Thus, probability of reaching the next state is the same for all states in the chain and likely independent on initial distribution. If a Markov chain is not stationary, then the transition probabilities change as the chain moves further from the initial state. As the length of the chain increases, the distribution of the states becomes less dependent on the initial state and more dependent on the transition probabilities (Geyser, 2011).

If an overall distribution defines each states transition probability matrix and if the chain becomes “long,” the distribution of the chain itself will approach the overall distribution (Smith, 2003). The most effective way to reach this desired distribution is by allowing the Markov chain array to grow large. The actual size necessary to properly represent a distribution will be dependent on the size of the data set, and the problem at hand. Cycle development extensively employs Markov chains, as shown in the next section. MCMC simulations expand on the idea of Markov chains, but flip the objective of the problem. Rather than attempting to find the distribution of a specific Markov chain, MCMC starts with a defined distribution and attempts to create a Markov chain that has this same distribution (Smith, 2003). When this is the case, one knows the distribution of interest and the Markov chain proceeds in the same stochastic manner until it has an array length, which is long enough to adequately represent the distribution. Researchers employ MCMC techniques to genomics, variable selection in regression, spatial statistics, longitudinal studies, mathematical chemistry, statistical physics, or when more simplistic techniques cannot be used.

The work performed by Smith is the most relevant example of engine test cycle development using a MCMC method for heavy-duty on-road vehicles (Smith, 1978). Researchers commonly use MCMC for generation of drive cycles. Ohio State University’s Center for Automotive Research generated driving cycles of plug-in hybrid electric vehicle (PHEV) using a Markov chain technique. Researchers collected data from nine PHEV passenger cars consisting of over 70,000 miles of data that were segmented into 530 “sessions” defined from key-on to key-off (Gong et al., 2011). Dembski et al., also of Ohio State University, used similar techniques and calculations of sequence occurrence probabilities to analyze drive cycles for passenger cars. Kinematic sequences defined by velocity and acceleration were developed and linked together based on the probability of each type of sequence occurring, essentially utilizing a Markov chain to create drive cycles (Dembski et al., 2002). Researchers at University of California–Davis used Markov chains to produce drive cycles for a number of research programs from light-duty vehicles and chase vehicles (Lin and Niemeier, 2002; Lin and Niemeier, 2003; Dai et al., 2008).

Data processing and second-by-second concatenation

OTR trucks

We used MATLAB (MathWorks, Natick, MA) for data processing and cycle development algorithms. Normalization of engine data is necessary to compare engine operation between different sizes, makes, and models. All normalized data points contain a value from 0 to 1, representing 0–100% of the maximum possible value. The parameter percent load at current speed is an example of a normalized parameter. This parameter is based on the lug curve of the engine and is broadcasted as a percentage of the maximum possible load at a given engine speed. Previous studies considered this value similar to percent engine torque, especially at high loads (Clark et al., 2007). We used this parameter as an equivalent value of normalized percent torque for this methodology. Engine speed normalization was based on the equation used for denormalizing speed from the CFR Title 40 Chapter I Subchapter C Part 86 Subpart N Section 86.1333. Normalizing engine speed was performed using eq 1. It is noted that these engines were all diesel fuel engines, which typically include a factor of 112% of rated speed. Any data within ±2% of 0% speed were set to 0, and the same was done for engine percent load. Any speed between rated and 102% of rated was set to rated speed or 1. Of the over 4 million data points collected, speeds between rated and governed only represented 0.29% of activity and were excluded. It is also noted that all OTR data were collected in the Appalachian Basin and OTR activity may differ by shale play.

(1)

Elimination of some data was necessary because the engines often broadcasted false data or trended to an infinite value if the size of the file became too large. Data were deemed invalid if the broadcasted value was outside of the defined limits of a specific parameter and the cause was not known. The second source of invalid data was due to hardware and software issues where an error in recording caused subsequent data to increase in value through the duration of the recording. After filtering of invalid data, speed and load were arranged in a two-column array and concatenated, for each engine.

Data binning used 5% increments for both speed and load. The bins contained values of speed and load greater than or equal to the governing value of the current bin and less than the governing value of the next highest bin. Figure 1 shows the distribution of data. Note that nearly 20% of the activity data occurred in a single bin. This was bin (1, 3) and it contained nearly all of the points that occurred at idle engine speed.

Figure 1. Truck data distribution.

Once the algorithm binned the current point in the distribution matrix, analysis of the next point occurred. The algorithm binned this data point into the transition matrix of the bin of the previously selected point and the overall distribution matrix. The program used a multilayered structure to allow each bin of the overall distribution matrix to also contain a secondary transition matrix of the same size. This was necessary because of the stochastic logic of the Markov chain technique. Once the entire length of the array was binned, a normalized distribution was created by dividing the size of each data bin by the total length of the array. Statistics of the data set were calculated. In addition, tabulated were average speed and load, and time at idle. The idle data points were excluded from a second calculation so that the non-idle speed and load averages could be determined (see Table 4).

Table 4. Trucking data statistics.

Parameter	Value
Average speed (%)	29.2
Average load (%)	23.1
Average non-idle speed (%)	46.9
Average non-idle load (%)	28.3
Normalized idle time (%)	37.6
Normalized non-idle time (%)	62.4
Total time (sec)	4,724,800

Drilling engines

Drilling engines operated at nearly constant rated speed and in this study operated at 1200 revolutions per minute (rpm). The binning logic used was similar to that used for OTR trucks except that it focused only on load. Figure 2 shows the distribution of engine load bins.

Figure 2. Drilling data distribution.

We further subdivided data into two categories of load, high load and low load, for the purpose of calculating statistics. High load was any value at or above 40% and low load below 40%. We selected 40% to ensure that all steady drilling operations were high load. Some of these events occurred at steady-state values below 50%, but none occurred below 40%. The idle time of the data was also calculated and included operation at 0% load—a small overall portion of the time. The average value and the normalized time of high- and low-load operations were calculated. Table 5 shows the statistics of these data.

Table 5. Drilling data statistics.

Parameter	Value
Average load (%)	43.0
Average high load (%)	60.2
Average low load (%)	19.0
Normalized high-load time (%)	58.4
Normalized low-load time (%)	41.2
Normalized idle time (%)	0.4
Total time (sec)	840,238

Hydraulic fracturing engines

Hydraulic fracturing data binning used the same methodology as drilling. Data bins were in 5% increments, forming an overall distribution matrix by stepping through each point of the array. Figure 3 shows the distribution of these data.

Figure 3. Fracturing data distribution.

Data were sorted into idle (0% load), low load (points below 40% load), and high load (points at or above 40% load). The average overall load, high load, and low load as well as normalized idle, low-load, and high-load times were calculated. Table 6 provides the tabular statistics of these data.

Table 6. Fracturing data statistics.

Parameter	Value
Average load (%)	35.5
Average high load (%)	61.5
Average low load (%)	20.3
Normalized high-load time (%)	46.8
Normalized low-load time (%)	33.2
Normalized idle time (%)	19.9
Total time (sec)	239,905

MCMC

The cycle began with a period of idle, and from there the Markov chain technique progressed freely based on the transition probability matrix of each of the respective data sets. Based on the nature of driving and observations made during data collection, it was determined that, at a minimum, the first 30 sec should be idle. Therefore, the algorithm sampled points between 1 and 30 from the zero-speed row of the distribution matrix and weighted on the transition probability matrix of the current bin.

Once 30 sec of idle were completed, the overall distribution matrix was open for sampling and the Marko chain proceeded. We employed a weighted sampling technique to pseudorandomly select the next bin of the cycle. Columns and rows had separate weights. The weights of each were a function of the ratio of the time in a given bin to the total time spent in the column or row. Then, the row of the next bin was selected by sampling 1–20 (the number of rows), based on the weights of the rows developed previously. Selection of the next column bin used the same manner. Once the row and column of the next bin were found, the bin was defined and the algorithm randomly selected from the bin’s population. This point represented the next point in the cycle and was added to the end of the current cycle array. The transition matrix of the selected bin then governed the weights used for selection of the next bin.

The Markov chain continued to sample using the defined method up to 1770 sec. It was determined that the vehicle should return to a state of idle 30 sec before the end of the test to simulate real-world activity and also mark the end of the cycle. If the speed bin of the cycle at this point was above 10, meaning normalized speed was above 50%, an intermediate point was required at point 1771, to create an acceptable derivative of speed. Sampling for that bin moved to the row that was five bins lower than the current bin. Once the current speed bin was below 10, the speed row was forced to 1, the idle speed row, and the process continued as it did in the first 30 sec of cycle construction, with selection of points based only on the load columns until the cycle reached a length of 1800 sec.

With the cycle completed, we imposed limitations on the cycles in two steps. The first was to search the cycle and determine if there were any single points of nonzero speed between points of zero speed. This type of activity made segmentation of cycles, discussed later, difficult and was not necessary to represent real-world engine operation. Once located, we forced these points to zero speed. As an example, in the creation of 100 cycles, each with a length of 1800 sec, only 29 of these points occurred (0.016%) and were forced to zero. The second step was that all points below 2% normal speed or load were forced to zero. The reduction of these points did not change the overall distribution of the cycle and had no effect on idle time, and we verified this by applying the change to a randomly generated sample population of 100 cycles.

We created a distribution of each cycle by stepping through the cycle and binning its points in the same manner as the overall in-use data. The same average statistics were determined from the cycles as from the data. These values included average speed, average load, average non-idle speed, average non-idle load, and normalized idle time. We compared the distribution and average statistics with those of the data and calculated values that defined the representativeness of the individual cycles. These values included the “statistical value” (SV), the “distribution value” (DV), and the “performance value” (PV). Absolute differences between the desired data set and the created cycles defined these values. The DV was defined by taking the absolute difference between the cycle and data of the time distribution of each bin. Each absolute difference was multiplied by 0.05 (equal weighting factor for each of the 20 bins, summing to unity), and the sum of those weighted absolute differences determined the value. Equation 2 gives the definition of DV.

(2)

The SV was similar except that it used the absolute differences of the statistics. Each of the five statistics previously defined were given an equal weighting of 0.2 (equal weighting factor for the five statistical categories, summing to unity), and the sum of the weighted differences was used to calculate the SV, as shown by eq 3.

(3)

The PV was calculated by combining the DV and SV with equivalent weights of 0.5 and summing them, as shown by eq 4.

(4)

We used the above approach for drilling and fracturing engines, with the only major difference being idle time at the beginning and end of each cycle. Based on in-use statistics, we determined that only 10 sec of low load were required at the beginning and end of the cycle.

The use of Markov chain logic should create cycles, which have a similar distribution to that of the overall data. Sometimes, however, a cycle can be stuck in a certain bin or “space,” or a cycle may move through the space while choosing low probability events—a side effect of the random nature of the Markov chain. MCMC simulation is the idea of creating a Markov chain that has an equivalent distribution to the one desired. This can theoretically be accomplished with the creation of one long chain or several smaller chains. Ideally, a cycle constructed under the logic of a Markov chain would exactly match the desired distribution if the cycle grew to infinite length. The other option available to best represent the desired distribution would be to create a “large” number of cycles until one matches the desired distribution. The first option could not be applied, as cycles must have a reasonable length. The second option raised a number of questions. Since an infinite number of cycles could not be created, the main issue was answering the question: “How large must the sample size be for a cycle of optimal representation to be found?” It was not practical to create hundreds of thousands of cycles to find an ideal candidate. In an initial attempt to find the optimal cycle for each of the three prime movers, we used Monte Carlo simulation to create populations of 10,000 cycles for each prime mover. The best PV from these populations of OTR truck, drilling, and hydraulic fracturing cycles were 0.0327, 0.0038, and 0.0323, respectively.

Optimization with a genetic algorithm

Review

The concept of the genetic algorithm (GA) stems from Darwin’s theory of evolution, developed in 1859. The governing principle of the theory is survival of the fittest. The key aspect that separates GAs from algorithms that are more traditional is the use of the population rather than analyzing one individual at a time (Sivanandam and Deepa, 2008). GAs are advantageous in solving complex problems for several reasons. These advantages include their ability to discover global optimum solutions, their ability to handle large and poorly understood search spaces, their performance with large-scale optimization, and their easy application to a large variety of problems (Sivandam and Deepa, 2008). These advantages make GAs a choice for difficult problems such as cycle development. Several works from West Virginia University (WVU) have used GAs for problems involving cycle development (Marlow, 2009; Perhinschi et al., 2011).

Design

The first issue encountered was determining how to define the chromosomes and genes of the GA population. The overall cycles constructed by the Markov chain technique defined the chromosomes or individuals. The individuals were composed of 1800-point arrays of speed and load pairs. These individuals all had SVs, DVs, and PVs associated with them as well as distribution matrices and average cycle statistics. We created genes from these individuals. Typically, when using GAs, individuals are converted to binary strings and bits in the strings represent the genes. It was impractical to convert the entire array to binary representation, and so another methodology defined genes. Each cycle was “segmented” into portions of different activity. We defined segments in a similar manner to “micro-trips” from previous studies (Dai et al., 2008; Clark et al., 2007; Perhinschi et al., 2011). Instead of using vehicle speed to define micro-trips, we selected engine speed. Periods of idle and non-idle represented different segments. An idle segment was any length of the cycle longer than 10 sec in which the speed remained at 0% of rated. A non-idle segment was any length of the cycle in which there were more than 10 consecutive seconds of non-idle operation. One hundred random cycles from the 10,000 cycles generated by the MCMC defined the initial population.

With genes defined, we then selected and implemented genetic operators. The GA applied mutation first. The cycles selected for mutation were randomly chosen based on the mutation rate. The GA applied a random number for each cycle in the population. If the number generated was less than the rate of mutation, then the cycle in question experienced mutation. Mutation was the replacement of a segment of the selected cycle by a segment randomly chosen from the pool of segments.

Once mutation occurred, the GA selected cycles for crossover. Determining cycles to experience crossover relied on random number generation for each of the cycles in the initial population of the current generation. Cycles added to the population due to mutation were not eligible for crossover, but the cycles from which they emanated were. If the random number applied to a given cycle was less than the defined crossover rate, then the cycle was selected for crossover. In order to perform crossover, we required another cycle from the population. The GA randomly selected this cycle from the remaining 99 cycles of the current generation.

The new population size was determined by adding the number of newly created cycles from mutation and crossover to the original population size of the generation, which we fixed at 100. We then calculated cycle statistics of the newly created cycles, since newly formed cycles did not contain any information other than the speed load array.

The next step in the GA was to sort the cycles based on their respective PVs in order to determine which individuals were most representative. The PV was desired to be as close to 0 as possible; however, it was not well scaled, between 0 and 1, over the entire population. Many of the cycles within a given population had PVs closer to 0 than 1. In order to more evenly distribute the population so that selection of the next population could be implemented, “evaluation” and “fitness” functions were defined. These functions were necessary to employ roulette-wheel selection to define the next generation. In order to map the PV over the range of [0, 1] to create selection probabilities, the largest and smallest PVs of the respective generation were defined as PV_max and PV_min, respectively. These values were then used to create coefficients a and b, and scale the PV into an “evaluation function” (EF), as shown in eqs 6, 7, and 8.

(6)

(7)

(8)

The use of these equations defined an evaluation function for every cycle. Since a lower PV was desirable, it was more valuable to define the opposite of the EF desired. Equation 9 shows the definition of the inverse EF.

(9)

A total fitness (TF) value (see eq 10) and a probability array were then constructed based on each individual’s evaluation function, as is shown in eq 11.

(10)

(11)

A “q” array was defined that acted as a number line with spacing equivalent to the probabilities defined from eq 12.

(12)

Equations 6–12 were all defined in Roeva et al. (2013). With the “q” array defined, roulette-wheel selection could take place to define the next population. In addition to roulette-wheel selection, the GA used elitist selection to ensure the cycle with the lowest PV advanced to the next generation. Elitist selection improves exploitation and reduces variability between generations (Roeva et al., 2013). Roulette-wheel selection defined the remaining population. To implement the selection process, we generated a random number between 0 and 1, and placed it on the “q” number line as described above. The index of the “q” value falling directly to the right of the random number was the index of the cycle that advanced. Note that we performed this without removal of the cycle from the pool. This meant that it was possible for the same cycle to have multiple copies in the next generation and that other cycles not selected became extinct. This process continued to define the entire new population of the next generation.

One issue in defining the GA was finding genetic operation rates that most efficiently found a solution. In other words, it was necessary to select the best form of the GA to use in cycle development. We applied several iterations of the GA in an attempt to find the best performing genetic operation rates. During these iterations, the GA ran for a life span of 30 generations with the fixed population of 100. The genetic operation rates varied during these iterations, and each of the iterations was ranked based on the value of the elite PV that was produced. The final GA used the operators that produced the lowest average PV amongst the three cycles. Table 7 shows the results of this analysis.

Table 7. Genetic algorithm iteration results.

GA	Rates of operators		Drilling		Fracking		Trucking		Total
Rank	Crossover rate	Mutation rate	Elite PV	Rank	Elite PV	Rank	Elite PV	Rank	PV
1	0.4	0.2	0.0025	1	0.0314	5	0.0301	1	0.0640
2	0.4	0.25	0.0031	2	0.0302	1	0.0336	5	0.0669
3	0.35	0.2	0.0031	2	0.0324	7	0.0322	2	0.0677
4	0.4	0.15	0.0032	3	0.0308	3	0.0340	6	0.0680
5	0.45	0.2	0.0032	3	0.0321	6	0.0335	4	0.0688
6	0.2	0.1	0.0044	4	0.0312	4	0.0334	3	0.0690
7	0.3	0.2	0.0044	4	0.0329	8	0.0322	2	0.0695

Based on the analysis, the final genetic operation rates of 0.4 or 40% for crossover and 0.2 or 20% for mutation yielded the best results. For each of the three types of cycles, an initial population of 100 cycles was applied to the final GA. Each final GA had a life span of 50 generations. We extended the life span of the GA with the goal of giving better results than those seen in the initial iterations, which were shorter to save computational time.

Results

OTR trucks

The initial population of 100 truck cycles had an average PV of 0.0910, and after 50 generations the average PV reduced by 61% to 0.0354. The elite PV of the initial population was 0.0388 and reduced by 18% to 0.0320 by the final generation. Figure 4 shows the improvements of the average and elite PVs as a function of the GA generation. At the 50th generation, the average PV improved by only 0.35% from the previous generation and the elite PV did not improve. Figure 5 shows the elite OTR truck cycle. Table 8 presents the statistics of the elite OTR truck cycles as compared with the in-use data.

Table 8. OTR truck cycle and data statistics.

Metric		Value
		Cycle	Data
Average	Speed	0.2911	0.2923
	Load	0.2269	0.2312
	Non-idle speed	0.4665	0.4685
	Non-idle load	0.2649	0.2826
Time	Idle	0.3761	0.3762
	Non-idle	0.6239	0.6238

Figure 4. OTR truck GA improvements.

Figure 5. Elite OTR truck cycle.

Drilling

The initial population of 100 created drilling cycles had an average PV of 0.0249. After 50 generations, the average PV reduced by 82% to 0.0044. The elite PV of the initial population was 0.0075 and reduced by 65% to 0.0026 by the end of the 50th generation. The average PV at the 50th generation actually was worse than it was in the 49th generation, and the elite PV showed no improvement. Figure 6 shows the elite drilling cycle, and Table 9 presents the statistics of the elite drilling cycle as compared with in-use data.

Table 9. Drilling cycle and data statistics.

Metric		Value
		Cycle	Data
Average	Load	0.4295	0.4298
	High load	0.5997	0.6024
	Low load	0.1934	0.1895
Time	High	0.5850	0.5837
	Low	0.4067	0.4123
	Idle	0.0083	0.0040

Figure 6. Elite drilling cycle.

Hydraulic fracturing

The initial population of 100 hydraulic fracturing cycles had an average PV of 0.0762. After 50 generations, the average PV reduced by 54% to 0.0347. The elite PV of the initial population was 0.0406 and reduced by 25% to 0.0305 by the end of the 50th generation. The 50th generation elite PV showed no improvement over the 49th generation, and the average PV only improved by 0.55%. Figure 7 shows the elite hydraulic fracturing cycle, and Table 10 presents the statistics of the elite cycle and the in-use data statistics.

Table 10. Fracturing cycle and data statistics.

Metric		Value
		Cycle	Data
Average	Load	0.3583	0.3554
	High load	0.6174	0.6149
	Low load	0.2167	0.2033
Time	High	0.4606	0.4681
	Low	0.3411	0.3324
	Idle	0.1983	0.1994

Figure 7. Elite fracturing cycle.

Implementing the dynamometer cycles

We exercised a test engine over the elite cycles to ensure their applicability to an engine dynamometer configuration. The engine used for the purpose of cycle verification was a 2008 Cummins 8.9-L ISL-G, rated at 280 hp. Table 11 provides engine specifications. A GE 800-hp dynamometer controlled engine speed, and load was measured via side-arm torque. Note that eq 1 was directly applicable to the ISL-G, which operates over the Otto cycle. Typically, denormalization for diesel engines includes values up to 112% of rated speed, but data beyond rated speed only represented 0.29% of over 4 million data points and were excluded from this analysis. In the case of diesel-fueled engines for drilling and fracturing, operation only occurred at rated speed.

Table 11. Cummins ISL-G engine specifications.

Specification	Value	Unit
Horsepower	280	hp
Peak torque	900	lb-ft
Governed speed	2200	rpm
Ignition	Spark
Arrangement	Inline—6 cylinder
Intake	Turbocharged
Displacement	540	in^3
Bore	4.49	in
Stroke	5.69	in
Oil capacity	7.3	US gal
Coolant capacity	13.1	US qts
System voltage	12	V
Aftertreatment	TWC
Fuel used	CNG

Note. CNG = compressed natural gas; TWC = three-way catalyst.

Regression and smoothing

Before use for research or measurement of emissions, a cycle regression analysis must occur. Regression statistics must be within a certain range according to the CFR for validation of the cycle. Table 12 provides the validation criteria as defined by CFR Title 40 Chapter I Subchapter U Part 1065 Subpart F Section 514.

Table 12. Cycle regression validation criteria for Cummins ISL-G.

Cycle validation criterion		Acceptable range
		Minimum	Maximum
Speed	Intercept	−70.1	70.1
	Slope	0.95	1.03
	R^2	0.97	1.00
	SEE	0.00	110.35
Torque	Intercept	−18.0	18.0
	Slope	0.83	1.03
	R^2	0.85	1.00
	SEE	0.00	90
Power	Intercept	−5.70	5.70
	Slope	0.83	1.03
	R^2	0.91	1.00
	SEE	0.00	28.50
Work	Actual bhp-hr	Within 5%

Note. SEE = standard error of the estimate.

Drilling and fracturing cycles

Tables 13 and 14 provide the passing results for the drilling and hydraulic fracturing cycles, respectively.

Table 13. Drilling cycle regression results.

Criterion		Min	Value	Max
Speed	Intercept	−70.1≤	−24.187≤	70.1
	Slope	0.95≤	1.011≤	1.03
	R^2	0.97≤	0.998≤	1
	SEE	0≤	5.994≤	110.35
Torque	Intercept	−18.0≤	2.786≤	18.0
	Slope	0.83≤	0.988≤	1.03
	R^2	0.85≤	0.941≤	1
	SEE	0≤	35.969≤	90
Power	Intercept	−5.7≤	1.584≤	5.7
	Slope	0.83≤	0.986≤	1.03
	R^2	0.91≤	0.948≤	1.03
	SEE	0≤	14.198≤	28.5
Work	Actual bhp-hr	57.85≤	60.816≤	63.94
	Reference		60.897

Note. SEE = standard error of the estimate.

Application of smoothing for the OTR truck cycle

The elite OTR truck cycle did not meet regression criteria in original form. The failure of regression for this cycle was most likely the result of the second-by-second concatenation technique. This technique allowed activity points to sometimes jump by intervals of as much as 10%. This did not allow the engine to respond to the dynamometers signals in time to meet the desired set points. We decided to implement a smoothing technique to reduce the rate of change between certain points of both speed and load. Several smoothing and filtering techniques were considered for this application, including moving average, Savitsky-Golay (SG), and local regression methods (loess and lowess). These techniques smoothed the speed and load array with the defined function and parameters and calculated the new statistics of the cycle. An analysis determined how the different smoothing techniques affected the cycle with different levels of span and degree. We performed a sensitivity analysis of the various techniques to examine their effects on the cycle PV. The goal was to apply a technique, which smoothed the cycle with minimal effect on the PV. Based on the analysis, the SG13-7 was determined to be the best method. The SG13-7 smoothing technique had an effect on 1799 of 1800 cycle points. The smoothing technique changed the affected speed points by an average of 4% and the load points by an average of 15%. Table 15 provides the statistics of the smoothed elite OTR truck cycle and in-use data. Table 16 presents the passing results of the regression test, and Figure 8 displays the final normalized OTR truck cycle.

Table 14. Hydraulic fracturing cycle regression results.

Criterion		Min	Value	Max
Speed	Intercept	−70.1≤	−22.419≤	70.1
	Slope	0.95≤	1.010≤	1.03
	R^2	0.97≤	0.998≤	1
	SEE	0≤	6.113≤	110.35
Torque	Intercept	−18.0≤	0.335≤	18.5
	Slope	0.83≤	1.000≤	1.03
	R^2	0.85≤	0.962≤	1
	SEE	0≤	34.967≤	90
Power	Intercept	−5.7≤	0.346≤	5.7
	Slope	0.83≤	0.999≤	1.03
	R^2	0.91≤	0.966≤	1.03
	SEE	0≤	13.738≤	28.5
Work	Actual bhp-hr	48.26≤	50.677≤	53.34
	Reference		50.803

Note. SEE = standard error of the estimate.

Table 15. Smooth OTR truck cycle and data statistics.

Metric		Value
		Cycle	Data
Average	Speed	0.291	0.292
	Load	0.227	0.231
	Non-idle speed	0.463	0.469
	Non-idle load	0.264	0.283
Time	Idle	0.372	0.376
	Non-idle	0.628	0.624
Values	Statistical	0.039
	Distribution	0.026
	Performance	0.032

Table 16. Smooth OTR truck cycle regression results.

Criterion		Min	Value	Max
Speed	Intercept	−70.1≤	0.666≤	70.1
	Slope	0.97≤	0.999≤	1.03
	R^2	0.97≤	0.999≤	1.00
	SEE	0.0≤	10.083≤	110.35
Torque	Intercept	−18.0≤	0.286≤	18.0
	Slope	0.83≤	0.973≤	1.03
	R^2	0.85≤	0.919≤	1.00
	SEE	0.0≤	54.134≤	90.0
Power	Intercept	−5.7≤	−0.052≤	5.7
	Slope	0.83≤	0.983≤	1.03
	R^2	0.91≤	0.934≤	1.00
	SEE	0.0≤	15.462≤	22.1
Work	Actual bhp-hr	20.819≤	22.186≤	23.011
	Reference		21.915

Note. SEE = standard error of the estimate.

Figure 8. Normalized OTR truck cycle with SG13-7 smoothing.

Comparison of the generated cycles with current cycles

OTR and FTP

As the 8.9-L ISL-G is a heavy-duty engine used in OTR trucks, we compared both the cycle statistics and emissions results from the OTR truck cycle and the standard FTP. Table 17 presents both the statistics and emissions comparisons with the FTP values as reference. Emissions data were collected using a 40 CFR 1065–compliant dilution tunnel and constant-volume sampling system (“Code of Federal Regulations”; Wu et al., 2009). Carbon monoxide (CO), oxide of nitrogen (NO_x), carbon dioxide (CO₂), and total hydrocarbon (THC) emissions were measured with a Horiba MEXA 7200 D analytical bench. Brake-specific emissions results are only shown for hot-start tests and are the average of three repeated tests. The OTR cycle has decreased time at idle in comparison with the FTP and therefore produces slightly higher CO₂ emissions (2.9% higher than FTP), which correlates to a slight increase in brake-specific fuel consumption. Of greater concern are the increases in CO, NO_x, and THC emissions over those produced from the FTP. The certification levels for CO and NO_x for this engine were 14.4 and 0.2 g/bhp-hr, respectively. We see that the FTP values are below or near the certification levels. However, when operated on the OTR cycle, CO, THC, and NO_x emissions increased by factors of 1.69, 2.85, and 1.79, respectively. Therefore, inventories that utilize FTP emissions values to represent this particular trucking sector could be significantly underestimating these regulated pollutants. Greater than 90% of the THC emissions from this engine was methane, a potent greenhouse gas, and its relative increase should be addressed for any greenhouse gas inventories or when reporting CO₂-equivalent emissions.

Table 17. Smooth OTR truck cycle regression results compared with FTP.

Metric		Value
		OTR	FTP	% Difference
Average	Speed	0.291	0.415	−42.6%
	Load	0.227	0.241	−6.17%
	Non-idle speed	0.463	0.719	−55.2%
	Non-idle load	0.264	0.397	−50.3%
Time	Idle	0.372	0.423	−13.8%
	Non-idle	0.628	0.577	8.1%
Values	Statistical	0.039	0.173	−342%
	Distribution	0.026	0.072	−183%
	Performance	0.032	0.123	−279%
Emissions (g/bhp-hr)	CO	3.19	1.89	40.8%
	CO2	460	447	2.83%
	HC	2.26	0.910	59.7%
	NOx	0.430	0.240	44.2%

NRSC and drilling and fracturing cycles

Both drilling and hydraulic fracturing engines are subject to certification when operated over the NRSC. The NSRC cycle consists of five modes, with torque set points of 100, 75, 50, 25, and 10%. The weighting factors for each mode are 0.05, 0.25, 0.3, 0.3, and 0.1, respectively. The average time-weighted load of the NSRC is 47.25%. The average load of the drilling cycle was 42.95%, and the average load of the hydraulic fracturing cycle was 35.83%. We set a high-load and a low-load threshold of 40% based on operating statistics from in-field data. The NSRC yields 60% time of high-load operation and 40% low-load operation. The time of high-load operation for drilling was 58.5%, whereas low-load was 41.5%. The time of high-load operation for hydraulic fracturing was 46.1%, whereas low-load was 53.9%. The in-use engines are not subjected to transient operation but do experience transient operation in the field and within the newly developed cycles. The transient operation of the drilling cycle reflects activities such as drilling, pipe connections, operation of draw works, circulation of fluids, and other intermittent activities. The transient operation of the hydraulic fracturing cycle reflects multiple stages of stimulation, which may occur at different loads due to varied pressures and gear selection of the pumps. Transient operation of turbocharged diesel engines can produce higher combustion temperatures, which lead to increased NO_x and particulate matter emissions (Rakopoulos and Giakoumis, 2006). Off-road spark ignited engines have also shown to increase THC plus NO_x emissions by a factor of 6 and CO emissions by a factor 5 when comparing real in-use transient operation with those emissions produced from the C2 NRSC (U.S. Environmental Protection Agency [EPA], 2001). Future work will examine NSRC emissions from the 8.9-L ISL-G engine as compared with those from this work.

Summary

To reduce fuel costs and emissions, the unconventional well completion industry has been attempting to utilize produced natural gas to meet power demands for cost savings. Several companies are doing this by displacing diesel fuel normally used for engines. The primary consumers, or prime movers, of diesel fuel in unconventional natural gas well development are OTR trucks, and high-horsepower drilling and hydraulic fracturing engines. In order to conduct cost- and time-effective research and evaluations, we developed a method to create engine test cycles for application to scaled engines in a laboratory setting.

We collected in-use engine activity data for the three prime movers. Data collected included engine speed and load from OTR trucks traveling to and from well sites hauling water, sand, and gravel, hydraulic fracturing engines used to power pumps for several jobs on two separate pads, and drilling engines used to drill two wells on a single pad. We employed a MCMC technique to create cycles. The goal was to replicate data by creating a large number of potential cycles. A Markov chain utilized second-by-second concatenation and a transition probability matrix. We combined the Markov chain technique with Monte Carlo simulations to create a population of 10,000 cycles for each prime mover. The elite PVs of the populations of OTR truck, drilling, and hydraulic fracturing cycles were 0.0327, 0.0042, and 0.0072, respectively. We then used a GA to optimize the cycle selection process. The GA employed elitist and roulette-wheel selection to determine the members of the population allowed to advance to the next generation. The final GA used genetic operation rates of 40% for crossover and 20% for mutation. It utilized an initial population of 100 cycles and had a life span of 50 generations.

After completion of the GA, the elite OTR truck cycle, drilling cycle, and hydraulic fracturing improved PVs by 18%, 62%, and 65%, respectively. The GA reduced PVs from 0.0388 to 0.0320 for OTR trucks, 0.0079 to 0.0030 for drilling, and from 0.0155 to 0.0054 for hydraulic fracturing. We exercised the cycles on an engine and dynamometer to confirm that they met regression criteria for transient cycles defined by the CFR.

The elite drilling and hydraulic fracturing cycles passed the regression criteria needed as initially defined. The OTR cycle, however, failed regression based on the R² of torque and power criteria. We applied a smoothing technique to the normalized cycle speed and load arrays. A sensitivity analysis determined the effect of smoothing techniques on the PV of the cycles. A SG filter with a span of 13 and polynomial degree of 7 was implemented. The smoothing resulted in a 1.25% increase in the PV from 0.0320 to 0.0324. We then exercised the smoothed OTR truck cycle on the same engine and dynamometer platform, and it passed all necessary criteria.

Based on our analyses, we feel these cycles are representative of actual engine operation in the field and will use them for research related to unconventional natural gas well development. We also provide the final cycles in tabular format for evaluation and use by other researchers. Cycle development will continue to be necessary for research as emissions standards become more stringent and different industries come under the scrutiny of regulatory institutions. For these reasons, it is important to look for new ways to apply and optimize the techniques of in-field data representation. Future research should also include catalyst temperature as a cycle metric. Older off-road engines are typically Tier 2 with no aftertreatment; however, newer engines that can be run on dual fuels have diesel oxidation catalysts (DOCs) and new Tier 4 engines may have both diesel particulate filters (DPFs) and selective catalytic reduction (SCR)—the efficiency of both are correlated with exhaust energy.

In addition, we have shown that from an emissions perspective, a dedicated natural gas engine operated over our OTR cycle produced more CO, NO_x, and THC emissions compared with the conventional FTP cycle. Although we did not compare the drilling and hydraulic fracturing cycle emissions with those from the NRSC, we plan to do so in the future. The drilling cycle high-load and low-load time distribution is similar to the NRSC but includes realistic transient operation. The hydraulic fracturing cycle yielded increased low-load operation compared with the NRSC. The transient nature of the in-use engines and our cycles will likely increase emissions over those collected during NRSC certification tests. As such, current test methods and inventory methods, which utilize current emissions factors, may be underreporting regulated emissions compared with those that occur in use.

Acknowledgment

The authors thank numerous industry collaborators for allowing them to collect these valuable data.

Funding

This work was made possible under a grant from the Department of Energy’s National Energy Technology Laboratory, DE-FE0013689, under direction of Mr. William Fincham. Some data were also collected during a collaborative effort under a separate grant, DE-FE0024297.

Additional information

Funding

This work was made possible under a grant from the Department of Energy’s National Energy Technology Laboratory, DE-FE0013689, under direction of Mr. William Fincham. Some data were also collected during a collaborative effort under a separate grant, DE-FE0024297.

Notes on contributors

Derek Johnson

Derek Johnson, Ph.D., is a research assistant professor at West Virginia University (WVU) conducting research focused on methane and other emissions associated with the use of natural gas.

Robert Heltzel

Robert Heltzel, M.S.M.E., is a research assistant at WVU’s Center for Alternative Fuels, Engines, and Emissions.

Andrew Nix

Andrew Nix, Ph.D., is an assistant professor at WVU conducting research on natural gas reciprocating engines and gas turbines.

Rebekah Barrow

Rebekah Barrow is an undergraduate research assistant pursuing bachelor's degrees in mechanical and aerospace engineering.