Validating the Model of a No-Till Coulter Assembly Equipped with a Magnetorheological Damping System

: Variability in soil conditions has a signiﬁcant inﬂuence on the performance of a no-till seeder in terms of an inconsistency in the depth of seeding. This occurs due to the inappropriate dynamic responses of the coulter to the variable soil conditions. In this work, the dynamics of a coulter assembly, designed with a magnetorheological (MR) damping system, were simulated, in terms of vertical movement and ground impact. The developed model used measured inputs from previously performed experiments, i.e., surface proﬁles and vertical forces. Subsequently, the actual coulter was reassembled with an MR damping system. Multiple sensors were attached to the developed coulter in order to capture its motion behavior together with the proﬁles, which were followed by the packer wheel. With the aim to validate the correctness of the simulation model, the simulation outputs, i.e., pitch angles and damper forces, were compared to the measured ones. The comparison was based on the root-mean-squared error (RMSE) in percentage, the root-mean-squared deviation (RMSD), and the correlation coe ﬃ cient. The average value of the RMSE for the pitch angle, for all currents applied on the MR damper, was below 10% and 8% for the speeds of 10 km h − 1 and 12 km h − 1 , respectively. For the damper force, these ﬁgures were 15% and 13%. The RMSD was below 0.5 deg and 1.3 N for the pitch angle and the damper force, respectively. The correlation coe ﬃ cient for all datasets was above 0.95 and 0.7 for the pitch angle and the damper force, respectively. Since the damper force indicated a comparatively lower correlation in the time domain, its frequency domain and coherence were investigated. The coherence value was above 0.9 for all datasets.


Introduction
In no-till seeding, a consistency in seeding depth plays a very important role in achieving proper seed germination and seedling emergence since both ultimately influence the crop growth, and in general the produced yield [1,2]. The inconsistency in the seeding depth mainly occurs due to the inappropriate dynamic behavior of the furrow components. Harsh soil conditions, such as rough soil surface, residue effects, variable soil density, etc., have a considerable effect on the dynamic performance of the seeding machines. The inability to regulate the vertical dynamics of the coulter assembly, in terms of reaction forces and displacements [3], results in variability of their dynamic response. The forces that are resulting from the interaction between the coulter tine and the soil could characterize the motion behavior of the assembly [4,5]. However, this depends on the components of the assembly that are excited by the soil reaction forces. In addition, the assembly might include the extra packing component, where the developed impact forces also need to be considered [6]. By considering the nature of the forces when regulating the vertical position of the coulter, its performance could be deviation (RMSD), and correlation coefficient. Furthermore, a correlation between the frequency content of the simulated and measured damper forces, which would based on the coherence value, should be carried out to indicate if the simulation model can produce a similar frequency range as the measured one.

Instrumentation and Implementation of the MR Damper
A metal frame with 2 wheels that could be attached on the 3-point hitch of the tractor was constructed to carry the coulter assembly (AMAZONEN-Werke H. Dreyer GmbH & Co. KG, Hasbergen, Germany) ( Figure 1). The coulter was fixed to a square rod using rubber rollers. The square rod had one degree of freedom and could rotate while a hydraulic cylinder was used that applied downward pressure to the coulter with a value of 6.89 MPa. To absorb the impact of the reaction forces, an RD-8040-1 MR damper (LORD, Baltimore, MD, USA) was added to the assembly using an extra shank and bearing system at the joint. During the field tests, a variable damping ratio was provided using a controlling device (Wonder box, LORD, Baltimore, MD, USA) by regulating the input voltage, which was applied on the magnetic field in the flow path of the MR fluid. The device regulated the output current, which was linearly proportional to the input voltage of the MR damper. More detailed information about the performance of the MR damper can be found in previous work by the authors of reference [27].
Multiple sensors were employed in order to acquire the field profiles that were followed by the packing wheel and the developed assembly dynamics (Figure 1). The in-field absolute geo-referenced position of the mainframe and the vertical movement of the developed coulter were obtained using an SPS930 total station (Trimble, Sunnyvale, CA, USA) and a DT50 range detector (SICK AG, Waldkirch, Germany), respectively. The total station detected the geo-referenced position of the frame by following an MT900 active prism powered by the tractor (Trimble, Sunnyvale, CA, USA), while the range finder recorded the vertical movement from a reflectance plate attached on the coulter assembly. The tilting information of the developed frame and the assembly were acquired using two VN-100 inertial measurement units (IMUs) (VectorNav, Dallas, TX, USA). The reaction and damper forces were extracted from the data recorded by 6 linear 350 Ohm DY41-1.5 strain gauges (HBM GmbH, Darmstadt, Germany) fixed on the coulter. operations. The comparison will be based on criteria such as root-mean-squared error (RMSE), rootmean-squared deviation (RMSD), and correlation coefficient. Furthermore, a correlation between the frequency content of the simulated and measured damper forces, which would based on the coherence value, should be carried out to indicate if the simulation model can produce a similar frequency range as the measured one.

