Therefore, if the costs of steps I and II are combined, it is possible to compare them with the price of cereal straw at origin. The total costs of these two steps ranged between 28.1 and 51 EUR t−1, which was in the range of half of the straw prices in 2019, which had an approximate value of 70 EUR t−1 [27]; however, in 2020, the price of cereal straw was about 42 EUR t−1 [28], which was less than the harvesting costs of small bales of rice straw.

Step III, the transport to destination points ranged from 5 EUR t−1 (for short distances of 25 km) to 26 EUR t−1 (for a 150 km one-way trip), using big trucks. Longer distances were not relevant, as the points of straw consumption were always to the West of Valencia, where we found the cereal growing areas, and it would not be possible to compete with either the price of rain-fed cereal or its transportation costs in years with normal straw yields.

#### **5. Conclusions**

The collection and packing processes of rice straw in the area of L'Albufera de València were found to not be very efficient, in comparison with the processes carried out for rain-fed cereal. This is due to the structural problems of the plots and the soft soil, which increased the associated costs.

Therefore, in years with low cereal straw prices, the process of collecting rice straw would run at a deficit, meaning that it would require some kind of official subsidy if farmers are to undertake the process; this also happens in other occidental rice producing countries, such as California in the USA, where straw-burning limitations have forced farmers to bale and transport straw, and the relevant authorities are subsiding these tasks [5].

However, in the years when straw reaches a high price, the costs of collecting rice straw from L'Albufera de València and transporting it to relatively close areas could be seen as a competitive activity. Furthermore, another advantage of rice straw is that, being an irrigated crop, production is regular and reliable, even in times of drought.

The packing and extraction of rice straw could, therefore, be viewed as a way to reduce the environmental impacts of burning or burying rice straw in the natural park.

**Author Contributions:** Conceptualization, A.T. and B.V.-M.; methodology, A.T. and J.M.G.; investigation, A.T. and J.M.G.; writing—original draft preparation, A.T.; writing—review and editing, A.T. and B.V.-M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by MINISTERIO DE CIENCIA, INNOVACIÓN Y UNIVERSI-DADES. AGENCIA ESTATAL DE INVESTIGACIÓN, grant number RTC-2017-6249-2 titled "Desarrollo de un nuevo insumo para la agricultura sostenible: mulch de paja de arroz con incorporación de bacterias promotoras del crecimiento de las plantas (PGPB) y mecanización integral de procesos (SMART MULCH)".

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Acknowledgments:** We gratefully acknowledge Gestión Integral de Residuos Sólidos S.A. (GIRSA) and, in particular, the technician Pere López; as well as the farmer José Luis Pecino for facilitating the field works, and the technicians José V. Seguí and Cristina Zarzuela from the UPV.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Review* **Current Knowledge and Future Directions for Improving Subsoiling Quality and Reducing Energy Consumption in Conservation Fields**

**Shangyi Lou 1, Jin He 1,\*, Hongwen Li 1, Qingjie Wang 1, Caiyun Lu 1, Wenzheng Liu 2, Peng Liu 1, Zhenguo Zhang 3 and Hui Li 4**

**Abstract:** Subsoiling has been acknowledged worldwide to break compacted hardpan, improve soil permeability and water storage capacity, and promote topsoil deepening and root growth. However, there exist certain factors which limit the wide in-field application of subsoiling machines. Of these factors, the main two are poor subsoiling quality and high energy consumption, especially the undesired tillage depth obtained in the field with cover crops. Based on the analysis of global adoption and benefits of subsoiling technology, and application status of subsoiling machines, this article reviewed the research methods, technical characteristics, and developing trends in five key aspects, including subsoiling shovel design, anti-drag technologies, technologies of tillage depth detection and control, and research on soil mechanical interaction. Combined with the research progress and application requirements of subsoiling machines across the globe, current problems and technical difficulties were analyzed and summarized. Aiming to solve these problems, improve subsoiling quality, and reduce energy consumption, this article proposed future directions for the development of subsoiling machines, including optimizing the soil model in computer simulation, strengthening research on the subsoiling mechanism and comprehensive effect, developing new tillage depth monitoring and control systems, and improving wear-resisting properties of subsoiling shovels.

**Keywords:** subsoiling machine; subsoiling shovel design; anti-drag; tillage depth detection; tillage depth control; tillage depth stability; soil mechanical interaction

#### **1. Introduction**

As the important factor in crop growth, soil is a valuable agricultural resource, a significant production factor, and an indispensable foundation supporting the sustainable development of agriculture [1]. The ideal farmland soil is composed of 50% soil particles with organic matter and 50% pores; additionally, with respect to pores, the moisture and air respectively account for 25% [2]. However, due to long-term conventional ploughing practices, hardpan has been shaped, which is a typical form of soil compaction [3,4]. Consequently, the failure layer is formed, resulting in the rearrangement of soil particles under external forces [5,6], and giving rise to a reduction in soil porosity and an increase in bulk density [7–9]. The root-system penetrability is affected by soil compaction, and especially the high-level compaction inhibits the root growth [10]. The hardpan binds root proliferation, reduces root penetration, and decreases root length and dry matter [11]. In addition, the key enzymes required for plant respiration show a downward trend with

**Citation:** Lou, S.; He, J.; Li, H.; Wang, Q.; Lu, C.; Liu, W.; Liu, P.; Zhang, Z.; Li, H. Current Knowledge and Future Directions for Improving Subsoiling Quality and Reducing Energy Consumption in Conservation Fields. *Agriculture* **2021**, *11*, 575. https:// doi.org/10.3390/agriculture11070575

Academic Editor: José Pérez-Alonso

Received: 21 May 2021 Accepted: 20 June 2021 Published: 23 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

increasing compaction stress, which leads to impaired root function and the loss of nutrient absorption [12–15].

Subsoiling can break the hardpan and alleviate the soil compaction without turning ploughing soil and disrupting the original topsoil structure [16]. Compared with conventional ploughing, the subsoiling tillage at 30 cm, 35 cm, and 40 cm decreased the mean bulk density by 4.59%, 7.13%, and 8.27% and reduced soil compactness by 17.62%, 23.63%, and 36.42%, respectively [17]. Meanwhile, subsoiling increased soil porosity, enhanced soil infiltration and water storage capacity, and decreased water inputs per growing day [18]. Additionally, in comparison with the rotary tillage, the subsoiling practice significantly increased the average root length density by 13% due to the better root penetration in both horizontal and vertical directions. Moreover, the average grain yield and dry matter weight of maize were respectively increased by 6.3% and 3.7% [19]. Furthermore, subsoiling improved the activity of plant-protective enzymes in maize root, reduced the degree of peroxide in cell membrane substances, and delayed the senescence of maize root. It is conducive to maintaining root system vitality in the later growth stage of maize [20]. Therefore, subsoiling is classified as a resource-saving and environmental-friendly technology, and is significant in promoting the sustainable development of agriculture [21–24].

In the early 1930s, some countries in North America and Europe began to use subsoiling technology to solve soil problems. Thanks to the significant effects of subsoiling on the repair and protection of arable land, the combination of tillage methods and comprehensive utilization of subsoiling technology has become increasingly close in these regions. The proportion of arable land that adopted subsoiling in Europe was 15.4% in 2005. Meanwhile, the proportion in North America exceeded 40% in 2006. Until 2015, the proportion in European and American areas had reached 60% [25]. In some Asian countries, subsoiling technology was adopted relatively late. To address the issue of the crop yield decline caused by soil compaction, China has promulgated a series of associated regulations to promote subsoiling and land preparation technologies since 2009 (Table 1). By the end of 2018, the subsoiling was practiced on 10.6 Mha, accounting for 7.41% of the total arable land [26]. At present, subsoiling technology is being widely used in scientific experiments and actual productions of various crops such as wheat, corn, cotton, sugar cane, tobacco, and soya beans.

**Table 1.** Relevant development regulations of subsoiling and land preparation technology in China.

[a] "Government Work Report" is a form of the official document of the Government of the People's Republic of China, which clearly points out the government's tasks for the current year; [b] "Central Document No.1" is the name originally given to the first document issued by the central authority every year, and focuses on the development of agriculture, rural areas, and farmers in 1982.

Subsoiling machines are indispensable to mechanized subsoiling, and their performance directly affects the hardpan's broken quality, topsoil structure, soil permeability, soil microbial quantity, root growth, and crop yield [36]. According to their working function, subsoiling machines can be divided into single-function subsoiling machines and multi-function subsoiling machines. The single subsoiling machines can be further divided into two categories: subsoiling ploughs and omni-directional subsoilers [37]. Some typically available subsoiling machines around the world are discussed in Table 2. Subsoiling machines in certain countries with large per-capita arable land areas (e.g., America, Canada, Australia, Russia) are mainly matched with high-power tractors which have several advantages, such as large tillage width, fast operation speed, high operation efficiency, and advanced machining technology. On the contrary, for those countries with small per-capita arable land areas, such as China and Japan, the per-capita arable land areas are respectively only about 0.08 ha and 0.03 ha, which are merely approximately one-sixth and one-twentieth of that in America, and about 0.038 and 0.014 times as large as Australia, respectively [38]. The subsoiling machines are mainly matched with small or medium-power tractors. The features of subsoiling machines in these countries include small size, light weight, flexible operation, and low cost. Nevertheless, current subsoiling machines have some problems with their performance, such as poor soil loosening quality, undesired soil disturbance, unstable tillage depth, high power consumption, and rapid wear of shovels.

**Table 2.** Comparison of different kinds of subsoiling machines.

In order to improve subsoiling quality and reduce energy consumption in the field with cover crops, several studies have been conducted on subsoiling shovel structure, anti-drag technologies, tillage depth detection and control, and soil-mechanical interaction [45–51]. Various subsoiling shovels are capable of meeting the needs of different regions in different countries. Reducing tillage resistance during subsoiling makes subsoiling machines more prominent with the advantages of resource saving and being environmentally friendly. Modern sensor detection technology and electronic-hydraulic control technology are used to accurately control the tillage depth of subsoiling machines, which benefits in terms of conserving moisture and preparing a good seedbed [19,52–55]. Additionally, soil-mechanical

interaction mechanisms presented in micro and macro ways, mainly clarified via theoretical analysis and simulation, are an important means to realize the efficient development of the abovementioned studies. The integration of abovementioned technologies contributes greatly to the improvement of subsoiling quality and the reduction of energy consumption in the field with cover crops.

The goal of this study was to comprehensively review the existing literature related to technologies for improving subsoiling quality and reducing energy consumption currently being used in subsoiling machines. This article was based on the following aspects and organized as follows. Section 2 states the methods of design and optimization of subsoiling shovels, commenting on the different types of available shovels. Then, Section 3 explores five ways to reduce tillage resistance. Sections 4 and 5 provide a review of tillage depth detection and control, applying sensor detection technology and electro hydraulic control technology. Section 6 delves into the application of theoretical analysis and discrete element analysis in research on soil-mechanical interaction. Finally, in Section 7, recommendations are given for the future development of subsoiling machines.

#### **2. Subsoiling Shovel Design**

Subsoiling shovels are the main working parts of a subsoiling machine [37]. Their design has a significant influence on subsoiling quality and tillage resistance [56,57]. The common types of commercial shovels are chisel, winged, "V" shaped, and side bended shovels. Through optimizing the penetrating angle, the opening angle of winged shovels, and the shovel-handle structure, some subsoiling shovels with better working performance have been developed. The comparisons among these shovels are shown in Table 3.

**Table 3.** Subsoiling shovel type, model, country, working width, working depth, matched power, number, and feature.

In order to improve quality, reduce resistance, and increase efficiency, a cambered subsoiler shank was designed by Chen [65]. The test showed that the subsoiling resistance was the smallest when the sweep angle was 50◦. Based on the orthogonal experiment's results, i.e., the penetrating angle had a significant effect on subsoiling resistance, Liu et al. [47] optimized the chisel shank's structure and found that the resistance was minimal when the rack angle was 21◦. Zhao et al. [66] designed a fitting curve-shaped shank and clarified the effect of the shovel tip on the soil particle's movement. The results of the comparison test between the fold-line-typed shank and the circular shank showed that the fitting curveshaped shank effectively reduced soil disturbance and tillage resistance. By simulation

analysis and digitized soil-bin experiments, Wang et al. [67] demonstrated that, on the circular shank, the optimum installing height of the winged shovel was 75 mm.

The shank type, penetrating angle, opening angle of the winged shovel, and shovelhandle structure are the main factors influencing the tillage resistance of subsoiling shovels. In order to improve subsoiling quality and reduce tillage resistance, these factors are used as the target parameters. Through subsoiling simulation and soil-bin and field trials, the effect of shovel structure on the tillage resistance as well as soil disturbance were analyzed and their impacts were clear. Then, these target parameters were optimized on the basis of the above results, and the shovel structure was accordingly improved to obtain better subsoiling quality and lower tillage resistance. Design and optimization of subsoiling shovels has always been one of the key technologies of subsoiling machines.

