Simulation of Crack Pattern Formation Due to Shrinkage in a Drying Material

Takuya Uehara

Department of Mechanical Systems Engineering, Yamagata University, Yonezawa, Japan.

DOI: 10.4236/ojmsi.2023.111001                

Abstract

Crack patterns observed in nature have attracted the interest of researchers in various fields, and the mechanism of the pattern formation has been investigated. However, the phenomenon is very complicated, and many factors affect the process. Therefore, we are motivated to construct a general simulation code with a simple algorithm. In this study, crack pattern formation due to shrinkage caused by the drying of a wet material was simulated. The process was simplified as follows: tensile force is generated in the model, and a crack is generated when the tension exceeds a critical value. The tensile forces in the x and y directions are independently evaluated. A crack propagates perpendicular to the tension until it reaches another crack or a boundary. Based on this modeling, simulations with a two-dimensional square domain were performed. Consequently, a cross-divided pattern was generated. Assuming zigzag crack propagation, more realistic patterns were obtained. The effects of the boundary and domain size were also considered, and various characteristic patterns were obtained. Furthermore, the orientation dependency was simulated, and 45˚ declined patterns and rectangularly divided patterns were generated. The model presented in this study is very simplified and is expected to be applicable to various objects.

Keywords

Desiccation Crack, Tensile Force, Crack Initiation, Crack Propagation, Pattern Formation, Computer Simulation

1. Introduction

Various crack patterns are observed in nature, for example, those that appear when wet clayey ground dries up, a fragile glass is broken, and some kinds of plants, such as a melon, are grown up. Each crack has a unique pattern, and although not two cracks are identical, they share many common features. For example, divided subregions have almost the same surface area regardless of the overall size of the region. Every subregion has mostly a polygonal shape; the edges are straight, zigzag, or slightly curved, but never circular or highly curved. In most cases, cracks in a material are caused by tensile force, except for cases caused by an impact load such as broken glass. Mechanically, the difference in the growth rates of the inner and outer surfaces, the external tensile force, and the shrinkage of the media in a fixed container are all potential sources of tension. However, these situations are relatively reduced to the shrinkage of the material without loss of generality, and the desiccation pattern is the representative case. The physical or chemical origin of cracking is the breaking of bonds in the material, which is sometimes caused by chemical bonding in the atomistic/molecular level and sometimes by macroscopic defects and continuum mechanical stress concentration. Another important factor is irregularity. These crack patterns are not typically observed in crystalline materials, which indicate that inhomogeneity plays a key role in cracking resistance.

In this regard, desiccation pattern formation is a very complicated phenomenon and has attracted the interest of researchers in various fields. Müller, for example, conducted experiments using a starch-water mixture to model the pattern formation in basalt columns and reported their similarity [1]. Nakahara and Matsuo conducted intensive experiments using calcium carbonite paste and found that these patterns are controllable by providing flow and vibration [2] [3] [4]. Geology is one of the most enthusiastic fields related to this problem, and the mechanism of the cracking process has been clarified through the measurements of the tensile stress and strain [5] [6]. The characterization of a crack pattern and divided cell morphology has also been investigated, and not only surface patterns but also three-dimensional morphology have been explored [7] [8]. In addition to laboratory-based studies, large-scale and long-term experiments have been conducted [9], and the introduction of information technology such as deep learning has been in progress [10]. Concerning the modeling and mechanism of crack formation, the effects of material properties and geometry of the specimen correlated with tensile stress and strain have been investigated [11] [12], and these studies have been performed not only in geology but also in other fields such as chemical engineering [13]. Based on this knowledge, numerical analyses and simulations have also been performed. Sima et al. applied a discrete-element method and successfully obtained the crack pattern in a clay layer [14]. Hasebe and Oguni also developed a simulation model using a particle discretization scheme finite element method by coupling the diffusion and fracture and effectively reproduced the hierarchical pattern in the desiccation cracking [15] [16] [17]. Vo et al. investigated cracking using a cohesive fracture method and demonstrated numerical simulations [18].