Instrumentation and Implementation of the MR Damper
A metal frame with 2 wheels that could be attached on the 3-point hitch of the tractor was constructed to carry the coulter assembly (AMAZONEN-Werke H. Dreyer GmbH & Co. KG, Hasbergen, Germany) ( Figure 1). The coulter was fixed to a square rod using rubber rollers. The square rod had one degree of freedom and could rotate while a hydraulic cylinder was used that applied downward pressure to the coulter with a value of 6.89 MPa. To absorb the impact of the reaction forces, an RD-8040-1 MR damper (LORD, Baltimore, MD, USA) was added to the assembly using an extra shank and bearing system at the joint. During the field tests, a variable damping ratio was provided using a controlling device (Wonder box, LORD, Baltimore, MD, USA) by regulating the input voltage, which was applied on the magnetic field in the flow path of the MR fluid. The device regulated the output current, which was linearly proportional to the input voltage of the MR damper. More detailed information about the performance of the MR damper can be found in previous work by the authors of reference [27].
Multiple sensors were employed in order to acquire the field profiles that were followed by the packing wheel and the developed assembly dynamics (Figure 1). The in-field absolute geo-referenced position of the mainframe and the vertical movement of the developed coulter were obtained using an SPS930 total station (Trimble, Sunnyvale, CA, USA) and a DT50 range detector (SICK AG, Waldkirch, Germany), respectively. The total station detected the geo-referenced position of the frame by following an MT900 active prism powered by the tractor (Trimble, Sunnyvale, CA, USA), while the range finder recorded the vertical movement from a reflectance plate attached on the coulter assembly. The tilting information of the developed frame and the assembly were acquired using two VN-100 inertial measurement units (IMUs) (VectorNav, Dallas, TX, USA). The reaction and damper forces were extracted from the data recorded by 6 linear 350 Ohm DY41-1.5 strain gauges (HBM GmbH, Darmstadt, Germany) fixed on the coulter.

Experiments
The experiments to measure the coulter dynamics and the field profiles were performed at the research farm of the University of Hohenheim located in 48 • 43 27.33" N, 9 • 11 07.69" E. Before performing seeding operation, experiments were conducted to determine the soil properties [28]. When capturing the dynamic performances, the MR damper was excited with different current input levels (0-1 A with an increment of 0.2 A). This was carried out for 2 driving speeds, i.e., 10 km h −1 and 12 km h −1 .
The spatial sampling frequency of the attached sensors for the speed of 10 km h −1 has been reported in reference [26]. The spatial sampling frequency for the speed of 12 km h −1 of strain gauges, IMUs and range finder was equal to 11 mm, 66 mm, and 53 mm, respectively. Since the various sampling rates of the sensors resulted in non-concurrent data, a linear interpolation approach was employed to obtain values at concurrent instances.

