Origin of nonlinear force distributions in a composite system.

Composite materials have been actively developed in recent years because they are highly functional such as lightweight, high yield strength, and superior load response. In spite of importance of the composite materials, mechanisms of the mechanical responses of composites have been unrevealed. Here, in order to understand the mechanical responses of composites, we investigated the origin and nature of the force distribution in heterogeneous materials using a soft particle model. We arranged particles with different softness in a lamellar structure and then we applied homogeneous pressure to the top surface of the system. It is found that the density in each region differently changes and then the density difference induces a nonlinear force distribution. In addition, it is found that the attractive interaction suppresses the density difference and then the force distribution is close to the theoretical prediction. Those findings may lead material designs for functional composite materials.

Introduction

Mechanical properties of single component materials are expressed by typical parameters such as Young’s modulus and Poisson’s ratio and those parameters are important for material designs, construction industries and so on. On the contrary of the single component materials, the mechanical properties of composites, which are substances of combining multiple materials, cannot be described by linear combination of those parameters of each component[1,2]. Macroscopic material properties depend on the internal structure or arrangement of the components. Thus, the composite materials with designed arrangement, which is sometimes called a metamaterial, have superior functions such as high strength despite lightness, toughness and highly fracture resistance exceeding the simple combinations[3-9]. Those composites can be well observed in nature and biology. For an example, pearl shells have lamella structure of hard and soft layers and the soft layers absorb the energy of the cracks and then the material has resistant to fracture[10]. Recently, the composite materials mimicked the biomaterial, which is called biomimetics, have been developed[11-13]. Then the composite materials have been put to practical uses in airplanes and cars by realizing weight reduction and improving heat resistance[14,15]. Therefore, development of novel composite materials is actively promoted and the elucidation of their physical properties is expected.

In order to design the composites, mechanical responses for the external force in the composites are crucial. The mechanical properties are strongly connected with the force distribution inside the material. For example, the external force locally propagates like a chain in granular materials and it is related with the unique feature of the granular materials[16]. The jamming systems, which are densely packed particles, such as colloidal dispersion systems[17], foams[18-20] and emulsions[21] also show the chain-like force propagation and then they show peculiar phenomena such as relaxation[22], shear thickening, and shear thinning[23-25]. Particle simulations are being actively performed on such systems, and it has been found that these phenomena are caused by slight changes such as the change of the number of contact particles and the rearrangement of the internal structure. Although the mechanical properties in the single component systems have been investigated, the mechanism of the distribution of the external force in the composites has been unclear.

Here, understanding the force distribution mechanism can become a base of composite material design such as mixing method and mixing ratio and it will greatly lead to industrial developments in future. In this study, we focus on the non-uniformity of elastic modulus inside the material, and perform a numerical simulation using repulsive soft particles with different softness. The purpose is to clarify the mechanism of the nonlinear distribution of the force.

Model and methods

In this study, we investigated a two-dimensional binary repulsive soft particle model where the softness of each particle is different. We consider this system as a simple model of a composite elastic body with binary components. We arranged 4096 soft disks in a 64 64 triangular grids. We set orientation of the triangular lattice so that the base of the triangle is parallel to the x axis in order that plastic rearrangement events hardly occur by the external force in y axis (see an inset of Fig. 1a). The periodic boundary condition is applied in the x direction and the wall is located at y = 0. The particles are arranged at 0.

Figure 1

(a) distribution. Dark color corresponds to = 3, while bright color corresponds to = 1. is a width of G region. The inset schematic shows the arrangement of the triangular lattice. (b) Magnitude of the force is represented by the color, as indicated by the color bar. The force is larger in G region.

Pairs of the disks i and j interact via the pairwise harmonic repulsive potential:where , where is a parameter corresponding to the softness of the disk i. is a center-to-center distance between disks i and j, and D is a diameter of disks. for and for , where is a positive constant and it corresponds to the strength of the attractive interaction. = 0 when the interaction is only repulsive. In this model, the overlap between pairs of particles is allowed and the friction between the disks is neglected for simplicity. We set and as a number of the layer from the origin in x direction and y direction, respectively. Here we normalized the length by the initial distance of the layers in y direction. Thus, before applying the external force, xy coordinate of the disk for and is equal to .