The modeling and simulation of the desiccation pattern formation have been performed, as mentioned above. However, these models are focused on a specific material, and their simulation algorithms are very complicated. Thus, the application field is limited to a specific area. Therefore, in this study, we are motivated to construct a general simulation model with a simple algorithm [19] [20], so that it can be applied to a wide range of fields. The modeling method is very simplified, but we accept it to ensure generality, and the theoretical validations are expected for individual study.

2. Simulation Method

2.1. Fundamental Equations

In this study, a simple algorithm is considered according to a conventional mechanics of materials. The model material shrinks as drying progresses. The degree of drying (dryness) of the material is represented by D, and the time variation is expressed as

$\stackrel{\dot{}}{D}=\frac{c\left(D_e-D\right)}{D_e}$, (1)

where D_e represents the equilibrium value in the dryness, and c represents a constant corresponding to the drying speed. The volumetric change is assumed to be proportional to the dryness, and the shrinkage strain is expressed as

$e=kD$, (2)

where k represents a constant. In this study, tensile force, rather than stress, is considered a trigger for cracking to account for the scale dependency, i.e., a larger segment is easier to crack. The amount of shrinkage d is represented as

$d=le$, (3)

where l represents the length of the divided subregion in which the crack is generated. Tensile force T is induced according to the shrinkage and is expressed by Hooke’s law as

$T=Ed=Ele$, (4)

where E represents a constant corresponding to Young’s modulus. Here, an additional assumption is introduced. The tension near the border of the cracking subregion is set lower and large tension is assumed in the central area of the subregion, using a function $T\left(x\right)=4x\left(l-x\right)T_0$, where x represents the relative position in the subregion with values between 0 and 1, l represents the length of the subregion, and T₀ represents the tension calculated using Equation (4).

The critical value in the tension T_c for the crack initiation is provided as a computational parameter. The tensile force is calculated in the x and y directions independently, and the generated crack progresses perpendicularly to the corresponding force direction until it reaches the boundary of the calculation area or other existing cracks. The crack is thought to propagate in two ways: in a straight line or randomly zigzag manner. When a crack is generated, the tensile force in the relevant direction is released in the entire subregion, while the perpendicular force is released only beside the crack in a thin width, as is explained later using the standard results in Section 3.

2.2. Simulation Conditions

A two-dimensional square region is considered. In a situation of lab-based experiment, wet material, such as a slurry starch dissolved in water, is placed in a shallow container. The generated cracks not only create a two-dimensional pattern but also progress from the surface to the bottom. However, in this model, the effect of thickness is neglected for simplicity, which is acceptable because the crack pattern generated on the surface can be characterized as a planar pattern.

Drying is assumed to proceed uniformly in the entire region following Equation (1), with coefficient c as a constant. The material properties are also considered homogeneous, and hence, the coefficients k and E in Equations (2) and (4) are also constant. The critical value in the tensile force, T_c, will be affected by the state of the surrounding boundary. In this study, this value is varied as a parameter in the following three cases: 1) a constant in the entire region (referred to as condition B1), 2) a high value near the boundary owing to the strong connection with the container wall (B2), and 3) a low value near the boundary because of the irregularity of the contact face (B3). The size dependency of the created pattern is expected to appear, because the generated tensile force depends on the size of the region. To achieve this, the lengths of the calculation domain are varied to be 200 × 200 (referred to as L200 model), 300 × 300 (L300), and 400 × 400 (L400). The values of the parameters are listed in Table 1. Note that all parameters have nondimensionalized values. Time increment Δt is used to calculate dryness, and lattice interval Δx multiplied by the number of lattices (200, 300 or 400) corresponds to the model size. In addition, slight variations are provided for the initial conditions in dryness and tension using random numbers, considering the microscopic inhomogeneity in the material.

The results using a standard model with straight or zigzag cracking and their dependency of the model size and effects of boundary are presented in the following section. Then, the criteria for the crack initiation and propagation are modified as extended models, and the results are presented in Section 4.