Dynamics of the Developed Assembly
The assembly consisted of the coulter, the wheel shank and the packer wheel. The coulter assembly, together with the joint part, was modelled as a forced pendulum. The pendulum was also considered as damped since the point of the packer wheel touching the ground comprised of a material made out of rubber. The motion equation of the coulter was developed for 2 cases. The first case was for the original coulter assembly, where the packer wheel was considered as a passively controlled system (see Figure 2a). Detailed information can be found in previous work by the authors in reference [24]. In the second case, the semi-active MR damper, as a controlled mass-spring-damper system, was designed to be located between the wheel shank and the packer wheel, having an extra shank ( Figure 2b). The extra shank was assumed to have freedom of rotation.

Experiments
The experiments to measure the coulter dynamics and the field profiles were performed at the research farm of the University of Hohenheim located in 48°43′27.33″ N, 9°11′07.69″ E. Before performing seeding operation, experiments were conducted to determine the soil properties [28]. When capturing the dynamic performances, the MR damper was excited with different current input levels (0-1 A with an increment of 0.2 A). This was carried out for 2 driving speeds, i.e., 10 km h −1 and 12 km h −1 .
The spatial sampling frequency of the attached sensors for the speed of 10 km h −1 has been reported in reference [26]. The spatial sampling frequency for the speed of 12 km h −1 of strain gauges, IMUs and range finder was equal to 11 mm, 66 mm, and 53 mm, respectively. Since the various sampling rates of the sensors resulted in non-concurrent data, a linear interpolation approach was employed to obtain values at concurrent instances.

Dynamics of the Developed Assembly
The assembly consisted of the coulter, the wheel shank and the packer wheel. The coulter assembly, together with the joint part, was modelled as a forced pendulum. The pendulum was also considered as damped since the point of the packer wheel touching the ground comprised of a material made out of rubber. The motion equation of the coulter was developed for 2 cases. The first case was for the original coulter assembly, where the packer wheel was considered as a passively controlled system (see Figure 2a). Detailed information can be found in previous work by the authors in reference [24]. In the second case, the semi-active MR damper, as a controlled mass-spring-damper system, was designed to be located between the wheel shank and the packer wheel, having an extra shank ( Figure 2b). The extra shank was assumed to have freedom of rotation. The weight of the wheel and the coulter were considered as an unsprung mass [kg], and a sprung mass [kg], respectively. The reaction force [N] acting on the wheel was calculated by the deformation and the damping forces of the rubber wheel. The stiffness and the damping forces were introduced by the vertical shift and velocity of the unsprung mass . This was formulated as follows: introduced by the vertical shift and velocity of the unsprung mass m w . This was formulated as follows: where x w is the shift in the position of the wheel in the vertical plane, in [m]; x p is the field profile that the wheel followed; K w is the stiffness coefficient of the rubber tyre [N m −1 ] and C w is the damping coefficient [N s m −1 ]. More detailed information on defining the stiffness and the damping coefficient of the wheel-tire can be found in reference [24]. In this model, the horizontal force F d was not considered in the modelling due to its minor effect on the vertical motion behavior of the assembly. The vertical movement of the coulter x c was also expressed with the combination of the pitch angle θ and the distance between the rotating rod and the impact force affecting point L a . Since the pitch angle of the assembly relative to the main frame was small (θ < 5 • ), the angle approximation (sin θ = θ and cos θ = 1) was used to simplify the equation of motion. The damping force F mr [N] resulting from the MR damping system was the affecting force to the assembly since the MR damper was attached on the fixed shank from one side. Therefore, the profile impact force F sp was substituted by F mr . By taking into account the moment of all the effecting forces and the moment of the constant force, which was applied to the square rod M f , the equations of motion behavior of the assembly, including the semi-active MR damping system, were the following: where