We investigated a lamella structure in which the region of soft disks (G region) and the region of hard disks (G region) are arranged parallel to x axis (Fig. 1a). We fix = 1 for the softer disk and we define for the harder disk. The width of G region is set to = 4, 6, or 8. The change of corresponds to that of the area ratio between G region and G region. We added in the y direction to the disks located at the top, that is, the external force direction is perpendicular to the lamella structure. The strength of corresponds to the force when the size of the softer disk with = 1 shrinks by 1 %. Each particle moves following the normalized overdamped equation.where is a position of disk i. The repulsive normal force between the disk i and the bottom wall is when . We set . When the maximum velocity in all particles becomes less than 1.010, we regard that the system reaches a steady state.

Force distribution in a repulsive system

Firstly, we investigated the force distribution in the repulsive system () with . Figure 1b shows the force distribution represented by the color at and then it is found that the force is larger in the G region than in the G region. It was also seen that the force distribution in the upper region is different from that in the lower region due to the influence of the boundary condition of the top surface. Here, we focus on the properties of the bulk and then we analyzed the lower region.

We also computed the displacement of the particles in x direction and in y direction at = 12 as a function of x with . Figure 2a,b show for = 4 (blue line), 6 (red line), and 8 (green line), respectively. It is found that shows a characteristic change. has a positive correlation with x in the G region, while has a negative correlation in the G region. It means that the particles at the interface between the G region and the G region are pushed toward the G region because is asymmetric at the interface. Therefore, the density in the G region becomes larger than in the G region. It is also found that the slopes of subtly decreases in the G region when increases (0.0039 for = 4 and 0.0031 for = 8). Meanwhile, is less dependent on the value of G although the constant external force is applied, not the external control of the displacement in y direction. We note here that the subtle difference of in x direction depends on , thus this comes from an effect of the top surface boundary. However, can be regarded as constant macroscopically. In addition, becomes larger with increasing .

Figure 2

(a) Displacement in the x axis direction and (b) displacement in the y axis direction at . The vertical dashed lines represent the interfaces between G region and G region. has a positive correlation with x in the G region. Meanwhile, is almost constant. (c) The normal force on the bottom surface. (d) The ratio of the normal stress at G region to that at G region . increases with increasing .

Then, we examined the normal force on the bottom surface. It is found that the normal force in G region becomes stronger than that in G region (Fig. 2c). We define as the ratio of the normal stress at = 32 to that at = 0. Figure 2d shows as a function of . It is found that slightly increases with increasing . We will discuss the mechanism for the force distribution later.

Next, we investigated the force distribution by changing between = 1.5 and 10. Figure 3a,b show at and = 12 for = 1.5 (black line), 3 (blue line), and 5 (red line), and 10 (green line), respectively. The positive relation of in the G region becomes stronger when increases. It is also found that increases with increasing . We also note that the difference of in x direction becomes larger with increasing . We confirmed that the effect of the top surface is also stronger when is larger. Therefore, we conclude that the subtle difference of in x direction is caused by the effect of the top surface boundary. However, the difference is still small and thus can be regarded as constant macroscopically. Figure 3c,d show on the bottom surface and , respectively. The dotted line in (d) corresponds to the theoretical prediction using the continuum approximation, which we discuss later. It is found that increases only from 1.2 to 1.8 and this is much less than the theoretical prediction.

Figure 3

(a) Displacement in the x axis direction and (b) displacement in the y axis direction at . (c) The normal force on the bottom surface. (d) The ratio of the normal stress at G region to that at G region . The dotted line in (d) corresponds to the theoretical prediction using the continuum approximation. is much less than the theoretical prediction. This deviation comes from the density change of G region.

Force distribution in an attractive system

Furthermore, we investigate the force distribution with the attractive interaction. Figure 4a shows with changing the spring constant at = 8 and = 3. When becomes larger, the gradient of in G region becomes smaller. It is natural that the attractive interaction suppresses the separation of the particles. Figure 4b shows with respect to . It is found that becomes larger for larger . Here we note that becomes close to the theoretical prediction (2.4 for = 3) at .

Figure 4