3. Results and Discussion for Basic Model

3.1. Basic Results with Straight Cracks

As the most simple and typical case, the simulation result for the L300 model with a straight crack without a boundary effect is presented in Figure 1. The dryness develops monotonously and saturates to 1.0, as depicted in Figure 1(a). The final crack pattern obtained is presented in Figure 1(b), and the progress in the domain division in the early stage is exhibited in Figure 1(c), where the color indicates the identification number of the subregions. The first crack is generated vertically, and the domain is divided into two subregions. Then, a horizontal crack appears in the lefthand-side subregion, followed by the righthand-side cracking. This process yields four square-like subregions, even though they are not completely square because the cracking position is not at the very center due to the random deviation in the initial condition. The blue curve in

Table 1. Simulation parameters.

Figure 1. Standard result. Calculation size is L300 with straight crack under condition B1: (a) variation of dryness D and number of subregions N_R, (b) final crack pattern obtained at the 1000^th time step, and (c) the early process of subregion division, where the color indicates the subregion number.

represents the variation of the number of subregions, N_R. At this point, the value reaches N_R = 4, and here, we refer to this as the Rank-1 division. Tension is released at this moment and increases again as drying progresses. When the value reaches the threshold, the second division occurs. Then, each of the four subregions is divided into four smaller subregions, and 16 subregions in total are generated as Rank 2. This series of processes, i.e., crack generation, tension release, and increase in tension, is repeated two more times, and the final pattern depicted in Figure 1(b) is Rank 4 with $N_R=4^4=256$. No additional cracks are generated, because the dryness at the Rank-4 division has reached near saturation, and the subsequent increases in tension do not raise the value to the criterion.

The crack formation process in relation to the tensile force is depicted in Figure 2. The second and third rows represent the distributions of lateral (x direction) and vertical (y) tension, respectively. As explained in Section 2, the tensile force is generated to enable the generation of large tension in the middle and small tension at both ends in each region. When the lateral tension reaches the critical value at the 12th time step, a vertical crack is generated. Then, the lateral tension is released to zero, whereas the vertical tension is released in only a small

Figure 2. Crack genertion and subregion division in relation to tension distribution. The second and third rows represent the lateral and veritical tension distributions, respectively, where red represents a high value.

area along the generated crack. Subsequently, the vertical tension reaches the threshold, and a lateral crack is generated in the left and righthand-side regions in succession. The order of these events (cracking), i.e., vertical/lateral and lefthand/righthand-side orders, leading to the generation of four subregions is ideally simultaneous but fluctuates because of the random variation in the initial condition. Similar processes are repeated in each subregion, resulting in the Rank-2 division and further pattern formation.

3.2. Zigzag Cracking

In the simulation of the previous section, the crack was assumed to progress straight, but the appearance of the result is rather unnatural because of the highly geometric regularity. In the observed pattern, both in nature and experiments, the cracks are not completely straight but curved or zigzag and sawtooth-shaped. Tracking the path with theoretical precision is too difficult and necessitates a very complicated formulation. Therefore, zigzag progress is implemented using a random number as follows. In the case of a vertical crack propagating upward from the crack origin, the choice of the path is either straight above, one lattice right, or one lattice left. This choice is determined by one-third probability using a uniform random number.

Figure 3 represents the results under the same conditions as in Figure 1, except that the cracking path is set to zigzag. Although the dryness varies in the same rate, as shown in Figure 3(a), the final pattern displayed in Figure 3(b) is drastically different from that in Figure 1(b). The fundamental process is the same, as shown in Figure 3(c). The first vertical crack and two subsequent lateral cracks create four crosswise subregions in the Rank-1 division and 16 subregions in the Rank-2 division. These are obvious in the variation in the number of

Figure 3. Calculation results with zigzag cracks: (a) variation of dryness D and number of subregions N_R, (b) crack pattern at the 1000^th step, and (c) the early process of subregion division.