Bouc-Wen Model
Many different hysteresis models could be applied for the MR damper model to address its hysteresis and nonlinear behavior [29,30]. The Bouc-Wen model (Figure 3) has been proven to be outperforming over the other models [24] when investigating the nonlinear behavior of the MR damper. In the developed model (Figure 2b), F c is the control force of the MR damper resulting from the pre-yield stress of the damper fluid.
where is the shift in the position of the wheel in the vertical plane, in [m]; is the field profile that the wheel followed; is the stiffness coefficient of the rubber tyre [N m −1 ] and is the damping coefficient [N s m −1 ]. More detailed information on defining the stiffness and the damping coefficient of the wheel-tire can be found in reference [24].
In this model, the horizontal force was not considered in the modelling due to its minor effect on the vertical motion behavior of the assembly. The vertical movement of the coulter was also expressed with the combination of the pitch angle and the distance between the rotating rod and the impact force affecting point . Since the pitch angle of the assembly relative to the main frame was small ( < 5°), the angle approximation (sin = and cos = 1) was used to simplify the equation of motion. The damping force [N] resulting from the MR damping system was the affecting force to the assembly since the MR damper was attached on the fixed shank from one side. Therefore, the profile impact force was substituted by . By taking into account the moment of all the effecting forces and the moment of the constant force, which was applied to the square rod , the equations of motion behavior of the assembly, including the semi-active MR damping system, were the following:

Bouc-Wen Model
Many different hysteresis models could be applied for the MR damper model to address its hysteresis and nonlinear behavior [29,30]. The Bouc-Wen model (Figure 3) has been proven to be outperforming over the other models [24] when investigating the nonlinear behavior of the MR damper. In the developed model (Figure 2b), is the control force of the MR damper resulting from the pre-yield stress of the damper fluid.
where the evolutionary variable is the following: where , , and are the time-dependent parameters that represent the control of the linear behavior; is a constant related to the accuracy of the model; is a response-shape characteristic The resulting force from the Bouc-Wen model was expressed by the following formula: where the evolutionary variable z is the following: . z = γz|x c ||z| n−1 − βx c |z| n + Ax c (5) where γ, β, and A are the time-dependent parameters that represent the control of the linear behavior; α is a constant related to the accuracy of the model; n is a response-shape characteristic parameter of the model; K s (u) and C s (u) are the voltage-dependent coefficients representing the stiffness and viscous damping, respectively. The model parameters K s (u), C s (u), and α can be defined as a linear function of the applied voltage and are obtained by the following formulas [31]: where K s1 , K s2 , C s1 , C s2 , α 1 , α 2 are the necessary parameters to characterize the rheological and magnetic behaviour of the MR damper. To define the model parameters of the Bouc-Wen model, its response was simulated using the Equations (4)- (8) with the different applied current values on the coil. The parameters K s1 , K s2 , C s1 , C s2 , α 1 , α 2 , were determined by leveraging an empirical formulation that used the experimental data [32,33].

