- Email: [email protected]
3D FE—DEM Simulation of a Thickening Growth Model for Japanese Radish Root
3D FE–DEM Simulation of a Thickening Growth Model for Japanese Radish Root Yasuhito Fujita ∗ Hiroshi Nakashima ∗∗ Hiroaki Tanaka ∗∗∗ Juro Miyasaka ∗∗∗∗ Katsuaki Ohdoi † Hiroshi Shimizu ‡ Graduate School of Agriculture, Kyoto University, Kyoto 606-8502, Japan (e-mail: [email protected]). ∗∗ Graduate School of Agriculture, Kyoto University, Kyoto 606-8502, Japan (e-mail: [email protected]) ∗∗∗ National Agricultural Research Center for Western Region, Kagawa 765-0053, Japan (e-mail: [email protected]) ∗∗∗∗ Graduate School of Agriculture, Kyoto University, Kyoto 606-8502, Japan (e-mail: [email protected]) † Graduate School of Agriculture, Kyoto University, Kyoto 606-8502, Japan (e-mail: [email protected]) ‡ Graduate School of Agriculture, Kyoto University, Kyoto 606-8502, Japan (e-mail: [email protected]) ∗
Abstract: The mechanical impedance of soil is well known to influence the shape formation of storage roots of root vegetables strongly. Contact interaction between the root and soil become important in the growth and shape formation of the storage root. This study was aimed at construction of a numerical model for thickening growth of radish tap roots based on the root–soil contact interaction. However in 2D simulation, the height of the soil Discrete Elements increased greatly compared with experimental observations. To solve this problem, 3D simulation was extended and improvements in the shape representation of soil Discrete Elements were made. Soil Discrete Elements with non-spherical shape were generated such that several spherical Discrete Elements were joined together to create a model particle. In this report, the results of preliminary simulation obtained for a small scale are discussed. For the 3D simulation, the shape of root Finite Elements was reproduced based on growth velocities. No distortion occurred in root Finite Elements. Additional examinations of the mechanical properties of radish roots and parameter settings for root Finite Elements are needed. Furthermore, results confirmed that the extension of FE–DEM for the root–soil system in 3D was able to decrease the final height of surface soil Discrete Elements from 18 mm to approximately 14 mm. Clumped Discrete Elements decrease the height to about 13 mm. Keywords: Radish, Root, Thickening growth, 3D FE–DEM, Non-spherical discrete elements, Plants, Numerical simulation 1. INTRODUCTION It is generally recognized that the growth of plant roots can be affected adversely by high mechanical impedance imposed by their ambient soil environment. Generally, highly compacted soil interrupts root growth. Especially for root vegetables, storage roots must move the surrounding soil particles to make space as they grow in a soil medium. Consequently, the hardness, density, and structure of the soil strongly influence the development and shape formation of the storage root. Radish roots used to be evaluated according to their shape and size from the view of efficiency in transport and processing. This might also affect competitiveness on the market. Usually, straight and cylindrical shapes are desired. Although the shape patterns of storage roots of radish plants are basically programmed in their genetic DNA sequence, physical soil conditions alter these traits dramatically through vegetative growth. For example, soil layers that are heavily
compacted by wheels and tracks of heavy agricultural machines and continuous tillage operations at the same depth would produce irregular shapes of storage roots in Japanese Radish. Therefore, root–soil contact interaction plays an important role in the radish root shape formation. This study specifically examines the effects of the mechanical impedance in the radish root shape formation in terms of the contact problem. The study was designed to develop a thickening growth model of a radish tap root based on the root–soil contact interaction. As for modeling of root growth, with the increasing speed of computers, numerical methods such as Discrete Element Method (DEM) and Finite Element Method (FEM), are gradually becoming more effective techniques for the analysis of granular materials and contact problems between soil and agricultural machine devices (Horner et al. 2001; Nakashima et al. 2009). Based on the idea of cone penetration, when the root elongates in soil, development of
Fig. 1. Results of the 2D simulation of thickening growth of radish root (left, initial growth stage; right, final growth stage, Ti (day), ti (s)). 160 155
Height(mm)
a stress field in the soil medium has been analyzed using FEM (Faure 1994). In this study, root deformation was not considered. Similarly, the effects of the radial division and enlargement in the meristematic zone on the root tip and the deformation and stress field in the vicinity of the root tip have been analyzed using FEM (Richards and Greacen 1986). Recently, modeling of the root elongation growth based on the dynamic contact interaction between root and ambient soil environment using DEM was proposed (Nakashima et al. 2008). Furthermore, in a previous study, the basis of the thickening growth model of radish root based on the root–soil contact interaction was constructed in 2D (Fujita 2010). In this developed model, the DEM and the FEM were applied, respectively, to the modeled soil and root. Furthermore, root–soil contact interaction was analyzed using the Finite Element – Discrete Element Method (FE–DEM). The calculated shape of root Finite Elements mostly agreed with the experimental results. However, the behavior of the soil Discrete Elements caused larger displacements. Particularly, the height measured at the surface of soil Discrete Elements was considerably greater. Figure 1 shows numerical results in 2D simulation obtained at the initial and final growth stages. Furthermore, Fig. 2 shows changes of the height in soil Discrete Elements with time. In Fig. 1 and Fig. 2, T (day) is the time elapsed in the experiment and t(s) is that which elapsed in the simulation. Because of the extremely long calculation time of DEM simulation, a virtual time scale was introduced. The correspondence between the time scale in the simulation and in experiment was explained in previous reports (Fujita 2010). The root Finite Elements thickened laterally; soil Discrete Elements were largely deformed. Also, the height of soil Discrete Elements increased from 130 mm to approximately 155 mm. Therefore, the soil Discrete Elements had almost no compressibility, which was against the observation that the soil used in the experiment was compacted. Furthermore, in the previous model, the elongation of the hypocotyls could not be simulated apparently. To solve this unnatural increase of soil Discrete Elements, extension of FE–DEM to model root–soil systems in 3D was done and improvements in the shape of soil Discrete Elements were made.
150 145 140 135 130 125 120
Soil Height 15
20
25
30
35
time(day) Fig. 2. Height of the soil Discrete Elements shown versus time.
The main objective of this study is to construct a 3D thickening growth model of radish root by application of 3D FE–DEM.
2. MATERIALS AND METHODS
Fig. 3. Setting growth rule for root Finite Elements in 3D.
In this section, the 3D thickening growth model is described. Then, numerical methods used in the previous 2D model, FEM, DEM, and FE–DEM were extended to a 3D model. Furthermore, improvement in the shape of Discrete Elements was made. As in the previous 2D analysis, to reproduce the growth of root Finite Elements, a simplified assumption for mechanism of the thickening growth of radish root was made under which root Finite Elements would obey the experimental growth rule of “Growth Velocity” data obtained by virtual marking method in this radish experiment. The radish outline was obtained by projection onto a plane. Furthermore, the root shape was assumed to be axisymmetric. The root shape was expressed as the radius and the distance from the proximal end. Furthermore, in the 3D case, the radish growth was
assumed to be the same for all radial directions (Fig. 3). Figure 3 shows the assumed direction and magnitude of growth in the radish root. 2.1 Root Model A radish root was modeled using FEM. When the storage root grows radially, the root exerts a force to displace the soil particles; it is simultaneously reacted upon by the soil particles. The deformation and development of the stress fields within the root were calculated using FEM. The basic shape of the radish root was assumed to be axisymmetric and a quarter of the target domain was modeled. The root Finite Elements consist of standard hexahedral eight-node isoparametric elements (Fig. 4). For
or
Fig. 4. Configuration of the 3D root shape (left, full scale; right, quarter of the target domain) the material properties of the radish root, a linearly elastic model was assumed. In this model, the elastic strains {!} are related to the stresses {σ} by Hook’s law.
The thickening growth of root Finite Elements was formulated based on the prescribed “Growth Velocity” data. The displacement for all FE nodes at each time step can be expressed as Xi (t + ∆t) = Xi (t) + Vi (t)∆t, (1) where Xi (t) represents the position vector of the ith FE node at time t, mm. Also, Vi (t) is the growth velocity vector of ith FE node at time t, mms−1 . This Vi (t) was extracted from experimental results obtained using a radish. Details of the root growth experiment using the radish and the calculation method of Vi (t) were described in a previous report (Fujita 2010). This displacement is not the imposed displacement of Finite Elements, but the growth of root Finite Elements, which indirectly expresses the sum of biological activities in radish root, and which did not cause the strain and stress in the root Finite Elements. 2.2 Soil Model Soil was modeled using DEM. The basic shape of the Discrete Elements was spherical. The contact model comprises a spring and a dashpot in parallel in both normal and tangential directions, with a frictional slider in the tangential direction known as the Voigt model (Fig. 5). When two spherical Discrete Elements are in contact, the contact force is calculated from the relative displacement vectors in normal and tangential directions. In this study, the rolling resistance was not implemented because it might be thought that the rotation of the soil particles with irregular shape would be simulated more accurately by the shape representation of Discrete Elements than by the rolling resistance. However, several studies have emphasized the effect of the rolling resistance (Belheine et al. 2009). The translational and rotational motions of ith Discrete Element are updated using the concept of a state vector Y (t) with explicit numerical integration (Euler’s method). Y (t + ∆t) = Y (t) +
dY (t) ×∆t dt
(2)
X(t + ∆t) X(t) V (t) ∗ q(t + ∆t) q(t) 0.5ω q(t) P (t + ∆t) = P (t) + F (t) ∆t L(t + ∆t) L(t) τ (t)
(3)
In those expressions, X(t) is the position vector of the ith Discrete Element at time t; q(t) is quaternion. Furthermore, the following variables are used: P (t) and L(t) respectively denote the linear and angular momentum; V (t) is the velocity vector; ω(t) represents the angular velocity; ω ∗ is a shorthand for [0,ω(t)], and F (t) and τ (t) are the total force and torque acting on the ith Discrete Element. The rotational motion of the rigid body can be expressed by the Euler equations of motion of a body in space in relation to the principal axes of inertia of the body, in this study, the angular velocity was simply calculated by ω(t) = [I]−1 L(t). This simple method using state vector Y (t) was popular in the area of rigid body simulation (Baraff1997). However, spherical Discrete Elements offer a great advantage in computational simplicity. Deviation from the actual particle shape creates an overestimation of the rotation. It cannot simulate interlocking among particles because the normal contact force does not contribute to the moment. In addition, the rotation of spherical Discrete Elements is constrained only by frictional forces because of the adjacent Discrete Elements. The overall behavior of granular materials is well known to be affected strongly by the geometrical properties of particles, such as size, shape, deformation and crushability (Foutz et al. 1993). In DEM analysis, the particle shape representation is critical to the accuracy of the compressibility or bulk mechanical properties of the soil. In this study, arbitrary irregular shape of Discrete Elements was generated in such a way that several spherical elements were joined together to create a model particle (Abbaspour-Fard 2004; Matsushima et al. 2009). This method has greater effectiveness and simplicity in the programming and stability of the computation than the use of polyhedral Discrete Elements. Such as the availability of the sphere–sphere contact detection and calculation of the contact force. In this method, any number of spherical elements with different radii is useful to produce a model particle with or without overlap. Because the relative distances between these spherical elements composing the model particle remain fixed, each group of joined spherical elements forms a rigid body. The calculated contact forces acting on the spherical elements are transferred to the centroid of the model particle to which they belong. Once the total force and torque acting on the particle centroid are determined, the motion of the particle is updated based on the Newton’s second law; then the positions and translational and rotational velocities of these spherical elements are updated. In this thickening growth model, water, which strongly influences the radish plant growth, is not considered. 2.3 Root–Soil Contact Interaction Root–soil contact interaction was analyzed using the Finite Element – Discrete Element Method (FE–DEM) as a contact problem between root Finite Elements and soil Discrete Elements. The object modeled by assemblies of
Fig. 5. Voigt model with frictional slider.
Fig. 7. Program flow. Element shapes of three types are compared: (i) circle in 2D, (ii) sphere in 3D, and (iii) a clumped shape in 3D. The quantitative strategy is the following: first, applying soil consolidation by weight of Discrete Elements for 0.2 s; secondly, cutting the soil higher than 12 mm. Finally, starting of the thickening growth of root Finite Elements.
Fig. 6. Discrete Element – Finite Element contact interaction. FE meshes creates a polyhedral surface consisting of triangular planes, edges, and vertices. The contact detection between Discrete Elements and these geometric components was conducted through hierarchical checking procedures with the principle that a contact with one entity precludes the contact detection with all overlapping entities of lower precedence. The contact detection between a sphere and these primitive components is straightforward because these shapes are described by few parameters. A Discrete Element and Finite Element in contact are presented in Fig. 6. The surfaces of Finite Elements are covered with triangle meshes. The contact model is also a Voigt model. Therefore, FE meshes are treated as a geometric boundary condition. Furthermore, a calculated contact force was distributed to Discrete Element and FE nodes. Moreover, these distributed contact reactions can be expressed by the idea of the shape function defined at a projected contact point on the target surface element. 2.4 Flow of the Numerical simulation The algorithm of this numerical model comprises the following three main procedures: (i) thickening growth of root Finite Elements; (ii) calculation of contact force between Discrete Elements, and Finite Elements and within Discrete Elements; and (iii) numerical integration and update DEM and FEM states. The program flow is shown in Fig. 7. Consequently, the calculated shape of root Finite Elements includes results of two combined and opposing actions, i.e., the radish root’s own growth and resistance of the surrounding soil particles. Because of the enormous calculation time, it is difficult to analyze the quarter of a plant pot shown in Fig. 4. For this study, a preliminary simulation in a small scale was examined. In this preliminary simulation, the top part of the radish root and quarter of the cylindrical plant pot with 15 mm radius were analyzed (Fig. 8). The root Finite Elements thickened approximately from 1 mm to 8 mm at maximum radius. In Table 1, main parameters for DEM and FEM are listed. The results obtained using Discrete
3. RESULTS AND DISCUSSION The configurations of the numerical result at the initial and final growth stage for Discrete Element shapes of three types are depicted in Fig. 9. The reason we did not simulate the quarter of the target domain shown in Fig. 4 while reducing the number of Discrete Elements and using Discrete Elements with a large radius is the following. First, generally, the particle shape would affect the physical properties and mechanical behavior of the granular medium, and the size of Discrete Elements changes the overall deformed shape of soil Discrete Elements when the root Finite Elements thickened laterally. The large radius Discrete Elements would not be deformed smoothly with large resistance. In addition, it does not fit the root Finite Elements surface, which is the second reason. In the experiment, although the size of soil particles was not recorded, the radii of the soil particles are almost all less than 1 mm. The contact interaction between the radish root and silt and clay soil will differ from that with sand or gravel soil. Therefore, the size of Discrete Elements is an important factor in the contact interaction between root Finite Elements and soil Discrete Elements. We would not use Discrete Elements with too large a radius relative to the size of meshes covering the root Finite Elements. Regarding the results of 3D simulation, (ii) and (iii), shown in Fig. 9, it was confirmed that the extension in 3D FE–DEM can implement the contact interaction between soil Discrete Elements and root Finite Elements sufficiently. The growth rule for root Finite Elements can reproduce the growing shape of a radish root adequately. Furthermore, contact interaction between the root and soil would scarcely impede the thickening growth of root Finite Elements in this preliminary simulation. Because differentiation in the shape of storage root and its development for individuals are so much more variable than the order of the deformation caused by the contact reaction from the Discrete Elements, the calculated shapes of root Finite Elements are mostly in accordance with those of experimental results. This might result from the fact that we assumed parameters for FEM which are much harder than adequate values. The values of Young’s modulus and Poisson’s ratio were determined from aluminum alloy. Sev-
eral studies have investigated the mechanical properties of biomaterials, such as roots and seeds (Paulsen 1978; Yamamoto et al. 2000). However, the changes and distributions of mechanical properties and FEM parameters of radish roots throughout their growth period have not been clarified. The changes of the height of soil Discrete Elements for the three shape types are compared in Fig. 10. In Fig. 10, the heights of soil Discrete Elements are shown versus the converted time (day) according to the correspondence between Ti (day) and ti (s). Figure 10 shows that the extension of FE–DEM for root–soil systems in 3D can decrease the height of a soil from 18 mm to approximately 14 mm. Consequently, the 3D FE–DEM simulation is more effective for the analysis of the apparent behavior of the soil used in experiment than 2D to some extent. The simple Voigt model generally can not allow a large overlap between Discrete Elements in contact. Therefore, soil Discrete Elements were almost not compacted by the thickening root growth. It is readily apparent that when the soil is compressed laterally, the change of the height is more intense in 2D than in 3D if it occurs in an isovolumetric process. Furthermore, the clumped shape decreases the height to 13 mm. The clumped shape has little effect on the problem of increased soil height. Changes of the porosity of circular and spherical soil Discrete Elements are shown in Fig. 11. Data depicted in Fig. 11 were obtained from the simulation of the expansion of the cylindrical root analogue, in which a cylindrical wall imitating the radish root expanded radially. Furthermore, porosity n was calculated as n = (Vt − Vs )/Vt ×100, where Vt represents the total volume of soil area (mm3 ) and Vs denotes the total volume of soil Discrete Elements (mm3 ). These data suggest that as the cylindrical root expanded, soil Discrete Elements were easily deformed and disturbed similarly to dry granular materials such as glass beads. The porosity and bulk density of soil Discrete Elements is an almost constant value. The slight compressibility and change of soil Discrete Elements’ porosity would result from the geometrical change of the positions of Discrete Elements. Furthermore, in particular in the 2D case, the motion of Discrete Elements is constrained horizontally, so large local and temporal overlaps between Discrete Elements occurred, to which the strong contact force was applied. Similarly to the preliminary simulation, as the root thickened laterally, the soil was deformed, preserving the initial volume when the soil consolidation and setting soil height were finished. Furthermore, the porosity and bulk density of soil Discrete Elements is an almost constant value. Efforts to improve the apparent unnatural increase of soil and to simulate the compressibility are needed. Moreover, it might be necessary that further experimental studies of the analysis of the distribution of the porosity and bulk density of the soil from vicinity of the radish root with its thickening growth and theoretical and numerical studies of the construction and improvement of the soil modeling. It is thought that several approaches exist for this problem. One is modification of the contact model to allow the large overlap between Discrete Elements in contact to camouflage the compaction and plastic deformation of the
Table 1. Calculation conditions time step Spring Constant (normal) Spring Constant (tangential) Frictional Coefficient Young’s Modulus Poisson Ratio
1.0×10−6 s 1000 N/m 250 N/m 0.5 69 GPa 0.3
Fig. 8. Preliminary simulation of thickening growth of radish root. soil. The other is improvement in the shape representation of Discrete Elements, as conducted in this study. As described in section 2.2, the sizes and shapes of particles are critical to compressibility. Clumped Discrete Elements differ greatly from the shape and size of actual soil particles. Matsushima et al. (2009) developed methods to imitate the particle shape of lunar regoliths by combining spherical elements together to simulate actual lunar regolith particles. 4. CONCLUSION In this study, to resolve the unnatural behavior of soil Discrete Elements described in a previous report, analyses were made of a thickening growth model of radish root. Extension to 3D FE–DEM and improvement in the shape representation of soil Discrete Elements was made. A preliminary simulation confirmed that the thickening growth of root Finite Elements was not impeded by contact with soil Discrete Elements, and that 3D simulation including non-spherical Discrete Elements was able to decrease the soil Discrete Element height from 18 mm to approximately 13 mm. These results show that further studies of the mechanical properties of radish root and the adequate setting of FEM parameters are needed for accurate analysis of the contact interaction between a root and soil, and analysis of the deformation of root Finite Elements. Furthermore, to simulate the bulk mechanical behavior of soil used in experiments, further improvements in the soil model are needed, including the modification of the contact model and other approaches discussed in Section 3. REFERENCES Abbaspour-Fard, M.H. (2004). Theoretical validation of a multi-sphere, discrete element model suitable for biomaterials. Biosystems Engineering, 88(2), 153–161. Baraff, D. (1997). An Introduction to Physically Based Modeling : Rigid Body Simulation I–Unconstrained Rigid Body Dynamics. Carnegie Mellon University, Available at:
(i) 2D simulation with circular Discrete Elements.
(ii) 3D simulation with spherical Discrete Elements.
(iii) 3D simulation with clumped Discrete Elements. Fig. 9. Configurations of preliminary simulation results (left, initial stage; right, final stage).
Height(mm)
25 20 15 10 Value in circle(2D) Value in sphere(3D) Value in clump(3D)
5 0
15
20
25
30
35
time(day) Fig. 10. Change of soil Discrete Element height.
Porosity(%)
1
value in 3D value in 2D
0.8 0.6 0.4 0.2 0
0
1
2
3
4
5
Radius of Cylinderical Root(mm) Fig. 11. Change of soil Discrete Element porosity.
http://www.cs.cmu.edu/~baraff/sigcourse/notesd1.pdf. [Accessed 26 March 2010]. Belheine, N., Plassiard, J.P., Donz´e, F.V., Darve, F., and Seridi, A. (2009). Numerical simulation of drained triaxial test using 3d discrete element modeling. Computers and Geotechnics, 36(1-2), 320–331. Faure, A.G. (1994). Stress field developed by root growth: theoretical approach. Journal of Agricultural Engineering Research, 58, 53–67. Foutz, T.L., Thompson, S.A., and Evans, M.D. (1993). Comparsion of loading response of packed grain and individual kernels. Transactions of the ASAE, American Society of Agricultural Engineers, 36, 508–576. Fujita, Y. (2010). Numerical Simulation of the Thickening Growth of Radish Root Based on the Contact Interaction with Soil Environment. Master’s thesis, Graduate School of Agriculture, Kyoto University. Horner, D., Peters, J.F., and Carrillo, A. (2001). Large scale discrete element modeling of vehicle-soil interaction. Journal of Engineering Mechanics, 127(10), 1027– 1032. Matsushima, T., Katagiri, J., Uesugi, K., Tsuchiyama, A., and Nakano, T. (2009). 3d shape characterization and image-based dem simulation of the lunar soil simulant fjs-1. Journal of Aerospace Engineering, 22, 15–23. Nakashima, H., Fujita, Y., Tanaka, H., and Miyasaka, J. (2008). A model of root elongation by dynamic contact interaction. Plant Root, 2, 58–66. Nakashima, H., Takatsu, Y., Shinone, H., Matsukawa, H., and Kasetani, T. (2009). Fe-dem analysis of the effect of tread pattern on the tractive performance of tires operating on sand. Journal of Mechanical Systems for Transportation and Logistics, 25(2), 55–65. Paulsen, M.R. (1978). Fracture resistance of soybeans to compressive loading. Transactions of the ASAE, American Society of Agricultural Engineers, 21, 1210– 1216. Richards, B.G. and Greacen, E.L. (1986). Mechanical stresses on an expanding cylindrical root analogue in granular media. Australian Journal of Soil Research, 24, 393–404. Yamamoto, R., Fujii, S., Tanimoto, E., and Nevins, D.J. (2000). Computer simulation of osmotic expansion and shrinkage in okra hypocotyls segement. Biorheology, 37, 213–223.
Abstract: The mechanical impedance of soil is well known to influence the shape formation of storage roots of root vegetables strongly. Contact interaction between the root and soil become important in the growth and shape formation of the storage root. This study was aimed at construction of a numerical model for thickening growth of radish tap roots based on the root–soil contact interaction. However in 2D simulation, the height of the soil Discrete Elements increased greatly compared with experimental observations. To solve this problem, 3D simulation was extended and improvements in the shape representation of soil Discrete Elements were made. Soil Discrete Elements with non-spherical shape were generated such that several spherical Discrete Elements were joined together to create a model particle. In this report, the results of preliminary simulation obtained for a small scale are discussed. For the 3D simulation, the shape of root Finite Elements was reproduced based on growth velocities. No distortion occurred in root Finite Elements. Additional examinations of the mechanical properties of radish roots and parameter settings for root Finite Elements are needed. Furthermore, results confirmed that the extension of FE–DEM for the root–soil system in 3D was able to decrease the final height of surface soil Discrete Elements from 18 mm to approximately 14 mm. Clumped Discrete Elements decrease the height to about 13 mm. Keywords: Radish, Root, Thickening growth, 3D FE–DEM, Non-spherical discrete elements, Plants, Numerical simulation 1. INTRODUCTION It is generally recognized that the growth of plant roots can be affected adversely by high mechanical impedance imposed by their ambient soil environment. Generally, highly compacted soil interrupts root growth. Especially for root vegetables, storage roots must move the surrounding soil particles to make space as they grow in a soil medium. Consequently, the hardness, density, and structure of the soil strongly influence the development and shape formation of the storage root. Radish roots used to be evaluated according to their shape and size from the view of efficiency in transport and processing. This might also affect competitiveness on the market. Usually, straight and cylindrical shapes are desired. Although the shape patterns of storage roots of radish plants are basically programmed in their genetic DNA sequence, physical soil conditions alter these traits dramatically through vegetative growth. For example, soil layers that are heavily
compacted by wheels and tracks of heavy agricultural machines and continuous tillage operations at the same depth would produce irregular shapes of storage roots in Japanese Radish. Therefore, root–soil contact interaction plays an important role in the radish root shape formation. This study specifically examines the effects of the mechanical impedance in the radish root shape formation in terms of the contact problem. The study was designed to develop a thickening growth model of a radish tap root based on the root–soil contact interaction. As for modeling of root growth, with the increasing speed of computers, numerical methods such as Discrete Element Method (DEM) and Finite Element Method (FEM), are gradually becoming more effective techniques for the analysis of granular materials and contact problems between soil and agricultural machine devices (Horner et al. 2001; Nakashima et al. 2009). Based on the idea of cone penetration, when the root elongates in soil, development of
Fig. 1. Results of the 2D simulation of thickening growth of radish root (left, initial growth stage; right, final growth stage, Ti (day), ti (s)). 160 155
Height(mm)
a stress field in the soil medium has been analyzed using FEM (Faure 1994). In this study, root deformation was not considered. Similarly, the effects of the radial division and enlargement in the meristematic zone on the root tip and the deformation and stress field in the vicinity of the root tip have been analyzed using FEM (Richards and Greacen 1986). Recently, modeling of the root elongation growth based on the dynamic contact interaction between root and ambient soil environment using DEM was proposed (Nakashima et al. 2008). Furthermore, in a previous study, the basis of the thickening growth model of radish root based on the root–soil contact interaction was constructed in 2D (Fujita 2010). In this developed model, the DEM and the FEM were applied, respectively, to the modeled soil and root. Furthermore, root–soil contact interaction was analyzed using the Finite Element – Discrete Element Method (FE–DEM). The calculated shape of root Finite Elements mostly agreed with the experimental results. However, the behavior of the soil Discrete Elements caused larger displacements. Particularly, the height measured at the surface of soil Discrete Elements was considerably greater. Figure 1 shows numerical results in 2D simulation obtained at the initial and final growth stages. Furthermore, Fig. 2 shows changes of the height in soil Discrete Elements with time. In Fig. 1 and Fig. 2, T (day) is the time elapsed in the experiment and t(s) is that which elapsed in the simulation. Because of the extremely long calculation time of DEM simulation, a virtual time scale was introduced. The correspondence between the time scale in the simulation and in experiment was explained in previous reports (Fujita 2010). The root Finite Elements thickened laterally; soil Discrete Elements were largely deformed. Also, the height of soil Discrete Elements increased from 130 mm to approximately 155 mm. Therefore, the soil Discrete Elements had almost no compressibility, which was against the observation that the soil used in the experiment was compacted. Furthermore, in the previous model, the elongation of the hypocotyls could not be simulated apparently. To solve this unnatural increase of soil Discrete Elements, extension of FE–DEM to model root–soil systems in 3D was done and improvements in the shape of soil Discrete Elements were made.
150 145 140 135 130 125 120
Soil Height 15
20
25
30
35
time(day) Fig. 2. Height of the soil Discrete Elements shown versus time.
The main objective of this study is to construct a 3D thickening growth model of radish root by application of 3D FE–DEM.
2. MATERIALS AND METHODS
Fig. 3. Setting growth rule for root Finite Elements in 3D.
In this section, the 3D thickening growth model is described. Then, numerical methods used in the previous 2D model, FEM, DEM, and FE–DEM were extended to a 3D model. Furthermore, improvement in the shape of Discrete Elements was made. As in the previous 2D analysis, to reproduce the growth of root Finite Elements, a simplified assumption for mechanism of the thickening growth of radish root was made under which root Finite Elements would obey the experimental growth rule of “Growth Velocity” data obtained by virtual marking method in this radish experiment. The radish outline was obtained by projection onto a plane. Furthermore, the root shape was assumed to be axisymmetric. The root shape was expressed as the radius and the distance from the proximal end. Furthermore, in the 3D case, the radish growth was
assumed to be the same for all radial directions (Fig. 3). Figure 3 shows the assumed direction and magnitude of growth in the radish root. 2.1 Root Model A radish root was modeled using FEM. When the storage root grows radially, the root exerts a force to displace the soil particles; it is simultaneously reacted upon by the soil particles. The deformation and development of the stress fields within the root were calculated using FEM. The basic shape of the radish root was assumed to be axisymmetric and a quarter of the target domain was modeled. The root Finite Elements consist of standard hexahedral eight-node isoparametric elements (Fig. 4). For
or
Fig. 4. Configuration of the 3D root shape (left, full scale; right, quarter of the target domain) the material properties of the radish root, a linearly elastic model was assumed. In this model, the elastic strains {!} are related to the stresses {σ} by Hook’s law.
The thickening growth of root Finite Elements was formulated based on the prescribed “Growth Velocity” data. The displacement for all FE nodes at each time step can be expressed as Xi (t + ∆t) = Xi (t) + Vi (t)∆t, (1) where Xi (t) represents the position vector of the ith FE node at time t, mm. Also, Vi (t) is the growth velocity vector of ith FE node at time t, mms−1 . This Vi (t) was extracted from experimental results obtained using a radish. Details of the root growth experiment using the radish and the calculation method of Vi (t) were described in a previous report (Fujita 2010). This displacement is not the imposed displacement of Finite Elements, but the growth of root Finite Elements, which indirectly expresses the sum of biological activities in radish root, and which did not cause the strain and stress in the root Finite Elements. 2.2 Soil Model Soil was modeled using DEM. The basic shape of the Discrete Elements was spherical. The contact model comprises a spring and a dashpot in parallel in both normal and tangential directions, with a frictional slider in the tangential direction known as the Voigt model (Fig. 5). When two spherical Discrete Elements are in contact, the contact force is calculated from the relative displacement vectors in normal and tangential directions. In this study, the rolling resistance was not implemented because it might be thought that the rotation of the soil particles with irregular shape would be simulated more accurately by the shape representation of Discrete Elements than by the rolling resistance. However, several studies have emphasized the effect of the rolling resistance (Belheine et al. 2009). The translational and rotational motions of ith Discrete Element are updated using the concept of a state vector Y (t) with explicit numerical integration (Euler’s method). Y (t + ∆t) = Y (t) +
dY (t) ×∆t dt
(2)
X(t + ∆t) X(t) V (t) ∗ q(t + ∆t) q(t) 0.5ω q(t) P (t + ∆t) = P (t) + F (t) ∆t L(t + ∆t) L(t) τ (t)
(3)
In those expressions, X(t) is the position vector of the ith Discrete Element at time t; q(t) is quaternion. Furthermore, the following variables are used: P (t) and L(t) respectively denote the linear and angular momentum; V (t) is the velocity vector; ω(t) represents the angular velocity; ω ∗ is a shorthand for [0,ω(t)], and F (t) and τ (t) are the total force and torque acting on the ith Discrete Element. The rotational motion of the rigid body can be expressed by the Euler equations of motion of a body in space in relation to the principal axes of inertia of the body, in this study, the angular velocity was simply calculated by ω(t) = [I]−1 L(t). This simple method using state vector Y (t) was popular in the area of rigid body simulation (Baraff1997). However, spherical Discrete Elements offer a great advantage in computational simplicity. Deviation from the actual particle shape creates an overestimation of the rotation. It cannot simulate interlocking among particles because the normal contact force does not contribute to the moment. In addition, the rotation of spherical Discrete Elements is constrained only by frictional forces because of the adjacent Discrete Elements. The overall behavior of granular materials is well known to be affected strongly by the geometrical properties of particles, such as size, shape, deformation and crushability (Foutz et al. 1993). In DEM analysis, the particle shape representation is critical to the accuracy of the compressibility or bulk mechanical properties of the soil. In this study, arbitrary irregular shape of Discrete Elements was generated in such a way that several spherical elements were joined together to create a model particle (Abbaspour-Fard 2004; Matsushima et al. 2009). This method has greater effectiveness and simplicity in the programming and stability of the computation than the use of polyhedral Discrete Elements. Such as the availability of the sphere–sphere contact detection and calculation of the contact force. In this method, any number of spherical elements with different radii is useful to produce a model particle with or without overlap. Because the relative distances between these spherical elements composing the model particle remain fixed, each group of joined spherical elements forms a rigid body. The calculated contact forces acting on the spherical elements are transferred to the centroid of the model particle to which they belong. Once the total force and torque acting on the particle centroid are determined, the motion of the particle is updated based on the Newton’s second law; then the positions and translational and rotational velocities of these spherical elements are updated. In this thickening growth model, water, which strongly influences the radish plant growth, is not considered. 2.3 Root–Soil Contact Interaction Root–soil contact interaction was analyzed using the Finite Element – Discrete Element Method (FE–DEM) as a contact problem between root Finite Elements and soil Discrete Elements. The object modeled by assemblies of
Fig. 5. Voigt model with frictional slider.
Fig. 7. Program flow. Element shapes of three types are compared: (i) circle in 2D, (ii) sphere in 3D, and (iii) a clumped shape in 3D. The quantitative strategy is the following: first, applying soil consolidation by weight of Discrete Elements for 0.2 s; secondly, cutting the soil higher than 12 mm. Finally, starting of the thickening growth of root Finite Elements.
Fig. 6. Discrete Element – Finite Element contact interaction. FE meshes creates a polyhedral surface consisting of triangular planes, edges, and vertices. The contact detection between Discrete Elements and these geometric components was conducted through hierarchical checking procedures with the principle that a contact with one entity precludes the contact detection with all overlapping entities of lower precedence. The contact detection between a sphere and these primitive components is straightforward because these shapes are described by few parameters. A Discrete Element and Finite Element in contact are presented in Fig. 6. The surfaces of Finite Elements are covered with triangle meshes. The contact model is also a Voigt model. Therefore, FE meshes are treated as a geometric boundary condition. Furthermore, a calculated contact force was distributed to Discrete Element and FE nodes. Moreover, these distributed contact reactions can be expressed by the idea of the shape function defined at a projected contact point on the target surface element. 2.4 Flow of the Numerical simulation The algorithm of this numerical model comprises the following three main procedures: (i) thickening growth of root Finite Elements; (ii) calculation of contact force between Discrete Elements, and Finite Elements and within Discrete Elements; and (iii) numerical integration and update DEM and FEM states. The program flow is shown in Fig. 7. Consequently, the calculated shape of root Finite Elements includes results of two combined and opposing actions, i.e., the radish root’s own growth and resistance of the surrounding soil particles. Because of the enormous calculation time, it is difficult to analyze the quarter of a plant pot shown in Fig. 4. For this study, a preliminary simulation in a small scale was examined. In this preliminary simulation, the top part of the radish root and quarter of the cylindrical plant pot with 15 mm radius were analyzed (Fig. 8). The root Finite Elements thickened approximately from 1 mm to 8 mm at maximum radius. In Table 1, main parameters for DEM and FEM are listed. The results obtained using Discrete
3. RESULTS AND DISCUSSION The configurations of the numerical result at the initial and final growth stage for Discrete Element shapes of three types are depicted in Fig. 9. The reason we did not simulate the quarter of the target domain shown in Fig. 4 while reducing the number of Discrete Elements and using Discrete Elements with a large radius is the following. First, generally, the particle shape would affect the physical properties and mechanical behavior of the granular medium, and the size of Discrete Elements changes the overall deformed shape of soil Discrete Elements when the root Finite Elements thickened laterally. The large radius Discrete Elements would not be deformed smoothly with large resistance. In addition, it does not fit the root Finite Elements surface, which is the second reason. In the experiment, although the size of soil particles was not recorded, the radii of the soil particles are almost all less than 1 mm. The contact interaction between the radish root and silt and clay soil will differ from that with sand or gravel soil. Therefore, the size of Discrete Elements is an important factor in the contact interaction between root Finite Elements and soil Discrete Elements. We would not use Discrete Elements with too large a radius relative to the size of meshes covering the root Finite Elements. Regarding the results of 3D simulation, (ii) and (iii), shown in Fig. 9, it was confirmed that the extension in 3D FE–DEM can implement the contact interaction between soil Discrete Elements and root Finite Elements sufficiently. The growth rule for root Finite Elements can reproduce the growing shape of a radish root adequately. Furthermore, contact interaction between the root and soil would scarcely impede the thickening growth of root Finite Elements in this preliminary simulation. Because differentiation in the shape of storage root and its development for individuals are so much more variable than the order of the deformation caused by the contact reaction from the Discrete Elements, the calculated shapes of root Finite Elements are mostly in accordance with those of experimental results. This might result from the fact that we assumed parameters for FEM which are much harder than adequate values. The values of Young’s modulus and Poisson’s ratio were determined from aluminum alloy. Sev-
eral studies have investigated the mechanical properties of biomaterials, such as roots and seeds (Paulsen 1978; Yamamoto et al. 2000). However, the changes and distributions of mechanical properties and FEM parameters of radish roots throughout their growth period have not been clarified. The changes of the height of soil Discrete Elements for the three shape types are compared in Fig. 10. In Fig. 10, the heights of soil Discrete Elements are shown versus the converted time (day) according to the correspondence between Ti (day) and ti (s). Figure 10 shows that the extension of FE–DEM for root–soil systems in 3D can decrease the height of a soil from 18 mm to approximately 14 mm. Consequently, the 3D FE–DEM simulation is more effective for the analysis of the apparent behavior of the soil used in experiment than 2D to some extent. The simple Voigt model generally can not allow a large overlap between Discrete Elements in contact. Therefore, soil Discrete Elements were almost not compacted by the thickening root growth. It is readily apparent that when the soil is compressed laterally, the change of the height is more intense in 2D than in 3D if it occurs in an isovolumetric process. Furthermore, the clumped shape decreases the height to 13 mm. The clumped shape has little effect on the problem of increased soil height. Changes of the porosity of circular and spherical soil Discrete Elements are shown in Fig. 11. Data depicted in Fig. 11 were obtained from the simulation of the expansion of the cylindrical root analogue, in which a cylindrical wall imitating the radish root expanded radially. Furthermore, porosity n was calculated as n = (Vt − Vs )/Vt ×100, where Vt represents the total volume of soil area (mm3 ) and Vs denotes the total volume of soil Discrete Elements (mm3 ). These data suggest that as the cylindrical root expanded, soil Discrete Elements were easily deformed and disturbed similarly to dry granular materials such as glass beads. The porosity and bulk density of soil Discrete Elements is an almost constant value. The slight compressibility and change of soil Discrete Elements’ porosity would result from the geometrical change of the positions of Discrete Elements. Furthermore, in particular in the 2D case, the motion of Discrete Elements is constrained horizontally, so large local and temporal overlaps between Discrete Elements occurred, to which the strong contact force was applied. Similarly to the preliminary simulation, as the root thickened laterally, the soil was deformed, preserving the initial volume when the soil consolidation and setting soil height were finished. Furthermore, the porosity and bulk density of soil Discrete Elements is an almost constant value. Efforts to improve the apparent unnatural increase of soil and to simulate the compressibility are needed. Moreover, it might be necessary that further experimental studies of the analysis of the distribution of the porosity and bulk density of the soil from vicinity of the radish root with its thickening growth and theoretical and numerical studies of the construction and improvement of the soil modeling. It is thought that several approaches exist for this problem. One is modification of the contact model to allow the large overlap between Discrete Elements in contact to camouflage the compaction and plastic deformation of the
Table 1. Calculation conditions time step Spring Constant (normal) Spring Constant (tangential) Frictional Coefficient Young’s Modulus Poisson Ratio
1.0×10−6 s 1000 N/m 250 N/m 0.5 69 GPa 0.3
Fig. 8. Preliminary simulation of thickening growth of radish root. soil. The other is improvement in the shape representation of Discrete Elements, as conducted in this study. As described in section 2.2, the sizes and shapes of particles are critical to compressibility. Clumped Discrete Elements differ greatly from the shape and size of actual soil particles. Matsushima et al. (2009) developed methods to imitate the particle shape of lunar regoliths by combining spherical elements together to simulate actual lunar regolith particles. 4. CONCLUSION In this study, to resolve the unnatural behavior of soil Discrete Elements described in a previous report, analyses were made of a thickening growth model of radish root. Extension to 3D FE–DEM and improvement in the shape representation of soil Discrete Elements was made. A preliminary simulation confirmed that the thickening growth of root Finite Elements was not impeded by contact with soil Discrete Elements, and that 3D simulation including non-spherical Discrete Elements was able to decrease the soil Discrete Element height from 18 mm to approximately 13 mm. These results show that further studies of the mechanical properties of radish root and the adequate setting of FEM parameters are needed for accurate analysis of the contact interaction between a root and soil, and analysis of the deformation of root Finite Elements. Furthermore, to simulate the bulk mechanical behavior of soil used in experiments, further improvements in the soil model are needed, including the modification of the contact model and other approaches discussed in Section 3. REFERENCES Abbaspour-Fard, M.H. (2004). Theoretical validation of a multi-sphere, discrete element model suitable for biomaterials. Biosystems Engineering, 88(2), 153–161. Baraff, D. (1997). An Introduction to Physically Based Modeling : Rigid Body Simulation I–Unconstrained Rigid Body Dynamics. Carnegie Mellon University, Available at:
(i) 2D simulation with circular Discrete Elements.
(ii) 3D simulation with spherical Discrete Elements.
(iii) 3D simulation with clumped Discrete Elements. Fig. 9. Configurations of preliminary simulation results (left, initial stage; right, final stage).
Height(mm)
25 20 15 10 Value in circle(2D) Value in sphere(3D) Value in clump(3D)
5 0
15
20
25
30
35
time(day) Fig. 10. Change of soil Discrete Element height.
Porosity(%)
1
value in 3D value in 2D
0.8 0.6 0.4 0.2 0
0
1
2
3
4
5
Radius of Cylinderical Root(mm) Fig. 11. Change of soil Discrete Element porosity.
http://www.cs.cmu.edu/~baraff/sigcourse/notesd1.pdf. [Accessed 26 March 2010]. Belheine, N., Plassiard, J.P., Donz´e, F.V., Darve, F., and Seridi, A. (2009). Numerical simulation of drained triaxial test using 3d discrete element modeling. Computers and Geotechnics, 36(1-2), 320–331. Faure, A.G. (1994). Stress field developed by root growth: theoretical approach. Journal of Agricultural Engineering Research, 58, 53–67. Foutz, T.L., Thompson, S.A., and Evans, M.D. (1993). Comparsion of loading response of packed grain and individual kernels. Transactions of the ASAE, American Society of Agricultural Engineers, 36, 508–576. Fujita, Y. (2010). Numerical Simulation of the Thickening Growth of Radish Root Based on the Contact Interaction with Soil Environment. Master’s thesis, Graduate School of Agriculture, Kyoto University. Horner, D., Peters, J.F., and Carrillo, A. (2001). Large scale discrete element modeling of vehicle-soil interaction. Journal of Engineering Mechanics, 127(10), 1027– 1032. Matsushima, T., Katagiri, J., Uesugi, K., Tsuchiyama, A., and Nakano, T. (2009). 3d shape characterization and image-based dem simulation of the lunar soil simulant fjs-1. Journal of Aerospace Engineering, 22, 15–23. Nakashima, H., Fujita, Y., Tanaka, H., and Miyasaka, J. (2008). A model of root elongation by dynamic contact interaction. Plant Root, 2, 58–66. Nakashima, H., Takatsu, Y., Shinone, H., Matsukawa, H., and Kasetani, T. (2009). Fe-dem analysis of the effect of tread pattern on the tractive performance of tires operating on sand. Journal of Mechanical Systems for Transportation and Logistics, 25(2), 55–65. Paulsen, M.R. (1978). Fracture resistance of soybeans to compressive loading. Transactions of the ASAE, American Society of Agricultural Engineers, 21, 1210– 1216. Richards, B.G. and Greacen, E.L. (1986). Mechanical stresses on an expanding cylindrical root analogue in granular media. Australian Journal of Soil Research, 24, 393–404. Yamamoto, R., Fujii, S., Tanimoto, E., and Nevins, D.J. (2000). Computer simulation of osmotic expansion and shrinkage in okra hypocotyls segement. Biorheology, 37, 213–223.