subregions, plotted in Figure 3(a). However, the Rank-3 division and subsequent divisions are not as obvious as that in Figure 1; the variation in the subregion number shown in Figure 3(a) is blurred compared with the drastic increase in Figure 1(a). This is because the zigzag cracks cause the variation in the size of the divided subregion, and significant difference in tension in each subregion is induced. As a result, the timing for the maximum value to reach the threshold does not coincide. The total number of subregions in the final pattern is 220, which is less than 256 in the ideal pattern. The shape of the subregions is mostly quadrangle but some of the regions seem triangular, pentagonal, or hexagonal, resulting in a realistic pattern appearance.

3.3. Size Dependency

In this model, the tensile force instead of stress is used to determine the crack initiation. This is based on the fact that a large segment is easier to fall apart than a small segment, and the threshold should be determined by a size-dependent value such as force instead of a value standardized in the length dimension. Figure 4(a) and Figure 4(b) show the final pattern calculated using three sizes for straight cracking and zigzag cracking, respectively. In both cases, the Rank-5 division occurred for L300 models, as already shown in Figure 1 and Figure 3. When the size is as small as L200, the cracking is stopped earlier at Rank 3, and finally, 64 subregions are generated. In contrast, even when the size is increased to L400, further division does not occur, and the domain division stops at Rank 5. Note that the actual size of the obtained subregion is larger in the L400 model than in the L300 model, whereas the distribution in both models is similar. According to Equations (2)-(4), the tensile force is calculated as $T=\left(kED\right)l$,

Figure 4. Crack patterns calculated with different domain sizes with (a) straight and (b) zigzag cracking by (i) 200 × 200, (ii) 300 × 300 and (iii) 400 × 400 models.

where k and E represent constants, and D is independent of size. Therefore, the tensile force is smaller in the smaller model for the same dryness, and the cracking is delayed. Consequently, in the L200 model, cracking occurs only by Rank 3. In contrast, in the larger model L400, larger tension is generated. This can induce a higher cracking rank, but the drying process finishes before the tension reaches the threshold value for the next-rank cracking. As a result, the same rank of cracking as that of the L300 model was observed.

3.4. Effect of Boundary

The container wall, corresponding to the boundary of the simulation area, would affect crack formation and propagation. In this section, the critical value was modified in the area near the boundary as $T_c^B$ against $T_c^{\text{0}}$ in the bulk, and simulations were carried out for Condition B2 for $T_c^B>T_c^0$, and B3 for $T_c^B<T_c^0$. Simulation conditions, except for T_c, are the same as those presented in Section 3.2.

Figures 5(a)-(c) show the results when $T_c^B$ is set higher than $T_c^0$ (Condition B2). As shown in Figure 5(a) for $T_c^B=60$ and $T_c^0=30$, subregions touching the boundary are slightly larger even though the overall pattern is not so much different. In Figure 5(b), for $T_c^B=40$ and $T_c^0=20$, the overall subregions become small because the critical tension is set lower. Then, the effect of the boundary is more obvious. The subregions in the boundary area become rectangular shapes, most of which have longer edges perpendicular to the boundary. When the difference between $T_c^B$ and $T_c^0$ is larger, as shown in Figure 5(c), the boundary effects spread farther; the second cells counted from the boundary are also enlarged, and fine subregions are concentrated in the central area.

Figure 5. Crack patterns calculated considering the effect of boundary; critical tension values in the boundary area $T_c^B$ and that in the bulk region $T_c^0$ were varied as condition B2 for (a)-(c), and B3 for (d)-(f).

In contrast, when $T_c^B$ is set lower than $T_c^0$ (Condition B3), the subregions near the boundary are smaller than those in the bulk area, as depicted in Figures 5(d)-(f). When the critical value in the boundary area is set to half of the bulk, the difference in cell size is moderate, as shown in Figure 5(d) and Figure 5(e), whereas the distribution is evident when the difference is set larger, as shown in Figure 5(f). This situation applies when the shrinkage of the material is significant and the contact between the wall and material is stiff, which causes strong constraints. In this study, the effect of mechanical constraints is not considered explicitly, and instead, these effects are expressed indirectly in the critical value.