(a) Displacement in the x axis direction at with the attractive interaction systems. = 3 and = 8. (b) as a function of . When the attractive interaction becomes stronger, becomes flat and then becomes close to the theoretical prediction.

From and dependence of (Figs. 2d, 4b), it is found that becomes close to the theoretical prediction when the gradient of in G region is small. In addition, the deviation of from the theoretical prediction is large for larger . Meanwhile, the slope of in the G region is larger for larger , smaller , and smaller . It suggests that the deviation of from the theoretical prediction should be correlated with the slope of in the G region.

Origin of the nonlinear force distribution

Here we discuss the reason why is almost constant even though the softness is inhomogeneous. At the interface between G and G regions, the particles in the G region and the particles in the G region are staggered. The particles at the interface are geometrically pinned with respect to the movement in y direction and then should be continuous at the interface. Since the structural deformation is small, is almost independent of . In order to check the pinned effect of the staggered arrangement, we examine the system where the orientation of the triangular lattice is rotated by 30 degrees, that is, the base of the triangle is parallel to the y axis. We investigate displacements and force distribution at = 8 and = 3. It is found that is almost constant and is large in the G region (see Fig. 5). The geometrical binding becomes much weaker in y direction and then the force propagates straightly in y direction. Then the particles mainly move in y direction rather than in x direction. Thus, we confirm that the staggered interface induced that is almost constant with respect to x.

Figure 5

(a) Schematic of the arrangement of the triangular lattice. (b) Displacement in the x axis direction and (c) displacement in the y axis direction at . (d) The normal force on the bottom surface.

Then we compare the continuum approximation in our system. We assume an isotropic elastic body where the height is H. The width of G region is and the width of G region is (). We also set the interface between G region and G region is bound. This assumption represents the geometric bound interface in our simulation. It is known that , where i, j, and l are x or y and the superscript a is s for G region or h for G region. and are the stress tensor and the strain tensor in region, respectively. The sign of and is defined as positive when the direction is from the inside to the outside of each region. and are the Lamé’s constant. Young’s modulus Y and Poisson’s ratio can be described as and , respectively. To simplify, we apply the constant deformation in y direction, not the constant pressure. Then the deformation in x direction is described as . Since the force in x direction is balanced at the interface, we obtained . Here we assume that , and the Poisson’s ratio is same. It leads and , where . Then we obtainwhere . ThenandFinally, we obtainMeanwhile, the density can be described as below;where is a density before the deformation. Here we assume since . We also set c is constant because the Poisson’s ratio weakly depends on the materials (). Then it is obtained that , , and . If k is independent of the density, k should be same as the ratio of G. Here is approximately 1/3 for metals and then we obtained c = 0.5. The dotted line in Fig. 3d is a prediction at c = 0.5. It is found that obtained by the simulation is much less than the prediction for larger .

Here, it is known that the elasticity dramatically increases with increasing the density near the close packing or jamming packing fraction. Thus k depends on the density after deformation, not constant. After the compression, is slightly smaller than () and then this subtle density difference induces the large difference of the elasticity, that is, the large decrease of k. As a result, decreases depended on , or . This is consistent with the simulation results where the deviation of from the theoretical prediction should be correlated with the slope of or in the G region.

Conclusion

In summary, a simulation was performed in which a force was applied from above in a soft-particle model arranged in a triangular lattice pattern. We considered a system in which the softness of the particles was dispersed and the hardness of the internal region was non-uniform. The force distribution is localized in the harder region but it is largely deviated from the prediction from the continuum approximation. It is found that the nonlinear distribution strongly depends on the expansion of the harder region in x due to the repulsive interaction. Meanwhile, the force distribution becomes close to the theoretical prediction when the strength of the attractive interaction is similar to that of the repulsive interaction, that is, the potential is close to the symmetric.

In this time, we have clarified the most basic state, but we consider that the future work is to introduce quantitative control of roughness and friction and expand it to more realistic problems. When there is friction between particles, it is expected that the spatial distribution of displacement will be suppressed and the force propagation will change significantly. Furthermore, it is expected to be applied to more industrial fields by clarifying the mechanical properties of composite material systems such as dynamic behaviors when a force is applied locally as well as uniform pressure.