**Figure 1.** (**a**) Chisel shovels of the Hedongxiongfeng 1S-150A subsoiler (http://www.hedongxiongfeng. cn/pd.jsp?id=52#\_pp=103\_366 (accessed on 15 September 2020)); (**b**) Winged shovels of the Zhongnongboyuan 1S-200 subsoiler (http://www.boyomac.com/product/17.html (accessed on 14 May 2021)); (**c**) V shaped shovels of the Woye 1SQ-340 subsoiler (https://www.nongjitong.com/product/3072.html (accessed on 15 September 2020)) (**d**) Side-bended shovels of the DaHua BaoLai 1S-250 subsoiler (https://www.dhbl.net/product/53.html (accessed on 15 September 2020)).

#### **3. Anti-Drag**

The energy consumption caused by subsoiling is three to five times as much as that of other planting procedures such as seeding, managing, and harvesting [68]. Decreasing tillage resistance during subsoiling is the main measure to reduce energy consumption. Several methods of decreasing resistance are available, such as vibration anti-drag technology, bionic anti-drag technology, layered subsoiling structure, shovel design with anti-drag structure, and surface coating technology of shovels.

#### *3.1. Vibration Anti-Drag*

Vibration anti-drag refers to adding an excitation source on the subsoiling machine to make subsoiling shovels constantly subject to the repeated positive and negative forces. The greatest advantage of the vibration subsoiling machine is its low traction resistance. According to the driving mode of excitation source, vibration machines can be divided into self-excited and forced-excited categories. Shovels of the forced-excited vibration subsoiler

are capable of vibrating via connecting an eccentric mechanism. Sun et al. [69] adopted two eccentric discs to drive connecting rods and shovels reciprocating vibration in the vertical direction. The minimum resistance and energy consumption were obtained under the optimum combination of forward speed, vibration frequency, amplitude, angle, and speed ratio [70–73]. A typical forced-excited vibration subsoiler is shown in Figure 2. Self-excited vibrating can be obtained by applying elastic elements, such as pressure springs (Figure 3a), leaf springs (Figure 3b), and hydraulics [74–76] (Figure 3c). This anti-drag method not only has little impact on soil, but also protects shovels and other components [77].

**Figure 2.** DaHua BaoLai 1SZ-180 forced-excited vibration subsoiler (http://www.nongjitong.com/ product/dahua\_1sz-270\_deep\_loosening\_machine.html (accessed on 12 November 2020)).

**a b c**

**Figure 3.** (**a**) Great Plains SS1300 subsoiler (http://www.greatplainsmfg.com.ua/en-gb/products/710/sub-soiler-narrowtillage (accessed on 8 December 2020)); (**b**) Kverneland GLG-II subsoiler (https://www.kverneland.cn/node\_80887/node\_ 80897/node\_80958# (accessed on 8 December 2020)); (**c**) AGRI-WELD subsoiler (https://www.plowmanbrothers.com/ plowman-agricultural/ (accessed on 8 December 2020)).

#### *3.2. Bionic Anti-Drag*

The primary task of bionic anti-drag is to extract biological information, which includes the function, structure, process, or behavioral characteristics and mechanisms of the biological system. On this basis, a technology similar to the function of biological-systems was developed and applied to shovels to reduce resistance [78]. The contour line of the claw toe from house mice (Mus musculus) was extracted and found to feature an exponential function curve, and Zhang et al. [79] converted the upper contour line of the longitudinal profile of their (Mus musculus) claws into the technical parameters of the shovel's cutting edge curve. Then, the contour curve was enlarged properly and applied to the structural design of the cutting edge of the shovel handle. The comparative field experiment between the shovel with the feature of an exponential curve and the conventional shovel showed that the tillage resistance of the bionic shovel was reduced by 8.5–39.5% compared with

the traditional shovel. Li et al. [80] obtained the accurate shape and dimension of a bear claw by a three-dimensional laser scanner and established a claw model. By the discrete element simulation, they found that the overall performance of the claw was the best when the rack angle was 30◦. A bionic vibratory subsoiler was designed by Zhang et al. [81], which combined vibration anti-drag with bionic anti-drag. It was reported that the draft force of the bionic vibratory subsoiler was respectively reduced by 13–18% and 8.5–39.5% compared with the 1SZ-460 lever-type subsoiler and bionic anti-drag subsoiler.

In addition to the above two methods, the layered subsoiling structure is capable of reducing tillage resistance and mainly depends on two shovels.  The shovel design with an anti-drag structure optimized the handle or tip via experiment or simulation, thereby obtaining better subsoiling quality and less tillage resistance [84,85]. The surface-coating technology aims at obtaining better soil-removal performance of shovels by changing the properties of surface material, hence decreasing tillage resistance [86].

The above mentioned anti-drag technologies can effectively reduce resistance, whereas in studies of reducing shovel resistance, less consideration is given to the thickness of the three soil layers (tillage layer, hardpan, subsoil layer) and the soil's physical properties, as well as the soil smashing process. Apart from reducing tillage resistance, the power consumption of subsoilers and straw blocking of shovels should also be reduced. Although the forced-excited vibration subsoiler is capable of decreasing the traction, the required fuel consumption is much higher. In order to balance the traction and energy requirements and to improve the overall efficiency, better methods should be explored through further studies so as to drive the vibrator.

#### **4. Tillage Depth Detection**

Initially, the tillage depth could only be measured manually, which had great labor intensity and low detection accuracy. Additionally, these measured depth data were limited and discrete, and could therefore not reflect continuous changes including the variation of tillage depth with surface and time; consequently, it is unable to study the dynamic characteristics of depth. To this purpose, sensors, which have higher detection accuracy and provide a basis for depth adjustment of subsoiling machines, are used for measuring tillage depth. At present, there are two main kinds of sensors used in tillage depth detection: ultrasonic sensors and inclination sensors.

#### *4.1. Tillage Depth Detection Based on Ultrasonic Sensor*

The detection method of ultrasonic sensors utilizes a constant speed transmission of acoustic waves. The signal generator generates acoustic waves with a certain oscillation frequency and which travels through the air at a constant speed. The ultrasonic echo is generated when the acoustic wave encounters impurities and interfaces, and then it is received by the signal-receiving end. Therefore, the distance is obtained by calculating the formula of time interval and sound speed [87].

Tests were carried out to evaluate detection accuracy of the ultrasonic sensor on different fields by Mouazen et al. [88], and the ultrasonic sensor was installed on the bottom of the frame. Compared with the manual measurement values, it was found that the detection accuracy of the ultrasonic sensor was better in the soft sandy loam, while it was lower in the soil with cover crops. A tillage depth measurement device based on the ultrasonic sensor was designed by Li [89], and proved to behave better than the resistancestain type measurement device in terms of measuring accuracy. Adamchuk et al. [90] developed a closed-loop automatic control system for tillage depth, and the ultrasonic sensor was applied to measure tillage depth. An online measuring device of subsoiling depth was proposed in the Chinese Academy of Agriculture Mechanization Sciences [91]. It used an ultrasonic sensor to detect the distance between the frame of the subsoiling machine and the surface, and thus the tillage depth was obtained. A tillage depth detection

system was designed by Suomi et al. [50], which used an ultrasonic sensor installed in front of the frame with an inclination sensor mounted on the connection arm of the gauge wheel and the frame to measure the height and the angle; the depth was obtained by calculating the formula related to these two physical properties. The field experiment showed that the detection error of this system was within 10 mm.

#### *4.2. Tillage Depth Detection Based on Inclination Sensor*

Inclination sensors determine the position of the object according to the three-dimensional angle change, which should be attached to the objects' surface to measure the inclination of the object from the horizontal plane [92]. When the inclination sensor is used to detect tillage depth, it is often on the hitch device or profiling mechanism to measure the angle change. Then, the tillage depth is obtained by formula calculation.

A detector, consisting of a frame, an inclination sensor, and a sliding plate which undulated with the surface, was designed by Zhao et al. [93] and used to measure tillage depth. By detecting the relative angle between the frame and the tractor, the depth was acquired through geometric modelling of the depth and angle. Li et al. [94] developed a measurement device, and the inclination sensor was on the linkage of the frame and gauge wheel. By angle measurement and calculation, the tillage depth was obtained. Results showed that the error of this measurement method was within 6%. Xie et al. [95] installed inclination sensors on the lift arm of a tractor to detect angle change from the horizontal plane. The actual tillage depth was calculated according to the inclination of the lift arm and the geometrical dimension relation of the linkage mechanism. In addition to the abovementioned detection methods based on ultrasonic sensors and inclination sensors, potentiometers and encoders, of which the detection principle is similar to that of the inclination sensor, are also used to detect tillage depth [96,97].

In existing research on the technology of tillage depth detection, ultrasonic sensors and inclination sensors are widely used because of their low cost, simple working principles, and convenient installation methods. However, there still exist some problems. The detection accuracy of ultrasonic sensors would be affected by weeds, straws, and clods, which could lead to unstable detection results. In addition, inclination sensors can be used to measure the horizontal dip angle of the rear-suspension lifting arm. Although this method avoids the negative influence of surface unevenness and other abovementioned factors on detection accuracy, other problems may appear. As there are many connecting rods between the tractor boom and the suspension unit, drivers may adjust the length of the lift rod and the upper pull rod in practice, which leads to a change in the geometric parameter. Thus, it is inconvenient to use this method due to the need to recalibrate the sensor. Besides, the method, applying inclination sensors to measure the angle change of the profiling mechanism and calculate the mathematical model to get the tillage depth, results in the mismatch between the adjustment completion point and detection point of tillage depth. This is because the profiling mechanism is generally placed behind the tillage component, which leads to the delay of profiling.

#### **5. Tillage Depth Control**

Because of the complexity and variability of the field-work environment, the tillage depth often needs to be adjusted many times to achieve stable and consistent outcomes. It is essential to develop the control technology of tillage depth to obtain precise and efficient adjustments. Currently, there are two major control modes of tillage depth: the method of unified adjustment of each row and the method of independent adjustment of a single row. The adjusting mechanism belonging to the first adjustment mode is based mainly on the attachment position of the tillage machine and its tractor. This method commonly combines the hydraulic device of three-point suspension with mechanical and electrical technology to achieve automatic control of tillage depth [98]. With regard to the method of independent adjustment of a single row, separate adjustment mechanisms, hydraulic mechanisms, or electrical mechanisms need to be designed and installed on the machinery, and the number of mechanisms is equal to the number of subsoiling shovels. This method is capable of avoiding undesired tillage depth of each subsoiling shovel.

#### *5.1. Method of Unified Adjustment of Each Row*

The method of unified adjustment of each row aims at adjusting the tillage depth of the whole machine (adjusted on the three-point suspension) or a group of tillage components of the combined land-preparation equipment. Therefore, the hydraulic suspension system with electronic control technology (EHC) is an important support for this mode. The method of controlling depth using EHC mainly involved three control methods, including force control, position control, and integrated control of the two abovementioned modes. Long-term research on the automatic control of tillage depth based on these three methods has been conducted. In order to improve control accuracy, Lee et al. [48] developed a tillage depth control system. Some sensors were used to detect inclination of the three-point suspension, pitch angle of the tractor, and distance from the transducer to the surface. The drive circuit of the electromagnetic valve was controlled by a controller to change the state of the hydraulic circuit, and thus the position of the three-point suspension was adjusted. For the electro-hydraulic suspension system composed of a cartridge valve, a suspension control scheme based on the Controller Area Network (CAN) bus was proposed by Xie et al. [99] to improve the control accuracy, which mainly included the suspension subsystem ECU and two intelligent nodes. The experiment demonstrated that the developed nodes were capable of meeting the working requirements of the suspension system. Nie et al. [100] developed a tillage depth automatic control system, which was based on the original hydraulic suspension system and applied sensors, microcomputers, and stepping motors to realize automatic control. Moreover, the reliability and stability of the system were verified by experiments.

In addition, combined tillage equipment with the function of monitoring and controlling depth was researched and developed by some agricultural machinery enterprises such as OPICO, John Deere, and CASE. The HE-VA combined land preparation machine consists of groups of subsoiling components and disc harrows, of which the tillage depth is controlled by adjusting the height of the suppress roller through the hydraulic system [101]. John Deere developed a TruSet monitoring system for the land preparation equipment. Depth detection sensors, which are mounted on the supporting wheels and frame, are capable of accurately measuring the depth of the whole machine. Additionally, the depth of each group of tillage components is adjusted via the matched AccuDpthTM hydraulic system [102]. The AFS (Advanced Farming System) tillage depth intelligent monitoring system was developed by CASE, and is capable of precisely adjusting the tillage depth of the whole machine according to different soil conditions [103]. Before operation, the prescription map containing the predetermined tillage depth information of the field is input into the system. Owing to the tractor being equipped with a GPS positioning system, when the equipment moves to the corresponding position, the tillage depth is automatically adjusted by the hydraulic system to the required tillage depth. A stroke-detection sensor is installed in the hydraulic cylinder to provide real-time access to length changes of the hydraulic cylinder, and this is the key to the precise control of tillage depth.

#### *5.2. Method of Independent Adjustment of Single Row*

The abovementioned technologies in Section 5.1 are capable of adjusting the tillage depth of the whole machine or a group of tillage components of the combined land preparation equipment. However, due to surface relief or soil resistance variation, there exists the problem of lateral instability of tillage depth of the machine with several subsoiling shovels, which leads to the inconsistency of soil conditions in the lateral distribution. Therefore, it is necessary to develop subsoilers with special subsoiling assemblies that are capable of independently adjusting tillage depth of each subsoiling shovel. The method of independent adjustment of a single row is capable of avoiding undesired tillage depth of each subsoiling shovel to improve the tillage depth stability between rows.

A universal-type device which automatically monitors and controls the tillage depth of each subsoiling shovel was developed [104]. Laser-ranging sensors were used to measure the distance between the stand and surface. Then, the corresponding signal was sent to the Programmable Logic Controller (PLC) to be analyzed and calculated, and the real-time tillage depth was recorded. Comparing the real depth with the set value, the PLC sent an execution signal to the driver circuit of the electromagnetic valve to change the solenoid valve switch, and the hydraulic cylinder was adjusted to extend or shorten. Therefore, the tillage depth, in real time, was controlled to be consistent. Wang et al. [105] designed an electric-hydraulic control system which was capable of adjusting the tillage depth of subsoiling assemblies in a timely manner (Figure 4a). Two inclination sensors were used to detect the dip angle in the horizontal direction, and were mounted horizontally on the subsoiler frame and inclined frame of the suppress roller, respectively. Through internal algorithm, the tillage depth was calculated. The hydraulic cylinder performed the action, elongating or shortening to obtain desired tillage depth. Moreover, the performance of the system was evaluated in the aspect of the tillage depth stability between rows and within rows. Wu [106] developed a subsoiling machine to automatically monitor and control tillage depth based on ultrasonic sensors and hydraulic drivers. It was capable of independently monitoring and controlling the depth of each subsoiling shovel.

**Figure 4.** (**a**) Subsoiling assembly with the electric-hydraulic control system (Wang et al., 2018); (**b**) Subsoiler equipped with four subsoiling assemblies with the electric-hydraulic control system (Wang, 2019).

> In conclusion, the control method of unified adjustment of each row adopts integrated adjustment and position adjustment combined with force adjustment. This mode not only ensures the tillage depth stability to a certain extent, but also takes the influence of tillage resistance into account. The control precision of tillage depth by this mode meets the operation requirements. However, relevant parameters in the proportional fuzzy control are difficult to determine, and the control process is more complicated. This control mode is suitable for the operations in flat fields. Furthermore, the single-row tillage depth cannot be independently adjusted according to undulation of local surfaces. By comparison, for the method of independent adjustment of a single row, the designed mechanism used to adjust depth has strong applicability. By adjusting the connecting parts of the control device, the adjustment mechanism can be installed on the different tillage machines. In addition, the main program of the automatic control system belonging to the separate adjustment mechanism is relatively simple. Moreover, tillage depth of subsoiling assemblies can be

independently adjusted. Therefore, this control mode is excels at obtaining consistent depths of each subsoiling shovel.

#### **6. Soil Mechanical Interaction**

The interaction mechanism analysis of subsoiling shovels and soil is the key to the study of the shovel force and the effect of shovels on the dynamic changes of soil disturbance [107]. Wheeler [108] developed a mechanical model and studied the influence of speed on tillage resistance. Manuwa [109] set up models for the resistance prediction of different shovels and analyzed the effect of loosening depth on draft force. Based on the mechanical principle of two-side wedge and three-side wedge in soil, Ma and Zhou [110,111] studied the straight-legged shank and the circular shank and built mathematical models. In recent years, the Discrete Element Method (DEM) has been widely used to study the tillage mechanism [112,113]. Shmulevich [114] used PFC2D discrete-element software to research the interaction between the wide shovel and soil, and four types of shanks were chosen for simulation. Results showed that soil accumulation in front of shovels increased horizontal resistance during movement, and the soil flow below the tip influenced the vertical force applied to the shovel. A soil-tool-residue interaction model was developed by Zeng et al. [115] using the DEM and was used to investigate and compare the effects of forward speed on soil displacement, residue displacement, and residue cover reduction resulting from four types of shovels (Figure 5). Considering the normal tension generated among elements, Momozu et al. [116] used coefficients to modify the traditional discrete element model, and then to represent the adhesion generated by moisture among soil particles. Tamas et al. [117] developed a discrete element model to analyze the effects of the sweep rack angle and speed on soil loosening and porosity. Results showed that the comparison between the measured and simulated draught of subsoiling shovels with a 30◦ sweep rack angle behaved as a good match. In the selected speed interval of 0.5–2.4 m s−1, the error range was 4–15%. A working model of subsoiling was established by Hang et al. [118], applying the DEM. Through simulation based on this model and an indoor soil-bin test, they compared and analyzed the micro-movement and macro-disturbance behavior of subsoil at different locations under the combined effect of two shovels (Figure 6).

**Figure 5.** The soil-tool-residue interaction model (Zeng et al., 2020).

**Figure 6.** Analysis of the soil disturbance process (Hang et al., 2017): (**a**) Schematic diagram of the disturbance process in the longitudinal-sectional; (**b**) Schematic diagram of the disturbance process in the transverse-sectional.

In the subsoiling process, the soil is sheared, extruded, and uplifted by the subsoiling shovel, which provides the reacting force from the soil. Although the soil-mechanical interaction is a complex process, the disruption and deformation of soil under the force of subsoiling components still have some regularity. According to the soil structure, soil mechanics, as well as the failure and deformation regularity, the virtual soil model can be built in the simulation software, which is capable of replacing actual soil within a certain range of error. Furthermore, the simulation test of subsoiling machines and soil models can be carried out to study machines' performance and optimize structural parameters, which is conducive to the shortening of the development cycle of subsoiling machines, improving research efficiency, and saving time and resources. Additionally, the soil smashing mechanism during subsoiling and interaction between the shovel and soil in the process of adjusting tillage depth can be revealed.

#### **7. Summaries and Recommendations**

A large amount of research on subsoiling shovel design, anti-drag, tillage depth detection and control, and soil mechanical interaction have been done to improve subsoiling quality and reduce energy consumption. However, most soil models established in simulation software lacking plant root or straw are inconsistent with the field soil condition. There is insufficient information regarding the broken soil characteristics, stubble-breaking process analysis, and comprehensive effects of subsoiling. Further research should focus on the interaction mechanisms between soil, subsoiling shovels, and plant root or straw. The combination of tillage depth detection and control can effectively improve tillage depth stability, whereas the influence of different soil conditions on the control system is unclear. It is necessary to develop monitoring and control systems of tillage depth with good applicability. Apart from using the abovementioned methods to save energy, the performance of subsoiling components should be considered with more attention. As shovels directly make contact with soil, they wear severely. For this reason, studies on the wear resistance of subsoiling components should be conducted. Therefore, to develop high-performance subsoiling machines with good subsoiling quality and low energy consumption, this article proposes the following four recommendations.

(1) It is very important to improve simulation accuracy by optimizing the soil model in computer simulation. The discrete element model of a root-soil composite model with straw should be established, and the key parameters of soil models need to be calibrated by simulation and laboratory tests. Furthermore, the accuracy and reliability of soil models should be verified by subsoiling simulation and field experiments.

**Author Contributions:** Conceptualization, S.L. and J.H.; writing—original draft preparation, S.L.; writing—review and editing, H.L. (Hongwen Li), Q.W., C.L., W.L., P.L., Z.Z., H.L. (Hui Li); funding acquisition, H.L. (Hongwen Li), J.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the China Agriculture Research System of MOF and MARA (Grant No.CARS-03); Innovative Research Team in University of China (Grant No. IRT13039).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** This research was supported by The 2115 Talent Development Program of China Agricultural University.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Article* **Development of Thermal Performance Metrics for Direct Gas-Fired Circulating Heaters**

**Benjamin C. Smith, Brett C. Ramirez \* and Steven J. Hoff**

Department of Agricultural and Biosystems Engineering, Iowa State University, Ames, IA 50011, USA<EMAIL_ADDRESS>(B.C.S<EMAIL_ADDRESS>(S.J.H.)

**\*** Correspondence<EMAIL_ADDRESS>Tel.: +1-515-294-0468

**Abstract:** Many climate-controlled agricultural buildings use direct gas-fired circulating heaters (DGFCH) for supplement heat. There is no standardized test to calculate thermal efficiency for these heaters. This study aimed to develop a measurement system and analytical analysis for thermal efficiency, quantify the measurement uncertainty, and assess economics of DGFCH efficiency. The measurement system developed was similar to the ASHRAE 103 standard test stand with adaptations to connect the apparatus to the DGFCH. Two different propane measurement systems were used: input ratings < 30 kW used a mass flow system and input ratings > 30 kW used a volumetric gas meter. Three DGFCHs (21.9, 29.3, 73.3 kW) were tested to evaluate the system. Thermal efficiencies ranged from 92.4% to 100.9%. The resulting uncertainty (coverage factor of 2; ~95% Confidence Interval) ranged from 13.1% to 30.7% for input ratings of 56.3 to 11.4 kW. Key sources of uncertainty were propane and mass flow of air measurement. The economic impact of 1% difference in thermal efficiency ranged from USD \$61.3 to \$72.0 per heating season. Refinement of the testing system and procedures are needed to reduce the uncertainty. The application of this system will aid building designers in selection of DGFCHs for various applications.

**Keywords:** heating cost; thermal efficiency; uncertainty

Performance Metrics for Direct Gas-Fired Circulating Heaters. *Agriculture* **2021**, *11*, 588. https:// doi.org/10.3390/agriculture11070588

Academic Editor: José Pérez-Alonso

**Citation:** Smith, B.C.; Ramirez, B.C.; Hoff, S.J. Development of Thermal

Received: 26 April 2021 Accepted: 22 June 2021 Published: 24 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

The use of Direct Gas-Fired Circulating Heaters (DGFCH) is common in climatecontrolled agricultural buildings (e.g., livestock, poultry, and greenhouses) as the heating component of the thermal environment modification system. Supplemental heat is needed to achieve the environmental control goal in an agricultural setting, that is, maintain an ideal environment for the growth of the agriculture commodity with minimal cost [1]. For both animal confinement buildings and greenhouses, heating during the winter months is necessary to maintain the optimal environment and presents a significant utility cost [1,2]. In most building types, a forced air, circulating style heater is used to provide circulation of the heated air to the occupied zone within the building. Gas-fired heaters are ubiquitous due to low operating costs and fuel availability to remote locations [2]. While DGFCH are popular for agricultural applications, there are downsides in implementing them.

DGFCH are simple heaters in comparison to residential units based on the construction and operation principles. DGFCHs use a propane burner manifold placed on the intake of a blower. This does not the vent the combustion gases out of the occupied space. The output of the blower can be equipped with various directional vanes to aid in the distribution of heated air. The placement of such heaters in agricultural buildings is highly variable with no standardized placement design currently available in literature [1,3]. Thermal efficiencies of DGFCHs are commonly reported with assumed values of approximately 99%; however, the water vapor and ventilation rates of a building can decrease the thermal efficiency to 80% [4]. In comparison to residential combustion heaters, which indirectly transfer heat to the air stream via a heat exchanger and combustion gases are vented outdoors, [5]. The use of the heater exchanger reduces the thermal efficiency to a range of 65% to 80%.

DGFCHs are direct-fired, resulting in combustion gases (i.e., CO and CO2) formed directly in the heated air. In agricultural applications, the occupied (animal or plant) spaces where DGFCHs are present are continuously ventilated, unlike residential applications. The use of DGFCHs in greenhouses could result in safety concerns for workers, if prolonged operation throughout a day is necessary due to lower ventilation rates per floor area compared to animal buildings [4]. ANSI/ASHRAE [5] define the residential furnace thermal efficiency standard; however, there is limited standard thermal efficiency testing standards and guidelines for DGFCHs. Standardization of thermal efficiency testing and calculation allow for accurate comparison of different DGFCHs designs and capacities as well as generate informative equipment design specifications to support improved supplement heat capacity sizing and layout for agricultural buildings.

There are two types (direct and indirect) of combustion tests used to determine thermal efficiency. The direct test is a measure of the enthalpy gain of the air heated by the heater. Thermal efficiency can then be determined by dividing enthalpy gain by the sum of the heating fuel gross energy and electrical energy of the blower. This ratio of energy output over energy input for a given heater is the most accurate definition of thermal efficiency as the unit must obey the conservation of energy. The indirect test is a measure of the heating value of the fuel to determine the heat output less the flue gas or combustion gas enthalpy gain [5]. There are several challenges associated with applying the indirect test methodology to evaluate DGFCHs as the flue gas is mixed with the output air. The thermal efficiency calculated with the indirect methods typically is less accurate due to measurement error. Thus, DFGCHs are most appropriately evaluated by direct test methodology.

The objectives of this study were to (1) develop and validate a testing system and analytical analysis for determining the steady-state thermal efficiency of a DGFCH through a direct testing method; (2) evaluate the testing system across a range of DGFCH units; and (3) evaluate the economics of various DGFCH for agricultural uses.

#### **2. Materials and Methods**

#### *2.1. Thermal Efficiency*

The following sections will describe the methodology used to develop a DGFCH thermal efficiency measurement system. The system was designed to accommodate a range of rated heater capacities. For development purposes, three different DGFCHs (21.9 kW (75,000 BTU h−1), 29.3 kW (100,000 BTU h−1), 73.3 kW (250,000 BTU h−1)) were tested, procedural differences for each capacity are discussed in the testing procedure section. All tests were performed in the Air Dispersion Laboratory (ADL). See [6] for specifics of the laboratory.

#### *2.2. Instrumentation*

A Data Acquisition and Control (DAQC) system was developed to measure all necessary parameters to calculate the thermal efficiency of a DGFCH (analytical analysis to be discussed in a later section). The DAQC system consisted of two sub-systems, one to measure the input energy into the DGFCHs and one to measure the output enthalpy of the output air.

The instrumentation system for measuring input energy consisted of sensors to measure the propane input and electrical energy used by the DGFCH. Two different set-ups were created for propane input measurement. For DGFCHs rated <30 kW, a mass flow system was used, and for DGFCHs rated >30 kW, a volumetric device was used. This was performed as the minimum flow rates for traditional style diaphragm gas flowmeters exhibit high uncertainty for 30 kW rated heaters and below based on manufacturer specifications due to the extremely low flowrates. The mass flow system was not used for all heaters as the uncertainty in propane input was high, discussed in later section.

The mass flow system consisted of a 45.4 kg (100 lb) propane tank placed on a steel platform with a single-point 150 kg load cell with a 2 mV V−1 linear response signal (see Table 1 for specifications). The load cell signal was conditioned by a programmable-gain instrumentation amplifier (see Table 1 for specifications). A dual-stage regulator rated for 58.6 kW at 2.74 kPa was used (see Table 1 for specifications) with a 7.9 mm diameter flexible hose to connect the gas supply to the heater connection manifold.

**Table 1.** Energy input instruments used in the study in order as described in the text. FS = Full Scale.

The gas flowmeter (see Table 1 for specifications) was connected to the propane utility supply (two 3.79 m3 propane tanks with a line pressure of 2.74 kPa) using a 12.7 mm diameter flexible hose. The flowmeter's output was connected to the heater connection manifold with a 12.7 mm Internal Diameter (ID) flexible hose. Flowrate was determined by dividing the difference between the recorded initial and final positions of the dial on the flowmeter by the duration of the test.

The heater connection manifold consisted of a sediment trap (following manufacturer's recommendation) and two taps that contained a silvered PT100 Resistance Temperature Detector (RTD; see Table 1 for specifications) and a pressure transducer (see Table 1 for specifications) to monitor supply propane temperature and pressure, respectively. A precision resistor (see Table 1 for specifications) was used to create a voltage divider circuit to condition the RTD output signal with a supply voltage of 2.5 Volt DC (VDC) and the 4 to 20 mA signal from the pressure transducer. The RTD temperature was determined using the alpha equation (Equations (1) and (2)) for RTDs, with alpha determined using the resistance at 0 ◦C and 100 ◦C. Heater electrical power consumption was measured using a portable power meter (see Table 1 for specifications).

$$T\_{\bar{i}} = \frac{\left(\frac{R\_T}{R\_0} - 1\right)}{\alpha} \tag{1}$$

where

*Ti* = temperature (◦C); where *i* is RTD, propane, or exiting air *RT* = RTD resistance at *Ti* (Ω) see Equation (2) *R*0 = RTD resistance at 0 ◦C (100 Ω)

*α* = temperature coefficient of resistance (3.851E−4 Ω ◦C−1)

$$R\_T = \,^lR\_{d1} \* \left(\frac{V\_s}{V\_{out}} - 1\right) \tag{2}$$

where

*Rd*1 = resistance of voltage divider resistor (1 k Ω) *Vs* = supply voltage (2.5 VDC) *Vout* = analog voltage output measured by DAQC board (VDC)

Heater combustion exhaust measurement system consisted of a 0.2 m diameter, 2.44 m long round, galvanized steel duct that contained a PT100 RTD and a velocity pressure pitot tube (Figures 1 and 2, Tables 2 and 3). The same conditioning circuit and temperature equations were used as previously discussed. The supply end of the duct featured a flange fitting that was sealed with foil tape and secured to the face of the heater around the exhaust outlet. For heaters with larger output openings, a reduction flare fitting, 0.3 m to 0.2 m diameter and 0.3 m in length, was included. All seams on the inside of the duct were sealed with mastic caulking, and the outside of the duct was sealed with foil tape. The duct was wrapped with 51 mm (2 in.) thick, foil-backed fiberglass insulation. Additional insulation was placed on the seams on the fiberglass insulation, and all seams were sealed with foil tape. Velocity pressure was measured using a multi-range, differential pressure transducer (see Table 3 for specifications).

**Figure 1.** Test duct with locations of temperature measurement taps (D) and velocity pressure pitot tube (E). Measurements for the specific labels are in Table 1. Note diagram is not to scale.

**Figure 2.** Test stand configuration and DAQC system in ADL for DGFCH test.

**Table 2.** Dimensions of the test duct in Figure 1.

**Table 3.** Energy output instruments used in the study in order as described in the text.

Heater supply air (room) Dry-Bulb Temperature (Tdb) was measured using a grouping of four plastic-coated thermistors (see Table 3 for specifications) installed 0.25 m below the intake of the heater on the side opposite of the Liquid Propane (LP) intake and burner manifold. The four thermistors were multiplexed (Model: CD4052BE, Texas Instruments, Dallas, TX, USA) and conditioned with a voltage divider circuit consisting of a precision resistor (see Table 3 for specifications). Temperature from the thermistors was calculated using the β equation (Equation (3)) and then converted to Celsius for analysis. The thermistors were calibrated in a variable environment chamber (Series 7064-3140, Parameter Generation & Control, Black Mountain, NC, USA) at three unique temperatures (25 ◦C, 35 ◦C, and 45 ◦C) using the exact same DAQC setup prior to testing. The calibration regression was developed by regressing the chamber temperature recoded from the control panel on the x-axis and the average thermistor temperature on the y-axis. The final calibration equation was determined by taking the inverse of the previous calibration equation [7].

$$T\_{TRi} = 1/\left[ (1/298.15) + (1/\beta) \* \ln((R\_{d10} \* (V\_{s1}/V\_{out1} - 1))/R\_{ref}) \right] \tag{3}$$

where

*TTR* = thermistor temperature *i* (K) *B* = thermistor coefficient (3435) *Rd*10 = resistance of voltage divider resistor (Ω) *Vs*1 = voltage supply (5.0 VDC) *Vout*1 = analog voltage output measured by DAQC (VDC) *Rref* = thermistor reference resistance at 25 ◦C (Ω)

Supply air Tdb and Relative Humidity (RH) were measured in the same location using a portable hygrometer with a Liquid Crystal Display (LCD) screen (see Table 3 for specifications). The remote hygrometer was recorded manually at the start and end of each test.

All sensors, excluding the remote sensing unit, were recorded with a custom software interface (Visual Basic v19, Microsoft, Redmond, WA, USA). The software allowed for user-selectable sampling frequency and a calibration procedure to calibrate the load cell.

#### *2.3. Testing Procedure*

Data were collected on the three nominally rated DGFCHs (21, 29, and 73 kW) at the full rated output and half output setting. For the 21 and 29 kW rated heaters, three replicates of each output were measured. For the 73 kW heater, two replicates of each output were completed. This was due to issues with the test duct retaining residual heat

between tests. Prior to test initiation, the load cell was verified with a 45.4 kg (100 lb) standard set of weights. During a test, the DAQC was configured to sample continuously at 1 Hz. For DGFCH ratings > 18 kW, the total sampling duration was 500 s with the DAQC system sampling at 1 Hz. For DGFCH ratings < 18 kW, the total sampling duration was 900 s, with the DAQC sampling at 1 Hz. The time duration was different for the heater and settings rated for <18 kW based on the anticipated propane consumption to have a minimum change in the propane tank mass double the load cells accuracy. Electrical energy consumption was recorded from the power meter at 3 min into the test and assumed constant for the duration of the test. The three-minute mark was selected based on ANSI/ASHRAE Standard 103, considering steady-state combustion at three minutes [5]. Between performance tests, the ADL was ventilated for 60 min to maintain similar room temperatures at the start of each test.

#### *2.4. Thermal Efficiency Development*

The energy balance of a DGFCH is the balance of the energy input and output of the heater (Equation (4)). Output DGFCH energy can be defined as the sum of the combustion heat in the output airflow and the radiative and convective heat loss of the heater jacket (Equation (5)). For agriculture buildings, the main objective is to warm the occupied zones of the house. Thus, the goal of a DGFCH thermal efficiency is to quantify the efficiency of heating the output air (Equation (6); [5]).

$$Q\_{out} = Q\_{in} \tag{4}$$

where

*Qout* = energy output of a direct gas-fired circulating heater (kW) *Qin* = energy input (propane and electrical) of a direct gas-fired circulating heater (kW)

$$Q\_{\rm out} = Q\_{\rm a} + Q\_{\rm j} \tag{5}$$

where

*Qa* = energy output via air (kW) *Qj* = energy output via heater jacket (kW)

$$
\dot{Q}\_{\text{out}} = \dot{Q}\_{\text{in}} \times n\_{\text{ss}} \tag{6}
$$

where

*Q*ˆ *out* = measured energy output of heater output air (full or half load; kW) *nss* = steady-state sensible energy output efficiency (%)

Energy input into a DGFCH can be quantified as the energy from the liquid propane plus the electrical energy (Equation (7)). Electrical energy is included as most of the DGFCH supply air is drawn in over the motor and electronic components for cooling purposes. For the 21 and 29 kW rated heaters, input propane was measured gravimetrically, while input propane for the 73 kW rated heater was measured volumetrically, with a temperature and pressure compensated gas meter.

$$Q\_{in} = \left(\dot{\mathbf{x}}\_{LPG} \times \ HHV\_{prop}\right) + \mathcal{W}\_{elec} \tag{7}$$

where

. *xLPG* = mass flow rate of propane to heater (kg s−1) or volume flow rate of propane (m3 s−1) *HHVprop* = higher heating value of propane (kJ kg or kJ m−3) *Welec* = electrical energy consumption of heater (kW)

$$\mathcal{Q}\_{\text{out}} = \dot{m}\_a \times (h\_{\text{out}} - h\_{\text{in}}) \tag{8}$$

where

. *ma* = mass flow rate of exhaust air from heater (kg s−1) *hout* = enthalpy of output air from heater (kJ kg−1) *hin* = enthalpy of oncoming air to heater (kJ kg−1)

$$
\dot{m}\_a = \stackrel{\rightarrow}{V} \times \left(\pi r^2\right) \times \rho \tag{9}
$$

where -

 *V* = fluid velocity (m s−1) *r* = radius of the duct (0.10 m)

*ρ* = density of the air exiting heater (kg m−3)

$$
\overrightarrow{V} = \sqrt{\frac{2\delta p}{\rho}}\tag{10}
$$

where

*δp* = velocity pressure measurement (Pa)

The output energy was quantified by the enthalpy change of the air moving through the heater (Equation (8)). The mass flow of the DGFCH exhaust air was calculated using the density measurements and the pitot tube measurement (Equations (9) and (10)). The specific volume and enthalpy were calculated using the humidity ratio of the starting room conditions [1,8]

By combining the input energy (Equation (4)) and the quantified energy output (Equation (5)), the steady-state thermal efficiency can be solved for the DGFCH (Equations (11) and (12)). The thermal efficiency equation is similar in nature to efficiency equations used to evaluate solar air heaters [9,10]

$$
\dot{m}\_a \times \begin{pmatrix} h\_{out} \ - \ h\_{in} \end{pmatrix} = \begin{pmatrix} \left( \dot{\mathbf{x}}\_{prop} \times \begin{pmatrix} HHV\_{prop} \end{pmatrix} + W\_{elec} \right) \times \begin{pmatrix} n\_{ss} \end{pmatrix} \end{pmatrix} \tag{11}
$$

$$m\_{ss} = \left(\dot{m}\_a \times \left(h\_{out} - h\_{in}\right)\right) \div \left(\left(\dot{\mathbf{x}}\_{prop} \times \,\,\, HHV\_{prop}\right) + \mathcal{W}\_{\text{clcc}}\right) \tag{12}$$

#### *2.5. Thermal Efficiency Calculations and Analysis*

From each replicate, the mean of six random sub-samples of 30 s (*n* = 30 per subsample) of consecutive data was used to calculate the steady-state thermal efficiency (Equation (12)). A moving average of five data points was used to reduce fluctuations in the output air temperature. This processing step was justified due to the random error caused by turbulence and non-constant combustion of the heater. The mass flow rate of the propane was calculated by determining the weight change of the propane cylinder divided by the experiment duration. For DGFCH's rated >18 kW, a 4 s average before the heater solenoid was opened and after the test was used to determine the pre- and post-weight. For DGFCH's rated <18 kW, a 3 s average was used. The averaging of the cylinder weight was necessary to account for variation in the weight caused by vibration of the cylinder due to propane flow.

A sensitivity analysis of the steady-state thermal efficiency was performed with respect to the Higher Heating Value (HHV) of propane. Four sources of the HHV were used in the analysis due to the slight differences among sources (Table 4).

**Table 4.** Higher Heating Value (HHV) values from different sources used in the sensitivity analysis of the thermal efficiency equation.

#### *2.6. Uncertainty Analysis*

The final standard uncertainty of the thermal efficiency (Δ*nss*) was calculated from the standard uncertainty obtained from all key measurement inputs propagated through (Equation (13) through Equation (20), Figure 3).

**Figure 3.** Fishbone diagram of the final standard uncertainty in the thermal efficiency calculation. Each line represents the integration of a measurement and its respective uncertainty into a standard uncertainty calculation.

> A zeroth-ordered uncertainty budget was created for each measured input: voltage measurement for exiting air temperature (*Tex*; Table 5), a Differential Pressure Transducer (DPT; Table 6), dry-bulb temperature, and relative humidity for humidity ratio calculation (Tdb and RH; Table 7), propane flow for mass and volumetric ( . *x*; Table 8), and electrical energy consumption (*Welec*; Table 9) and included the manufacturer's accuracy and long-term stability, quantization error from the 14-bit Analog-to-Digital Converter (ADC), and the Standard Error (SE) from the experimental data (SE was removed from the budget as it changes with each experiment; [13]).

[a] ±0.5 Manufacturer's stated accuracy with ±5. VDC input range. [b] ±0.5 ADC resolution = 14-bit ADC resolution, ±5 VDC reference range = 3.1E−4 V BL−1. BL = Binary Level.

[a] Root sum square (at constant tdb), ±1.0% full scale (*δp*1) 0 to 248.8 Pa; [b] ±0.5 sensor resolution = (14-bit ADC resolution, ±10 VDC reference range = 6.1E−4 V BL−1) (sensor sensitivity)−1.

**Table 7.** Uncertainty budget for temperature and relative humidity sensor used to calculate humidity ratio.

[a] Manufacturer stated accuracy; [b] ±0.5 smallest display value = 0.01; [c] long-term stability is included in accuracy term from manufacturer specifications.

**Table 8.** Uncertainty budget for propane for both mass and flow basis.

[a] Manufacturer stated accuracy; [b] ±0.5 smallest display value = 0.01; [c] long-term stability is included in accuracy term from manufacturer specifications; [d] ±0.5 sensor resolution = (14-bit ADC resolution, 1 VDC reference range = 1.2E−4 V BL−1) (sensor sensitivity)−1.

**Table 9.** Uncertainty budget for electrical energy consumption of the heater.

[a] ±0.5 smallest display value = 50.

The uncertainty budget for the Tdb and RH includes a display resolution component as it was manually recorded (Table 7). The propane volumetric measurement also had a display resolution component from the dial on the meter display from manually recording measurements.

The standard uncertainty associated with the area of the duct was determined by propagating the standard uncertainty of the duct radius (Δr = 2.70 × 10−3 m, one-half the reading resolution set by the manufacturer) through the area equation of a circle resulting in a standard uncertainty of the area Δa = 3.39 × 10−3 m2.

$$
\Delta VDC\_i^2 = \left(\frac{\Delta VDC}{n}\right)^2 + SE^2 \tag{13}
$$

where

Δ*VDCi* = combined standard uncertainty in analog voltage measurement (*Tex* VDC) *SE* = standard error of the voltage measurement (VDC)

$$
\Delta T\_{\rm ex}^2 = \left(\frac{\partial T\_{\rm ex}}{\partial VDC} \, \Delta VDC\_i\right)^2 + \left(\frac{\partial T\_{\rm ex}}{\partial R\_{d1}} \, \Delta R\_{d1}\right)^2 + \left(\frac{\partial T\_{\rm ex}}{\partial R\_{\rm ref}} \, \Delta R\_{\rm ref}\right)^2 + RMSE^2 \tag{14}
$$

where

Δ*Tent* = combined standard uncertainty in entering air temperature (◦C) Δ*Rd*1 =1kΩ resistor in divider circuit (±0.05%; Ω; rectangular distribution) Δ*Rref* = 100 Ω at 0 ◦C (±0.12%; Ω; rectangular distribution)

*RMSE* = root mean square error from nonlinear regression (◦C)

The standard uncertainty of the entering temperature measured by the thermistor did not include propagation of uncertainty through (Equation (3)) as the DAQC system was identical in both calibration and testing. Thus, it was assumed that the RMSE of the calibration encompasses the uncertainty associated with the measurements propagated through the equation.

$$
\Delta T\_{\rm ent}^2 = \, A \text{CC}^2 + \, RMSE^2 + \, SE^2 \tag{15}
$$

where

Δ*Tent* = combined standard uncertainty in entering air temperature (◦C) *ACC* = manufacturer's accuracy (±0.2 ◦C; rectangular distribution) *RMSE* = root mean square error from linear calibration regression (◦C) *SE* = standard error of the mean (◦C)

$$
\Delta \mathcal{W}^2 = \left(\frac{\partial \mathcal{W}}{\partial t\_{db}} \Delta t\_{db}\right)^2 + \left(\frac{\partial \rho}{\partial r h} \Delta r h\right)^2 \tag{16}
$$

where

Δ*W* = combined standard uncertainty in humidity ratio (kgH2O kgda−1)

$$
\Delta \rho^2 = \left(\frac{\partial \rho}{\partial t\_{\varepsilon \chi}} \Delta t\_{\varepsilon \chi}\right)^2 + \left(\frac{\partial \rho}{\partial W} \Delta W\right)^2 \tag{17}
$$

where

Δ*ρ* = combined standard uncertainty in moist air density (kg m−3)

The standard uncertainty in the mass flow (Equation (18)) was determined from the propagation of the velocity pressure uncertainty, density uncertainty, area uncertainty, and accuracy of the pitot tube system.

$$
\Delta \dot{m} \dot{a}^2 = \left(\frac{\partial \dot{m} a}{\partial dp} \, \Delta dp\right)^2 + \left(\frac{\partial \dot{m} a}{\partial \rho} \, \Delta \rho\right)^2 + \left(\frac{\partial \dot{m} a}{\partial a} \Delta a\right)^2 + A \, \text{CC}\_v^{-2} \tag{18}
$$

where

Δ . *ma* = combined standard uncertainty in mass flow moist air (kg s−1) *ACCv* = accuracy of the velocity measurement (1.0%, [5])

$$\left(\Delta h\_{\rm ex-cnt}\right)^{2} = \left(\frac{\partial h\_{\rm out-in}}{\partial W}\,\Delta W\right)^{2} + \left(\frac{\partial h\_{\rm out-in}}{\partial T\_{\rm cut}}\,\Delta T\_{\rm cut}\right)^{2} + \left(\frac{\partial h\_{\rm out-in}}{\partial T\_{\rm ex}}\,\Delta T\_{\rm ex}\right)^{2} \tag{19}$$

where

Δ*hout*−*in* = combined standard uncertainty in change in enthalpy (kJ kg−1)

The standard uncertainty of . *xLPG* included in (Equation (20)) included the SE term for the initial and final measurement from the gravimetric system as a multi-second average was used. It was assumed that there was no uncertainty in the run time of each test factored into . *xLPG*.

$$
\Delta n\_{\rm ss}{}^2 = \left(\frac{\partial n\_{\rm ss}}{\partial \dot{m} \dot{a}} \Delta \dot{m}\right)^2 + \left(\frac{\partial n\_{\rm ss}}{\partial h\_{\rm crx-cnt}} \Delta h\_{\rm crx-cnt}\right)^2 + \left(\frac{\partial n\_{\rm ss}}{\partial \dot{x}} \Delta \dot{x}\_{LPG}\right)^2 + \left(\frac{\partial n\_{\rm ss}}{\partial W\_{\rm clxc}} \Delta W\_{\rm clxc}\right)^2 \tag{20}
$$

$$
\text{where}
$$

Δ*nss* = combined standard uncertainty in thermal efficiency (%)

The combined standard uncertainty was calculated using the data from each replicate of the heater model and output setting. Thus, for each replicate and output, a standard deviation was calculated for the uncertainty due to variations in measurements between each replicate.

#### *2.7. Economics of Heater Efficiency*

Economic implications of heating agricultural buildings were assessed by modifying the degree-day heating method (Equation (21); [3]). The method is simpler than other heating cost methods that use binned ambient weather data, with the consequence of overestimating the heating season costs. The efficiency represents the Annual Fuel Energy Use Efficiency (AFEUE). The AFEUE is based on the ASHRAE standard 103 [5]. For DGFCH without a pilot light, the AFEUE is the same as the steady-state thermal efficiency. For a DGFCH with a pilot light, the calculation of AFEUE considers the proportion of time with the pilot light on and its efficiency over the duration of the heating season.

$$E = \mathcal{C}\_D \times \left[ (\left( H\_L \times D \times \mathcal{S} \right) \div \left( \Delta t \times n\_{\rm ss} / A \rm LI} \times HHV\_{prop} \right) \right] \tag{21}$$

where

*E* = propane usage per season (kg) *CD* = correction factor for degree-days (0.6) *HL* = design heating load at minimum ventilation (kW) *D* = degree-days (◦C day) *S* = seconds per day (86,400 s d−1) Δ*t* = design temperature difference, inside/outside (◦C) *HHVprop* = higher heating value for propane (kJ kg−1)

The calculation of the heating load for a building (heat gain/heat loss) was developed to account for the wide variety of building specifications (Equation (22)). This was performed by accounting for the heat gain on the basis of floor area. The heating load is to be calculated at the wintertime 99% design temperature for the given location.

$$HL = \left[ (UL \times \Delta T) + \left( \dot{V}\_b \div \upsilon \times \mathbb{C}\_p \times \Delta T \right) \right] + \left[ q\_\% \times A\_f \right] \tag{22}$$

where

*UA* = overall building heat loss coefficient (W ◦C−1) . *Vb* = minimum ventilation rate (m3 s−1) υ = specific volume of air/water mixture (kg m−3) *Cp* = specific heat of air (J kg−1) *Qg* = heat production per unit of floor area (W m−2) *Af* = floor area of building (m2)

A sensitivity analysis of the energy usage (Equation (21)) was performed based on recommended thermal resistance (R values) for various building components, internal design temperature, two different minimum ventilation rates, and heater efficiency. This analysis used a high heat gain per unit floor area of 1.6 W m−2. It was assumed that the outside design temperature was −29 ◦C (99% heating value for Minnesota; [8]). For this economic analysis, the building size was 45.7 m (L), 15.2 m (W) with a 2.5 m ceiling (150 ft L, 50 ft W with 8 ft ceiling). The building's thermal analysis assumed no windows and doors and did not include any perimeter or floor heat loss. The R values, internal building temperature, and minimum ventilation rates (Table 10). For the analysis of the first three components (level of insulation, design internal Tdb, and minimum ventilation rate), the heater efficiency was 99%, and the cost of propane was assumed at USD \$0.53 kg−1 (\$1 gallon−1).

**Table 10.** Sensitivity analysis inputs for R values, internal building temperature, and minimum ventilation rate. Recommended levels of insulation for homes based on climate [3].

[a] ACH = air changes per hour.

#### **3. Results and Discussion**

#### *3.1. Thermal Efficiency*

The average values of the operational parameters and results for the three DGFCHs tested at each output performed as expected (Table 11). The procedures developed were capable of detecting the relatively small differences in input and output energy, although the 73 kW did not reach rated operation. The time to steady-state assumed from ASHRAE 103 was represented in the results of this study (Figure 4), with all heaters and replicates showing similar trends in temperature change over time. Future tests should involve a high accuracy CO and O2 sensor to verify the thermal efficiency through stoichiometry (indirect combustion test). A unique trend was noted in the comparison of the half and full output thermal efficiencies. Heater models A and B exhibited higher thermal efficiencies at full output, while heater model C exhibited a higher thermal efficiency at half output. The thermal efficiency of over 100% is related to the system's uncertainty, discussed in the next section. The general trend of efficiency increasing with heater output could be influenced by the test stand gaining heat during the test and affecting measurements. This was the reason behind two replicates with heater model C. The heater models used in the study had limited thermal efficiency values reported, model A had no reported values, and models B and C reported 99.99% thermal efficiency. The reported thermal efficiencies in this study are higher than what is assumed for DGFCH in buildings 80%, with high humidity and low ventilation [4]. However, it was not within the scope of this study to evaluate the heaters under such conditions that are believed to drastically reduce the performance of these heaters in buildings. The thermal efficiency of the DGFCH used in this study is significantly higher (>90%) compared to the thermal efficiency of common gas-fired residential units (63% to 44%; [13]). This drastic difference is mainly contributed to the difference in the heaters as the residential heaters.

**Table 11.** Average (standard deviation) of test conditions and measured performance values for all three DGFCH models and outputs.

[a] Measured input calculated using HHV from [5]. [b] Thermal efficiency calculated using HHV from [5].

**Figure 4.** Temperature gain of the exiting air over the course of a test for heater model A at full output.

The measured input energy to heater model C at full output was significantly lower than its nominal rating of 73.3 kW. Upon investigation of the operational parameters of the tests, it was noted that manifold pressure exhibits a negative trend in relation to the measured input (Figure 5). The temperature of the propane did not vary between any heaters and replicates. This is most likely caused by an improper design of the propane supply of ADL that limits the maximum output due to propane flow and not the heater reaching its rated output.

**Figure 5.** Measured input of propane (x-axis) versus the manifold propane pressure (y-axis).

The sensitivity analysis of the thermal efficiency equation in relation to the HHV was performed for heater model A at full output. The maximum difference noted in the analysis is between HHV source A (97.3%) and D/C (96.0%), with source B (96.5%) falling in-between. While all four sources are relatively similar in value, this highlights the differences in reported values among groups (sources A, B, C) and the federal government (source D). A unified HHV between all organizations is imperative to avoid discrepancies in design and in the rating of all types of propane combustion heaters.

#### *3.2. Uncertainty Analysis*

The expanded (coverage factor = 2; ~95% confidence interval) uncertainty ranged from 13.1% (Standard Deviation; SD = 0.1%) to 30.7% (SD = 8.0%) for the input energy rating 56.3 to 11.4 kW; the SD is from the multiple replicates of each heater output tested. There is a distinct difference in the variability of the uncertainty between the replications for a heater based on the propane measurement system. The mass flow system had significantly greater variability (Figure 6). The breakdown of the uncertainty contributions (Table 12) demonstrates that for the mass flow system, the propane mass flow and air mass flow are the largest contributors while using volumetric flow measurement, the air mass flow, and change in enthalpy are the largest contributors. Overcoming the variability in the propane mass flow for future tests could require a hardware-based signal conditioning filter to account for the signal noise due to vibrations. Alternatively, extended test duration could be used to reduce variability.

**Figure 6.** Expanded thermal efficiency (y-axis) versus heater input rating (x-axis). The blue circles represent the measurements using the mass flow system, and the red squares represent the volumetric flow measurements.

**Table 12.** Expanded uncertainty component contributions for Equation (20) based on the propane measurement system. The values in parentheses represent 1 standard deviation of the mean.

The expanded uncertainty component contributions for the volumetric propane measurement system suggest that if a high accuracy, low flow rate option was feasible, it would reduce the uncertainty considerably. To reduce uncertainty in future experiments, a higheraccuracy pressure transducer and RTD would be warranted.

#### *3.3. Economics of Heater Efficiency*

The economic estimates for heating the specified building using design conditions show a wide heating cost for a heating season (Table 13). There are a few limitations to the degree-day method for estimating this cost. First, it assumes a constant balance temperature for a building. Balance temperature is the ambient temperature that separates the need for heating or cooling, which is different for the various types of agricultural buildings (livestock versus greenhouse). For livestock buildings, this approach cannot account for the increase in the animal's heat production through a growth cycle and the change in sensible heat production with indoor temperature. This economic assessment confirms that the greatest heating cost factor is the minimum ventilation rate, with building internal temperature and insulation level being the other factors impacting cost. It should be cautioned against decreasing the ventilation rate of a building to reduce the seasonal heating cost as air quality and moisture levels will be directly impacted and could negatively impact production within the building. The costs and corresponding energy inputs are higher than reported values from field studies for animal agricultural buildings with similar heat production values per floor area [14–16]. These differences are most likely attributed to the method used and the local weather data differences from the reported study's locations.

**Table 13.** Economic analysis (values reported in USD) for a general agricultural building for a heating season with a heat production of 1.6 W m−2.

The sensitivity analysis for the scenario regarding heater efficiency (Table 14) shows the economic impact of reducing the DGFCH thermal efficiency by 1% for the combinations of internal temperature and building envelope insulation level. The range of costs is caused by the different range of supplemental heat needs based on internal temperature and insulation value. The magnitude of the costs would be significantly higher if an AEU efficiency for a heater with a pilot light.

**Table 14.** Economic cost (values reported in USD) of reducing the DGFCH thermal efficiency by 1% for a heating season.

The methods used to perform this economic analysis could be adapted to a range of scenarios to evaluate utility costs for building types, climates, and construction for a cost estimate. Other methods for estimating heating costs are binning methods, where data is binned by day for heating demand (time and temperature difference). This method has advantages in being able to incorporate livestock heat production curves and variable ventilation rates. The downside is that this method requires extensive calculations for each time period within each day and is not a simple calculation for agriculture professions to estimate heating costs for a building.

#### **4. Conclusions**

A custom testing apparatus and instrumentation system were developed and evaluated for Direct Gas-Fired Circulating Heaters (DGFCH) commonly used in agricultural buildings. Three DGFCH were used for the evaluation testing of the system. The steadystate thermal efficiency of the tests was determined using data from the testing system and based on calculations from a standard for residential systems [5]. An uncertainty analysis was performed on the thermal efficiency calculation for the given measured inputs. An economic analysis for heating common agricultural buildings was performed. An inadequate propane service limited the systems testing ability for 73 kW rated heaters. A

sensitivity analysis of the thermal efficiency in respect to the Higher Heating Value (HHV) results in a maximum difference between sources of HHV of 1.3% efficiency. This highlights the need for standardization between all organizations involved. The uncertainty of the system ranged from 13.1% to 30.7% across the input range of 56.2 to 11.4 kW, respectively. The propane mass was the key source of uncertainty at the lower ratings. A more robust hardware system and longer run times would reduce this key source of uncertainty. At higher input ratings, the mass flow of air was the key source of uncertainty. Reducing this would require the use of higher-accuracy sensors. The economic analysis showed that for a 1% difference in heater efficiency across the combination of inputs (insulation and inside temperature), the costs per heating season ranges from USD \$61.3 to \$72.0. Future evolutions of this system are needed to reduce the overall uncertainty in thermal efficiency. It is suggested that a flowmeter with lower minimum flow ratings be used or longer test durations to reduce the main source of uncertainty. Other instrumentation improvements are needed for the output air and velocity measurement in the duct. Overall this system shows a great prospect for testing DGFCH commonly used in agricultural buildings. This study has highlighted the need for standardized procedures for quantifying DGFCH thermal efficiency specific to agricultural applications. The unique nature of these heaters presents challenges for testing and has a large economic impact on the agricultural operation during heating seasons.

**Author Contributions:** Conceptualization, B.C.S., B.C.R., and S.J.H.; methodology, B.C.S., B.C.R., and S.J.H.; software, B.C.S. and S.J.H.; validation, B.C.S.; formal analysis, B.C.S.; resources, B.C.R.; data curation, B.C.S.; writing—original draft preparation, B.C.S.; writing—review and editing, B.C.S., B.C.R., and S.J.H.; supervision, B.C.R.; project administration, B.C.R.; funding acquisition, B.C.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research is a product of the Iowa Agriculture and Home Economics Experiment Station, Ames, Iowa. Project Number IOW04100 is sponsored by the Hatch Act and State of Iowa funds. The content of this article is, however, solely the responsibility of the authors and does not represent the official views of the USDA.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** The authors would like to thank the contributions of undergraduate assistants, Devri Cliburn and Julia Bowman, for their assistance with the calibration and construction of the testing apparatus.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

#### **References**

## *Article* **Design and Experimental Research on Soil Covering Device with Linkage and Differential Adjustment of Potato Planter**

**Zhiqi Zheng, Zuoli Fu \*, Chenyang Wang, Yuxiang Huang and Jinpu He**

College of Mechanical and Electronic Engineering, Northwest A&F University, Yangling 712100, China<EMAIL_ADDRESS>(Z.Z<EMAIL_ADDRESS>(C.W<EMAIL_ADDRESS>(Y.H<EMAIL_ADDRESS>(J.H.)

**\*** Correspondence<EMAIL_ADDRESS>
**Abstract:** When the potato planter works on sloping field, it will cause problems such as poor film mulching quality due to the difference in volume of soil covering both sides of the discs and the inconvenient adjustment of the soil covering disc. The soil covering device with linkage and differential adjustment was designed to improve the mulching quality. The main research content includes explaining the structure and principle of the soil covering device and analyzing the structural parameters of the adjustment mechanism. The field experiment was completed to verify the performance of soil covering device, which takes the stability coefficient and uniformity coefficient of the volume of covering soil as factors. The result shows the following: (1) The volume of covering soil changes exponentially with the angle of the disc through data fitting, which can standardize the angle of covering disc; and (2) when the angle of disc is 30◦ and 60◦, respectively, the uniformity coefficient of volume of covering soil is lower than 1.4, which has premium soil covering quality. When the angles of the discs on both sides differ greatly, the stability coefficient of volume of covering soil is 0.41, which can meet the requirements of the mulching quality of potato planter. This research provides the technical support for high-quality potato planting.

**Keywords:** potato planter; soil covering device; design; experiment

#### **1. Introduction**

In Ningxia, Gansu, Shaanxi and other regions of China, the plastic film mulching potato planter are widely used. The cultivated land in these areas is mainly hills and is mountainous and there are a lot of planters working on sloping fields. The main function of mulching is to maintain warmth and moisture levels, to promote the rapid germination of seeds and to achieve the purpose of increasing yield and reducing the growth period. The effect of covering soil is an important evaluation of the quality of mulching and there is a large impact on the volume of covering soil on both sides of mulch in sloping fields. Furthermore, it is necessary to make tentative adjustments based on experience and it is difficult to achieve best the covering quality. Therefore, it is of great significance to develop a soil covering device for potato planters that achieve the same volume of covering soil on both sides in sloping fields [1,2].

The research on the soil covering device of potato planter has been conducted earlier in foreign countries with advanced technology and high reliability. The typical equipment includes GL series potato planter produced by GRIMME in Germany and the LK series potato planter produced by STRUIK in Netherlands. These machines have high efficiency and quality, which are suitable for areas with large-scale fields rather than hills and mountainous areas. Recently, the development of potato planting machinery technology and equipment in China has become rapid. The potato planter was designed and tested and includes the sowing depth control device and wide ridge pattern suitable planter [3–5]. Liu and other scholars designed the film mulching potato planter, which applied to rice,

**Citation:** Zheng, Z.; Fu, Z.; Wang, C.; Huang, Y.; He, J. Design and Experimental Research on Soil Covering Device with Linkage and Differential Adjustment of Potato Planter. *Agriculture* **2021**, *11*, 665. https://doi.org/10.3390/ agriculture11070665

Academic Editor: José Pérez-Alonso

Received: 30 May 2021 Accepted: 12 July 2021 Published: 14 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

quinoa, potato and so on [6–8]. Chen and other scholars put forward film mulching technology, soil conveying technology, etc., which provide technical support for the research and development of potato planters [9–13]. The study on the optimization of working parameters and improvement of the design structure of the potato planter and mulching device was conducted by Xin and others [14–18]. These technologies and equipment have been experimented and verified to basically suit the requirements of field working. When working on flat fields, the quality of covering soil is higher; however, in the sloping field, there will be inconsistencies in the volume of covering soil on both sides and the stability of cover will decrease. The covering quality is reduced, which is not conducive to heat and moisture preservation and will affect the high yield and premature delivery of potatoes [19–21].

Therefore, a soil covering device was designed to improve seeding quality by linkage and differential adjustments of the disc separately and the field experiment was completed to verify the performance of soil covering device, which takes the stability coefficient and uniformity coefficient of the volume of covering soil as factors.

#### **2. Material and Methods**

#### *2.1. The Device of Experiment*

As shown in Figure 1, the potato planter mainly includes fertilizer box, seeder, seed box, ridging machine, laminating machine and soil covering device. The prototype takes the Qingdao HONGZHU potato planter and the planting method utilizes double rows on a ridge. The fertilizer from the fertilizer box is applied to the middle position of the ridge and planted on both sides of the ridge.

**Figure 1.** The structure scheme of potato planter.

A spoon-chain seeder metering device is used to take out the lumpy potato seeds from the seed box and to shake the excess potato seeds back into the seed box under the action of the vibrating seed cleaning device. The ridging machine trims the ridge shape according to the set value. The laminating machine covers the film on the ridge surface and the packer wheel and soil covering machine takes the soil covering both sides of film.

Figure 2 shows the structure diagram of soil covering device and the device mainly includes covered soil discs, adjusting mechanism, frame and so on. The distance between the selection centers of the covered soil discs on both sides is adjustable from 1100 to 1500 mm, which can suit the work of 0.8–1.2 m width film and the maximum adjustment of the rotation angle of the covering disc is 75◦.

**Figure 2.** The structure scheme of soil covering device.

The tractor pulls the device forward and the press roller presses the film on the ridge surface. The packer wheel then presses down the film on the ridge side and the covered soil disc pushes the loose soil of the furrow to the film on the ridge side to complete the entire film.

The uniformity and stability of the volume of covering soil is an important factor for the quality of the soil covering device. When working on flat ground, the angle of the covered soil discs on both sides is generally required to be consistent and the covering quality can basically meet the operation requirements. However, as shown in Figure 3, when working on sloping field, the volume of soil on both sides will be affected by the slope angle and the amount of soil on the upper side is much larger than that on the lower side. Therefore, it is necessary to adjust the angle of the covered soil disc to keep the volume of covering soil on both sides consistent.

**Figure 3.** The soil covering at the sloping field.

#### *2.2. The Data Analysis*

The main function of the adjustment mechanism is to adjust the angle of the soil covering discs on both sides to improve the mulching quality.

As shown in Figure 4, the adjustment mechanism mainly includes the covered soil disks on both sides, the rocker arm, the four-bar mechanism composing the frame, the connecting rod and the swing arm. The connecting rod adopts a two-way screw structure and, by rotating the middle nut, the length of the connecting rod can be changed and the angle of the discs on both sides can be adjusted in linkage under the action of the four-bar mechanism. The two sides of the connecting rod and the swing arm are correspondingly connected by hinged shafts. On the other hand, the rocker arm can rotate around the end point. During the rotation, the connecting rod will rotate synchronously to realize the differential adjustment of the covered soil disc. The angles α and β of both side discs are correspondingly related to the mechanism parameters of the four-bar mechanism.

(**a**) The linkage adjustment (**b**) The differential adjustment

**Figure 4.** The schematic diagram of the adjusting mechanism.

Firstly, the model equations of the angle and related parameters are established through a theoretical calculation and then the corresponding relationships between different slope angles and adjustments are obtained based on the test data.

(1) The linkage adjustment

The intercept of the two-way nut screw is *p*, the length of the swing arm is *Ls* and *LLs*, respectively. When the function of the hinge axis is not considered, the Equation (1) exists and is described as follows:

$$\begin{cases} \frac{np}{\Delta a} = \frac{Ls}{LLs} \\ a - a\_0 = 2arcsin\frac{np}{2Ls} \end{cases} \tag{1}$$

where *n* is number of rotations and Δ*a* is the turning distance of the end of the swing arm in mm.

(2) The differential adjustment

When the four-bar mechanism is composed of the frame, the connecting rod and the swing arm is arranged in a parallel. If the rocker arm is deflected counterclockwise and the angle is *γ*, then there exists the following equation:

$$
\alpha - \mathfrak{a}\_0 = \beta - \beta\_0 = \gamma
$$

where *a*0 and *β*0 are the initial angles of the discs on both sides.

When the four-bar mechanism is not arranged in a parallel, the calculation method of steering trapezoid is referred to and described in the following.

$$\begin{cases} \begin{aligned} \Delta a &= \gamma \\ \Delta \beta &= \gamma - \arcsin \frac{\sin((\gamma + \gamma\_0)}{\sqrt{\left(\frac{L\_s}{L\_s}\right)^2 + 1 - 2\frac{L\_s}{L\_s}\cos(\gamma + \gamma\_0)}} - \arccos \frac{\frac{L\_s}{L\_s} [2\cos \gamma\_0 - \cos(\gamma + \gamma\_0) - \cos 2\gamma\_0]}{\sqrt{\left(\frac{L\_s}{L\_s}\right)^2 + 1 - 2\frac{L\_s}{L\_s}\cos(\gamma + \gamma\_0)}} \end{aligned} \end{cases} \end{cases}$$

where *γ*0 is the initial bottom angle of trapezoidal four-bar mechanism in degrees and *L* is the length of the connecting rod (short side) in mm.

The size of adjusting mechanism in potato planter is shown in Figure 5. The initial bottom angle of the trapezoidal four-bar mechanism is determined by the length of the swing arm, the intercept of the nut screw and the number of rotations. By using the Equation (3), the rotation angle of the swing arm can be quantified. It is also necessary to calibrate the inclination requirements of the upper and lower soil covering discs according to different slope angles. So far, the design parameters of the four-bar mechanism can be obtained.

**Figure 5.** The size of adjusting mechanism.

The purpose of the experiment is to calibrate the volume of covered soil under different angles of the disc. On the other hand, when the angles of the discs on both sides are inconsistent, the experiment detects the difference and stability of the volume of covered soil on both sides. The flat fields have been used to complete the field experiment by using the stability coefficient of the difference in the volume of covered soil on both sides to evaluate the quality of the covering soil.

The potato planter used in the experiment adopts a self-made soil covering device with adjusting mechanism, which based on the HONGZHU-1.2 large ridge double-row planter. The adjusting mechanism was shown in Figure 6 and the mulch uses ordinary plastic film with a thickness of 0.0012 mm and a width of 1.2 m.

**Figure 6.** The adjusting mechanism.

The tractor uses DONGFANGHONG-504 from LUOTUO (about 36.8 kW) and the operating speed is about 3.6 km/h.

#### *2.3. The Method of Experiment*

On 9 April 2020 at the Northwest A&F University test farm, the operation performance test of soil covering device was carried out. The field of experiment has a length of 110 m and a width of 10 m. It is a typical columbine soil condition. Moreover, corn was planted in the first stubble and two rotary tillage treatments were carried out after harvest. The soil moisture content of the cultivated layer (5–10 cm) is 18.5%(w.b.) and the average firmness of soil is 2.05 MPa. The soil covering disc adopts a concave disc with a diameter of 45 cm and the depth of the disc is adaptively adjusted by spring preloading and force adjustment, the measured data shows that the average depth of the disc into the soil is 14.27 mm.

In order to calibrate the volume of soil covering under the adjustment of different disc angles, five points were randomly selected to measure in each level and 30◦, 45◦ and 60◦, respectively, were used as linkage adjustment angles of the covered soil disc. Moreover, the angles of the covering discs on both sides are selected to be 30◦ and 60◦, respectively, as the difference adjustment. The Table 1 shows the method for adjusting the angle of the covered soil discs on both sides. The sampling points were selected according to the five-point method and randomly selected on a ridge in stable working area as the sampling points. Repeat three ridge operations under the same angle of discs and take the average for calculation.

As shown in Figure 7, the volume of covered soil is calculated by the area of the section and the length of the sampling (take 1 m), as shown in Equation (3) as follows:

$$N = \frac{a+b-h}{2}hl\tag{3}$$

where V is the volume of covered soil, m3; a is the length of the section, m; *b* is the width of the section, m; *h* is the height of the section, m; *l* is the length of the sampling, m.

**Figure 7.** The diagram of the volume of covering soil and section.

Therefore, the main measurement factors include the width, height and length of the covered soil section. Considering that the ridge structure is relatively stable (the slope angle is about 45◦), the area of the covering soil section can be formed by subtracting the area of an isosceles right triangle from a right-angled trapezoid for which the side lengths are a, b and h.

The area of the covering soil section can be formed by a right-angled trapezoid with width, length and height as three sides and minus the area of the isosceles right triangle with the height as the right side.

The stability coefficient and uniformity of the volume of covered soil are used as indicators to investigate the feasibility of the differential adjustment and linkage adjustment. Refer to the "uniformity determination" and "stability determination" methods in the test method of the quality evaluation technical specification for fertilization machinery (NY/T 1003-2006) for calculation. The uniformity mainly measures the uniformity coefficient of the volume of covered soil. When the angles of the soil covering discs on both sides are different, the uniformity coefficient of each side is measured separately, which can characterize the feasibility of the differential adjustment. The stability mainly measures the stability coefficient of the volume of covered soil. When the angle of the soil covering discs on both sides are same, the stability coefficient of each side is measured separately, which can characterize the stability of the linkage adjustment.

#### 2.3.1. The Calculation of Uniformity Coefficient

The uniformity coefficient of the volume of covered soil includes standard deviation and the variation coefficient. Before the experiment, the parameters of the soil covering plate were adjusted and the tractor was moving forward steadily at an average speed of 3.6 km/h. Five points were randomly selected in the test area and the section parameters of the covered soil on both sides were correspondingly measured. The volume of covered soil was calculated according to the above method and, according to Equation (4), the standard deviation and the variation coefficient of the uniformity on both sides were calculated as follows:

$$\begin{cases} \quad \mathbf{x} = \frac{\sum\_{i=1}^{n} x\_i}{n} \\ \quad \mathbf{S} = \frac{\sum\_{i=1}^{n} (x\_i - \mathbf{x})^2}{n-1} \\ \quad V = \frac{\mathbf{S}}{\mathbf{x}} \times 100\% \end{cases} \tag{4}$$

where *x* is average volume of covered soil, m3; *xi* is volume of covered soil of the point, m3; *S* is the standard deviation of the uniformity; *V* is the variation coefficient of the uniformity, %; *n* is the number of measuring points, (*n* = 5).

#### 2.3.2. The Calculation of Stability Coefficient

The stability coefficient of the volume of covered soil includes the standard deviation and the variation coefficient. The test method is the same as the calculation of the uniformity coefficient. Refer to the calculation in Equation (3), where, *x* is average value of the difference in the volume of covered soil on both sides, m3; *xi* is difference of volume of covered soil on both sides of the point, m3; *S* is the standard deviation of the stability; *V* is the variation coefficient of the stability, %; *n* is the number of measuring points, (*n* = 5).

#### **3. Results and Discussion**

As shown in Figure 8, when the angles of the soil covering discs on both sides are 30◦ and 45◦, the volume of covered soil at the measuring points is 0.28–0.73 × 10−4 m3 and 1.07–1.48 × 10−4 m3, respectively. The volume of covered soil is quite different between both the angles of discs, but the difference is relatively close at different points. The maximum difference of 0.79 × 10−4 m3 appears at the first point and the minimum difference of 0.61 × 10−4 m3 appears at the fourth point. When the angle of the covering disc is 30◦, the average volume of covered soil is 0.55 × 10−4 m3 and when the angle of the covering disc is 45◦, the average volume is 1.12 × 10−4 m3.

**Figure 8.** The volume of covered soil at different sampling points (Both the angles of left and right disc are 45◦and 30◦, respectively).

It shows that the volume of covered soil is different at each point but, at the same point, there is a relatively fixed difference of the soil covering discs on both sides.

As shown in Figure 9, when the angles of the soil covering discs on both sides are 45◦ and 60◦, the volume of covered soil at the measuring points is 0.67–1.11 × 10−4 m3 and 1.24–1.52 × 10−4 m3, respectively. The maximum difference of 0.60 × 10−4 m3 appears at the first point and the minimum difference of 0.40 × 10−4 m3 appears at the second point. When the angle of the covering disc is 45◦, the volume of covered soil is 0.85 × 10−4 m3 and when the angle of the covering disc is 60◦, the average volume is 1.37 × 10−4 m3.

**Figure 9.** The volume of covered soil at different sampling points (Both the angles of left and right disc are 60◦and 45◦, respectively).

When the angle of covering disc is 45◦, the volume of covered soil is about 31% smaller than that of the previous group. The reason may be that there will be some interference between adjacent rows during the operation. Specifically, the engaged previous row of soil covering will reduce the amount of soil in the ditch and, as a result, the amount of floating soil on the adjacent side becomes less while the volume of covered soil is reduced.

As shown in Figure 10, when the angles of the soil covering discs on both sides are 30◦ and 60◦ the volume of covered soil at the measuring points is 0.61–0.88 × 10−4 m3 and 1.19–1.49 × 10−4 m3, respectively. The maximum difference of 0.70 × 10−4 m3 appears at the first point and the minimum difference of 0.58 × 10−4 m3 appears at the fourth point. When the angle of the covering disc is 30◦, the volume of covered soil is 0.77 × 10−4 m3 and when the angle of the covering disc is 60◦, the average volume is 1.40 × 10−4 m3.

**Figure 10.** The volume of covered soil at different sampling point (Both the angles of left and right disc are 30◦and 60◦, respectively).

When the angle of the covering disc is 30◦, the amount of soil covering is about 29% compared to that of the first group. The reason may be that the first row is close to the edge of the plot during operation and the soil is not as loose as the middle position and so the covering amount is lesser. When the inclination angle is 60◦, the amount of soil covered is close to that of the second group with a difference of less than 2.2%.

Therefore, the reason for the larger change in the volume of covered soil may be caused by soil conditions and the curve can be explained as shown in Figures 8–10.

When the angles of covering discs 30◦, 45◦ and 60◦ are, respectively, fitted and the relationship between the volume of covered soil and the angle of disc is obtained as a quadratic function, as shown in Equation (5). The angle of the disc can be calibrated through the following equation:

$$\mathbf{Y} = 0.001\mathbf{X}^2 - 0.0701\mathbf{x} + 1.9597\tag{5}$$

where X is the angle of disc, degrees; is the volume of covering soil, 10−4 m3.

According to the test data, the uniformity coefficient and stability coefficient of the volume of covered soil under different angles of the soil covering discs are calculated, as shown in Table 2.

**Table 2.** Coefficient of uniformity and stability of covering soil at different angles.

Note: 1 The factor column indicates the angle of the soil covering disc and the brackets indicate that the angles of the soil covering disc are different on the left and right sides; 2 the relevant dimensions of the covering soil section during the calculation of the standard deviation in mm.

According to the standard, the variation coefficient of uniformity should be less than 40% and the variation coefficient of stability should be less than 7.8%. When the coefficient of variation is less than the specified value, it indicates that the machine possesses good quality covering soil.

As shown in Table 2, when the angle of the soil covering disc is 30◦ and 60◦, respectively, the uniformity coefficient is relatively small (less than 40%), indicating that the

volume of covered soil is more consistent and that the quality is good. When the angle of the soil covering disc is 45◦, the uniformity coefficient is relatively large (more than 40%), indicating that the uniformity of the volume of covered soil is not good. When the difference between the angles of the discs on both sides is large (30◦ and 60◦), the standard deviation and variation coefficient of stability coefficient are the lowest (less than 7.8%), indicating that the greater the difference between the discs on the both sides, the easier it is to achieve differential overburden, which is beneficial under sloping fields.

#### **4. Conclusions**

This paper designs a soil covering device that can perform linkage and differential adjustments separately and the field experiment was completed to verify the performance of soil covering device. The results show that this device can obtain higher quality potato planting. However, due to the condition of test sloping fields being hard to obtain, in this paper, the flat fields were used to complete the field experiment, which needs verification in further research. Moreover, the Soil-bin Testing Platform for Sloping Field will be installed in the Northwest A&F University, which can be used to test the performance of device in sloping field.

In order to improve the mulching quality, a soil covering device that can realize linkage and differential adjustment was designed and field experiments were carried out to obtain the fitting equation of the volume of covered soil and the angle of the covering disc; the uniformity coefficient and stability coefficient are measured and calculated.

(1) The structure and principle of the soil covering device are put forward and the design calculation and analysis are focused on the four-bar adjustment mechanism and reasonable structural parameters are obtained. The linkage adjustment is carried out through the two-way screw and nut mechanism and the differential adjustment is carried out through the steering mechanism.

(2) By using data fitting, it was found that the volume of covered soil changes in a binomial manner with the angle of the discs. Field test results show the following: when the angles of discs on both sides are the same (30◦ or 60◦), the uniformity coefficient of volume of covering soil is lower than 14% and when the angles of the discs are different (respectively, 30◦ and 60◦), the stability coefficient of volume of covering soil is 4.1% (less than 7.8%), which indicates that the linkage and difference adjustment device can meet the quality requirements of the potato film side covering soil.

**Author Contributions:** Conceptualization, Z.F.; methodology, Z.Z. and Z.F.; validation, J.H.; investigation, Y.H. and Z.Z.; resources, Z.F.; writing—original draft preparation, C.W.; writing—review and editing, C.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Key research and development projects of Shaanxi Province (2019NY-169 and 2018ZDCXL-NY-03-06).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available upon requrest from the correspondence author at (fuzuoli@nwsuaf.edu.cn).

**Acknowledgments:** This work was financially supported by the Northwest A&F University. All supports and assistance are sincerely appreciated.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Article* **An Investigation of Pull-Out Force of Semi-Buried Lotus Roots after Hydraulic Scouring**

**Rong Zeng 1,2, Yitao Lin 1,2, Zhihua Wan 1,2, Ming Tu 1,2, Jun Jiao 1,2 and Guozhong Zhang 1,2,\***

**Abstract:** Hydraulic scouring is the most effective approach to harvest lotus roots, but its application is largely restricted by the low harvesting rate. After hydraulic scouring, some mature lotus roots are still partially buried in the soil. Therefore, it is highly necessary to develop an appropriate approach to harvest semi-buried lotus roots. In this work, we for the first time studied the interaction between semi-buried lotus roots and soil, as well as analyzing the pull-out process and the factors affecting the pull-out force of semi-buried lotus roots. Firstly, a simple testing platform was designed based on the virtual prototype technology and the tests on pull-out force were conducted on five lotus roots with similar shape but different weights, with the burial depth and pull-out speed as the experimental factors. The results revealed that the maximum pull-out force is not significantly affected by the pull-out speed, whereas it is significantly influenced by the burial depth and the surface area of lotus roots. The maximum pull-out force increased with increasing lotus root surface area and burial depth. In addition, the discrete element method was employed to simulate the pull-out process of lotus root at different pull-out speeds. The simulation results indicated that a higher pull-out speed would result in a greater pull-out force at the same displacement of the lotus root from the soil. Both experimental and simulation results revealed that soil adhesion contributes the most to the pull-out resistance. It was also observed that a slight loosening of semi-buried lotus roots could drastically reduce the pull-out force. These results suggest that some kind of mechanical structure or improvement of water flow can be applied to the existing lotus root harvester to reduce the adhesion between lotus roots and soil. Overall, our findings provide a novel direction for optimizing hydraulic harvesting machines of lotus roots.

**Keywords:** semi-buried lotus root; pull-out force; hydraulic scouring; soil adhesion

#### **1. Introduction**

Lotus root is one of the largest aquatic vegetables in China. According to statistics, the planting area of lotus root in China in 2014 was about 500,000 to 700,000 hectares [1]. With crispy and refreshing flavor and rich nutrition, lotus roots are popular among Chinese people. Due to the huge economic benefits, the number of lotus root growers and the planting area have been increasing year by year. Lotus roots are mostly planted in low-lying and rotten fields with deep silt [2]. The special planting environment and growth conditions result in great difficulties in mechanized harvesting. As a result, manual harvest is still the main harvesting method. With the rapid decrease in labor resources and increase in labor cost in rural China, the production cost of lotus roots has increased dramatically, which seriously restricts the scale development of the lotus root industry. Hence, some scholars have developed several simple harvesting machines for lotus roots based on the planting characteristics, and the mechanized harvesting of lotus roots has been partly realized.

**Citation:** Zeng, R.; Lin, Y.; Wan, Z.; Tu, M.; Jiao, J.; Zhang, G. An Investigation of Pull-Out Force of Semi-Buried Lotus Roots after Hydraulic Scouring. *Agriculture* **2021**, *11*, 706. https://doi.org/10.3390/ agriculture11080706

Academic Editor: José Pérez-Alonso

Received: 17 May 2021 Accepted: 26 July 2021 Published: 27 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Currently, shovel digging and hydraulic scouring are the two most commonly used methods for mechanized harvesting of lotus roots. Due to the fragile texture and slender shape of lotus roots, harvesting with shovel digging causes high damage and loss. The harvesting machine with hydraulic scouring uses a high-pressure water jet to wash the mud around lotus roots. When the lotus roots are free from the mud, they will float to the water's surface so that the operators can easily pick them up. The pressure of the water jet can be adjusted to reduce damage to the lotus roots. Therefore, hydraulic harvesting machines are more widely applied to the harvesting of lotus roots [3–5]. The shape of lotus root is similar to a drum. Sludge attached to the upper surface of lotus root can be washed away by the water jet, whereas the sludge attached to the lower surface is obstructed and cannot be washed away by the water jet. Due to high viscosity of soil, some lotus roots are still partially bound by sludge, and will be in a semi-buried state and thus cannot float up. These lotus roots need to be harvested manually. To develop a reasonable harvesting method, it is essential to analyze the interaction between the semi-buried lotus roots and soil, which may be best reflected by the pull-out force of semi-buried lotus roots. Aiming to improve the working performance of harvesting machines for crops from soil, research has been focused on the pull-out force of crops from soil. For example, Li et al. designed a device to measure the pull-out force of cotton stalk and analyzed the relationship between the pull-out force and cotton root diameter, pull-out direction and pull-out speed [6]. Similarly, Song et al. studied the factors affecting the pull-out resistance of cotton stalk and proposed that the diameter of the root of cotton stalk and the soil properties had significant influence on the pull-out resistance [7]. Moreover, Xie et al. simulated the pull-out force of leaf vegetables from soil by EDEM software and determined the optimal pull-out force [8]. Based on the simulation results, they also proposed the optimal working parameters of the harvesting machine. In addition, according to the grip strength, Seongin et al. developed a harvesting robot for lettuce [9].

Understanding the pull-out force of semi-buried lotus roots from soil may help to develop effective methods for the subsequent harvesting of semi-buried lotus roots. At present, there has been no research on the pull-out force of lotus roots. Hence, in order to explore the factors affecting the pull-out force and analyze the pull-out process of lotus roots, this paper constructed a simple platform to test the pull-out force of semi-buried lotus roots. The pull-out force of different lotus roots was investigated at different pull-out speeds and burial depths, and the effect of these factors on the pull-out force was evaluated. In addition, the pull-out process of a testing lotus root was simulated by the discrete element method via the EDEM software. Experimental and simulation analysis results suggest the importance of finding an effective approach to reduce the pull-out resistance, so as to improve the harvesting machine with hydraulic scouring for lotus roots.

#### **2. Materials and Methods**

#### *2.1. Analysis of Pull-Out Force of Semi-Buried Lotus Roots after Hydraulic Scouring*

After washing by water jet, the lotus roots may naturally float out of the water surface. Taking one lotus root as an example, the floating process is demonstrated in Figure 1. Firstly, the lotus root is buried in a certain depth under the mud (a). Then, the hydraulic harvesting machine starts to wash the mud layer (b), and the mud surrounding the lotus root is washed away (c). Finally, the lotus root is separated from the mud layer and floats out of the water surface (d). Since the burial depth of each lotus root in mud layer is not consistent, the distance between the nozzle of the water jet and each lotus root is uncertain, resulting in the failure of deeper lotus roots to float out. In addition, the water jet of hydraulic scouring will generate a downward impact force. If the mud layer is soft, the impact force may press the lotus roots into a deeper layer, preventing the lotus roots from floating. These lotus roots will be semi-buried in the mud and require further measures to be harvested.

**Figure 1.** Floating process of a lotus root. (**a**) Initial state of the lotus root; (**b**) Start of floating; (**c**) Continuous floating; (**d**) Total floating.

In this section, the forces acting on the semi-buried lotus roots are analyzed. When the semi-buried lotus root tends to move upward, it is under the buoyancy of water, gravity, the soil adhesion force, the normal force of the mud layer and the pull-out force. Taking the cross-section center of the lotus root as the origin to establish the coordinates, the free-body diagram of the lotus root is depicted in Figure 2.

**Figure 2.** Free-body diagram of the semi-buried lotus root.

Figure 2 shows the vertical force acting on the lotus root when it is pulled out, and the resultant force can be expressed as:

$$
\sum \stackrel{\rightarrow}{F} = \stackrel{\rightarrow}{F}\_f + \stackrel{\rightarrow}{F}\_p + \stackrel{\rightarrow}{F}\_N + \stackrel{\rightarrow}{G} + \stackrel{\rightarrow}{F}\_{a\prime} \tag{1}
$$

where ∑*F* represents the resultant force; *Ff* and *Fp* denote the buoyancy of water and pull-out force, respectively; *FN* is the normal force of the mud layer; *G* is the gravitational force; and *Fa* describes the equivalent soil adhesion force in the vertical direction. When the pull-out force is lower than the vertical resistance, the lotus root will stay static in the mud layer, and thus, ∑*F* is 0. Therefore, the equilibrium equation is given by:

$$0 = \stackrel{\rightarrow}{F}\_p + \stackrel{\rightarrow}{F}\_{pz\_{\prime}} \tag{2}$$

and the vertical resistance is:

$$
\stackrel{\rightarrow}{F}\_{pz} = \stackrel{\rightarrow}{F}\_f + \stackrel{\rightarrow}{F}\_N + \stackrel{\rightarrow}{F}\_a + \stackrel{\rightarrow}{G}\_\prime \tag{3}
$$

where *Fpz* is the pull-out resistance. If the pull-out resistance is overcome, the lotus root will start to float.

#### *2.2. Pull-Out Force Test*

#### 2.2.1. Testing Platform

A testing platform for the measurement of pull-out force of semi-buried lotus roots is designed in this section. Figure 3 depicts the working principle.

**Figure 3.** Sketch map of the testing platform.

A single lotus root is buried in the mud layer in the collection box. A clamping claw is employed to grasp the lotus root. The wire rope on the clamping claw is driven upward by the motor to provide a pull-out force. A S-type tension sensor is connected to the wire rope to measure the pull-out force. The output voltage of the tension sensor is acquired by a USB type data acquisition (DAQ) card and converted to pull-out force by a computer with the LabVIEW software. To establish the relationship between the output voltage and the pull-out force, the sensor was calibrated first, and the linear relationship was obtained as:

$$F\_p = K\mathcal{U} \tag{4}$$

where *K* denotes the linear coefficient, which was calibrated to be 872.15, and *U* is the output voltage of the S-type tension sensor.

#### 2.2.2. Experimental Materials

Lotus roots generally present a drum shape, and there is certain regularity in the diameter distribution. In this experiment, five single lotus roots of the No. 5 Hubei variety were selected as the test objects according to the diameter distribution. The testing lotus roots showed similar shapes but different sizes without obvious damage in appearance. The average length of the five single lotus roots was 22.58 cm, and the minimum and

maximum diameters were 55.1 mm and 94.5 mm, respectively. First, their weights were measured and the values are presented in Table 1.

**Table 1.** Weights of the testing lotus roots.

The mud layer in the collection box was set as 30 cm in depth, and the water layer was set as 10 cm in depth. Soil samples were collected from a lotus root field located at Binhu fishing ground in Dongsheng Town, Shishou City, Hubei Province, China. The soil hardness and moisture content were measured and presented in Tables 2 and 3, respectively. Since the distribution of soil layer was not uniform, the soil hardness was measured five times at different positions to obtain the mean value (Table 2). Table 3 shows the moisture content in different layers of soil. As the surface layer is mainly solid with silt, it has higher fluidity and moisture content. With increasing soil depth, the soil hardness increases, accompanied by a decrease in moisture content.

**Table 2.** Soil hardness.

**Table 3.** Soil moisture content (%).

#### 2.2.3. Experimental Methods

Since lotus roots are usually buried in a disordered state, different burial depths of lotus roots will require different pull-out force. Therefore, the burial depth was taken as one experimental factor. During the test, the burial depth was defined as shown in Figure 4. Before the test, the testing lotus root was slowly pressed into the mud layer, and the distance between the upmost point of the lotus root and the water layer surface was measured by a tape measure, which was recorded as the burial depth. Then, the mud attached to the upper surface was removed to form a semi-buried state. As the soil hardness increases with increasing soil depth, it is difficult to bury lotus roots to depths greater than 20 cm. Therefore, the burial depths of 15 cm and 20 cm were chosen to measure the pull-out force of the semi-buried lotus roots.

**Figure 4.** Definition of the burial depth.