4. Results and Discussion for Extended Models

4.1. Critical Force Component

In the simulations thus far, the initiation of cracks was determined by the normal components of the tensile force, F_x and F_y, independently. Here, the criterion was modified to the mean stress (F_x + F_y)/2. Then, the direction of crack propagation should be modified to the perpendicular direction to the force. For simplicity, the crack direction was fixed to ±45˚ directions. Their sign, corresponding to top-right to bottom-left or top-left to bottom-right, is selected randomly.

Simulation results are shown in Figure 6, in which Figure 6(a) and Figure 6(b) show the results for T_c = 20 and 10, respectively, without boundary effects. Oblique patterns are obtained as expected, and a finer division is generated for T_c = 10 because the critical tension is set smaller. Figure 6(c) shows the result when the boundary effect was considered as $T_c^B=30$ and $T_c^0=10$. Large regions are generated in the boundary area, and the effects are validly reflected. Theoretical support from the mechanical viewpoint is necessary, but the characteristic pattern formation will broaden the application field of this model.

Figure 6. Crack patterns calculated by criterion based on the mean stress without boundary effect (a), (b), and with the boundary effect (c).

4.2. Anisotropic Properties

In this section, based on the previous model with an independent evaluation by the x and y components, the critical values in tension are set separately in the x and y directions, as $T_c^x$ and $T_c^y$, respectively. Figure 7(a) shows the crack patterns obtained by $T_c^x=20$ and $T_c^y=40$. Divided subregions have a vertically long shape. This is because the threshold for the lateral tension in the x direction is set smaller than that in the vertical direction, and the cracking occurs more frequently in the lateral direction. This tendency is more obvious when the difference is significant, as shown in Figure 7(b). In contrast, laterally long subregions are generated for $T_c^x=40$ and $T_c^y=20$, as shown in Figure 7(c), following the same manner, and the anisotropic pattern formation was successfully simulated. These patterns have been experimentally observed and reported by Nakahara et al. [2]: When an external vibrative motion is applied in one direction, rectangular divisions similar to those shown in Figure 7 are generated. The external vibration creates an anisotropic distribution in the density and induces a difference in the strength or resistance to cracking along and across the vibration direction. This tendency is reflected in the difference in the critical values $T_c^x$ and $T_c^y$ in this study, and the effects on the cracking are validly reflected.

4.3. Shortened Cracks

In crack patterns observed in nature, some cracks do not always reach the boundary or other cracks and sometimes terminate in the middle. To regenerate this feature, an algorithm for crack propagation was modified. The basic procedure from the crack origin is the same as zigzag progress, but an additional criterion was imposed to determine whether the crack continues to propagate or terminates in the middle: when the tension of the preceding point is higher than a critical value $T_c^P$ , the crack propagates, but it stops if the tension is below $T_c^P$. This criterion becomes considerable when the crack tip comes close to the edge because of the reducing tensile force distribution, as shown in Figure 2. It may seem worrisome that if the crack does not reach the existing crack or boundary, the original region is not completely divided into subregions. However, in this algorithm, the subregion is not included in the main flow of the program but is only used to explain the results. Therefore, this modification of the condition is not crucial in the simulation scheme.

Figure 7. Crack patterns calculated with anisotropic criterion by applying different critical values $T_c^x$ and $T_c^y$ in the lateral and vertical directions, respectively.

Figure 8. Crack patterns considering the termination of crack propagation when tension on the front is below $T_c^P$. In (c), the cracks in (b) are colored by the generation order, i.e., light blue represents early stage cracks and red represents later stage cracks.