Evaluation of the Model Performance
All the Equations from (1) to (8) for the coulter dynamics and the hysteresis model of the MR damper were applied and simulated in MATLAB and Simulink. The acquired field surface profiles ( Figure 4) and the vertical forces ( Figure 5) from the dynamic measurements at the speeds of 10 km h −1 and 12 km h −1 , when the MR damper was excited with 0.3 A, 0.5 A, and 0.7 A, were used as inputs to the simulation model. The acquired strains were utilized to determine the vertical forces and the damper forces that were acting on the coulter tine and the fixed shank. To define the field profiles, the geo-referenced coordinates of the active prism were transformed to the ground touching point of the wheel. The detailed formulas on calculating all the forces and transformations, in terms of rotations and translations for the profile determination, are given in reference [27].
The outputs of the model were the pitch angle of the assembly movement and the damper force. All parameters and constants of the Bouc-Wen model, as well as the coulter model, can be found in reference [24]. The real-life test of the assembly performance was performed with different current input values loaded on the MR damper (0.0, 0.1, 0.3, 0.5, 0.7, 0.9, 1.0 A). Two different traveling speeds of 10 km h −1 and 12 km h −1 were set for all performed experiments. The field tests with the current values 0.3 A, 0.5 A, and 0.7 A had the best performance, in terms of coulter dynamics [27], thus they were chosen to be further analyzed in the present work.
In order to validate the correctness of the simulation model, its outputs such as the simulated pitch angles of the assembly vertical motion and the damper forces were compared with the measured ones. The comparison analyses were carried out by calculating the RMSE in percentage, absolute RMSD, and the correlation coefficient between the simulated and the measured data. The RMSE estimated the relative error in percentage between the measured and the simulated datasets. Taking into account the effect of the neglected lateral forces on the tire dynamics, the value of 20% for the RMSE in percentage was considered as a threshold that should not be exceeded [34,35]. The formula of the RMSE in percentage that was used is the following: (z s (t)) 2 /N; z m (t) and z s (t) are the measured and the simulated data, respectively, within the time t that was needed to travel the specified distance. N is the total number of samples. Subsequently, the average difference between the measured and the simulated outputs was evaluated by determining the RMSD value. This was defined by the following expression: In addition to the error estimation, the fitting accuracy, in terms of similarity relationship, of the simulated data to the measured data were also addressed with a correlation coefficient. The values higher than 0.6 could be considered as practically acceptable based on the same reasons as above. The correlation coefficient ρ was based on the following equation: where µ z m and σ z m are the mean value and the standard deviation of z m (t), respectively, and µ z s and σ z s are the mean value and the standard deviation of z s (t), respectively.      Due to the wide range of signal wavelengths of the damper forces, the correlation of the measured and the simulated forces in the frequency domain was carried out based on the coherence analysis. In this case, the relationship between the measured and the simulated forces can be properly evaluated by analyzing the coherence value that was resulted from the cross-spectrum of the wavelengths in the frequency domain. The power spectral density (PSD) of the measured and the simulated datasets were computed to determine the coherence based on the following equation: where P F m and P F s are the resulted PSD datasets for the measured and the simulated forces, respectively, and P F m ,F s is the cross-spectrum density of these 2 data sets. A coherence value equal to 1 indicates the maximum possible agreement while a value of 0 denotes no correlation [36]. Figure 6 presents both measured and simulated pitch angles of the assembly movement for the driving speeds of 10 km h −1 and 12 km h −1 . Three different cases are highlighted, i.e., when the MR damper was supplied with 0.3 A, 0.5 A, and 0.7 A.

