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.

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][4][5][6][7][8][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][12][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][19][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][24][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. www.nature.com/scientificreports/ 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 y > 0. Pairs of the disks i and j interact via the pairwise harmonic repulsive potential: where k ij = (G i + G j )/2 , where G i is a parameter corresponding to the softness of the disk i. r ij is a center-tocenter distance between disks i and j, and D is a diameter of disks. f (r) = k att for r < 0 and f (r) = 1 for r > 0 , where k att is a positive constant and it corresponds to the strength of the attractive interaction. k att = 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 n x and n y 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 n x and n y is equal to (x, y) = (n x / √ 3, n y ). We investigated a lamella structure in which the region of soft disks (G s region) and the region of hard disks (G h region) are arranged parallel to x axis (Fig. 1a). We fix G i = 1 for the softer disk and we define G i = G h for the harder disk. The width of G h region W n is set to W n = 4, 6, or 8. The change of W n corresponds to that of the area ratio between G s region and G h region. We added F ex = −0.01 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 F ex corresponds to the force when the size of the softer disk with G i = 1 shrinks by 1 %. Each particle moves following the normalized overdamped equation.
where r i is a position of disk i. The repulsive normal force between the disk i and the bottom wall is F N = (D/2 − y i )(G w + G i )/2 when y i < D/2 . We set G w = 1000 . When the maximum velocity in all particles becomes less than 1.0×10 −6 , 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 ( k att = 0 ) with G h = 3 . Figure 1b shows the force distribution represented by the color at W n = 8 and then it is found that the force is larger in the G h region than in the G s 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 u x and in y direction u y at n y = 12 as a function of x with G h = 3 . Figure 2a,b show u x u y for W n = 4 (blue line), 6 (red line), and 8 (green line), respectively. It is found that u x shows a characteristic change. u x has a positive correlation with x in the G h region, while u x has a negative correlation in the G s region. It means that the particles at the interface between the G s region and the G h region are pushed toward the G s region because k ij is asymmetric at the interface. Therefore, the density in the G s region becomes larger than in the G h region. It is also found that the slopes of u x subtly decreases in the G h region when W n increases (0.0039 for W n = 4 and 0.0031 for W n = 8). Meanwhile, u y 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 u y in x direction depends on n y , thus this comes from an effect of the top surface boundary. However, u y can be regarded as constant macroscopically. In addition, u y becomes larger with increasing W n .
Then, we examined the normal force F N on the bottom surface. It is found that the normal force in G h region becomes stronger than that in G s region (Fig. 2c). We define α as the ratio of the normal stress at n x = 32 to that www.nature.com/scientificreports/ at n x = 0. Figure 2d shows α as a function of W n . It is found that α slightly increases with increasing W n . We will discuss the mechanism for the force distribution later. Next, we investigated the force distribution by changing G h between G h = 1.5 and 10. Figure 3a,b show u x u y at W n = 8 and n y = 12 for G h = 1.5 (black line), 3 (blue line), and 5 (red line), and 10 (green line), respectively. The positive relation of u x in the G h region becomes stronger when G h increases. It is also found that u y increases with increasing G h . We also note that the difference of u y in x direction becomes larger with increasing G h . We confirmed that the effect of the top surface is also stronger when G h is larger. Therefore, we conclude that the subtle difference of u y in x direction is caused by the effect of the top surface boundary. However, the difference is still small and thus u y can be regarded as constant macroscopically. Figure 3c,d show F N 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.

Force distribution in an attractive system
Furthermore, we investigate the force distribution with the attractive interaction. Figure 4a shows u x (x) with changing the spring constant k att at W n = 8 and G h = 3. When k att becomes larger, the gradient of u x in G h region becomes smaller. It is natural that the attractive interaction suppresses the separation of the particles. Figure 4b shows α with respect to k att . It is found that α becomes larger for larger k att . Here we note that α becomes close to the theoretical prediction (2.4 for G h = 3) at k att = 1.
From W n and k att dependence of α (Figs. 2d, 4b), it is found that α becomes close to the theoretical prediction when the gradient of u x in G h region is small. In addition, the deviation of α from the theoretical prediction is large for larger G h . Meanwhile, the slope of u x in the G h region is larger for larger G h , smaller W n , and smaller k att . It suggests that the deviation of α from the theoretical prediction should be correlated with the slope of u x in the G h region.

Origin of the nonlinear force distribution
Here we discuss the reason why u y is almost constant even though the softness is inhomogeneous. At the interface between G s and G h regions, the particles in the G s region and the particles in the G h region are staggered. The particles at the interface are geometrically pinned with respect to the movement in y direction and then u y should be continuous at the interface. Since the structural deformation is small, u y is almost independent of G i . 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 www.nature.com/scientificreports/ displacements and force distribution at W n = 8 and G h = 3. It is found that u x is almost constant and u y is large in the G h 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 u y is almost constant with respect to x. Then we compare the continuum approximation in our system. We assume an isotropic elastic body where the height is H. The width of G s region is W (s) and the width of G h region is W (h) ( W (h) ≪ W (s) ). We also set the interface between G s region and G h region is bound. This assumption represents the geometric bound interface in our simulation. It is known that P   ij 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 Y = µ(3 + 2µ)/( + µ) and σ = /2( + µ) , respectively. To simplify, we apply the constant deformation HE yy in y direction, not the constant pressure. Then the deformation in x direction is described as xx . Since the force in x direction is balanced at the interface, we obtained P Meanwhile, the density can be described as below; . www.nature.com/scientificreports/ where ρ 0 is a density before the deformation. Here we assume γ ≈ 0 since W (h) ≪ W (s) . We also set c is constant because the Poisson's ratio weakly depends on the materials ( 0 < σ < 0.5 ).
Then it is obt aine d t hat E . 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 G h .
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, ρ (h) is slightly smaller than ρ (s) ( �ρ > 0 ) 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 E (h) xx . This is consistent with the simulation results where the deviation of α from the theoretical prediction should be correlated with the slope of u x or E (h) xx in the G h 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.