Figure 8(a) and Figure 8(b) show the results for $T_c^P=0.4T_c^0$ and $T_c^P=0.6T_c^0$, respectively. In both cases, most of the cracks terminate before reaching the other cracks. Figure 8(c) shows the same result as Figure 8(b), but each crack is colored by the generation order, i.e., the cracks in light blue are generated in the early stage, and those in orange or red are generated in the later stage. The early cracks propagate a long distance, whereas those generated later terminate within a short distance. The higher $T_c^P$ values raise the possibility of the termination of crack propagation and result in creation of short cracks. These patterns are frequently observed in nature, and the results also broaden the application area of this model.

5. Conclusions

In this study, crack pattern formation caused by shrinkage during the drying of a wet material was simulated. To simulate the complicated phenomena using a simple algorithm, detailed material behavior and mechanics were omitted, and the following main points were focused on. Tensile force is generated according to the dryness of the model material, and a crack is generated when the tension exceeds the threshold value. The forces in the vertical and lateral directions are treated independently, and they are released when cracking occurs. Cracks propagate straightly perpendicular to the tensile force; moreover, zigzag crack propagation is also modeled. The influence of the container wall is reflected in the threshold tension value. As a result, realistic pattern formation was successfully simulated. Furthermore, the model is extended to more scenarios: oblique crack patterns caused by other stress components, anisotropic patterns, and short crack formation were simulated, and wider variations of patterns were generated. In this way, the model has great flexibility to adjust to the objective patterns. These results will broaden the applicability of this algorithm to various pattens observed in nature.

For further improvements of the model, the following points are considered as future works. The drying phenomenon in the model material was assumed to be uniform, in this study; however, in a real material, it proceeds from the surface and borders to the inner region, in real material. Therefore, the calculation of the drying process should be coupled with this model. Furthermore, the effect of the thickness or depth of the drying material is not considered, and three-dimensional modeling is necessary for more precise regeneration of natural crack patterns. Mechanical consideration is of course necessary to theoretically validate this model; nevertheless, the simplified algorithm in this study will offer excellent potential for a wide range of real phenomena.

Acknowledgements

