A second-degree micromorphic continuum model for 1D granular structures

Abstract

This paper introduces a second-degree micromorphic continuum model for one-dimensional (1D) granular structures. The proposed framework aims to capture the mechanical behavior of materials with a distinct microstructure, such as granular rods or layered composites, by incorporating internal degrees of freedom. The kinematics are described by both macro-scale displacement and micro-scale deformation fields, where the micro-displacement within a representative line element (LE) is approximated by a third-order polynomial expansion. This approach introduces higher-order kinematic variables (micro-stretch and micro-curvature) and their gradients, accounting for non-local effects and size-dependent responses. The governing equilibrium equations and the corresponding boundary conditions (BCs) are obtained using the principle of minimum potential energy. The resulting system of sixth-order differential equations is solved using the state-space approach. A systematic investigation of 12 BC scenarios (Dirichlet and Neumann types) reveals the critical role of microstructural constraints in energy distribution. It is demonstrated that Dirichlet-type geometric clamping leads to pronounced boundary layers and energy localization peaks, whereas Neumann-type conditions promote field homogenization. The results highlight that micro-curvature and higher-order gradients are indispensable for predicting internal redistribution and potential failure zones in microstructured materials, providing a foundation for the design of architected metamaterials with tailored properties.

Keywords

micromorphic continuum model second-degree micromorphic theory granular structures one-dimensional (1D) materials micro-displacement polynomial expansion micro-stretch and micro-curvature

1. Introduction

Materials with pronounced microstructure – such as granular media, geomaterials, layered composites, and architected metamaterials – exhibit mechanical responses that cannot be adequately captured by classical Cauchy elasticity. The classical theory assumes that a material point is structureless and that stresses depend solely on the first gradient of displacement. As a result, it fails to account for size effects, internal length scales, and microstructural interactions that become significant when the characteristic size of grains or layers is comparable to the macroscopic dimensions of the body. These limitations have been widely recognized in the literature on generalized continua and nonlocal theories [1 –4]. In particular, experimental and theoretical studies have repeatedly shown that classical elasticity systematically underestimates stiffness at small scales, cannot reproduce dispersion curves observed in structured media, and fails to predict stress concentrations near interfaces in heterogeneous materials.

To overcome these deficiencies, several generalized continuum theories have been developed, including micropolar, micro-stretch, strain-gradient, and micromorphic models [1,3,5,6]. All these approaches share a common idea: the kinematics of a material point must be enriched to reflect the presence of an underlying microstructure. Instead of a single displacement field, additional internal variables – such as micro-rotations, micro-deformations, or higher-order gradients – are introduced to represent the relative motion and deformation of sub-structures within a representative volume element. Micropolar elasticity, in particular, provides an enriched kinematic description by introducing independent rotations and couple-stresses. The contributions of Nikabadze and Ulukhanyan to the theory of thin bodies and multilayer structures [7,8], as well as recent advances in deriving dispersion equations and wave velocities for micropolar media [9], demonstrated that these extended kinematic descriptions and proper formulation of interlayer contact conditions are essential for capturing boundary effects and scale-dependent behavior that are absent in classical elasticity. These findings highlight the necessity of continuum models that explicitly account for microstructural degrees of freedom.

A parallel and highly influential line of research is the granular micromechanics approach developed by Misra and co-authors, which derives micromorphic continuum models directly from grain-scale interactions [10 –14]. In this framework, the internal variables of a micromorphic continuum are not introduced ad hoc but emerge as systematic averages of discrete particle motions and contact forces. As a result, micro-deformation modes, micro-rotations, and couple-stress effects acquire a clear physical interpretation and can be directly related to the geometry and topology of the granular assembly.

The general micromorphic theory of degree n proposed by Nejadsadeghi and Misra [15] provides a flexible continuum framework in which the micro-deformation field is approximated by a polynomial expansion of arbitrary order. In the present work, we adopt this general formulation as a conceptual basis but focus specifically on the second-degree micromorphic case. This second-degree specialization introduces micro-curvature and higher-order deformation modes, which are essential for capturing boundary-layer effects and internal bending mechanisms in granular rods.

The applicability of this higher-order micromorphic framework to granular structures was further demonstrated in the subsequent work of Nejadsadeghi and Misra [16], where the statics and dynamics of granular-microstructured rods were analyzed. Their results showed that internal length scales and higher-order micro-deformation modes lead to non-classical distributions of strain energy and a strong dependence of the response on BCs, even in nominally one-dimensional (1D) settings.

Most existing micromorphic models, however, are limited to first-degree approximations, which account only for linear micro-stretch. Such models cannot capture micro-curvature, micro-bending, or nonlinear deformation patterns within individual grains. This limitation becomes especially severe when the microstructure exhibits pronounced bending, twisting, or higher-order deformation modes, as is the case in granular chains, beam-like metamaterials, and pantographic structures. Higher-order micromorphic and generalized continua have also been extensively developed by Francesco dell’Isola and collaborators [17,18], who demonstrated that enriched kinematics with multiple internal variables and characteristic lengths are essential for capturing the behavior of metamaterials, pantographic structures, and architected lattices. Their studies on higher-order continua and internally constrained systems provide numerous examples where classical and first-order generalized models fail to reproduce experimentally observed stiffness, dispersion, and localization patterns, whereas higher-order micromorphic descriptions succeed. Their work provides strong evidence that higher-order gradients and micro-curvature terms are not optional refinements but necessary ingredients for accurately modeling materials with complex internal geometry [19 –23]. In this sense, second-degree micromorphic models can be viewed as a minimal extension capable of capturing micro-curvature and internal bending modes while remaining analytically and computationally tractable.

Boundary layer phenomena are central to generalized continua. Near boundaries or interfaces, sharp variations in stresses and strains arise that classical elasticity cannot predict. These boundary layers are closely linked to the presence of internal length scales: when the characteristic size of the microstructure is not negligible compared to the structural dimensions, the influence of the boundary penetrates into the bulk over a finite distance, leading to non-uniform distributions of micro-deformation and energy. In micropolar, strain-gradient, and micromorphic models, the presence of additional kinematic fields and higher-order stresses naturally leads to localized deformation patterns near supports, interfaces, or material discontinuities [6,24]. Similar behavior appears in higher-order micromorphic models, where micro-curvature and higher gradients generate localized energy peaks that classical Cauchy elasticity cannot predict. Granular micromechanics-based micromorphic models likewise demonstrate boundary-layer formation and energy localization in granular rods and microstructured solids [25]. These results collectively indicate that the correct description of boundary layers in microstructured media requires both internal degrees of freedom and higher-order spatial derivatives in the constitutive framework.

Granular rods serve as a fundamental model system for studying the interplay between macro- and micro-deformations. They provide insight into the behavior of more complex two-dimensional (2D) and three-dimensional (3D) structures, including geomaterials, concrete, layered composites, and architected metamaterials. The works by Tan and Poh demonstrate that even in periodic layered structures, microstructural effects strongly influence wave propagation and static stress distributions [26,27]. From a modeling perspective, 1D granular rods offer a clean setting in which the impact of internal length scales, higher-order gradients, and BCs can be analyzed in detail without the additional complications of full 3D geometry [23,28 –32]. This makes them an ideal testbed for assessing the capabilities and limitations of second-degree micromorphic theories and for identifying the regimes in which classical or lower-order models break down. Taken together, these developments across micropolar, strain-gradient, and higher-order micromorphic theories consistently indicate that internal rotations, micro-curvature, and higher-order gradients are indispensable for describing scale-dependent behavior, wave dispersion, and boundary-layer phenomena in microstructured solids.

In parallel to these developments, a complementary line of research has focused on the continuum modeling of 1D chain-like systems with internal resonators and architected microstructures. In particular, variational approaches to finite-length Maxwell–Rayleigh-type lattices have demonstrated how microstructural design and BCs influence dispersion properties and the emergence of frequency band gaps [33]. Furthermore, higher-order continuum formulations, including second-gradient models, have been successfully employed to describe wave dispersion and attenuation phenomena in microstructured materials [34]. Although these approaches are primarily developed in the context of wave propagation, they emphasize the fundamental role of internal degrees of freedom and higher-order gradients, which are also central to the present micromorphic framework.

In this work, we develop a second-degree micromorphic theory for 1D granular rods, extending the granular micromechanics-based framework of Nejadsadeghi and Misra [15,16]. Our primary focus is the analysis of boundary layers and energy localization under various types of BCs, including both classical (Dirichlet) and generalized (Neumann) constraints. We demonstrate that the second-degree micromorphic model accurately captures the redistribution of strain energy within the microstructure, including the emergence of localized deformation modes that cannot be described by lower-order theories. Special attention is paid to the role of BCs imposed on the micro-fields: by comparing different combinations of Dirichlet- and Neumann-type constraints at the macro- and micro-levels, we show how the choice of BCs controls the thickness of boundary layers, the intensity of energy localization, and the effective stiffness of the rod. The results highlight the importance of micro-curvature and higher-order gradients for the correct modeling of granular materials and provide a foundation for the design of new metamaterials with tailored mechanical properties. Moreover, the proposed framework offers a systematic way to connect discrete granular models with continuum descriptions of higher order, thereby contributing to the broader effort of developing multiscale models that remain predictive across a wide range of length scales.

2. Continuum framework: a micromorphic model of degree 2

2.1. Kinematic variables

Let us consider a 1D object of length L that possesses an underlying granular microstructure. This microstructure is envisioned as being composed of numerous grains, each characterized by its own random mechanical and inertial properties (henceforth referred to as a 1D granular rod), as depicted in Figure 1.

Figure 1.

Schematic representation of a 1D granular structure in its reference and current configurations, illustrating macro- and micro-scale kinematics.

At the larger, macroscopic scale, the object can be treated as a continuum. Here, a material point P is identified by the macro-scale coordinate system X. We define X as the initial position of point P, and $x = Ξ (X)$ as its current position, where Ξ is the macro-scale placement function. The overall macro-scale displacement is then given by $u = x - X$ .