Results and Discussion
For the same current input levels applied on the MR damper, i.e., 0.3 A, 0.5 A, and 0.7 A, both measured and simulated damper forces for the driving speeds of 10 km h −1 and 12 km h −1 are shown in Figure 7.
The results of the correlation, in terms of the RMSE, RMSD and correlation coefficient, which were calculated using Equations (9)-(11), respectively, are given in Table 1.
The analyses of the RMSE percentage and RMSD showed good agreement between the simulated and the measured outputs. The average values of the RMSE and the RMSD for the pitch angle at the speed of 10 km h −1 were equal to around 10% and 0.4 deg, respectively. These figures were below 8% and around 0.5 deg for the speed of 12 km h −1 . Based on the above a high correlation coefficient was determined with an average value of 0.96 and 0.97 for the traveling speed of 10 km h −1 and 12 km h −1 , respectively. The correlation analyses of the damper forces indicated slightly different results, in terms of a higher RMSE percentage and RMSD value, and smaller correlation coefficient compared to the ones resulted from the pitch angles. The average value of the RMSE percentage was 15% and 13% for the operation speed of 10 km h −1 and 12 km h −1 , respectively. The average of the RMSD, indicated a value of 1.2 N at the speed of 10 km h −1 and 1.3 N at the speed of 12 km h −1 . Consequently, the higher error resulted in a lower correlation coefficient with an average value of around 0.7 and 0.69 for these speeds.  The analyses of the RMSE percentage and RMSD showed good agreement between the simulated and the measured outputs. The average values of the RMSE and the RMSD for the pitch angle at the speed of 10 km h −1 were equal to around 10% and 0.4 deg, respectively. These figures were below 8% and around 0.5 deg for the speed of 12 km h −1 . Based on the above a high correlation coefficient was determined with an average value of 0.96 and 0.97 for the traveling speed of 10 km h −1 and 12 km h −1 , respectively. The correlation analyses of the damper forces indicated slightly Since the correlation analyses in the time domain indicated lower agreement between the measured and the simulated damper forces, a correlation was carried out based on the coherence analyses in the frequency domain (Equation (12)). Figures 8 and 9 present the frequency content of the measured and the simulated damper forces, for the traveling speed of 10 km h −1 and 12 km h −1 , respectively, together with the corresponding coherence, with current inputs of 0.3 A, 0.5 A, and 0.7 A applied on the MR damper. Since the correlation analyses in the time domain indicated lower agreement between the measured and the simulated damper forces, a correlation was carried out based on the coherence analyses in the frequency domain (Equation (12)). Figures 8 and 9 present the frequency content of the measured and the simulated damper forces, for the traveling speed of 10 km h −1 and 12 km h −1 , respectively, together with the corresponding coherence, with current inputs of 0.3 A, 0.5 A, and 0.7 A applied on the MR damper.
Despite the fact that there was a wide range of wavelengths of the damper forces in the time domain, which resulted in lower correlation, the analyses of their frequency content showed good agreement according to the resulted coherence values. The resulting coherence values were above 0.96 at the speed of 10 km h −1 for all the current inputs applied on the MR damper. This figure was slightly lower for the speed of 12 km h −1 with a coherence value of 0.91.
The analyses in the time domain denoted that all the RMSE percentage, RMSD, and correlation coefficient values were within an acceptable range in terms of the predefined thresholds. In addition, the pattern as the correlation coefficient was high when the RMSE indicated a low error, which confirmed the correctness of the performed correlation. Considering all the above-given results in both time and frequency domain, it could be stated that the simulation model that described the dynamics of the assembly equipped with the MR damper was capable of producing output values within acceptable deviation ranges compared to the observed ones. Appl. Sci. 2019, 9, x FOR PEER REVIEW 11 of 14

Conclusions
A mathematical model of the vertical motion dynamics of the coulter assembly with the MR damper was developed. The necessary instrumentation to measure the vertical dynamics of the assembly and the field profiles were presented. The dynamic response, in terms of the damper forces and the pitch angles of the assembly vertical motion, were simulated using the measured surface profiles and vertical forces. The performance of the mathematical model for two different traveling speeds of the seeder (10 km h −1 and 12 km h −1 ) was validated by comparing the RMSE percentage, the RMSD, and the correlation coefficient of the simulated outputs and the measured ones. The comparison analyses depicted a good agreement, with an RMSE percentage below 15% and a correlation coefficient above 0.7, between the simulated and measured datasets for all current inputs on the MR damper. The average difference, in terms of RMSD, also indicated acceptable values for both pitch angle and damper force. This analysis confirmed that the mathematical model is capable of producing similar results compared to the observed ones. Thus, the model could be used as a basis for finding the optimal settings for the MR damper but also to further develop and optimize the entire seeding assembly under various field and technical conditions. Future steps should include the implementation of different control techniques, such as the real-time adaptive PID and fuzzy hybrid controller [37], in order to optimize the dynamic performance of the MR damper and thus improve the assembly performance under different soil conditions.

Funding:
The authors would like to thank GA nr 213-2723/001-001-EM Action2 TIMUR project that financially covered the study visit for the data acquisition.