The author acknowledges Mr. Kai Watanabe, a former undergraduate student at Yamagata University, for his preliminary work on this model.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Müller, G. (1998) Starch Columns: Analog Model for Basalt Columns. Journal of Geophysical Research, 103, 15,239-15,253. https://doi.org/10.1029/98JB00389
[2]	Nakahara, A. and Matsuo, Y. (2005) Imprinting Memory into Paste and Its Visualization as Crack Patterns in Drying Process. Journal of the Physical Society of Japan, 74, 1362-1365. https://doi.org/10.1143/JPSJ.74.1362
[3]	Nakahara, A., Shinohara, Y. and Matsuo, Y. (2011) Control of Crack Pattern Using Memory Effect of Paste. Journal of Physics: Conference Series, 319, Article ID: 012014. https://doi.org/10.1088/1742-6596/319/1/012014
[4]	Matsuo, Y. and Nakahara, A. (2012) Effect of Interaction on the Formation of Memories in Paste. Journal of the Physical Society of Japan, 81, Article ID: 024801. https://doi.org/10.1143/JPSJ.81.024801
[5]	Wang, L.-L., Tang, C.-S., Shi, B., Cui, Y.-J., Zhang, G.-Q. and Hilary, I. (2018) Nucleation and Propagation Mechanisms of Soil Desiccation Cracks. Engineering Geology, 238, 27-35. https://doi.org/10.1016/j.enggeo.2018.03.004
[6]	Sawada, M., Sumi, Y. and Mimura, M. (2021) Measuring Desiccation-Induced Tensile Stress during Cracking Process. Soils and Foundations, 61, 915-928. https://doi.org/10.1016/j.sandf.2021.03.006
[7]	An, N., Tang, C.-S., Cheng, Q., Wang, D.-Y. and Shi, B. (2020) Application of Electrical Resistivity Method in the Characterization of 2D Desiccation Cracking Process of Clayey Soil. Engineering Geology, 265, Article ID: 105416. https://doi.org/10.1016/j.enggeo.2019.105416
[8]	Zhuo, Z., Zhu, C., Tang, C.-S., Xu, H., Shi, X. and Mark, V. (2022) 3D Characterization of Desiccation Cracking in Clayey Soils Using a Structured Light Scanner. Engineering Geology, 299, Article ID: 106566. https://doi.org/10.1016/j.enggeo.2022.106566
[9]	Cordero, J.A., Prat, P.C. and Ledesma, A. (2021) Experimental Analysis of Desiccation Cracks on a Clayey Silt from a Large-Scale Test in Natural Conditions. Engineering Geology, 292, Article ID: 106256. https://doi.org/10.1016/j.enggeo.2021.106256
[10]	Han, X.-L., Jiang, N.-J., Yang, Y.-F., Choi, J., Singh, D.N., Beta, P., Du, Y.-J. and Wang, Y.-J. (2022) Deep Learning Based Approach for the Instance Segmentation of Clayey Soil Desiccation Cracks. Computers and Geotechnics, 146, Article ID: 104733. https://doi.org/10.1016/j.compgeo.2022.104733
[11]	Zeng, H., Tang, C-S., Cheng, Q., Inyang, H.I., Rong, D.-Z., Lin, L. and Shi, B. (2019) Coupling Effects of Interfacial Friction and Layer Thickness on Soil Desiccation Cracking Behavior, Engineering Geology, 260, Article ID: 105220. https://doi.org/10.1016/j.enggeo.2019.105220
[12]	Tang, C.-S., Cheng, Q., Lin, L., Tian, B.-G., Zeng, H. and Shi, B. (2022) Study on the Dynamic Mechanism of Soil Desiccation Cracking by Surface Strain/Displacement Analysis. Computers and Geotechnics, 152, Article ID: 104998. https://doi.org/10.1016/j.compgeo.2022.104998
[13]	Shepidchenko, T., Zhang, J., Tang, X., Liu, T., Dong, Z., Zheng, G. and Yang, L. (2020) Experimental Study of the Main Controlling Factors of Desiccation Crack Formation from Mud to Shale. Journal of Petroleum Science and Engineering, 194, Article ID: 107414. https://doi.org/10.1016/j.petrol.2020.107414
[14]	Sima, J. Jiang, M. and Zhou, C. (2014) Numerical Simulation of Desiccation Cracking in a Thin Clay Layer Using 3D Discrete Element Modeling. Computers and Geotechnics, 56, 168-180. https://doi.org/10.1016/j.compgeo.2013.12.003
[15]	Hirobe, S. and Oguni, K. (2016) Coupling Analysis of Pattern Formation in Desiccation Cracks. Computer Methods in Applied Mechanics and Engineering, 307, 470-488. https://doi.org/10.1016/j.cma.2016.04.029
[16]	Hirobe, S. and Oguni, K. (2017) Numerical Simulation of Desiccation Cracking Process by Weak Coupling of Desiccation and Fracture. International Journal of GEOMATE, 12-33, 8-13. https://doi.org/10.21660/2017.33.2535
[17]	Hirobe, S. and Oguni, K. (2017) Modeling and Numerical Investigations for Hierarchical Pattern Formation in Desiccation Cracking. Physica D, 359, 29-38. https://doi.org/10.1016/j.physd.2017.08.002
[18]	Vo, T.D., Pouya, A., Hemmati, S. and Tang, A.M. (2017) Numerical Modelling of Desiccation Cracking of Clayey Soil Using a Cohesive Fracture Method. Computers and Geotechnics, 85, 15-27. https://doi.org/10.1016/j.compgeo.2016.12.010
[19]	Watanabe, K. (2021) Simulation of Crack Pattern Formation Using Computer Simulation. Graduation Thesis, Yamagata University, Yonezawa. (In Japanese)
[20]	Uehara, T. and Watanabe, K. (2022) Simulation of Crack-Pattern Formation Due to Drying. In: Proceedings of the 7th Symposium on Multi-Scale Mechanics of Materials, The Society of Materials Science, Japan, Online, P43. (In Japanese)