Moving to a finer spatial resolution, the material point P is understood to represent a collection of grains, forming a line element (LE) with a characteristic length $L^{'} = dX$ . For granular structures that exhibit periodicity, such as a composite material with repeating grains, this LE is analogous to a unit cell. In the case of non-periodic granular arrangements, the LE is chosen as the domain over which local (micro-scale) deformation is homogenized. This choice is guided by the mechanical problem’s requirements, ensuring the LE contains a sufficiently large number of grains to uphold the continuum assumption.

The position of each individual grain within the LE is tracked using the micro-scale coordinate system $X^{'}$ . This local coordinate system is attached to the center of mass (COM) of the material point P, remains aligned with the macro-scale coordinate system X, and translates with the macro-scale displacement u. The micro-scale displacement of a grain $u^{'}$ is expressed as $u^{'} = x^{'} - X^{'}$ , where $X^{'}$ and $x^{'} = Ξ^{'} (X, X^{'})$ represent the position vectors of a grain centroid at the initial and current configurations, respectively, and $Ξ^{'}$ is the micro-scale placement function. For the present static analysis, both micro- and macro-scale deformations are assumed to be infinitesimal, continuous, and differentiable functions of their respective coordinates up to the third order for micro-displacement. Thus, we can articulate the displacements as:

\begin{matrix} u = u (x), u^{'} = u^{'} (x, x^{'}) . \end{matrix}

(1)

Central to the proposed framework is a micromorphic theory of degree 2. In this context, the micro-scale displacement $u^{'}$ within the LE is approximated by a polynomial expansion in the micro-coordinate $x^{'}$ , specifically around the COM where $x^{'} = 0$ . Drawing upon established approaches [15], we incorporate terms up to the third order in this expansion:

\begin{matrix} u^{'} (x, x^{'}) = ψ_{11} (x) x^{'} + ψ_{111} (x) {(x^{'})}^{2} + ψ_{1111} (x) {(x^{'})}^{3}, \end{matrix}

(2)

here, $ψ_{11}$ represents a micro-stretch, $ψ_{111}$ represents a micro-curvature, and $ψ_{1111}$ represents the gradient of micro-curvature.

In (2), the coefficients $ψ_{11}$ , $ψ_{111}$ , and $ψ_{1111}$ are the functions of x only, capturing the local deformation within the LE at a given macro-position and ${(x^{'})}^{2}$ and ${(x^{'})}^{3}$ are the square and cube of micro-coordinate $x^{'}$ . It is important to acknowledge that for systems with a limited number of particles (or layers, for instance, in 1D composite models), a more precise displacement approximation might be achieved by partitioning the LE into distinct regions with varying strain regimes (e.g., as discussed in [26,27]). However, for scenarios involving a large number of particles or layers, linear, quadratic, or cubic approximations spanning the entire LE often prove to be the most practical [35,36]. We also note the ongoing research in formal homogenization and continualization of mass-spring systems, such as in [31,32], which frequently explore multi-scale decompositions of the displacement field to develop continuum models.

By coherently combining the macro-scale displacement u with the micro-scale displacement $u^{'}$ (from (2)), the total displacement vector for individual grains within the LE can be formally written as:

\begin{matrix} ϕ (x, x^{'}) = u + u^{'} = \bar{ϕ} (x) + ψ_{11} (x) x^{'} + ψ_{111} (x) {(x^{'})}^{2} + ψ_{1111} (x) {(x^{'})}^{3}, \end{matrix}

(3)

here, $\bar{ϕ} = u$ is adopted to maintain consistency in variable nomenclature with prior publications [12,16].

2.2. Deformation measures and intergranular kinematics

To quantify the deformation state within the granular rod and subsequently formulate the constitutive laws, we introduce several relative deformation measures. These measures are designed to capture the intrinsic differences between the macro-scale displacement gradients and the micro-scale kinematic variables. For a micromorphic theory of degree 2, these are systematically defined as follows [1,2,12]:

\begin{matrix} γ_{11} = {\bar{ϕ}}_{, x} - ψ_{11}, γ_{111} = ψ_{11, x} - ψ_{111}, γ_{1111} = ψ_{111, x} - ψ_{1111} . \end{matrix}

(4)

In equation (4), differentiation is consistently performed with respect to the macro-scale coordinate system. The subscript $, x$ henceforth denotes this spatial differentiation.

The terms introduced can be interpreted in the standard 1D setting as follows: $γ_{11}$ represents a relative stretch (or macro–micro displacement gradient), $γ_{111}$ corresponds to a relative curvature measure (or micro-rotation gradient), and $γ_{1111}$ characterizes a higher-order curvature gradient (or relative micro-strain gradient).

A key constitutive assumption for a micromorphic theory of degree 2 is made at this point: we assume that the highest-order relative deformation measure, $γ_{1111}$ , vanishes. This assumption effectively establishes the relationship $ψ_{1111} = ψ_{111, x}$ , thereby simplifying the overall model. This constraint implies that the highest-order micro-deformation mode is not treated as an independent kinematic variable but is instead fully determined by the spatial gradient of the micro-curvature. From a physical standpoint, this corresponds to neglecting independent third-order micro-fluctuations within the representative LE while retaining the dominant contributions associated with micro-curvature and its gradient. Such a simplification is commonly adopted in higher-order micromorphic and granular micromechanics-based continua to ensure a minimal yet well-posed model.

To construct the constitutive model from underlying grain interactions, it is essential to define the relative displacement between any two neighboring grains n and p, denoted $δ^{np}$ . Utilizing the total displacement field given in (3) and the relative deformation measures from (4), this objective intergranular relative displacement for the 1D case can be expressed as:

\begin{matrix} δ^{np} = ϕ^{p} - ϕ^{n} = δ^{M} - δ^{m_{1}} - δ^{m_{2}} + δ^{g} + J_{2}^{np} ψ_{11, x} . \end{matrix}

(5)

The terms on the right-hand side of (5) represent distinct physical contributions to $δ^{np}$ , defined by the following micro-scale kinematic measures:

\begin{matrix} \begin{matrix} δ^{M} = J_{1}^{np} {\bar{ϕ}}_{, x}, δ^{m_{1}} = J_{1}^{np} γ_{11}, δ^{m_{2}} = J_{2}^{np} γ_{111}, δ^{g} = J_{3}^{np} ψ_{111, x} . \end{matrix} \end{matrix}

(6)

Equation (6) meticulously decomposes the intergranular relative displacement $δ^{np}$ into its constituent physical components. These components are: $δ^{M}$ is the macro-scale deformation component, representing the portion of the relative displacement driven by the macro-scale displacement gradient ${\bar{ϕ}}_{, x}$ . $δ^{m_{1}}$ is the first-order micro-fluctuation, which quantifies the contribution stemming from the relative macro-micro displacement gradient $γ_{11}$ (the difference between ${\bar{ϕ}}_{, x}$ and $ψ_{11}$ ). $δ^{m_{2}}$ is the second-order micro-fluctuation, accounting for the contribution arising from the relative micro-stretch gradient $γ_{111}$ (the difference between $ψ_{11, x}$ and $ψ_{111}$ ). $δ^{g}$ is the higher-order micro-deformation, which specifically denotes the portion of the relative displacement attributed to the gradients of the micro-scale kinematic measures, often collectively termed higher-order terms.

Furthermore, we define the geometry moment measures $J_{1}^{np}$ , $J_{2}^{np}$ , and $J_{3}^{np}$ . These moments are given by $J_{1}^{np} = l^{p} - l^{n}$ , $J_{2}^{np} = l^{p} l^{p} - l^{n} l^{n}$ , and $J_{3}^{np} = l^{p} l^{p} l^{p} - l^{n} l^{n} l^{n}$ , where $l^{q}$ represents the scalar distance of the grain q centroid from the COM of the LE. These geometry moments are crucial as they characterize the relative spatial positions of grains n and p within the LE. More explicitly, we can write $J_{2}^{np} = (l^{p} - l^{n}) (l^{p} + l^{n}) = J_{1}^{np} (l^{p} + l^{n})$ , and $J_{3}^{np} = (l^{p} - l^{n}) ({(l^{p})}^{2} + l^{p} l^{n} + {(l^{n})}^{2}) = J_{1}^{np} ({(l^{p})}^{2} + l^{p} l^{n} + {(l^{n})}^{2})$ . It is noteworthy that the higher-order moments, $J_{2}^{np}$ and $J_{3}^{np}$ , directly imply that the greater the distance of grains n and p from the LE’s COM, the more significant their contribution to the micro-fluctuations $δ^{m_{2}}$ and the higher-order term $δ^{g}$ .

2.3. Constitutive equations

We now proceed to establish the correspondence between the kinematic measures and the various stress quantities. We postulate that the macro-scale deformation energy density, denoted by W, is a function of the continuum kinematic measures: ${\bar{ϕ}}_{, x}, γ_{11}$ , $γ_{111}$ and $ψ_{111, x}$ . Therefore, its functional form is $W = W ({\bar{ϕ}}_{, x}, γ_{11}, γ_{111}, ψ_{111, x})$ .

The macro-scale stress measures are defined as work-conjugates to these continuum kinematic measures. Specifically, they are expressed as:

\begin{matrix} \begin{matrix} τ_{11} = \frac{∂W}{\partial {\bar{ϕ}}_{, x}} \\ σ_{11} = \frac{∂W}{\partial γ_{11}}, \\ σ_{111} = \frac{∂W}{\partial γ_{111}}, \\ μ_{1111} = \frac{∂W}{\partial ψ_{111, x}} \end{matrix} \end{matrix}

(7)

This macro-scale energy density W can be systematically derived through the homogenization of the micro-scale deformation energy. If $W^{α}$ represents the deformation energy associated with the α-th interacting grain pair within an LE of length $L^{'}$ , then W can be expressed as:

\begin{matrix} W = \frac{1}{L^{'}} \sum_{α} W^{α} (δ^{αM}, δ^{α m_{1}}, δ^{α m_{2}}, δ^{αg}) \end{matrix}

(8)

here, α represents all possible grain-pair interactions (and not only nearest neighbors) and is assumed sufficiently large to ensure the validity of the continuum approximation. The numerical value of α depends on the nature of grain-pair interactions and the grain arrangement.

In turn, the intergranular forces, which are the fundamental elements of the micro-scale interactions, can be defined as work-conjugates to the micro-scale kinematic measures introduced in (6). These forces are expressed as:

\begin{matrix} \frac{\partial W^{α}}{\partial δ^{αM}} = f^{αM}, \frac{\partial W^{α}}{\partial δ^{α m_{1}}} = f^{α m_{1}}, \frac{\partial W^{α}}{\partial δ^{α m_{2}}} = f^{α m_{2}}, \frac{\partial W^{α}}{\partial δ^{αg}} = f^{αg} . \end{matrix}

(9)

By substituting the relationship given in (8) into the stress definitions (7), and concurrently employing the definitions from equations (6) and (9), the macro-scale stress measures can be explicitly formulated. These expressions relate the macro-scale stresses to the averaged micro-scale forces and their corresponding geometry moments:

\begin{matrix} \begin{matrix} τ_{11} = \frac{1}{L^{'}} \sum_{α} f^{αM} J_{1}^{α}, σ_{11} = \frac{1}{L^{'}} \sum_{α} f^{α m_{1}} J_{1}^{α}, σ_{111} = \frac{1}{L^{'}} \sum_{α} f^{α m_{2}} J_{2}^{α}, \\ μ_{1111} = \frac{1}{L^{'}} \sum_{α} f^{αg} J_{3}^{α} . \end{matrix} \end{matrix}

(10)

Equation (10) thus defines the macro-scale stress measures in terms of intergranular force measures and their associated geometry moments, where $J_{1}^{n}$ , $J_{2}^{n}$ , and $J_{3}^{n}$ correspond to the geometry moments for the n and p grain pair interacting within the LE.

To formulate the micro-scale constitutive equations, which link micro-scale kinematic measures to their conjugate intergranular force measures, we propose the following quadratic form for the micro-scale deformation energy $W^{α}$ of the α-th grain pair. This formulation implies a linear elastic response at the micro-level:

\begin{matrix} \begin{matrix} W^{α} = \frac{1}{2} (K^{αM} {(δ^{αM})}^{2} + K^{α m_{1}} {(δ^{α m_{1}})}^{2} + K^{α m_{2}} {(δ^{α m_{2}})}^{2} + K^{αg} {(δ^{αg})}^{2}) + \\ + 2 (K^{αM m_{1}} δ^{αM} δ^{α m_{1}} + K^{αM m_{2}} δ^{αM} δ^{α m_{2}} + K^{αMg} δ^{αM} δ^{αg} + \\ + K^{α m_{1} m_{2}} δ^{α m_{1}} δ^{α m_{2}} + K^{α m_{1} g} δ^{α m_{1}} δ^{αg} + K^{α m_{2} g} δ^{α m_{2}} δ^{αg}) \end{matrix} \end{matrix}

(11)

Based on equation (11), we identify sixteen distinct linear mechanisms contributing to the deformation of a grain pair in contact. Each mechanism is quadratic in form, where $K^{αij}$ , with $i, j = M, m_{1}, m_{2}, g$ , represent the intergranular stiffnesses associated with their corresponding interaction mechanisms. The intergranular forces introduced in equation (9) are then obtained by differentiating equation (11) with respect to the corresponding kinematic measures:

\begin{matrix} \begin{matrix} f^{αM} = \frac{\partial W^{α}}{\partial δ^{αM}} = K^{αM} δ^{αM} + K^{αM m_{1}} δ^{α m_{1}} + K^{αM m_{2}} δ^{α m_{2}} + K^{αMg} δ^{αg}, \\ f^{α m_{1}} = \frac{\partial W^{α}}{\partial δ^{α m_{1}}} = K^{α m_{1}} δ^{α m_{1}} + K^{αM m_{1}} δ^{αM} + K^{α m_{1} m_{2}} δ^{α m_{2}} + K^{α m_{1} g} δ^{αg}, \\ f^{α m_{2}} = \frac{\partial W^{α}}{\partial δ^{α m_{2}}} = K^{α m_{2}} δ^{α m_{2}} + K^{αM m_{2}} δ^{αM} + K^{α m_{1} m_{2}} δ^{α m_{1}} + K^{α m_{2} g} δ^{αg}, \\ f^{αg} = \frac{\partial W^{α}}{\partial δ^{αg}} = K^{αg} δ^{αg} + K^{αMg} δ^{αM} + K^{α m_{1} g} δ^{α m_{1}} + K^{α m_{2} g} δ^{α m_{2}}, \end{matrix} \end{matrix}

(12)

Finally, by substituting these expressions for the micro-forces (12) into the averaged definitions of the macro-stresses (10), we derive the macroscopic constitutive relations. These relations explicitly describe the linear elastic relationship between the macro-scale stresses and the kinematic variables for a degree 2 micromorphic rod:

\begin{matrix} τ_{11} = C^{M} {\bar{ϕ}}_{, x} + C^{M m_{1}} ({\bar{ϕ}}_{, x} - ψ_{11}) + C^{M m_{2}} (- ψ_{111} + ψ_{11, x}) + C^{Mg} ψ_{111, x} \\ σ_{11} = C^{M m_{1}} {\bar{ϕ}}_{, x} + C^{m_{1}} ({\bar{ϕ}}_{, x} - ψ_{11}) + C^{m_{1} m_{2}} (- ψ_{111} + ψ_{11, x}) + C^{m_{1} g} ψ_{111, x} \\ σ_{111} = C^{M m_{2}} {\bar{ϕ}}_{, x} + C^{m_{1} m_{2}} ({\bar{ϕ}}_{, x} - ψ_{11}) + C^{m_{2}} (- ψ_{111} + ψ_{11, x}) + C^{m_{2} g} ψ_{111, x} \\ μ_{1111} = C^{Mg} {\bar{ϕ}}_{, x} + C^{m_{1} g} ({\bar{ϕ}}_{, x} - ψ_{11}) + C^{m_{2} g} (- ψ_{111} + ψ_{11, x}) + C^{g} ψ_{111, x} \end{matrix}

(13)

where the macro-scale stiffnesses $C^{ij}$ , $i, j = M, m_{1}, m_{2}, g$ , are expressed as

\begin{matrix} \begin{matrix} C^{M} = \frac{1}{L^{'}} \sum_{α} K^{αM} J_{1} J_{1} C^{M m_{1}} = \frac{1}{L^{'}} \sum_{α} K^{αM m_{1}} J_{1} J_{1} \\ C^{M m_{2}} = \frac{1}{L^{'}} \sum_{α} K^{αM m_{2}} J_{1} J_{2} C^{Mg} = \frac{1}{L^{'}} \sum_{α} K^{αMg} J_{1} J_{3} C^{m_{1} m_{2}} = \frac{1}{L^{'}} \sum_{α} K^{α m_{1} m_{2}} J_{1} J_{2} \\ C^{m_{1} g} = \frac{1}{L^{'}} \sum_{α} K^{α m_{1} g} J_{1} J_{3} C^{m_{2} g} = \frac{1}{L^{'}} \sum_{α} K^{α m_{2} g} J_{2} J_{3} \\ C^{m_{2}} = \frac{1}{L^{'}} \sum_{α} K^{α m_{2}} J_{2} J_{2} C^{m_{1}} = \frac{1}{L^{'}} \sum_{α} K^{α m_{1}} J_{1} J_{1} C^{g} = \frac{1}{L^{'}} \sum_{α} K^{αg} J_{3} J_{3} \end{matrix} \end{matrix}

(14)

The superscripts in the macro-scale stiffnesses from (14) carry specific physical meanings. The superscript M denotes the stiffnesses related to macro-scale deformation, while $m_{1}$ and $m_{2}$ denote the stiffnesses associated with first-order (relative displacement) and second-order micro-fluctuations, respectively. The superscript g denotes the stiffnesses corresponding to higher-order micro-deformation. Combined superscripts such as $M m_{1}$ , $M m_{2}$ , Mg, $m_{1} m_{2}$ , $m_{1} g$ , and $m_{2} g$ represent the coupling (cross-linking) stiffnesses between the respective macro- and micro-scale deformation modes.

It is important to note that all stiffness measures defined in equation (14) naturally incorporate characteristic lengths. These multiple length scales are a direct consequence of the assumed kinematic field of grain motion within the LE. Consequently, for the analysis of static problems in this second-degree micromorphic model, the material behavior is characterized by several intrinsic static length scales, which will be explicitly defined as dimensionless parameters in subsequent sections dedicated to the dimensionless form of the equations.

3. Governing equations and BCs

3.1. Variational formulation

The principle of minimum potential energy is used to obtain the equilibrium equations for the 1D granular rod. This principle requires that the total potential energy of the system be minimized in equilibrium and is expressed as

\begin{matrix} \begin{matrix} δŴ - δ Ŵ^{ext} = 0, \end{matrix} \end{matrix}

(15)

where δ is the variation symbol, $δŴ = \int_{L} δWdx$ represents the variation of total macro-scale deformation energy, which is expressed as

\begin{matrix} \begin{matrix} δŴ = - \int_{L} {(τ_{11} + σ_{11})}_{, x} δ \bar{ϕ} dx - \int_{L} (σ_{111, x} + σ_{11}) δ ψ_{11} dx - \int_{L} (σ_{111} + μ_{1111, x}) δ ψ_{111} dx + \\ + (τ_{11} + σ_{11}) δ \bar{ϕ} |_{0}^{L} + σ_{111} δ ψ_{11} |_{0}^{L} + μ_{1111} δ ψ_{111} |_{0}^{L} . \end{matrix} \end{matrix}

(16)

Finally, the term $δ Ŵ^{ext}$ in equation (15) corresponds to the variation of total external energy, which is defined as

\begin{matrix} δ Ŵ^{ext} = \int_{L} f_{1} δ \bar{ϕ} dx + \int_{L} f_{2} δ ψ_{11} dx + \int_{L} f_{3} δ ψ_{111} dx + h_{1} δ \bar{ϕ} |_{0}^{L} + h_{2} δ ψ_{11} |_{0}^{L} + h_{3} δ ψ_{111} |_{0}^{L} . \end{matrix}

(17)

In (17), we define the external contributions: $f_{1}$ as the non-contact body force per unit length, $f_{2}$ as the non-contact body double force per unit length, and $f_{3}$ as the non-contact body triple force per unit length. Correspondingly, $h_{1}$ represents the contact traction, $h_{2}$ represents the contact double traction, and $h_{3}$ represents the contact triple traction. By integrating (16) and (17) into the variational principle expressed in (15), the full variational statement is obtained as:

\begin{matrix} \begin{matrix} \int_{L} (- {(σ_{11} + τ_{11})}_{x} δ \bar{ϕ} - (σ_{111, x} + σ_{11}) δ ψ_{11} - (μ_{1111, x} + σ_{111}) δ ψ_{111}) dx + ((σ_{11} + τ_{11}) δ \bar{ϕ} + σ_{111} δ ψ_{11} + μ_{1111} δ ψ_{111}) |_{0}^{L} \\ - \int_{L} (f_{1} δ \bar{ϕ} + f_{2} δ ψ_{11} + f_{3} δ ψ_{111}) dx + (h_{1} δ {\bar{ϕ}}_{1} + h_{2} δ ψ_{11} + h_{3} δ ψ_{111}) |_{0}^{L} = 0 \end{matrix} \end{matrix}

(18)

3.2. Equilibrium equations in strong form

From the comprehensive variational statement in (18), the full system of equilibrium equations and the corresponding BCs can be rigorously derived.

Specifically, by applying the fundamental lemma of variational calculus, we obtain the strong form of the equilibrium equations. This is achieved by setting the coefficients of the arbitrary variations $\bar{ϕ}$ , $ψ_{11}$ , and $ψ_{111}$ within the integral term to zero. Concurrently, setting the boundary terms to zero yields the necessary BCs.

The resulting strong form of the equilibrium equations for the kinematic variables $\bar{ϕ}$ , $ψ_{11}$ , and $ψ_{111}$ is:

\begin{matrix} {\begin{matrix} {(τ_{11} + σ_{11})}_{, x} + f_{1} = 0 \\ σ_{11} + σ_{111, x} + f_{2} = 0 \\ σ_{111} + μ_{1111, x} + f_{3} = 0 \end{matrix} \end{matrix}

(19)

And the corresponding natural BCs, derived from the variational principle, are:

\begin{matrix} {\begin{matrix} (τ_{11} + σ_{11} - h_{1}) δ \bar{ϕ} |_{0}^{L} = 0 \\ (σ_{111} - h_{2}) δ ψ_{11} |_{0}^{L} = 0 \\ (μ_{1111} - h_{3}) δ ψ_{111} |_{0}^{L} = 0 \end{matrix} \end{matrix}

(20)

Substituting the constitutive equations, as detailed in (13), into this equilibrium system (19), and assuming that the macro-scale stiffnesses are spatially invariant (i.e., independent of the coordinate x), we obtain a coupled system of differential equations for the primary kinematic variables $\bar{ϕ}$ , $ψ_{11}$ , and $ψ_{111}$ . This system is presented as:

\begin{matrix} \begin{matrix} (C^{M} + 2 C^{M m_{1}} + C^{m_{1}}) {\bar{ϕ}}_{, xx} + (C^{m_{1} m_{2}} + C^{M m_{2}}) ψ_{11, xx} + (C^{m_{1} g} + C^{Mg}) ψ_{111, xx} \\ - (C^{m_{1}} + C^{M m_{1}}) ψ_{11, x} - (C^{m_{1} m_{2}} + C^{M m_{2}}) ψ_{111, x} + f_{1} = 0 \\ (C^{M m_{2}} + C^{m_{1} m_{2}}) {\bar{ϕ}}_{, xx} + C^{m_{2}} ψ_{11, xx} + C^{m_{2} g} ψ_{111, xx} + (C^{M m_{1}} + C^{m_{1}}) {\bar{ϕ}}_{, x} + \\ + (C^{m_{1} g} - C^{m_{2}}) ψ_{111, x} - C^{m_{1}} ψ_{11} - C^{m_{1} m_{2}} ψ_{111} + f_{2} = 0 \\ (C^{Mg} + C^{m_{1} g}) {\bar{ϕ}}_{, xx} + C^{m_{2} g} ψ_{11, xx} + C^{g} ψ_{111, xx} + (C^{M m_{2}} + C^{m_{1} m_{2}}) {\bar{ϕ}}_{, x} \\ + (- C^{m_{1} g} + C^{m_{2}}) ψ_{11, x} + (C^{m_{2} g} - C^{m_{2} g}) ψ_{111, x} - C^{m_{1} m_{2}} ψ_{11} - C^{m_{2}} ψ_{111} + f_{3} = 0 \end{matrix} \end{matrix}

(21)

Similarly, the BCs in (20) are transformed by substituting the constitutive relations, yielding the explicit form:

\begin{matrix} \begin{matrix} [(C^{M} + 2 C^{M m_{1}} + C^{m_{1}}) {\bar{ϕ}}_{, x} + (C^{m_{1} m_{2}} + C^{M m_{2}}) ψ_{11, x} + (C^{m_{1} g} + C^{Mg}) ψ_{111, x} \\ - (C^{m_{1}} + C^{M m_{1}}) ψ_{11} - (C^{m_{1} m_{2}} + C^{M m_{2}}) ψ_{111} - h_{1}] δ \bar{ϕ} |_{0}^{L} = 0 \\ [(C^{M m_{2}} + C^{m_{1} m_{2}}) {\bar{ϕ}}_{, x} + C^{m_{2}} ψ_{11, x} + C^{m_{2} g} ψ_{111, x} \\ - C^{m_{1} m_{2}} ψ_{11} - C^{m_{2}} ψ_{111} - h_{2}] δ ψ_{11} |_{0}^{L} = 0 \\ [(C^{Mg} + C^{m_{1} g}) {\bar{ϕ}}_{, x} + C^{m_{2} g} ψ_{11, x} + C^{g} ψ_{111, x} - C^{m_{1} g} ψ_{11} - C^{m_{2} g} ψ_{111} - h_{3}] δ ψ_{111} |_{0}^{L} = 0 \end{matrix} \end{matrix}

(22)

To simplify the representation of this complex system, we introduce a set of coefficients based on the macro-scale stiffnesses:

\begin{matrix} \begin{matrix} D_{1} = C^{M} + C^{m_{1}} + 2 C^{M m_{1}}, & D_{5} = C^{m_{2}}, & D_{9} = - C^{m_{1} m_{2}} \\ D_{2} = C^{m_{1} m_{2}} + C^{M m_{2}}, & D_{6} = C^{m_{2} g} & D_{10} = C^{g} \\ D_{3} = C^{m_{1} g} + C^{Mg} & D_{7} = C^{m_{1} g} - C^{m_{2}} & D_{11} = - C^{m_{1} g} \\ D_{4} = - C^{m_{1}} - C^{M m_{1}} & D_{8} = - C^{m_{1}} \end{matrix} \end{matrix}

(23)

Utilizing these newly defined coefficients, the system of differential equations (21) for the problem domain $0 \leq x \leq L$ and the associated BCs (22) can be expressed in a more compact form.

The compact form of the equilibrium equations is:

\begin{matrix} \begin{matrix} D_{1} {\bar{ϕ}}_{, xx} + D_{2} ψ_{11, xx} + D_{3} ψ_{111, xx} + D_{4} ψ_{11, x} - D_{2} ψ_{111, x} + f_{1} = 0 \\ D_{2} {\bar{ϕ}}_{, xx} + D_{5} ψ_{11, xx} + D_{6} ψ_{111, xx} - D_{4} {\bar{ϕ}}_{, x} + D_{7} ψ_{111, x} + D_{8} ψ_{11} + D_{9} ψ_{111} + f_{2} = 0 \\ D_{3} {\bar{ϕ}}_{, xx} + D_{6} ψ_{11, xx} + D_{10} ψ_{111, xx} + D_{2} {\bar{ϕ}}_{, x} - D_{7} ψ_{11, x} + D_{9} ψ_{11} - D_{5} ψ_{111} + f_{3} = 0 \end{matrix} \end{matrix}

(24)

And the compact form of the BCs is:

\begin{matrix} \begin{matrix} [D_{1} {\bar{ϕ}}_{, x} + D_{2} ψ_{11, x} + D_{3} ψ_{111, x} + D_{4} ψ_{11} - D_{2} ψ_{111} - h_{1}] δ \bar{ϕ} |_{0}^{L} = 0 \\ [D_{2} {\bar{ϕ}}_{, x} + D_{5} ψ_{11, x} + D_{6} ψ_{111, x} + D_{9} ψ_{11} - D_{5} ψ_{111} - h_{2}] δ ψ_{11} |_{0}^{L} = 0 \\ [D_{3} {\bar{ϕ}}_{, x} + D_{6} ψ_{11, x} + D_{10} ψ_{111, x} + D_{11} ψ_{11} - D_{6} ψ_{111} - h_{3}] δ ψ_{111} |_{0}^{L} = 0 \end{matrix} \end{matrix}

(25)

4. Solution for the static problem

4.1. Operator formulation

For a more compact and illustrative representation of the system of differential equations and BCs derived in the previous section, an operator form is employed. We introduce column vectors for the kinematic variables $u (x)$ , non-contact body forces $f (x)$ , and contact tractions

\begin{matrix} u (x) = (\begin{matrix} \bar{ϕ} \\ ψ_{11} \\ ψ_{111} \end{matrix}), f = (\begin{matrix} f_{1} \\ f_{2} \\ f_{3} \end{matrix}), h = (\begin{matrix} h_{1} \\ h_{2} \\ h_{3} \end{matrix}) \end{matrix}

(26)

Next, we define the differential operators $L_{1} = (L_{1 (ij)})$ for the governing equations and $L_{2} = (L_{2 (ij)})$ for the BCs, where $i, j = 1, 2, 3$ :

\begin{matrix} \begin{matrix} L_{1} = (\begin{matrix} D_{1} \partial_{x}^{2} & D_{4} \partial_{x} + D_{2} \partial_{x}^{2} & - D_{2} \partial_{x} + D_{3} \partial_{x}^{2} \\ - D_{4} \partial_{x} + D_{2} \partial_{x}^{2} & D_{8} + D_{5} \partial_{x}^{2} & D_{9} + D_{7} \partial_{x} + D_{6} \partial_{x}^{2} \\ D_{2} \partial_{x} + D_{3} \partial_{x}^{2} & D_{9} - D_{7} \partial_{x} + D_{6} \partial_{x}^{2} & - D_{5} + D_{10} \partial_{x}^{2} \end{matrix}) \end{matrix} \end{matrix}

(27)

\begin{matrix} L_{2} = (\begin{matrix} D_{1} \partial_{x} & D_{4} + D_{2} \partial_{x} & - D_{2} + D_{3} \partial_{x} \\ D_{2} \partial_{x} & D_{9} + D_{5} \partial_{x} & - D_{5} + D_{6} \partial_{x} \\ D_{3} \partial_{x} & D_{11} + D_{6} \partial_{x} & - D_{6} + D_{10} \partial_{x} \end{matrix}) \end{matrix}

(28)

Using these introduced operators, the original system of differential equations (24) and BCs (25) for the problem domain $0 \leq x \leq L$ can be written in a compact matrix form as follows:

\begin{matrix} \begin{matrix} L_{1} u + f = 0 \end{matrix} \end{matrix}

(29)

\begin{matrix} (L_{2} u - h) |_{0}^{L} = 0 \end{matrix}

(30)

4.2. Dimensionless form of the governing equations

To simplify the analysis, reduce the number of independent physical parameters, and highlight size effects, it is advantageous to transform the governing equations into a dimensionless form. For this purpose, we introduce the following dimensionless variables:

\begin{matrix} \begin{matrix} \tilde{x} = \frac{x}{L}, \tilde{\bar{ϕ}} = \frac{\bar{ϕ}}{L}, {\tilde{ψ}}_{11} = ψ_{11}, {\tilde{ψ}}_{111} = L ψ_{111} \end{matrix} \end{matrix}

(31)

The dimensionless spatial domain of the problem is, therefore, defined as $0 \leq \tilde{x} \leq 1$ . Derivatives with respect to the physical coordinate x can be transformed into derivatives with respect to the dimensionless coordinate $\tilde{x}$ as follows:

\begin{matrix} {\bar{ϕ}}_{, x} = {\tilde{\bar{ϕ}}}_{, \tilde{x}}, {\bar{ϕ}}_{, xx} = \frac{1}{L} {\tilde{\bar{ϕ}}}_{, \tilde{x} \tilde{x}}, ψ_{11, x} = \frac{1}{L} {\tilde{ψ}}_{11, \tilde{x}}, ψ_{11, xx} = \frac{1}{L^{2}} {\tilde{ψ}}_{11, \tilde{x} \tilde{x}} \\ ψ_{111, x} = \frac{1}{L^{2}} {\tilde{ψ}}_{111, \tilde{x}}, ψ_{111, xx} = \frac{1}{L^{3}} {\tilde{ψ}}_{111, \tilde{x} \tilde{x}} \end{matrix}

Using these transformations, we also define a set of dimensionless material parameters and characteristic length scales that inherently reflect the influence of the granular microstructure:

\begin{matrix} \begin{matrix} l_{1} = \frac{1}{L} \frac{C^{M m_{2}}}{C^{M}}, l_{2} = \frac{1}{L} \sqrt{\frac{C^{m_{2}}}{C^{M}}}, l_{3} = \frac{1}{L} \sqrt[3]{\frac{C^{m_{2} g}}{C^{M}}}, l_{4} = \frac{1}{L} \sqrt[4]{\frac{C^{g}}{C^{M}}}, \\ α_{m_{1}} = \frac{C^{m_{1}}}{C^{M}}, α_{M m_{1}} = \frac{C^{M m_{1}}}{C^{M}}, β_{m_{1} m_{2}} = \frac{C^{m_{1} m_{2}}}{C^{M m_{2}}}, γ_{m_{1} g} = \frac{C^{m_{1} g}}{C^{m_{2}}}, γ_{Mg} = \frac{C^{Mg}}{C^{m_{2}}}, \\ \frac{C^{m_{1} m_{2}}}{C^{M}} = l_{1} L β_{m_{1} m_{2}}, \frac{C^{m_{1} g}}{C^{M}} = \frac{C^{m_{1} g}}{C^{m_{2}}} \frac{C^{m_{2}}}{C^{M}} = l_{2}^{2} L^{2} γ_{m_{1} g}, \frac{C^{Mg}}{C^{M}} = l_{2}^{2} L^{2} γ_{Mg} \end{matrix} \end{matrix}

(32)

Furthermore, we introduce the dimensionless parameter n, defined as the ratio of the macroscopic rod length L to the characteristic length of the LE $L^{'}$ : $n = L ∕ L^{'}$ . This parameter, previously introduced in [16], quantifies the scale separation between the macro- and microstructure. For numerical demonstrations, a representative value for n is typically chosen, as exemplified in [16]. The dimensionless micro-coordinate ${\tilde{x}}^{'}$ , representing the position within a LE relative to its COM, is defined as ${\tilde{x}}^{'} = x^{'} ∕ L^{'}$ , varying from $- 0.5$ to 0.5.

These dimensionless parameters are crucial for understanding the system’s behavior. The parameters $l_{1}$ , $l_{2}$ , $l_{3}$ , and $l_{4}$ represent the dimensionless static length scales that naturally arise from various stiffness components. Larger values of these parameters indicate more pronounced higher-gradient effects, suggesting a more significant microstructural influence relative to the overall rod length. These parameters also reflect the size effect, where the relative scale of the microstructure to the macroscopic rod length impacts the material response. The parameter $α_{m_{1}}$ represents the ratio of micro-scale (relative deformation) stiffness to macro-scale stiffness, indicating the relative importance of micro-deformations. Similarly, $α_{M m_{1}}$ represents the ratio of the coupling (cross-linking) stiffness between macro- and first-order micro-scale deformations to the macro-scale stiffness, quantifying the interaction strength between these scales. The parameter $β_{m_{1} m_{2}}$ is the ratio of the coupling stiffness between first- and second-order micro-scale fluctuations, highlighting the interrelation between different levels of micro-deformation. Finally, $γ_{m_{1} g}$ and $γ_{Mg}$ characterize further couplings involving gradient stiffnesses, providing insight into how higher-order microstructural effects are linked across different kinematic quantities.

By applying these dimensionless variables and parameters, the system of governing equations (24) and BCs (25) can be reformulated in a dimensionless form. The dimensionless equilibrium equations are given as:

\begin{matrix} \begin{matrix} (1 + 2 α_{M m_{1}} + α_{m_{1}}) {\tilde{\bar{ϕ}}}_{, \tilde{x} \tilde{x}} + l_{1} (1 + β_{m_{1} m_{2}}) {\tilde{ψ}}_{11, \tilde{x} \tilde{x}} + l_{2}^{2} (γ_{m_{1} g} + γ_{Mg}) {\tilde{ψ}}_{111, \tilde{x} \tilde{x}} \\ - (α_{m_{1}} + α_{M m_{1}}) {\tilde{ψ}}_{11, \tilde{x}} - l_{1} (1 + β_{m_{1} m_{2}}) {\tilde{ψ}}_{111, \tilde{x}} + \underset{̲}{f_{1}} = 0 \\ l_{1} (1 + β_{m_{1} m_{2}}) {\tilde{\bar{ϕ}}}_{, \tilde{x} \tilde{x}} + l_{2}^{2} {\tilde{ψ}}_{11, \tilde{x} \tilde{x}} + l_{3}^{3} {\tilde{ψ}}_{111, \tilde{x} \tilde{x}} + (α_{M m_{1}} + α_{m_{1}}) {\tilde{\bar{ϕ}}}_{, \tilde{x}} \\ + l_{2}^{2} (γ_{m_{1} g} - 1) {\tilde{ψ}}_{111, \tilde{x}} - α_{m_{1}} {\tilde{ψ}}_{11} - l_{1} β_{m_{1} m_{2}} {\tilde{ψ}}_{111} + \underset{̲}{f_{2}} = 0 \\ l_{2}^{2} (γ_{Mg} + γ_{m_{1} g}) {\tilde{\bar{ϕ}}}_{, \tilde{x} \tilde{x}} + l_{3}^{3} {\tilde{ψ}}_{11, \tilde{x} \tilde{x}} + l_{4}^{4} {\tilde{ψ}}_{111, \tilde{x} \tilde{x}} + l_{1} (1 + β_{m_{1} m_{2}}) {\tilde{\bar{ϕ}}}_{, \tilde{x}} \\ + l_{2}^{2} (- γ_{m_{1} g} + 1) {\tilde{ψ}}_{11, \tilde{x}} - l_{1} β_{m_{1} m_{2}} {\tilde{ψ}}_{11} - l_{2}^{2} {\tilde{ψ}}_{111} + \underset{̲}{f_{3}} = 0 \end{matrix} \end{matrix}

(33)

where $\underset{̲}{f_{1}} = \frac{f_{1} L}{C^{M}}$ is the dimensionless non-contact body force per unit length, $\underset{̲}{f_{2}} = \frac{f_{2}}{C^{M}}$ is the dimensionless non-contact double force per unit length, and $\underset{̲}{f_{3}} = \frac{f_{3}}{C^{M} L}$ is the dimensionless non-contact triple force per unit area. The dimensionless spatial domain of the problem is $0 \leq \tilde{x} \leq 1$ .

The dimensionless form of the BCs (34) is:

\begin{matrix} \begin{matrix} (\underset{̲}{h_{1}} - [(1 + 2 α_{M m_{1}} + α_{m_{1}}) {\tilde{\bar{ϕ}}}_{, \tilde{x}} + l_{1} (1 + β_{m_{1} m_{2}}) {\tilde{ψ}}_{11, \tilde{x}} \\ + l_{2}^{2} (γ_{m_{1} g} + γ_{Mg}) {\tilde{ψ}}_{111, \tilde{x}} - (α_{m_{1}} + α_{M m_{1}}) {\tilde{ψ}}_{11} - l_{1} (1 + β_{m_{1} m_{2}}) {\tilde{ψ}}_{111}]) |_{0}^{L} = 0 \\ (\underset{̲}{h_{2}} - [l_{1} (1 + β_{m_{1} m_{2}}) {\tilde{\bar{ϕ}}}_{, \tilde{x}} + l_{2}^{2} {\tilde{ψ}}_{11, \tilde{x}} + l_{3}^{3} {\tilde{ψ}}_{111, \tilde{x}} - l_{1} β_{m_{1} m_{2}} {\tilde{ψ}}_{11} - l_{2}^{2} {\tilde{ψ}}_{111}]) |_{0}^{L} = 0 \\ (\underset{̲}{h_{3}} - [l_{2}^{2} (γ_{Mg} + γ_{m_{1} g}) {\tilde{\bar{ϕ}}}_{, \tilde{x}} + l_{3}^{3} {\tilde{ψ}}_{11, \tilde{x}} + l_{4}^{4} {\tilde{ψ}}_{111, \tilde{x}} - l_{2}^{2} γ_{m_{1} g} {\tilde{ψ}}_{11} - l_{3}^{3} {\tilde{ψ}}_{111}]) |_{0}^{L} = 0 \end{matrix} \end{matrix}

(34)

where $\underset{̲}{h_{1}} = \frac{h_{1}}{C^{M}}$ is the dimensionless contact traction, $\underset{̲}{h_{2}} = \frac{h_{2} L}{C^{M}}$ is the dimensionless contact double traction per unit length, and $\underset{̲}{h_{3}} = \frac{h_{3} L^{2}}{C^{M}}$ is the dimensionless contact triple traction per unit area, respectively.

Finally, to fully characterize the material’s deformation, the total dimensionless displacement ${\tilde{ϕ}}_{total}$ of a point within the microstructure must be considered. Building upon the kinematic expansion in equation (3) and applying the dimensionless variables and parameters, ${\tilde{ϕ}}_{total}$ explicitly depends on both the dimensionless macro-coordinate $\tilde{x}$ and the dimensionless micro-coordinate ${\tilde{x}}^{'}$ :

\begin{matrix} {\tilde{ϕ}}_{total} (\tilde{x}, {\tilde{x}}^{'}) = \tilde{\bar{ϕ}} (\tilde{x}) + \frac{1}{n} {\tilde{ψ}}_{11} (\tilde{x}) {\tilde{x}}^{'} + \frac{1}{n^{2}} {\tilde{ψ}}_{111} (\tilde{x}) {({\tilde{x}}^{'})}^{2} + \frac{1}{n^{3}} {\tilde{ψ}}_{111, \tilde{x}} (\tilde{x}) {({\tilde{x}}^{'})}^{3} \end{matrix}

(35)

here, ${\tilde{x}}^{'} = x^{'} ∕ L^{'}$ varies from $- 0.5$ to 0.5 within each LE. It follows from (35) that the micro-kinematic contributions enter the total displacement through terms scaled by the factors $1 ∕ n$ , $1 ∕ n^{2}$ , and $1 ∕ n^{3}$ . Since $n = L ∕ L^{'}$ denotes the ratio between the macroscopic length and the characteristic size of the LE, these scaling factors naturally reduce the relative influence of the micro-kinematic variables ${\tilde{ψ}}_{11}$ , ${\tilde{ψ}}_{111}$ , and ${\tilde{ψ}}_{111, \tilde{x}}$ on the total displacement as n increases. This attenuation is a direct consequence of the adopted non-dimensionalization and does not imply that the micro-kinematic fields themselves become small; rather, their contribution to the macroscopic displacement field becomes less pronounced when the macroscopic length scale is much larger than the microstructural one. In the limit $n \to \infty$ , the continuum recovers classical (non-microstructural) behavior, with ${\tilde{ϕ}}_{total} \to \bar{ϕ}$ .

An important observation often overlooked in classical approaches concerns the potential size dependence of inclusion stiffness. In many traditional formulations, the material stiffness of inclusions is assumed to remain constant regardless of their size. However, in realistic scenarios, the stiffness properties of inclusions may depend on their characteristic dimensions and, consequently, on the scale parameter $n = L ∕ L^{'}$ . Such size-dependence would manifest through variations in the material constants ( $C^{m_{1}}$ , $C^{m_{2}}$ , $C^{M m_{1}}$ , $C^{m_{1} m_{2}}$ , etc.) as n changes. Since the dimensionless parameters ( $l_{1}$ , $l_{2}$ , $l_{3}$ , $l_{4}$ , $α_{m_{1}}$ , $α_{M m_{1}}$ , $β_{m_{1} m_{2}}$ , etc.) are derived from these material constants, they too would vary with n. Consequently, the micro-kinematic fields $\tilde{\bar{ϕ}}$ , ${\tilde{ψ}}_{11}$ , and ${\tilde{ψ}}_{111}$ themselves could change with n, not only through the geometric scaling factors $1 ∕ n$ , $1 ∕ n^{2}$ , $1 ∕ n^{3}$ appearing explicitly in (35) but also through the implicit n-dependence of the material parameters governing their evolution. While this size-dependent stiffness effect represents a physically relevant phenomenon worthy of investigation, the exploration of this aspect is beyond the scope of the present work.

5. Solution for the static problem

This section is dedicated to developing the analytical solution for the static deformation of the 1D granular rod, described by our second-degree micromorphic theory. While previous studies, such as [16], focused on a degree 1 micromorphic theory and employed methods like reducing the system to a single higher-order ordinary differential equation (ODE), our work presents a more generalized framework.

For the coupled system arising from our degree 2 micromorphic model, we utilize the state-space representation as a robust and systematic analytical approach. We will outline this method, discuss the fundamental conditions for the existence of a physically meaningful solution, and present the general solution. Specific case studies with various BCs will be explored in detail in the subsequent section.

5.1. Analytical solution using state-space method

The dimensionless equilibrium equations, previously derived in equation (33), constitute a set of three coupled second-order ODEs for $\tilde{\bar{ϕ}}, {\tilde{ψ}}_{11}, {\tilde{ψ}}_{111}$ . To facilitate a comprehensive analytical solution, especially for complex coupled systems, the state-space representation method is adopted. This approach transforms the higher-order coupled system into a larger system of first-order differential equations.

To simplify algebraic manipulations, we introduce a compact notation for the coefficients, which will be used in the subsequent representation of the system:

\begin{matrix} a_{1} & = (1 + 2 α_{M m_{1}} + α_{m_{1}}) & b_{1} & = l_{1} (1 + β_{m_{1} m_{2}}) & c_{1} & = l_{2}^{2} (γ_{m_{1} g} + γ_{Mg}) \\ a_{2} & = l_{1} (1 + β_{m_{1} m_{2}}) & b_{2} & = l_{2}^{2} & c_{2} & = l_{3}^{3} \\ a_{3} & = l_{2}^{2} (γ_{Mg} + γ_{m_{1} g}) & b_{3} & = l_{3}^{3} & c_{3} & = l_{4}^{4} \\ d_{1} & = (α_{m_{1}} + α_{M m_{1}}) & e_{1} & = l_{1} (1 + β_{m_{1} m_{2}}) \\ d_{2} & = (α_{M m_{1}} + α_{m_{1}}) & e_{2} & = l_{2}^{2} (γ_{m_{1} g} - 1) & g_{2} & = α_{m_{1}} & h_{2} & = l_{1} β_{m_{1} m_{2}} \\ d_{3} & = l_{1} (1 + β_{m_{1} m_{2}}) & e_{3} & = l_{2}^{2} (- γ_{m_{1} g} + 1) & g_{3} & = l_{1} β_{m_{1} m_{2}} & h_{3} & = l_{2}^{2} \end{matrix}

Using this simplified notation, the system of dimensionless governing equations can be succinctly rewritten as:

\begin{matrix} \begin{matrix} a_{1} {\tilde{\bar{ϕ}}}_{, \tilde{x} \tilde{x}} + b_{1} {\tilde{ψ}}_{11, \tilde{x} \tilde{x}} + c_{1} {\tilde{ψ}}_{111, \tilde{x} \tilde{x}} - d_{1} {\tilde{ψ}}_{11, \tilde{x}} - e_{1} {\tilde{ψ}}_{111, \tilde{x}} + f_{1} = 0 \\ a_{2} {\tilde{\bar{ϕ}}}_{, \tilde{x} \tilde{x}} + b_{2} {\tilde{ψ}}_{11, \tilde{x} \tilde{x}} + c_{2} {\tilde{ψ}}_{111, \tilde{x} \tilde{x}} + d_{2} {\tilde{\bar{ϕ}}}_{, \tilde{x}} + e_{2} {\tilde{ψ}}_{111, \tilde{x}} - g_{2} ψ_{11} - h_{2} ψ_{111} + f_{2} = 0 \\ a_{3} {\tilde{\bar{ϕ}}}_{, \tilde{x} \tilde{x}} + b_{3} {\tilde{ψ}}_{11, \tilde{x} \tilde{x}} + c_{3} {\tilde{ψ}}_{111, \tilde{x} \tilde{x}} + d_{3} {\tilde{\bar{ϕ}}}_{, \tilde{x}} + e_{3} {\tilde{ψ}}_{11, \tilde{x}} - g_{3} ψ_{11} - h_{3} ψ_{111} + f_{3} = 0 \end{matrix} \end{matrix}

(36)

5.1.1. State-space formulation

For the analytical solution of this coupled system, we transform it into a system of first-order ODEs, known as the state-space form. We define six state variables, based on the kinematic fields and their first derivatives:

\begin{matrix} y_{1} = \tilde{\bar{ϕ}} & y_{2} = {\tilde{\bar{ϕ}}}_{, \tilde{x}} & y_{3} = {\tilde{ψ}}_{11} & y_{4} = {\tilde{ψ}}_{11, \tilde{x}} & y_{5} = {\tilde{ψ}}_{111} & y_{6} = {\tilde{ψ}}_{111, \tilde{x}} \end{matrix}

The state vector is $Y (\tilde{x}) = {(y_{1}, y_{2}, y_{3}, y_{4}, y_{5}, y_{6})}^{T}$ . Its derivative with respect to $\tilde{x}$ is

Y^{'} (\tilde{x}) = {(y_{1}^{'}, y_{2}^{'}, y_{3}^{'}, y_{4}^{'}, y_{5}^{'}, y_{6}^{'})}^{T} = {({\tilde{\bar{ϕ}}}_{, \tilde{x}}, {\tilde{\bar{ϕ}}}_{, \tilde{x} \tilde{x}}, {\tilde{ψ}}_{11, \tilde{x}}, {\tilde{ψ}}_{11, \tilde{x} \tilde{x}}, {\tilde{ψ}}_{111, \tilde{x}}, {\tilde{ψ}}_{111, \tilde{x} \tilde{x}})}^{T}

The system (36) can then be rearranged and expressed in the standard state-space form $Y^{'} (\tilde{x}) = A Y (\tilde{x}) + G (\tilde{x})$ . First, we write the system as $S Y^{'} = K Y + F$ , where the matrices are constructed using the defined coefficients:

S = (\begin{matrix} 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & a_{1} & 0 & b_{1} & 0 & c_{1} \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & a_{2} & 0 & b_{2} & 0 & c_{2} \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & a_{3} & 0 & b_{3} & 0 & c_{3} \end{matrix}), K = (\begin{matrix} 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & d_{1} & 0 & e_{1} \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & - d_{2} & g_{2} & 0 & h_{2} & - e_{2} \\ 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & - d_{3} & g_{3} & - e_{3} & h_{3} & 0 \end{matrix}), F = (\begin{matrix} 0 \\ - f_{1} \\ 0 \\ - f_{2} \\ 0 \\ - f_{3} \end{matrix})

Multiplying by $S^{- 1}$ , we obtain $Y^{'} = A Y + G$ , where $A = S^{- 1} K$ and $G = S^{- 1} F$ . The inverse matrix $S^{- 1}$ is constructed from the inverse of the $3 \times 3$ submatrix

S_{sub} = (\begin{matrix} a_{1} & b_{1} & c_{1} \\ a_{2} & b_{2} & c_{2} \\ a_{3} & b_{3} & c_{3} \end{matrix})

whose elements are denoted $S_{ij}^{'}$ :

S^{- 1} = (\begin{matrix} 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & S_{22}^{'} & 0 & S_{24}^{'} & 0 & S_{26}^{'} \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & S_{42}^{'} & 0 & S_{44}^{'} & 0 & S_{46}^{'} \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & S_{62}^{'} & 0 & S_{64}^{'} & 0 & S_{66}^{'} \end{matrix})

The system matrix $A = S^{- 1} K$ is then explicitly derived as:

A = (\begin{matrix} 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & (- S_{24}^{'} d_{2} - S_{26}^{'} d_{3}) & (S_{24}^{'} g_{2} + S_{26}^{'} g_{3}) & (S_{22}^{'} d_{1} - S_{26}^{'} e_{3}) & (S_{24}^{'} h_{2} + S_{26}^{'} h_{3}) & (S_{22}^{'} e_{1} - S_{24}^{'} e_{2}) \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & (- S_{44}^{'} d_{2} - S_{46}^{'} d_{3}) & (S_{44}^{'} g_{2} + S_{46}^{'} g_{3}) & (S_{42}^{'} d_{1} - S_{46}^{'} e_{3}) & (S_{44}^{'} h_{2} + S_{46}^{'} h_{3}) & (S_{42}^{'} e_{1} - S_{44}^{'} e_{2}) \\ 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & (- S_{64}^{'} d_{2} - S_{66}^{'} d_{3}) & (S_{64}^{'} g_{2} + S_{66}^{'} g_{3}) & (S_{62}^{'} d_{1} - S_{66}^{'} e_{3}) & (S_{64}^{'} h_{2} + S_{66}^{'} h_{3}) & (S_{62}^{'} e_{1} - S_{64}^{'} e_{2}) \end{matrix})

5.1.2. General solution and conditions for existence

The general solution for the homogeneous system ( $G = 0$ ) is obtained by finding the eigenvalues and eigenvectors of the matrix $A$ .

The eigenvalues $λ_{i}$ ) are found by solving the characteristic equation $\det (A - λ I) = 0$ . This yields a polynomial of sixth degree in λ. These eigenvalues $λ_{1}, \dots, λ_{6}$ represent the fundamental exponential modes of the system. For each eigenvalue $λ_{i}$ , the corresponding eigenvector $e_{i}$ is determined by solving the linear system $(A - λ_{i} I) e_{i} = 0$ . In cases of repeated eigenvalues, generalized eigenvectors may be required to form a complete basis of solutions.

The homogeneous solution for the state vector is then given by $Y_{h} (\tilde{x}) = C_{1} e_{1} e^{λ_{1} \tilde{x}} + C_{2} e_{2} e^{λ_{2} \tilde{x}} + \dots + C_{6} e_{6} e^{λ_{6} \tilde{x}}$ , where $C_{1}, \dots, C_{6}$ are the arbitrary integration constants. A particular solution, if non-zero body forces are present ( $F \neq 0$ ), can be found using methods such as variation of parameters, or by assuming a form analogous to $G$ . For constant body forces, the particular solution often simplifies to a constant vector. The total solution to the system is the sum of the homogeneous and particular solutions: $Y (\tilde{x}) = Y_{h} (\tilde{x}) + Y_{p} (\tilde{x})$ .

The existence of a physically meaningful and stable solution is fundamentally dependent on the material parameters, which are embedded within the stiffness coefficients $C_{ij}$ and consequently in the dimensionless parameters ( $a_{k}, b_{k}, \dots$ ) and the eigenvalues ( $λ_{i}$ ). For a solution to be physically meaningful, the strain energy density of the material must be positive definite. This condition ensures that the material accumulates energy upon deformation and is inherently stable. The dimensionless strain energy density $\tilde{W}$ , derived from its dimensional counterpart, can be expressed in quadratic form as:

\begin{matrix} \tilde{W} = \frac{1}{2} {\tilde{q}}^{T} \tilde{C} \tilde{q} \end{matrix}

(37)

where

\tilde{q} = (\begin{matrix} {\tilde{\bar{ϕ}}}_{, \tilde{x}} \\ {\tilde{\bar{ϕ}}}_{, \tilde{x}} - {\tilde{ψ}}_{11} \\ {\tilde{ψ}}_{11, \tilde{x}} - {\tilde{ψ}}_{111} \\ {\tilde{ψ}}_{111, \tilde{x}} \end{matrix})

is the vector of dimensionless kinematic measures, and $\tilde{C}$ is the dimensionless stiffness matrix:

\tilde{C} = (\begin{matrix} 1 & α_{M m_{1}} & l_{1} & γ_{Mg} l_{2}^{2} \\ α_{M m_{1}} & α_{m_{1}} & β_{m_{1} m_{2}} l_{1} & γ_{m_{1} g} l_{2}^{2} \\ l_{1} & β_{m_{1} m_{2}} l_{1} & l_{2}^{2} & l_{3}^{3} \\ γ_{Mg} l_{2}^{2} & γ_{m_{1} g} l_{2}^{2} & l_{3}^{3} & l_{4}^{4} \end{matrix})

For $\tilde{W}$ to be positive definite, the stiffness matrix $\tilde{C}$ must be positive definite. According to Sylvester’s criterion, this requires all principal minors of $\tilde{C}$ to be strictly positive. The first principal minor $M_{1} = \det (1) = 1$ is always positive.

The second principal minor

M_{2} = \det (\begin{matrix} 1 & α_{M m_{1}} \\ α_{M m_{1}} & α_{m_{1}} \end{matrix}) = α_{m_{1}} - α_{M m_{1}}^{2}

must be positive. This condition ensures the stability of macro-scale and first-order micro-scale deformations and notably matches the positive definiteness criterion for the degree 1 micromorphic model discussed in Nejadsadeghi and Misra [16].

The third principal minor

M_{3} = \det (\begin{matrix} 1 & α_{M m_{1}} & l_{1} \\ α_{M m_{1}} & α_{m_{1}} & β_{m_{1} m_{2}} l_{1} \\ l_{1} & β_{m_{1} m_{2}} l_{1} & l_{2}^{2} \end{matrix})

must also be positive. Expanding this determinant yields: $l_{2}^{2} (α_{m_{1}} - α_{M m_{1}}^{2}) - l_{1}^{2} (α_{m_{1}} + β_{m_{1} m_{2}}^{2}) + 2 l_{1} α_{M m_{1}} β_{m_{1} m_{2}} l_{1} > 0$ . This condition extends stability requirements to include the influence of the second-order micro-fluctuation ( ${\tilde{ψ}}_{11, \tilde{x}} - {\tilde{ψ}}_{111}$ ) and associated parameters like $l_{1}$ and $β_{m_{1} m_{2}}$ .

The fourth principal minor $M_{4} = \det (\tilde{C})$ provides the most comprehensive condition, ensuring the overall stability of the system by considering all interactions between the four kinematic measures. This is the most comprehensive condition, ensuring the overall stability of the system. Expanding this determinant is algebraically complex but essential for a full stability assessment. The full expansion of this $4 \times 4$ determinant is algebraically complex and extensive, making its explicit representation here impractical for readability. However, it must be strictly positive to ensure the overall stability of the system, considering all interactions between the four kinematic measures, including all length scales ( $l_{1}$ , $l_{2}$ , $l_{3}$ , $l_{4}$ ) and dimensionless parameters. The full explicit expression for $M_{4}$ is:

\begin{matrix} M_{4} & = l_{4}^{4} [l_{2}^{2} (α_{m_{1}} - α_{M m_{1}}^{2}) - l_{1}^{2} (α_{m_{1}} + β_{m_{1} m_{2}}^{2}) + 2 l_{1} α_{M m_{1}} β_{m_{1} m_{2}} l_{1}] \\ - l_{3}^{3} [l_{1} (α_{m_{1}} γ_{Mg} l_{2}^{2} - α_{M m_{1}} γ_{m_{1} g} l_{2}^{2}) - β_{m_{1} m_{2}} l_{1} (γ_{Mg} l_{2}^{2} - α_{M m_{1}} l_{1}) \\ + l_{2}^{2} (α_{M m_{1}} β_{m_{1} m_{2}} l_{1} - α_{m_{1}} l_{1})] + γ_{Mg} l_{2}^{2} [l_{1} (l_{1} γ_{m_{1} g} l_{2}^{2} - α_{m_{1}} l_{3}^{3}) \\ - α_{M m_{1}} (β_{m_{1} m_{2}} l_{1} l_{3}^{3} - l_{2}^{2} γ_{m_{1} g} l_{2}^{2}) + l_{1} (α_{m_{1}} l_{2}^{2} - β_{m_{1} m_{2}} l_{1}^{2})] - \\ - γ_{m_{1} g} l_{2}^{2} [l_{1} (α_{m_{1}} l_{3}^{3} - β_{m_{1} m_{2}} l_{1} γ_{Mg} l_{2}^{2}) - α_{M m_{1}} (l_{1} l_{3}^{3} - γ_{Mg} l_{2}^{2} l_{2}^{2}) \\ + γ_{m_{1} g} l_{2}^{2} (α_{M m_{1}} β_{m_{1} m_{2}} l_{1} - α_{m_{1}} l_{1})] > 0 \end{matrix}

This explicit expression, while algebraically intricate, must hold true for the system to maintain overall stability and ensure physically meaningful results across all scales of deformation.

The material parameters must satisfy the positivity conditions associated with the micro-scale deformation energy (see (11)). These conditions ensure that the energy functional is bounded from below and that the resulting static boundary-value problem admits a physically meaningful solution. Violation of these requirements leads to loss of ellipticity of the governing operator, which may result in non-existence or non-uniqueness of the static solution. Furthermore, for the state-space representation to be well-defined, the coefficient matrix associated with the second derivatives in the governing equations (the submatrix $S_{sub}$ ) must be invertible. This condition guarantees that the system can be rewritten in first-order form. It is not a physical requirement but an algebraic one. When $S_{sub}$ is invertible, the system matrix $A$ is well-defined, and the characteristic equation $\det (A - λ I) = 0$ yields six (generally distinct) eigenvalues, which form the basis of the exponential modes in the analytical solution. In the case of repeated eigenvalues, generalized eigenvectors are employed to construct a complete solution basis.

5.1.3. BCs and determination of constants

To uniquely determine the six integration constants $C_{1}, \dots, C_{6}$ , six independent BCs are required. These conditions are applied at the boundaries of the domain ( $\tilde{x} = 0$ and $\tilde{x} = 1$ ) and directly involve the components of the state vector $Y (\tilde{x})$ , i.e., $\tilde{\bar{ϕ}}, {\tilde{\bar{ϕ}}}_{, \tilde{x}}, {\tilde{ψ}}_{11}, {\tilde{ψ}}_{11, \tilde{x}}, {\tilde{ψ}}_{111}, {\tilde{ψ}}_{111, \tilde{x}}$ at specified points in $\tilde{x}$ . The dimensionless BCs from equation (25) (i.e., for $\underset{̲}{h_{1}}, \underset{̲}{h_{2}}, \underset{̲}{h_{3}}$ ) can be directly translated into expressions involving the state variables and their derivatives. For example, a BC involving a second derivative (e.g., ${\tilde{\bar{ϕ}}}_{, \tilde{x} \tilde{x}} ({\tilde{x}}_{0}) = K$ ) can be converted into an expression involving the state variables by utilizing the governing equations from $Y^{'} = A Y + G$ .

6. Static behavior analysis: numerical scenarios

In this section, a systematic investigation is conducted into the impact of various $BCs$ on the distribution of kinematic variables and the resulting strain energy density in a 1D granular rod. Twelve distinct scenarios, labeled $B C_{1}$ – $B C_{12}$ , are analyzed. These are categorized into two primary groups: Dirichlet-type conditions (representing geometric clamping) and Neumann-type conditions (prescribing microstructural gradients).

To maintain a consistent baseline for comparison, all numerical simulations utilize a uniform set of material and geometric parameters: $α_{m_{1}} = 0.5$ , $α_{M m_{1}} = - 0.4$ , $l_{1} = 0$ , $l_{2} = 0.05$ , $l_{3} = 0$ , $l_{4} = {(0.01)}^{1 ∕ 4} \approx 0.316$ , $β_{m_{1} m_{2}} = 0$ , $γ_{m_{1} g} = 0$ , $γ_{Mg} = 0.02$ , and $n = 100$ . This ensures that any observed variations in the system’s response are exclusively attributable to the changing BCs applied at the rod’s ends, as detailed in Table 1.

Table 1.

Boundary conditions for scenarios $B C_{1}$ – $B C_{12}$ .

	Left end ( $\tilde{x} = 0$ )			Right end ( $\tilde{x} = 1$ )
Scenario	$\tilde{\bar{ϕ}}$	${\tilde{ψ}}_{11}$	${\tilde{ψ}}_{111}$	$\tilde{\bar{ϕ}}$	${\tilde{ψ}}_{11}$	${\tilde{ψ}}_{111}$
$B C_{1}$	0	0	0	0.01	0	0
$B C_{2}$	0	0	0	0.01	0.01	0
$B C_{3}$	0	0	0	0.01	0	0.01
$B C_{4}$	0	0	0	0.01	0.01	0.01
$B C_{5}$	0.01	0	0.01	0	0	0
$B C_{6}$	0.01	0.01	0	0	0	0
$B C_{7}$	0.01	0	0	0	0	0.01
$B C_{8}$	0.01	0	0.01	0	0.01	0
Scenario	$\tilde{\bar{ϕ}}$	${\tilde{ψ}}_{11, x}$	${\tilde{ψ}}_{111}$	$\tilde{\bar{ϕ}}$	${\tilde{ψ}}_{11, x}$	${\tilde{ψ}}_{111}$
$B C_{9}$	0	0	0	0.01	0.01	0
$B C_{10}$	0	0	0.01	0.01	0	0
$B C_{11}$	0	0	0	0.01	0	0.01
$B C_{12}$	0	0	0	0	0	0.01

6.1. Response to geometric BCs (Dirichlet-type)

In scenarios $B C_{1}$ – $B C_{8}$ , BCs are imposed directly as values of the kinematic fields at the boundaries ( $\tilde{x} = 0$ and $\tilde{x} = 1$ ) (see Figures 2 and 3).

Figure 2.

Profiles of dimensionless kinematic variables ( $\tilde{\bar{ϕ}} (\tilde{x})$ , ${\tilde{ψ}}_{11} (\tilde{x})$ , ${\tilde{ψ}}_{111} (\tilde{x})$ ) for scenarios $B C_{1}$ – $B C_{4}$ under various geometric boundary conditions: (a) $\tilde{\bar{ϕ}} (1) = 0.01$ , ${\tilde{ψ}}_{11} (1) = 0$ , ${\tilde{ψ}}_{111} (1) = 0$ ; (b) $\tilde{\bar{ϕ}} (1) = 0.01$ , ${\tilde{ψ}}_{11} (1) = 0.01$ , ${\tilde{ψ}}_{111} (1) = 0$ ; (c) $\tilde{\bar{ϕ}} (1) = 0.01$ , ${\tilde{ψ}}_{11} (1) = 0$ , ${\tilde{ψ}}_{111} (1) = 0.01$ ; and (d) $\tilde{\bar{ϕ}} (1) = 0.01$ , ${\tilde{ψ}}_{11} (1) = 0.01$ , ${\tilde{ψ}}_{111} (1) = 0.01$ .

Figure 3.

Profiles of dimensionless kinematic variables ( $\tilde{\bar{ϕ}} (\tilde{x})$ , ${\tilde{ψ}}_{11} (\tilde{x})$ , ${\tilde{ψ}}_{111} (\tilde{x})$ ) for scenarios $B C_{5}$ – $B C_{8}$ under various geometric boundary conditions: (a) ${\tilde{ψ}}_{11} (0) = 0$ , ${\tilde{ψ}}_{11} (1) = 0$ , ${\tilde{ψ}}_{111} (0) = 0.01$ , ${\tilde{ψ}}_{111} (1) = 0$ ; (b) ${\tilde{ψ}}_{11} (0) = 0.01$ , ${\tilde{ψ}}_{11} (1) = 0$ , ${\tilde{ψ}}_{111} (0) = 0$ , ${\tilde{ψ}}_{111} (1) = 0$ ; (c) ${\tilde{ψ}}_{11} (0) = 0$ , ${\tilde{ψ}}_{11} (1) = 0$ , ${\tilde{ψ}}_{111} (0) = 0$ , ${\tilde{ψ}}_{111} (1) = 0.01$ ; and (d) ${\tilde{ψ}}_{11} (0) = 0$ , ${\tilde{ψ}}_{11} (1) = 0.01$ , ${\tilde{ψ}}_{111} (0) = 0.01$ , ${\tilde{ψ}}_{111} (1) = 0$ .

6.1.1. Symmetry and boundary layer formation $(B C_{1}$ , $B C_{3}$ , $B C_{5}$ , $B C_{7})$

Scenarios $B C_{1}$ , $B C_{3}$ , $B C_{5}$ , and $B C_{7}$ exhibit a high degree of symmetry. The macro-scale field $\tilde{\bar{ϕ}} (\tilde{x})$ under conditions of pure extension or compression, the macro-displacement field remains perfectly linear, as predicted by classical elasticity. The micro-stretch ${\tilde{ψ}}_{11}$ under $(0, 0)$ BCs forms pronounced boundary layers near the ends, with the function rapidly transitioning to a plateau in the central region of the rod. The width of these layers is inherently linked to the internal characteristic length scales $l_{i}$ . Scenarios $B C_{3}$ and $B C_{5}$ act as mutual antipodes, where reversing the sign of the boundary values results in a mirrored profile of the ${\tilde{ψ}}_{11}$ plateau (transitioning from a positive dome to a negative trough), confirming the model’s linearity and stability.

6.1.2. Asymmetry and energy concentration $(B C_{2}$ , $B C_{4}$ , $B C_{6}$ , $B C_{8})$

When non-zero micro-deformation values are prescribed at a single boundary (e.g., ${\tilde{ψ}}_{11} (1) = 0.01$ in $B C_{2}$ and $B C_{4}$ ), the symmetry is broken. The extrema of the micro-kinematic fields shift significantly toward the active boundary, while the second-degree variable ${\tilde{ψ}}_{111}$ exhibits an amplitude several orders of magnitude lower than the primary fields (order of $1 0^{- 5}$ ), serving as a higher-order adjustment to the internal strain redistribution.

6.2. Transition to Neumann-type BCs $(B C_{9}$ – $B C_{12})$

The transition from fixing field values to fixing gradients ( ${\tilde{ψ}}_{11, x}$ ) fundamentally alters the internal energetics of the granular rod, as illustrated in Figure 4. In scenarios involving $(0, 0)$ Neumann BCs, the micro-stretch ${\tilde{ψ}}_{11}$ is no longer constrained to zero at the boundaries but floats to a stable plateau level throughout the entire domain. This represents a state where the internal degrees of freedom of the grains adapt uniformly to the macro-extension. It is observed that when the macro-field $\tilde{\bar{ϕ}}$ and the higher-order micro-curvature ${\tilde{ψ}}_{111}$ are co-directional (as seen in $B C_{11}$ ), the system achieves maximum microstructural homogeneity.

Figure 4.

Profiles of dimensionless kinematic variables ( $\tilde{\bar{ϕ}} (\tilde{x})$ , ${\tilde{ψ}}_{11} (\tilde{x})$ , ${\tilde{ψ}}_{111} (\tilde{x})$ ) for scenarios $B C_{9}$ – $B C_{12}$ under various geometric boundary conditions: (a) $\tilde{\bar{ϕ}} (0) = 0$ , $\tilde{\bar{ϕ}} (1) = 0.01$ , ${\tilde{ψ}}_{11, x} (0) = 0$ , ${\tilde{ψ}}_{11, x} (1) = 0.01$ , ${\tilde{ψ}}_{111} (0) = 0$ , ${\tilde{ψ}}_{111} (1) = 0$ ; (b) $\tilde{\bar{ϕ}} (0) = 0$ , $\tilde{\bar{ϕ}} (1) = 0.01$ , ${\tilde{ψ}}_{11, x} (0) = 0$ , ${\tilde{ψ}}_{11, x} (1) = 0$ , ${\tilde{ψ}}_{111} (0) = 0.01$ , ${\tilde{ψ}}_{111} (1) = 0$ ; (c) $\tilde{\bar{ϕ}} (0) = 0$ , $\tilde{\bar{ϕ}} (1) = 0.01$ , ${\tilde{ψ}}_{11, x} (0) = 0$ , ${\tilde{ψ}}_{11, x} (1) = 0$ , ${\tilde{ψ}}_{111} (0) = 0$ , ${\tilde{ψ}}_{111} (1) = 0.01$ ; and (d) $\tilde{\bar{ϕ}} (0) = 0$ , $\tilde{\bar{ϕ}} (1) = 0$ , ${\tilde{ψ}}_{11, x} (0) = 0$ , ${\tilde{ψ}}_{11, x} (1) = 0$ , ${\tilde{ψ}}_{111} (0) = 0$ , ${\tilde{ψ}}_{111} (1) = 0.01$ .

6.3. The limit case $B C_{12}$ : microstructural dominance

Scenario $B C_{12}$ represents a critical limit in the study. In this case, the macro-displacement $\tilde{\bar{ϕ}}$ is set to zero at both boundaries, causing its amplitude to collapse to the order of $1 0^{- 8}$ . In this regime, ${\tilde{ψ}}_{111}$ emerges as the dominant driver of the system’s behavior. Physically, this describes a state where external macro-scale loading is absent, and the rod’s response is governed exclusively by internal micro-stresses. This confirms the model’s ability to capture effects such as residual micro-strains or self-organization in granular media.

6.4. Analysis of dimensionless strain energy density $\tilde{W} (\tilde{x})$

The visualization of the strain energy density $\tilde{W}$ through heatmaps, shown in Figure 5, provides a final synthesis of the physical phenomena. In Dirichlet scenarios ( $B C_{2}$ , $B C_{4}$ , $B C_{6}$ , $B C_{8}$ ), energy is strictly localized at the boundaries, with bright zones on the heatmaps corresponding to regions of intense conflict between macro-extension and the rigid micro-clamping.

Figure 5.

Distribution of dimensionless strain energy density $\tilde{W} (\tilde{x})$ for scenarios: (a) $B C_{1}$ – $B C_{4}$ , (b) $B C_{5}$ – $B C_{8}$ , and (c) $B C_{9}$ – $B C_{12}$ .

In contrast, Neumann scenarios ( $B C_{9}$ – $B C_{11}$ ) exhibit uniform energy distribution throughout the volume, with bright intensity across the entire heatmap indicating that internal degrees of freedom are engaged throughout the volume. The heatmap for $B C_{12}$ appears entirely black, representing the BC12 energy threshold. As energy is a quadratic form, the collapse of field amplitudes to $1 0^{- 8}$ reduces the energy density to negligible values ( $1 0^{- 13}$ – $1 0^{- 16}$ ), falling below the macroscopic sensitivity threshold. This proves that without a macro-scale driver ( $\tilde{\bar{ϕ}}$ ), microstructural fluctuations alone cannot perform significant macroscopic work.

6.5. Parametric analysis of macro- and micro-field nonlinearity

In the primary scenarios ( $B C_{1}$ – $B C_{12}$ ) presented earlier, the non-linear effects in the macro-displacement field $\tilde{\bar{ϕ}} (\tilde{x})$ were subtle and visually indistinguishable within the scale of total displacement (as seen in Figures 2 to 4). To explicitly highlight and analyze these inherent non-linear capabilities of our micromorphic model, and to demonstrate its sensitivity to material coupling, we conducted a parametric study by varying the dimensionless coupling coefficient $α_{m_{1}}$ (which relates micro-scale and macro-scale deformations).

The results are presented in Figure 6. The left panels illustrate the profiles of ( $\tilde{\bar{ϕ}} (\tilde{x})$ , ${\tilde{ψ}}_{11} (\tilde{x})$ , and ${\tilde{ψ}}_{111} (\tilde{x})$ ) for four different values of $α_{m_{1}} = 0.5$ (0.5, 0.2, 0.9, and 1.5). The top-left panel, labeled ’Non-linear $\tilde{\bar{ϕ}} (\tilde{x})$ , demonstrates the deviation of $δ \tilde{\bar{ϕ}} (\tilde{x})$ from a linear response. While this deviation is minimal for $α_{m_{1}} = 0.5$ (blue line), other $α_{m_{1}}$ values (e.g., $α_{m_{1}} = 0.2, 0.9, 1.5$ ) lead to significantly more pronounced, complex S-shaped non-linearities in the macro-displacement profile. This clearly shows that the model captures strong non-linear interactions between the micro- and macro-scales, which become apparent under varying material characteristics. Concurrently, ${\tilde{ψ}}_{11} (\tilde{x})$ and ${\tilde{ψ}}_{111} (\tilde{x})$ also exhibit profound changes in their profiles (middle and bottom-left panels), highlighting the dynamic adaptation of the microstructure in response to varying coupling strengths.

Figure 6.

Parametric study of kinematic fields and strain energy density for varying coupling coefficient $α_{m_{1}}$ .

The corresponding distributions of strain energy density $\tilde{W} (\tilde{x})$ for these $α_{m_{1}}$ values are shown in the right panels of Figure 6. The energy localization patterns are highly sensitive to $α_{m_{1}}$ , transitioning from spatially uniform distribution to distinct boundary layer formations, or more complex non-symmetric localization, directly reflecting the interplay of macro- and micro-field non-linearities. This further underscores the model’s ability to capture intricate macro–micro interactions and their impact on the energy landscape of granular structures.

7. Conclusion

In this study, a detailed investigation into the static behavior of a 1D granular rod was conducted using a second-degree micromorphic continuum model. The inclusion of higher-order kinematic variables, micro-stretch ( ${\tilde{ψ}}_{11}$ ) and micro-curvature ( ${\tilde{ψ}}_{111}$ ), allowed for the detection of microstructural effects that are elusive in classical theories.

The analysis of 12 numerical scenarios ( $B C_{1}$ – $B C_{12}$ ) led to the following conclusions:

BC Sensitivity: The distribution of internal micro-deformations is critically dependent on the BC type. Dirichlet conditions facilitate the formation of sharp boundary layers, while Neumann conditions promote the field homogenization.

Scale Hierarchy: The model successfully separates scales of influence. The macro-displacement $\tilde{\bar{ϕ}}$ defines the overall energetic baseline ( $1 0^{- 3}$ ), while higher-order micro-variables govern subtle force redistributions at scales of $1 0^{- 5}$ to $1 0^{- 8}$ .

Failure Prediction via Energy: Analysis of $\tilde{W} (\tilde{x})$ visualized stress concentration zones. Rigid geometric clamping localizes energy at interfaces, serving as a potential predictor for material failure.

Mathematical Robustness: Scenario $B C_{12}$ demonstrated that the model remains stable even in extreme limit regimes, accurately transitioning between macroscopic and microstructural dominance.

The scientific novelty of this work lies in the systematization of transition processes from macro- to micro-scales. The developed framework and energy maps provide a robust tool for the design of advanced composites and metamaterials where properties are tuned through purposeful microstructural modification.

The scientific novelty of this work lies in the systematization of transition processes from macro- to micro-scales. Crucially, the model also demonstrates the capacity to capture significant non-linearity in the macro-displacement field, directly influenced by microstructural coupling parameters, thus offering a richer description than classical continuum theories. The developed framework and energy maps provide a robust tool for the design of advanced composites and metamaterials where properties are tuned through purposeful microstructural modification.

Extension to dynamic loading and multi-dimensional settings will be addressed in future work.

Footnotes

ORCID iD

Armine Ulukhanyan

Funding

The authors received no financial support for the research, authorship, and/or publication of this article.

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

References

Mindlin

. Micro-structure in linear elasticity. Arch Ration Mech Anal 1964; 16: 51–78.

Germain

. The method of virtual power in continuum mechanics. Part 2: microstructure. SIAM J Appl Math 1973; 25(3): 556–575.

Eringen

. Microcontinuum field theories I: foundations and solids. Springer, 1999. https://doi.org/10.1007/978-1-4612-0555-5.

Aifantis

. Update on a class of gradient theories. Mech Mater 2003; 35(3–6): 259–280.

Cosserat

. Théorie des corps déformables. A. Hermann et Fils, 1909.

Jirásek

. Nonlocal theories in continuum mechanics. Acta Polytechnica 2004; 44: 16–34.

Nikabadze

Ulukhanyan

. Modeling of multilayer thin bodies. Contin Mech Thermodyn 2020; 32(3): 817–842.

Nikabadze

Ulukhanyan

. On the interlayer contact conditions in multilayer thin body theory and some issues of splitting initial-boundary value problems. Lobachevskii J Math 2022; 43(7): 1945–1961.

Ulukhanyan

. Dispersion equations and wave velocities in micropolar theory of elasticity. Lobachevskii J Math 2024; 45(8): 3773–3782.

10.

Chang

Misra

. Theoretical and experimental study of regular packings of granules. J Eng Mech (ASCE) 1989; 115: 704–720.

11.

Chang

Misra

. Packing structure and mechanical properties of granulates. J Eng Mech (ASCE) 1990; 116: 1077–1093.

12.

Misra

Poorsolhjouy

. Granular micromechanics based micromorphic model predicts frequency band gaps. Contin Mech Thermodyn 2016; 28: 215–234.

13.

Tordesillas

Walsh

SDC

. Incorporating rolling resistance and contact anisotropy in micromechanical models of granular media. Powder Technol 2002; 124: 106–111.

14.

Kruyt

Rothenburg

. A micromechanical study of dilatancy of granular materials. J Mech Phys Solids 2016; 95: 411–427.

15.

Nejadsadeghi

Misra

. Extended granular micromechanics approach: a micromorphic theory of degree n. Math Mech Solids 2020; 25(3): 407–429.

16.

Nejadsadeghi

Misra

. On the statics and dynamics of granular-microstructured rods with higher order effects. Math Mech Solids 2021; 26(12): 1815–1837.

17.

Alibert

Seppecher

Dell’Isola

. Truss modular beams with deformation energy depending on higher displacement gradients. Math Mech Solids 2003; 8(1): 51–73.

18.

dell’Isola

Seppecher

Alibert

, et al. Pantographic metamaterials: an example of mathematically driven design and of its technological challenges. Contin Mech Thermodyn 2019; 31(4): 851–884.

19.

dell’Isola

Seppecher

Spagnuolo

, et al. Advances in pantographic structures: design, manufacturing, models, experiments and image analyses. Contin Mech Thermodyn 2019; 31: 1231–1282.

20.

Forest

Sievert

. Elastoviscoplastic constitutive frameworks for generalized continua. Int J Solids Struct 2002; 39(19): 4835–4856.

21.

Forest

. Micromorphic approach for gradient elasticity, viscoplasticity, and damage. J Eng Mech 2009; 135(3): 117–131.

22.

Neff

Ghiba

Madeo

, et al. A unifying perspective: the relaxed linear micromorphic continuum. Contin Mech Thermodyn 2014; 26: 639–681.

23.

Ulukhanyan

Placidi

Fedele

, et al. Embedding 1D Euler beam in 2D second-gradient continua. Math Mech Solids. Epub ahead of print 1 October 2025. https://doi.org/10.1177/10812865251364514

24.

Maugin

Metrikine

. Mechanics of generalized continua: one hundred years after the Cosserats, Advances in mechanics and mathematics, vol. 21. Springer, 2010. https://doi.org/10.1007/978-1-4419-5695-8.

25.

Misra

Nejadsadeghi

. Longitudinal and transverse elastic waves in 1D granular materials modeled as micromorphic continua. Wave Motion 2019; 90: 175–195.

26.

Tan

Poh

. Homogenized gradient elasticity model for plane wave propagation in bilaminate composites. J Eng Mech 2018; 144(7): 04018075.

27.

Tan

Poh

. Enriched homogenized model for viscoelastic plane wave propagation in periodic layered composites. Adv Model Simul Eng Sci 2020; 7: 4.

28.

Reddy

. Microstructure-dependent couple stress theories of functionally graded beams. J Mech Phys Solids 2011; 59: 2382–2399.

29.

Papargyri-Beskou

Tsepoura

Polyzos

, et al. Bending and stability analysis of gradient elastic beams. Int J Solids Struct 2003; 40: 385–400.

30.

Wellmann

Lillie

Wriggers

. Homogenization of granular material modeled by a three-dimensional discrete element method. Comput Geotech 2008; 35: 394–405.

31.

Charlotte

Truskinovsky

. Towards multi-scale continuum elasticity theory. Contin Mech Thermodyn 2008; 20: 133–161.

32.

Charlotte

Truskinovsky

. Lattice dynamics from a continuum viewpoint. J Mech Phys Solids 2012; 60(8): 1508–1544.

33.

Placidi

Di Girolamo

Fedele

. Variational study of a Maxwell–Rayleigh-type finite length model for the preliminary design of a tensegrity chain with a tunable band gap. Mech Res Commun 2024; 136: 104255.

34.

De Fazio

Placidi

Fabbrocino

, et al. Analysis of ultrasonic wave dispersion in presence of attenuation and second-gradient contributions. CivilEng 2025; 6(3): 37.

35.

Kunin

. Elastic media with microstructure I: one-dimensional models, Springer series in solid-state sciences, vol. 26. Springer, 1982. https://doi.org/10.1007/978-3-642-81748-9.

36.

Askes

Metrikine

. Higher-order continua derived from discrete media: continualisation aspects and boundary conditions. Int J Solids Struct 2005; 42(1): 187–202.

A second-degree micromorphic continuum model for 1D granular structures

Abstract

Keywords

1. Introduction

2. Continuum framework: a micromorphic model of degree 2

2.1. Kinematic variables

2.2. Deformation measures and intergranular kinematics

2.3. Constitutive equations

3. Governing equations and BCs

3.1. Variational formulation

3.2. Equilibrium equations in strong form

4. Solution for the static problem

4.1. Operator formulation

4.2. Dimensionless form of the governing equations

5. Solution for the static problem

5.1. Analytical solution using state-space method

5.1.1. State-space formulation

5.1.2. General solution and conditions for existence

5.1.3. BCs and determination of constants

6. Static behavior analysis: numerical scenarios

6.1. Response to geometric BCs (Dirichlet-type)

6.1.1. Symmetry and boundary layer formation ( B C 1 , B C 3 , B C 5 , B C 7 )

6.1.2. Asymmetry and energy concentration ( B C 2 , B C 4 , B C 6 , B C 8 )

6.2. Transition to Neumann-type BCs ( B C 9 – B C 12 )

6.3. The limit case B C 12 : microstructural dominance

6.4. Analysis of dimensionless strain energy density W ~ ( x ~ )

6.5. Parametric analysis of macro- and micro-field nonlinearity

7. Conclusion

Footnotes

ORCID iD

Funding

Declaration of conflicting interests

References

6.1.1. Symmetry and boundary layer formation $(B C_{1}$ , $B C_{3}$ , $B C_{5}$ , $B C_{7})$

6.1.2. Asymmetry and energy concentration $(B C_{2}$ , $B C_{4}$ , $B C_{6}$ , $B C_{8})$

6.2. Transition to Neumann-type BCs $(B C_{9}$ – $B C_{12})$

6.3. The limit case $B C_{12}$ : microstructural dominance

6.4. Analysis of dimensionless strain energy density $\tilde{W} (\tilde{x})$