Statistical mechanics and phantom network modelling of disordered biopolymer gels

Abstract

Disordered biopolymer gels, such as alginate gels, contain disordered amorphous regions as well as ordered regions known as junction zones, the latter playing a paramount role in maintaining the integrity of the gel network. The considerable dimensions of these junction zones and their physical interactions are responsible for distinctive phenomena, such as self-healing, observed in biopolymer gels. Compared to rubber-like materials, where polymer chains are interconnected by chemical cross-links, statistical mechanics modeling of disordered biopolymer gels poses new challenges due to the need to capture the contributions from both the disordered and ordered regions. In the literature, the disordered domains are usually modeled by a collection of random coils, while each junction zone is represented by a rigid rod. Attempt has been made to introduce the two-node coil-rod structure, consisting of a random coil connected to a rigid rod, into classical polymer network models (e.g., the Arruda–Boyce eight-chain model) to predict the mechanical response of the biopolymer gels. Although this approach provides valuable insights and reasonable predictions that agree with experimental data, the interactions between the junction zones and the amorphous regions are significantly simplified. This study aims to extend the two-node coil-rod structure to a four-node coil-rod structure, in which two polymer chains share a junction zone, more explicitly representing chain association induced by physical forces in biopolymer gels. By incorporating statistical mechanics and the phantom network model, the entropy of the system is derived, from which the stress–stretch relationships for the network are obtained. Comparisons are made with the model based on the two-node coil-rod structure, and the conditions for establishing equivalency between the two methods are examined. This work establishes a foundation for developing more advanced network models that incorporates the simultaneous interaction of multiple junction zones, a phenomenon often present in biopolymer gels.

Nomenclature

${\left[\mathbf{0}\right]}_{3\times 3}$ The $3\times 3$ zero matrix ${\left[\mathbf{I}\right]}_{3\times 3}$ The $3\times 3$ identity matrix a Length of the rod A Matrix defined in equations (27) and (95) $A_0$ Normalization factor in the probability distribution (45) b Kuhn length c An arbitrary constant D Matrix defined in equation (64) ${\widehat{\mathbf{e}}}_x$ , ${\widehat{\mathbf{e}}}_y$ , ${\widehat{\mathbf{e}}}_z$ Unit vectors of the Cartesian coordinate system $f_1$ , $f_2$ Labels of nodes at the rod ends in the general network with one rod; see Figure 10(a) G The modulus equals to $\varrho k_{B}T$ i Imaginary unit number, $\sqrt{-1}$ $I_0$ Expression defined in equation (76) ${\left(k_{\mathrm{ij}}\right)}_x$ , ${\left(k_{\mathrm{ij}}\right)}_y$ , ${\left(k_{\mathrm{ij}}\right)}_z$ Components of ${\mathbf{k}}_{\mathrm{ij}}$ in Cartesian coordinates K Matrix defined in equations (21a) and (92a) k Position vector in the Fourier space corresponding to r ${\mathbf{K}}_a$ , ${\mathbf{K}}_b$ Block matrices defined in (69) with dimensions $(3p-3)\times 1$ and $3\times 1$ , respectively $k_B$ The Boltzmann constant ${\mathbf{k}}_{\mathrm{ij}}$ Position vector in the Fourier space corresponding to ${\mathbf{R}}_{\mathrm{ij}}$ L Matrix defined in equations (73) and (83) $\stackrel{\sim }{\mathbf{L}}$ Matrix defined in equation (102) $l_1$ , $l_2$ Dimensions of the two-dimensional rectangle cell along horizontal and vertical directions in the stress-free state; see Figure 4(a) m Extra number of segments 46 and 63 compared to 15 and 52 M Orthogonal matrix for diagonalizing Γ n Number of Kuhn segments in the coil ${\overline{n}}_i$ The number of Kuhn segments of the ith chain in the equivalent affine network; see Figure 9(b) $n_f$ Count of free nodes $n_r$ Number of rods $n_{\mathrm{ij}}$ Number of Kuhn segments of a coil between nodes i and j p Number of independent end-to-end vectors of the $p+1$ -node structure; see Figures 9(a) and 10(a) $P\left({\mathbf{x}}_1,\cdots ,{\mathbf{x}}_{t-1}\right)$ Probability that the t nodes of the structure is fixed with vectors ${\mathbf{x}}_1,\cdots ,{\mathbf{x}}_{t-1}$ , fixed temperature and number of Kuhn segments n $P_1$ , $P_2$ , $P_3$ Principal components of first Piola–Kirchhoff stress tensor $\mathbf{R}$ Matrix defined in equation (67) ${\overline{\mathbf{R}}}_i$ The end-to-end vector of the ith chain in the equivalent affine network; see Figure 9(b) R Matrix defined in equations (21b) and (92b) r End-to-end vector of the two-node chain ${\mathbf{R}}_a$ , ${\mathbf{R}}_b$ Block matrices defined in (70) with dimensions $(3p-3)\times 1$ and $3\times 1$ , respectively ${\mathbf{R}}_{\mathrm{ij}}$ Chain vector from node i to node j ${\mathbf{R}}_i$ Position vector of node i ${\left(R_{\mathrm{ij}}\right)}_x$ , ${\left(R_{\mathrm{ij}}\right)}_y$ , ${\left(R_{\mathrm{ij}}\right)}_z$ Components of ${\mathbf{R}}_{\mathrm{ij}}$ in Cartesian coordinates $\stackrel{\sim }{\mathbf{R}}$ Matrix defined in equations (31) and (98) ${\stackrel{\sim }{\mathbf{R}}}_a$ , ${\stackrel{\sim }{\mathbf{R}}}_b$ Block matrices defined in (99) with dimensions $(3p-3)\times 1$ and $3\times 1$ , respectively S The entropy of the deformed state per unit reference volume $s_1,s_2,\cdots ,s_{p+1}$ Labels of fixed nodes in the general network with $p+1$ fixed nodes; see Figures 9(a) and 10(a) T Absolute temperature t Number of nodes on one side of the general network with one rod; see Figure 10(a) u Variable vector defined in equation (9b) v Variable vector defined in equation (9a) w Variable vector defined in equation (29) $\overline{W}\left({\mathbf{k}}_{\mathrm{ij}}\right)$ Fourier transform of the probability distribution $W\left({\mathbf{R}}_{\mathrm{ij}}\right)$ $W^{\mathrm{R}}\left(\left|\mathbf{r}\right|\right)$ The probability that the end-to-end distance of the free rod is $\left|\mathbf{r}\right|$ , at fixed temperature $W_n^{\mathrm{CR}}(|\mathbf{r}|)$ Probability that the end-to-end distance of the free two-node coil-rod structure is $\left|\mathbf{r}\right|$ , with fixed temperature and number of Kuhn segments n in the coil $W_n^{\mathrm{C}}(|\mathbf{r}|)$ Probability that the end-to-end distance of the free coil structure is $\left|\mathbf{r}\right|$ , with fixed temperature and number of Kuhn segments n $w_x$ , $w_y$ , $w_z$ Components of w in Cartesian coordinates based on equation (29) defined in (26b) X Vector defined in equation (60) $X_i$ The ith element of vector X Y Matrix defined in equations (25), (26a), (93), and (94) β A dimensionless parameter defined in equation (46) $\delta \left(\mathbf{x}\right)$ Three-dimensional Dirac delta function $\delta \left(x\right)$ One-dimensional Dirac delta function ${\eta }_i$ The ith eigenvalue of matrix Γ ${\gamma }_{\mathrm{ij}}$ Element ij in matrix Γ defined in equation (71) $\mathbf{\Gamma }$ Matrix defined in equations (23), (71), and (90) ${\mathbf{\Gamma }}_{\mathrm{aa}}$ , ${\mathbf{\Gamma }}_{\mathrm{ab}}$ , ${\mathbf{\Gamma }}_{\mathrm{bb}}$ Block matrices defined in equation (77) $\stackrel{\sim }{\mathbf{\Gamma }}$ Matrix defined in equations (30) and (97) ${\stackrel{\sim }{\mathbf{\Gamma }}}_{\mathrm{aa}}$ , ${\stackrel{\sim }{\mathbf{\Gamma }}}_{\mathrm{ab}}$ , ${\stackrel{\sim }{\mathbf{\Gamma }}}_{\mathrm{bb}}$ Block matrices defined in equation (100) $\mathbf{\Lambda }$ Matrix defined in equations (24) and (91) ${\lambda }_1$ , ${\lambda }_2$ , ${\lambda }_3$ Principal stretches of the bulk material ν Total number of Gaussian chains in the general network with fixed $p+1$ nodes ϕ and θ Spherical coordinate variables Ψ The Helmholtz free energy of the deformed state per unit reference volume ${\mathrm{\Psi }}_0$ The Helmholtz free energy of the undeformed state per unit reference volume ϱ Chain density per unit reference volume ${\tilde{\zeta }}_{\mathrm{ij}}$ Element ij in matrix ${\stackrel{\sim }{\mathbf{\Gamma }}}^{-1}$ ${\zeta }_{\mathrm{ij}}$ Element ij in matrix ${\mathbf{\Gamma }}^{-1}$ defined in equation (75)

1. Introduction

Polymer networks are established through the copious interpenetration of chains. The complex interactions between them as well as structural features such as topological loops [1] make studying their statistical thermodynamics formidable. In the method proposed by James and Guth [2,3], known as the affine network model, the chain vectors are situated randomly within the network, and their extension is assumed to follow the macroscopic deformation controlled by external loading. In this framework, since the interaction of chains are only considered via cross-links, the partition function of each chain provides sufficient information to obtain the partition function of the entire network. Building upon this pioneering work, subsequent literature endeavors to develop more advanced models that describe observed phenomena in polymer networks while being numerically affordable. To name but a few, these includes the three-chain model [4], the full-network model [5,6], the eight-chain model [7,8], and the micro-sphere model [9]. In the method proposed by Flory [10], known as the phantom network model, the constraints of macroscopic deformation are relaxed, allowing the chains to maintain their topology while freely passing through each other. In this framework, the network is treated as a collection of chains where only a subset of cross-links is restricted by macroscopic deformation, while the remaining cross-links fluctuate around their average positions.

In disordered biopolymer gels, amorphous regions known as coils are interconnected through ordered regions called junction zones [11,12]. The process of junction zone formation varies depending on the specific biopolymers involved. For instance, in agar gels, the chains intertwine to create a double helix structure, while in alginate gels, the junction zone is formed through the mediation of external cation agents, resulting in egg-box structures. Disordered biopolymer gels can be decomposed into two fundamental building blocks: coils and junction zones. Coils are often modeled as freely jointed chains, while junction zones are represented by rigid rods. Within the framework of affine network, the end-to-end extension of each component (coil or rod) is assumed to be proportional to the corresponding macroscopic displacement. Since the rigid rod has infinite stiffness and therefore cannot undergo extension, the macroscopic behavior encompassing the rod reflects an unrealistically high stiffness under affinity assumptions.

To mitigate this issue, Higgs and Ball [13], and later Moosavian and Tang [14], considered the network as a collection of “coil-rod” structures, in each structure the freely jointed chain is connected to a rigid rod (see Figure 1(a)). The end-to-end extension of this combined structure (rather than treating the coil and rod separately) still follows the macroscopic deformation akin to an affine network. However, the node connecting the coil and rod can freely move resembling the phantom network. Since the macroscopic deformation can be related to the distance between two end nodes of the structure, it is hereafter referred to as the “two-node coil-rod structure”. This structure can be directly integrated into well-known affine network models in the literature. In addition, it allows us to capture the phenomenon of zipping/unzipping, where the junction zones expand or shrink due to applied loading. Recently, it has been demonstrated that implementing the two-node coil-rod structure with zipping/unzipping into the eight-chain model effectively captures both the dissipation caused by unzipping and the phenomenon of permanent set observed during cyclic loading and unloading [15]. Despite these advances, the two-node coil-rod structure in the aforementioned model is treated as a single entity. The decomposition of the complicated network, involving different regions of coils and rods, into a collection of two-node coil-rod structures is a simplification worth further examination. This prompts us to advance the model into a more explicit network where two or more coils are connected by sharing a single junction zone in the middle (see Figure 1(b)). Such a representation takes steps toward a more accurate depiction of chain interactions in biopolymer gels. By providing the new model, the ranges of applicability and limitations of the two-node coil-rod structure can also be explored, addressing whether it can adequately model the complex interactions between coils and rods.

Figure 1.

(a) The two-node coil-rod structure with nodes 1 and 2 attached to the macroscopic network and fluctuating node 3. The zigzag drawing is a schematic representation of the egg-box structure, with hollow circles indicating the gelling agents. Segment between 1 and 3 is modelled by a rigid rod and the coil between 2 and 3 is modelled by freely jointed chain. (b) The topology of the four-node coil-rod structure with nodes 1 to 4 attached to the macroscopic network and fluctuating nodes 5 and 6. The zigzag structure between 5 and 6 represents the rigid rod and the remaining four coils follow the freely jointed chain model.

In Section 2, the model of the two-node coil-rod structure is briefly reviewed, and a new four-node coil-rod structure is formulated within the framework of statistical mechanics in Section 2.1. Two formalisms for deriving the probability distribution are explained in Sections 2.2 and 2.3. Subsequently, the results are simplified by invoking Gaussian statistics for the coils in Section 2.4. It is shown that under this assumption, the four-node coil-rod structure can be decomposed into a collection of one two-node coil-rod structure and two pure coils. Section 3 integrates this finding into network models, and the results are compared with the eight-chain network model containing two-node coil-rod structures. Finally, Section 4 discusses the limitations and advantages of the proposed model for further applications.

2. Statistical mechanics of a four-node coil-rod structure

Prior to developing the four-node coil-rod structure, the formulation of a two-node coil-rod structure is reviewed. For a freely jointed chain with Kuhn length b and number of Kuhn segments n, if the end-to-end distance $\left|\mathbf{r}\right|$ is sufficiently small and n is sufficiently large, one can find the probability distribution $W_n^{\mathrm{C}}(\left|\mathbf{r}\right|)$ as [16]

$$\begin{array}{ccc}\multicolumn{1}{c}{W_n^{\mathrm{C}}(\left|\mathbf{r}\right|)={\left(\frac{3}{2\mathrm{\pi n}{b}^2}\right)}^{3∕2}\exp \left[-\frac{3{\left|\mathbf{r}\right|}^2}{2n{b}^2}\right],}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(1)

where r represents the end-to-end vector of the coil and |⋯| denotes the norm of the vector. This probability is commonly referred to as Gaussian distribution. In addition, for a single rigid rod with length a, the probability distribution can be described through [14]

$$\begin{array}{ccc}\multicolumn{1}{c}{W^{\mathrm{R}}(\left|\mathbf{r}\right|)=\frac{1}{4\pi a^2}\delta (\left|\mathbf{r}\right|-a),}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(2)

where $\delta (\cdots )$ denotes the one-dimensional Dirac delta function applied to a scalar argument. In the work by Moosavian and Tang [14], it was shown that if the Gaussian coil is followed by a rigid rod (see Figure 1(a)), the probability distribution of the resultant structure is modified to

$$\begin{array}{ccc}\multicolumn{1}{c}{W_n^{\mathrm{CR}}(\left|\mathbf{r}\right|)=\sqrt{\frac{3}{{\left(2\pi \right)}^{3}n{a}^2{b}^2}}\frac{1}{\left|\mathbf{r}\right|}\exp \left[-\frac{3}{2n{b}^2}({\left|\mathbf{r}\right|}^2+a^2)\right]\sinh \left(\frac{3a\left|\mathbf{r}\right|}{n{b}^2}\right).}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(3)

The superscripts CR, R, and C, respectively, signifies the two-node coil-rod, rod, and coil structure. For $a\to 0$ , the rod vanishes, and equation (3) converts back to that of Gaussian coil (1). The result given by (3) can be applied to the affine network models by simply replacing the coils in the traditional models such as eight-chain network [7] with the two-node coil-rod structures. As elaborated in the forthcoming sections, information on the probability distribution of structures and their placement within the network model is sufficient to derive the constitutive relations of the material. Specifically, the entropy of each structure is obtained by taking the natural logarithm of the probability distribution. In the context of entropic elasticity, the change in Helmholtz free energy under a constant temperature is solely determined by the change in entropy. Once a connection is established between the deformation of the microstructures and that of the macroscopic network, the Helmholtz free energy of the network can be calculated in terms of the deformation gradient tensor, which in turn enables the derivation of stress–stretch relations.

2.1. Four-node coil-rod structure

The simplest scenario illustrating the explicit interaction of coils and rods involves two coils sharing a single junction zone between them. Let us consider two coils, labeled 13 and 24 in Figure 1(b), which are laterally associated, thereby creating the junction zone 56 and regenerating four new coils 15, 25, 36, and 46. In three-dimensional space, the four nodes can form a cell as demonstrated by a tetrahedral in Figure 1(b). Hereafter, such a structure is referred to as the “four-node coil-rod structure”.

By defining ${\mathbf{R}}_{\mathrm{ij}}={\mathbf{R}}_j-{\mathbf{R}}_i$ , with ${\mathbf{R}}_i$ being the position vector of node i and setting node 1 fixed, the probability of locating node i within $\mathrm{d}{\mathbf{R}}_{1i}$ for $i=2,\dots ,6$ is expressed as $P\left({\mathbf{R}}_{12},{\mathbf{R}}_{13},{\mathbf{R}}_{14},{\mathbf{R}}_{15},{\mathbf{R}}_{16}\right){\mathrm{\prod }}_{i=2}^6\mathrm{d}{\mathbf{R}}_{1i}$ . To derive the probability distribution function of such system, all coils are represented by freely jointed chain models, while the junction zone is envisioned as a rigid rod of length a. In the affine network model, all nodes 1 to 6 are tethered to the macroscopic element, ensuring that the extension of vectors ${\mathbf{R}}_{\mathrm{ij}}$ follows the applied macroscopic deformation. Although deriving the probability distribution based on this viewpoint is relatively straightforward, the presence of a rigid rod imparts significant stiffness to the macroscopic element, rendering the overall response unrealistic. To circumvent this issue, the phantom network is adopted, enabling nodes 1 to 4 to undergo manipulation by macroscopic deformation, while nodes 5 and 6 are capable of fluctuating around their mean positions. Accordingly, the probability of finding the nodes i inside the volume element $\mathrm{d}{\mathbf{R}}_{1i}$ ( $i=2,3,4$ ) is denoted by $P\left({\mathbf{R}}_{12},{\mathbf{R}}_{13},{\mathbf{R}}_{14}\right)\mathrm{d}{\mathbf{R}}_{12}\mathrm{d}{\mathbf{R}}_{13}\mathrm{d}{\mathbf{R}}_{14}$ . The repetition of such structures results in a network forming a body in the undeformed state. Therefore, the entropy per unit volume of the undeformed state (i.e., per unit reference volume) is presented by:

$$\begin{array}{ccc}\multicolumn{1}{c}{S=\varrho k_B\ln \left[\mathrm{d}{\mathbf{R}}_{12}\mathrm{d}{\mathbf{R}}_{13}\mathrm{d}{\mathbf{R}}_{14}\right]+\varrho k_B\ln P({\mathbf{R}}_{12},{\mathbf{R}}_{13},{\mathbf{R}}_{14})+c,}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(4)

where ϱ is the chain density per unit reference volume, $k_B$ is the Boltzmann constant, and c is an arbitrary constant. The entropy change from the undeformed to the deformed state is readily obtained from the difference between their respective values. Once $\mathrm{\Delta }S$ is at hand, the Helmholtz free energy per unit reference volume, Ψ, is obtained as

$$\begin{array}{ccc}\multicolumn{1}{c}{\mathrm{\Psi }=-T\mathrm{\Delta }S,}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(5)

where T is the absolute temperature. Depending on the geometry of four-node coil-rod structure, the dimensions $\mathrm{d}{\mathbf{R}}_{1i}$ ( $i=2,3,4$ ) should be related to the principal stretches of the bulk material, ${\lambda }_1$ , ${\lambda }_2$ , and ${\lambda }_3$ . For compressible material, the principal components of first Piola–Kirchhoff stress $P_1$ , $P_2$ , and $P_3$ can be obtained as [17]

$$\begin{array}{ccc}\multicolumn{1}{c}{P_i=\frac{\partial \mathrm{\Psi }}{\partial {\lambda }_i},i=1,2,3.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(6)

It can be seen that the constitutive equation directly depends on the probability distribution $P({\mathbf{R}}_{12},{\mathbf{R}}_{13},{\mathbf{R}}_{14})$ . In general, there are two approaches to find this probability: the real space method and the Fourier space method, which will be elaborated upon in the next sections.

2.2. Real space representation of the probability

Suppose the probability distribution for a coil with $n_{\mathrm{pq}}=n_{\mathrm{qp}}$ Kuhn segments trapped between nodes p and q is denoted as $W_{n_{\mathrm{pq}}}^{\mathrm{C}}(|\mathbf{r}|)$ . By setting node 1 fixed, the probability of locating node i within $\mathrm{d}{\mathbf{R}}_{1i}$ for $i=2,\dots ,6$ is expressed as:

$$\begin{array}{cccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{P\left({\mathbf{R}}_{12},{\mathbf{R}}_{13},{\mathbf{R}}_{14},{\mathbf{R}}_{15},{\mathbf{R}}_{16}\right)\prod _{i=2}^6\mathrm{d}{\mathbf{R}}_{1i}}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{=W_{n_{15}}^{\mathrm{C}}\left(\left|{\mathbf{R}}_{15}\right|\right)W_{n_{52}}^{\mathrm{C}}\left(\left|{\mathbf{R}}_{52}\right|\right)W_{n_{63}}^{\mathrm{C}}\left(\left|{\mathbf{R}}_{63}\right|\right)W_{n_{46}}^{\mathrm{C}}\left(\left|{\mathbf{R}}_{46}\right|\right)W^{\mathrm{R}}\left(\left|{\mathbf{R}}_{56}\right|\right)\prod _{i=2}^6\mathrm{d}{\mathbf{R}}_{1i}.}& \multicolumn{1}{c}{}\end{array}$$

(7)

Now, the probability of finding nodes 2 to 4 within their own volume elements $\mathrm{d}{\mathbf{R}}_i$ ( $i=2,3,4$ ) is given by integration of the aforementioned probability over all possible locations of nodes 5 and 6, i.e.,

$$\begin{array}{cccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{P\left({\mathbf{R}}_{12},{\mathbf{R}}_{13},{\mathbf{R}}_{14}\right)=\frac{1}{4\pi a^2}\times }& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{\int \int W_{n_{15}}^{\mathrm{C}}\left(\left|{\mathbf{R}}_{15}\right|\right)W_{n_{52}}^{\mathrm{C}}\left(\left|{\mathbf{R}}_{52}\right|\right)W_{n_{63}}^{\mathrm{C}}\left(\left|{\mathbf{R}}_{63}\right|\right)W_{n_{46}}^{\mathrm{C}}\left(\left|{\mathbf{R}}_{46}\right|\right)\delta (\left|{\mathbf{R}}_{56}\right|-a)\mathrm{d}{\mathbf{R}}_{15}\mathrm{d}{\mathbf{R}}_{16}.}& \multicolumn{1}{c}{}\end{array}$$

(8)

By employing the following change of variables

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{\mathbf{v}=\frac{1}{2}\left({\mathbf{R}}_{15}-{\mathbf{R}}_{16}\right),}& \multicolumn{1}{c}{}\end{array}$$

(9a)

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{\mathbf{u}=\frac{1}{2}\left({\mathbf{R}}_{15}+{\mathbf{R}}_{16}\right),}& \multicolumn{1}{c}{}\end{array}$$

(9b)

equation (8) is rephrased as

$$\begin{array}{cccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{P\left({\mathbf{R}}_{12},{\mathbf{R}}_{13},{\mathbf{R}}_{14}\right)=\frac{2}{\pi a^2}\int \int W_{n_{15}}^{\mathrm{C}}\left(\left|\mathbf{v}+\mathbf{u}\right|\right)W_{n_{52}}^{\mathrm{C}}\left(\left|{\mathbf{R}}_{12}-\mathbf{v}-\mathbf{u}\right|\right)}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{\times W_{n_{63}}^{\mathrm{C}}\left(\left|{\mathbf{R}}_{13}+\mathbf{v}-\mathbf{u}\right|\right)W_{n_{46}}^{\mathrm{C}}\left(\left|{\mathbf{R}}_{14}+\mathbf{v}-\mathbf{u}\right|\right)\delta (2|\mathbf{v}|-a)\mathrm{d}\mathbf{u}\mathrm{d}\mathbf{v}.}& \multicolumn{1}{c}{}\end{array}$$

(10)

Utilizing $\delta \left(2\left|\mathbf{v}\right|-a\right)=\frac{1}{2}\delta \left(\left|\mathbf{v}\right|-\frac{a}{2}\right)$ , and applying the sifting property of the Dirac delta function, one layer of integration is eliminated, resulting in

$$\begin{array}{cccc}\multicolumn{1}{c}{P\left({\mathbf{R}}_{12},{\mathbf{R}}_{13},{\mathbf{R}}_{14}\right)=\frac{1}{4\pi }}& \multicolumn{1}{c}{\int [{\int }_{\varphi =0}^{\pi }{\int }_{\theta =0}^{2\pi }{W}_{n_{15}}^{\mathrm{C}}\left(\left|\mathbf{v}+\mathbf{u}\right|\right)W_{n_{52}}^{\mathrm{C}}\left(\left|{\mathbf{R}}_{12}-\mathbf{v}-\mathbf{u}\right|\right)}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{\times }& \multicolumn{1}{c}{W_{n_{63}}^{\mathrm{C}}\left(\left|{\mathbf{R}}_{13}+\mathbf{v}-\mathbf{u}\right|\right)W_{n_{46}}^{\mathrm{C}}\left(\left|{\mathbf{R}}_{14}+\mathbf{v}-\mathbf{u}\right|\right)\sin \varphi \mathrm{d}\varphi \mathrm{d}\theta ]\mathrm{d}\mathbf{u},}& \multicolumn{1}{c}{}\end{array}$$

(11)

where v is given in spherical coordinate as

$$\begin{array}{ccc}\multicolumn{1}{c}{\mathbf{v}=\frac{a}{2}\left(\sin \varphi \cos \theta ,\sin \varphi \sin \theta ,\cos \varphi \right).}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(12)

The schematics of the above variables are illustrated in Figure 2. It should be recalled that in freely jointed chain models, coils cannot extend beyond their fully extended state, where the associated probability is zero. Analogously, the maximum value of $\left|\mathbf{u}\right|$ is the fully extended length of the coil trapped between nodes 1 and 5 plus half of the rod length, $n_{15}b+a∕2$ . Any probability larger than this critical value is zero. However, constraints from other ends can impose stricter conditions on the probability, leading to smaller critical values of $\left|\mathbf{u}\right|$ , which may not be immediately apparent. This critical value helps us in performing numerical integration without encountering infinity in the bounds of integrals.

Figure 2.

The five layers of integration in equation (11) include three-dimensional integration over all possible vectors u and the surface integral with variable ϕ and θ over the sphere which has its center located at the end of u and radius $a∕2$ .

In the context of more complex system of coils and rods, the number of integration layers is given by $3n_f-n_r$ , where $n_f$ represents the count of free nodes (in this case, two nodes 5 and 6), and $n_r$ denotes the number of rigid rods. In the scenario shown in Figure 2, the number of integration layers is $3\times 2-1=5$ .

2.3. Fourier space representation of the probability

Fourier space representation sometimes offers a simpler approach to calculating the probability distribution. Since nodes 1 to 4 are attached to the macroscopic element, the following constraints between the end-to-end vectors hold:

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{{\mathbf{R}}_{12}={\mathbf{R}}_{15}+{\mathbf{R}}_{52},}& \multicolumn{1}{c}{}\end{array}$$

(13a)

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{{\mathbf{R}}_{13}={\mathbf{R}}_{15}+{\mathbf{R}}_{56}+{\mathbf{R}}_{63},}& \multicolumn{1}{c}{}\end{array}$$

(13b)

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{{\mathbf{R}}_{14}={\mathbf{R}}_{15}+{\mathbf{R}}_{56}+{\mathbf{R}}_{64}.}& \multicolumn{1}{c}{}\end{array}$$

(13c)

Now, the probability (8) is restated as follows:

$$\begin{array}{cccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{P\left({\mathbf{R}}_{12},{\mathbf{R}}_{13},{\mathbf{R}}_{14}\right)=\int W_{n_{15}}^{\mathrm{C}}(|{\mathbf{R}}_{15}|)W_{n_{52}}^{\mathrm{C}}(|{\mathbf{R}}_{52}|)W^{\mathrm{R}}(|{\mathbf{R}}_{56}|)W_{n_{63}}^{\mathrm{C}}(|{\mathbf{R}}_{63}|)W_{n_{64}}^{\mathrm{C}}(|{\mathbf{R}}_{64}|)}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{\times \delta ({\mathbf{R}}_{12}-{\mathbf{R}}_{15}-{\mathbf{R}}_{52})\delta ({\mathbf{R}}_{13}-{\mathbf{R}}_{15}-{\mathbf{R}}_{56}-{\mathbf{R}}_{63})\delta ({\mathbf{R}}_{14}-{\mathbf{R}}_{15}-{\mathbf{R}}_{56}-{\mathbf{R}}_{64})}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{\times \mathrm{d}{\mathbf{R}}_{15}\mathrm{d}{\mathbf{R}}_{52}\mathrm{d}{\mathbf{R}}_{56}\mathrm{d}{\mathbf{R}}_{63}\mathrm{d}{\mathbf{R}}_{64},}& \multicolumn{1}{c}{}\end{array}$$

(14)

where the constraints (13) are enforced through three-dimensional Dirac delta functions $\delta (\mathbf{r})$ . It can be readily verified that the above integration can be transformed back into equation (8) by applying the sifting property of the three-dimensional Dirac delta function, together with a change of variables. By utilizing the following property of the three-dimensional Dirac delta function,

$$\begin{array}{ccc}\multicolumn{1}{c}{\delta ({\mathbf{R}}_{\mathrm{ij}})={(2\pi )}^{-3}\int \exp \left[-\mathit{i}{\mathbf{k}}_{\mathrm{ij}}\cdot {\mathbf{R}}_{\mathrm{ij}}\right]\mathrm{d}{\mathbf{k}}_{\mathrm{ij}},}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(15)

where $\mathit{i}=\sqrt{-1}$ is the imaginary unit, and rearranging the integration in (14), we obtain the following expression:

$$\begin{array}{cccc}\multicolumn{1}{c}{P}& \multicolumn{1}{c}{\left({\mathbf{R}}_{12},{\mathbf{R}}_{13},{\mathbf{R}}_{14}\right)={(2\pi )}^{-9}}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{\int \overline{W_{n_{15}}^{\mathrm{C}}}(|{\mathbf{k}}_{12}+{\mathbf{k}}_{13}+{\mathbf{k}}_{14}|)\overline{W_{n_{52}}^{\mathrm{C}}}(|{\mathbf{k}}_{12}|)\overline{W^{\mathrm{R}}}(|{\mathbf{k}}_{13}+{\mathbf{k}}_{14}|)\overline{W_{n_{63}}^{\mathrm{C}}}(|{\mathbf{k}}_{13}|)\overline{W_{n_{64}}^{\mathrm{C}}}(|{\mathbf{k}}_{14}|)}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{\times \exp \left[-\mathit{i}{\mathbf{k}}_{12}\cdot {\mathbf{R}}_{12}\right]\exp \left[-\mathit{i}{\mathbf{k}}_{13}\cdot {\mathbf{R}}_{13}\right]\exp \left[-\mathit{i}{\mathbf{k}}_{14}\cdot {\mathbf{R}}_{14}\right]\mathrm{d}{\mathbf{k}}_{12}\mathrm{d}{\mathbf{k}}_{13}\mathrm{d}{\mathbf{k}}_{14}.}& \multicolumn{1}{c}{}\end{array}$$

(16)

In (16), the bar over a symbol denotes the Fourier transform of the function defined as:

$$\begin{array}{ccc}\multicolumn{1}{c}{\overline{W}({\mathbf{k}}_{\mathrm{ij}})=\int W({\mathbf{R}}_{\mathrm{ij}})\exp \left(\mathit{i}{\mathbf{k}}_{\mathrm{ij}}\cdot {\mathbf{R}}_{\mathrm{ij}}\right)\mathrm{d}{\mathbf{R}}_{\mathrm{ij}}.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(17)

For the rigid rod with length a, the Fourier transform of equation (2) gives rise to

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{\overline{W^{\mathrm{R}}}(|\mathbf{k}|)=\frac{\sin \left(\left|\mathbf{k}\right|a\right)}{\left|\mathbf{k}\right|a}.}& \multicolumn{1}{c}{}\end{array}$$

(18)

2.4. Gaussian coil approximation

Thus far, the probability distribution has been expressed in the form of multi-variable integrations, involving a substantial numerical burden. In this section, it will be shown that Gaussian approximation of coils can significantly simplify the calculations.

The probability distribution function of the Gaussian phantom network (with no rods) has been obtained by Flory [10] using real space representation. In Appendix 1, the same network is formulated with the aid of Fourier transform and it is shown that the result has an analytical closed form. Likewise, the integration described in (16) can be simplified by assuming a Gaussian distribution for coils 15, 52, 63, and 46. Suppose that $W_n^{\mathrm{C}}(r)$ follows a Gaussian distribution (1). Then, the Fourier transform of equation (1) leads to

$$\begin{array}{ccc}\multicolumn{1}{c}{\overline{W_n^{\mathrm{C}}}(|\mathbf{k}|)=\exp \left[\frac{-n{b}^2{\left|\mathbf{k}\right|}^2}{6}\right].}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(19)

Let us represent the vectors ${\mathbf{R}}_{\mathrm{ij}}$ and ${\mathbf{k}}_{\mathrm{ij}}$ in the Cartesian coordinate system with unit vectors $({\widehat{\mathbf{e}}}_x,{\widehat{\mathbf{e}}}_y,{\widehat{\mathbf{e}}}_z)$ as

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{{\mathbf{R}}_{\mathrm{ij}}={\left(R_{\mathrm{ij}}\right)}_x{\widehat{\mathbf{e}}}_x+{\left(R_{\mathrm{ij}}\right)}_y{\widehat{\mathbf{e}}}_y+{\left(R_{\mathrm{ij}}\right)}_z{\widehat{\mathbf{e}}}_z,}& \multicolumn{1}{c}{}\end{array}$$

(20a)

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{{\mathbf{k}}_{\mathrm{ij}}={\left(k_{\mathrm{ij}}\right)}_x{\widehat{\mathbf{e}}}_x+{\left(k_{\mathrm{ij}}\right)}_y{\widehat{\mathbf{e}}}_y+{\left(k_{\mathrm{ij}}\right)}_z{\widehat{\mathbf{e}}}_z.}& \multicolumn{1}{c}{}\end{array}$$

(20b)

Now, by defining the matrices

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{{\mathbf{K}}^{\mathrm{T}}=\{{(k_{12})}_x{(k_{12})}_y{(k_{12})}_z{(k_{13})}_x{(k_{13})}_y{(k_{13})}_z{(k_{14})}_x{(k_{14})}_y{(k_{14})}_z\},}& \multicolumn{1}{c}{}\end{array}$$

(21a)

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{{\mathbf{R}}^{\mathrm{T}}=\{{(R_{12})}_x{(R_{12})}_y{(R_{12})}_z{(R_{13})}_x{(R_{13})}_y{(R_{13})}_z{(R_{14})}_x{(R_{14})}_y{(R_{14})}_z\},}& \multicolumn{1}{c}{}\end{array}$$

(21b)

with label T denotes the transpose of the matrix, the integration (16) is rephrased as the following compact form:

$$\begin{array}{ccc}\multicolumn{1}{c}{P\left({\mathbf{R}}_{12},{\mathbf{R}}_{13},{\mathbf{R}}_{14}\right)={(2\pi )}^{-9}}& \multicolumn{1}{c}{\int \overline{W^{\mathrm{R}}}(\sqrt{{\mathbf{K}}^{\mathrm{T}}\mathbf{\Lambda }\mathbf{K}})\exp [-\frac{1}{2}{\mathbf{K}}^{\mathrm{T}}\mathbf{\Gamma }\mathbf{K}-\mathit{i}{\mathbf{R}}^{\mathrm{T}}\mathbf{K}]\mathrm{d}\mathbf{K},}& \multicolumn{1}{c}{}\end{array}$$

(22)

where

$$\begin{array}{ccc}\multicolumn{1}{c}{\mathbf{\Gamma }=\frac{b^2}{3}\left(\begin{array}{ccc}\multicolumn{1}{c}{\left(n_{15}+n_{52}\right){\left[\mathbf{I}\right]}_{3\times 3}}& \multicolumn{1}{c}{n_{15}{\left[\mathbf{I}\right]}_{3\times 3}}& \multicolumn{1}{c}{n_{15}{\left[\mathbf{I}\right]}_{3\times 3}}\\ \multicolumn{1}{c}{n_{15}{\left[\mathbf{I}\right]}_{3\times 3}}& \multicolumn{1}{c}{\left(n_{15}+n_{63}\right){\left[\mathbf{I}\right]}_{3\times 3}}& \multicolumn{1}{c}{n_{15}{\left[\mathbf{I}\right]}_{3\times 3}}\\ \multicolumn{1}{c}{n_{15}{\left[\mathbf{I}\right]}_{3\times 3}}& \multicolumn{1}{c}{n_{15}{\left[\mathbf{I}\right]}_{3\times 3}}& \multicolumn{1}{c}{\left(n_{15}+n_{64}\right){\left[\mathbf{I}\right]}_{3\times 3}}\end{array}\right),}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(23)

and

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{\mathbf{\Lambda }=\left(\begin{array}{ccc}\multicolumn{1}{c}{{\left[\mathbf{0}\right]}_{3\times 3}}& \multicolumn{1}{c}{{\left[\mathbf{0}\right]}_{3\times 3}}& \multicolumn{1}{c}{{\left[\mathbf{0}\right]}_{3\times 3}}\\ \multicolumn{1}{c}{{\left[\mathbf{0}\right]}_{3\times 3}}& \multicolumn{1}{c}{{\left[\mathbf{I}\right]}_{3\times 3}}& \multicolumn{1}{c}{{\left[\mathbf{I}\right]}_{3\times 3}}\\ \multicolumn{1}{c}{{\left[\mathbf{0}\right]}_{3\times 3}}& \multicolumn{1}{c}{{\left[\mathbf{I}\right]}_{3\times 3}}& \multicolumn{1}{c}{{\left[\mathbf{I}\right]}_{3\times 3}}\end{array}\right).}& \multicolumn{1}{c}{}\end{array}$$

(24)

In the above relations, ${\left[\mathbf{I}\right]}_{3\times 3}$ and ${\left[\mathbf{0}\right]}_{3\times 3}$ , respectively, denote the $3\times 3$ identity and zero matrices. Furthermore, let us define the change of variable

$$\begin{array}{ccc}\multicolumn{1}{c}{\mathbf{Y}=\mathbf{A}\mathbf{K}}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(25)

where

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{{\mathbf{Y}}^{\mathrm{T}}=\{{(k_{12})}_x{(k_{12})}_y{(k_{12})}_z{(k_{13})}_x{(k_{13})}_y{(k_{13})}_z{w}_x{w}_y{w}_z\},}& \multicolumn{1}{c}{}\end{array}$$

(26a)

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{w_x={(k_{13})}_x+{(k_{14})}_x,w_y={(k_{13})}_y+{(k_{14})}_y,w_z={(k_{13})}_z+{(k_{14})}_z,}& \multicolumn{1}{c}{}\end{array}$$

(26b)

and

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{\mathbf{A}=\left(\begin{array}{ccc}\multicolumn{1}{c}{{\left[\mathbf{I}\right]}_{3\times 3}}& \multicolumn{1}{c}{{\left[\mathbf{0}\right]}_{3\times 3}}& \multicolumn{1}{c}{{\left[\mathbf{0}\right]}_{3\times 3}}\\ \multicolumn{1}{c}{{\left[\mathbf{0}\right]}_{3\times 3}}& \multicolumn{1}{c}{{\left[\mathbf{I}\right]}_{3\times 3}}& \multicolumn{1}{c}{{\left[\mathbf{0}\right]}_{3\times 3}}\\ \multicolumn{1}{c}{{\left[\mathbf{0}\right]}_{3\times 3}}& \multicolumn{1}{c}{{\left[\mathbf{I}\right]}_{3\times 3}}& \multicolumn{1}{c}{{\left[\mathbf{I}\right]}_{3\times 3}}\end{array}\right).}& \multicolumn{1}{c}{}\end{array}$$

(27)

Since $\det \mathbf{A}=1$ , the integration (22) is rewritten as:

$$\begin{array}{ccc}\multicolumn{1}{c}{P\left({\mathbf{R}}_{12},{\mathbf{R}}_{13},{\mathbf{R}}_{14}\right)={(2\pi )}^{-9}}& \multicolumn{1}{c}{\int \frac{\sin (\sqrt{{\mathbf{w}}^{\mathrm{T}}\mathbf{w}}a)}{\sqrt{{\mathbf{w}}^{\mathrm{T}}\mathbf{w}}a}\exp [-\frac{1}{2}{\mathbf{Y}}^{\mathrm{T}}\stackrel{\sim }{\mathbf{\Gamma }}\mathbf{Y}-\mathit{i}{\stackrel{\sim }{\mathbf{R}}}^{\mathrm{T}}\mathbf{Y}]\mathrm{d}\mathbf{Y},}& \multicolumn{1}{c}{}\end{array}$$

(28)

where

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{\mathbf{w}=w_x{\widehat{\mathbf{e}}}_x+w_y{\widehat{\mathbf{e}}}_y+w_z{\widehat{\mathbf{e}}}_z,}& \multicolumn{1}{c}{}\end{array}$$

(29)

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{\stackrel{\sim }{\mathbf{\Gamma }}={\mathbf{A}}^{-\mathrm{T}}\mathbf{\Gamma }{\mathbf{A}}^{-1}=\frac{b^2}{3}\left(\begin{array}{ccc}\multicolumn{1}{c}{\left(n_{15}+n_{52}\right){\left[\mathbf{I}\right]}_{3\times 3}}& \multicolumn{1}{c}{{\left[\mathbf{0}\right]}_{3\times 3}}& \multicolumn{1}{c}{n_{15}{\left[\mathbf{I}\right]}_{3\times 3}}\\ \multicolumn{1}{c}{{\left[\mathbf{0}\right]}_{3\times 3}}& \multicolumn{1}{c}{\left(n_{63}+n_{64}\right){\left[\mathbf{I}\right]}_{3\times 3}}& \multicolumn{1}{c}{-n_{64}{\left[\mathbf{I}\right]}_{3\times 3}}\\ \multicolumn{1}{c}{n_{15}{\left[\mathbf{I}\right]}_{3\times 3}}& \multicolumn{1}{c}{-n_{64}{\left[\mathbf{I}\right]}_{3\times 3}}& \multicolumn{1}{c}{\left(n_{15}+n_{64}\right){\left[\mathbf{I}\right]}_{3\times 3}}\end{array}\right),}& \multicolumn{1}{c}{}\end{array}$$

(30)

and $\stackrel{\sim }{\mathbf{R}}={\mathbf{A}}^{-\mathrm{T}}\mathbf{R}$ , i.e.,

$$\begin{array}{ccc}\multicolumn{1}{c}{{\stackrel{\sim }{\mathbf{R}}}^{\mathrm{T}}=\{{(R_{12})}_x{(R_{12})}_y{(R_{12})}_z{(R_{43})}_x{(R_{43})}_y{(R_{43})}_z{(R_{14})}_x{(R_{14})}_y{(R_{14})}_z\}.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(31)

Label $-\mathrm{T}$ in the superscript represents the inverse transpose of the matrix. Based on the results provided in Appendix 2 (equation (72) with $p=3$ ), equation (28) can be written in a closed form:

$$\begin{array}{cccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{P\left({\mathbf{R}}_{12},{\mathbf{R}}_{13},{\mathbf{R}}_{14}\right)=}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{\sqrt{\frac{{\tilde{\zeta }}_{33}}{{(2\pi )}^9\det {\stackrel{\sim }{\mathbf{\Gamma }}}_{\mathrm{aa}}}}\frac{1}{a\left|\stackrel{\sim }{\mathbf{L}}\right|}\exp \left[-\frac{1}{2}\left({\stackrel{\sim }{\mathbf{R}}}_{\mathrm{a}}^{\mathrm{T}}{\stackrel{\sim }{\mathbf{\Gamma }}}_{\mathrm{aa}}^{-1}{\stackrel{\sim }{\mathbf{R}}}_{\mathrm{a}}\right)-\frac{{\tilde{\zeta }}_{33}}{2}\left(a^2+{\left|\stackrel{\sim }{\mathbf{L}}\right|}^2\right)\right]\sinh \left[a\left|\stackrel{\sim }{\mathbf{L}}\right|{\stackrel{\mathrm{~}}{\zeta }}_{33}\right],}& \multicolumn{1}{c}{}\end{array}$$

(32)

in which

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{{\stackrel{\sim }{\mathbf{\Gamma }}}_{\mathrm{aa}}^{-1}=\frac{3}{b^2}\left(\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{{\left(n_{15}+n_{52}\right)}^{-1}{\left[\mathbf{I}\right]}_{3\times 3}}& \multicolumn{1}{c}{{\left[\mathbf{0}\right]}_{3\times 3}}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{{\left[\mathbf{0}\right]}_{3\times 3}}& \multicolumn{1}{c}{{\left(n_{63}+n_{64}\right)}^{-1}{\left[\mathbf{I}\right]}_{3\times 3}}\end{array}\right),}& \multicolumn{1}{c}{}\end{array}$$

(33a)

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{{\stackrel{\sim }{\mathbf{R}}}_{\mathrm{a}}^{\mathrm{T}}=\{{(R_{12})}_x{(R_{12})}_y{(R_{12})}_z{(R_{43})}_x{(R_{43})}_y{(R_{43})}_z\},}& \multicolumn{1}{c}{}\end{array}$$

(33b)

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{\stackrel{\sim }{\mathbf{L}}=\frac{n_{63}}{n_{63}+n_{64}}{\mathbf{R}}_{14}+\frac{n_{64}}{n_{63}+n_{64}}{\mathbf{R}}_{13}-\frac{n_{15}}{n_{15}+n_{52}}{\mathbf{R}}_{12},}& \multicolumn{1}{c}{}\end{array}$$

(33c)

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{{\tilde{\zeta }}_{33}=\frac{3}{b^2}\frac{\left(n_{15}+n_{52}\right)\left(n_{63}+n_{64}\right)}{n_{52}{n}_{63}{n}_{64}+n_{15}\left(n_{63}{n}_{64}+n_{52}\left(n_{63}+n_{64}\right)\right)},}& \multicolumn{1}{c}{}\end{array}$$

(33d)

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{\det {\stackrel{\sim }{\mathbf{\Gamma }}}_{\mathrm{aa}}=\frac{b^{12}}{3^6}{\left(n_{15}+n_{52}\right)}^3{\left(n_{63}+n_{64}\right)}^3.}& \multicolumn{1}{c}{}\end{array}$$

(33e)

${\stackrel{\sim }{\mathbf{\Gamma }}}_{\mathrm{aa}}$ is the $6\times 6$ matrix defined as the upper-left block of matrix $\stackrel{\sim }{\mathbf{\Gamma }}$ (see equation (100)). ${\stackrel{\sim }{\mathbf{R}}}_{\mathrm{a}}$ is defined as the $6\times 1$ block vector of $\stackrel{\sim }{\mathbf{R}}$ , as shown in equation (99). ${\tilde{\zeta }}_{33}$ is the the coefficient of the last block matrix ${\left[\mathbf{I}\right]}_{3\times 3}$ in ${\stackrel{\sim }{\mathbf{\Gamma }}}^{-1}$ analogous to equation (75). The position vectors ${\mathbf{R}}_{12}$ , ${\mathbf{R}}_{13}$ , and ${\mathbf{R}}_{14}$ along with the geometric interpretation of $\stackrel{\sim }{\mathbf{L}}$ are shown in Figure 3(a) and (b).

Figure 3.

(a) The topology of the four-node coil-rod structure with rod length a and corresponding end-to-end vectors ${\mathbf{R}}_{12}$ , ${\mathbf{R}}_{13}$ , and ${\mathbf{R}}_{14}$ . (b) The geometric interpretation of $\stackrel{\sim }{\mathbf{L}}$ . This vector is constructed based on (33c) composed of vectors ${\mathbf{R}}_{10}$ and ${\mathbf{R}}_{10^{\prime }}$ . (c) The equivalent network model with three chains: one two-node coil-rod and two pure coils based on equation (34).

The probability distribution of the four-node coil-rod structure can be generalized to ( $p+1$ )-node structure with one rod as shown in Appendix 3. Based on the interpretation given in Appendix 3, a ( $p+1$ )-node structure with one rod can be decomposed to $p-1$ coils and one two-node coil-rod structure. For the simple case of four-node coil-rod structure, as shown in Figure 3(c), the associated probability distribution can be decomposed to three chains: one two-node coil-rod structure with rod length a, $\frac{3}{b^2{\tilde{\zeta }}_{33}}$ Kuhn segments in the coil and end-to-end vector $\stackrel{\sim }{\mathbf{L}}$ ; one Gaussian coil with $n_{15}+n_{52}$ Kuhn segments and end-to-end vector ${\mathbf{R}}_{12}$ ; and the other Gaussian coil with $n_{63}+n_{64}$ Kuhn segments and end-to-end vector ${\mathbf{R}}_{43}$ . With the aid of (3) and (1), the probability distribution function (32) is rewritten as:

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{P\left({\mathbf{R}}_{12},{\mathbf{R}}_{13},{\mathbf{R}}_{14}\right)=W_{n_{15}+n_{52}}^{\mathrm{C}}\left(\left|{\mathbf{R}}_{12}\right|\right)W_{n_{63}+n_{64}}^{\mathrm{C}}\left(\left|{\mathbf{R}}_{43}\right|\right)W_{3b^{-2}{\tilde{\zeta }}_{33}^{-1}}^{\mathrm{CR}}\left(\left|\stackrel{\sim }{\mathbf{L}}\right|\right).}& \multicolumn{1}{c}{}\end{array}$$

(34)

The above representation not only makes the notation simpler but also lays the solid base for comparing the present models with models in the literature based on two-node coil-rod structure and pure coils. Now suppose that the chains 15 and 52, as well as chains 63 and 46 are identical such that $n_{15}=n_{52}=n$ and $n_{63}=n_{46}=n+2m$ , where n is the number of Kuhn segments in chains 15 and 52, while $2m$ indicates the extra number of segments 46 and 63 compared to 15 and 52. Then, equations (33c) and (33d) can be further simplified to

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{\stackrel{\sim }{\mathbf{L}}=\frac{1}{2}\left({\mathbf{R}}_{14}+{\mathbf{R}}_{13}-{\mathbf{R}}_{12}\right),}& \multicolumn{1}{c}{}\end{array}$$

(35a)

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{{\tilde{\zeta }}_{33}=\frac{3}{\left(n+m\right)b^2}.}& \multicolumn{1}{c}{}\end{array}$$

(35b)

Therefore, equation (34) is rewritten as

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{P\left({\mathbf{R}}_{12},{\mathbf{R}}_{13},{\mathbf{R}}_{14}\right)=W_{2n}^{\mathrm{C}}\left(\left|{\mathbf{R}}_{12}\right|\right)W_{2n+4m}^{\mathrm{C}}\left(\left|{\mathbf{R}}_{43}\right|\right)W_{n+m}^{\mathrm{CR}}\left(\frac{1}{2}\left|{\mathbf{R}}_{14}+{\mathbf{R}}_{13}-{\mathbf{R}}_{12}\right|\right).}& \multicolumn{1}{c}{}\end{array}$$

(36)

3. Implementation into network

As mentioned earlier, once the probability distribution is established, the mechanical response of the gel can be determined by placing the four-node coil-rod structure within a macroscopic network with specific geometry. In this manner, the macroscopic deformation can be properly translated to the constituent structure. To examine the features of the four-node coil-rod structure, a simple scenario is first considered where the four nodes are placed in the same plane in the form of a rectangular cell as seen in Figure 4(a). It should be noted that while the nodes form a two-dimensional cell, the chains can occupy the three-dimensional space.

Figure 4.

(a) A four-node coil-rod structure with the four nodes placed in the same plane in a rectangular manner with deformed dimensions ${\lambda }_1{l}_1$ and ${\lambda }_2{l}_2$ . (b) The Helmholtz free energy contour $\left(\mathrm{\Psi }+{\mathrm{\Psi }}_0\right)∕G$ of the four-node coil-rod structure with $n=25$ , $a=5b$ , showing the minimum at $l_1{\lambda }_1=4.09b$ and $l_2{\lambda }_2=6.33b$ ,

3.1. In-plane placement of four nodes

Referring to Figure 4(a), suppose the nodes 1 to 4 are initially positioned at the vertices of a rectangle with edge lengths of $l_1$ and $l_2$ . To avoid any unnecessary complexity, the number of Kuhn segments of all coils is set to be equal to n. Deformations are assumed to be applied such that the principal stretches ${\lambda }_1$ and ${\lambda }_2$ align with the horizontal and vertical edges of the rectangle. By setting the Cartesian coordinate system at node 1 with coordinate axes along the rectangle edges, the position vectors of the nodes of the deformed cell can be expressed as:

$$\begin{array}{ccccccccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{{\mathbf{R}}_1=(0,0,0),}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{{\mathbf{R}}_2=({\lambda }_1{l}_1,0,0),}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{{\mathbf{R}}_3=({\lambda }_1{l}_1,{\lambda }_2{l}_2,0),}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{{\mathbf{R}}_4=(0,{\lambda }_2{l}_2,0),}& \multicolumn{1}{c}{}\end{array}$$

(37)

where setting ${\lambda }_1={\lambda }_2=1$ reverts the cell to its undeformed state. Such representation is capable of capturing any in-plane deformations including shear and extension within the xy plane. By substituting (37) in equation (36), and considering an equal number of Kuhn segments n for the coils ( $m=0$ ), one can conclude that

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{P({\mathbf{R}}_{12},{\mathbf{R}}_{13},{\mathbf{R}}_{14})=W_n^{\mathrm{CR}}({\lambda }_2{l}_2)W_{2n}^{\mathrm{C}}({\lambda }_1{l}_1)W_{2n}^{\mathrm{C}}({\lambda }_1{l}_1).}& \multicolumn{1}{c}{}\end{array}$$

(38)

Now, by substitution of equation (38) into equations (4) and (5), the Helmholtz free energy can be calculated as follows:

$$\begin{array}{ccc}\multicolumn{1}{c}{\mathrm{\Psi }+{\mathrm{\Psi }}_0=-G\ln \left[{\lambda }_1^2{\lambda }_2\right]-G\ln \left[\sinh \left(\frac{3a{\lambda }_2{l}_2}{n{b}^2}\right)\right]+\frac{3G}{2n{b}^2}\left(l_1^2{\lambda }_1^2+l_2^2{\lambda }_2^2\right),}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(39)

in which $G=\varrho k_{B}T$ and $\varrho =1∕(l_1{l}_2)$ . The first term on the right hand side is obtained after combining the first term of (4) and the denominator of $W_n^{\mathrm{CR}}({\lambda }_2{l}_2)$ in (3). ${\mathrm{\Psi }}_0$ is the Helmholtz free energy per unit reference volume where ${\lambda }_1={\lambda }_2=1$ so that the condition $\mathrm{\Psi }({\lambda }_1={\lambda }_2=1)=0$ is imposed. Based on the free energy (39), equation (6) leads to

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{P_1=\frac{\partial \mathrm{\Psi }}{\partial {\lambda }_1}=G\left(\frac{3l_1^2{\lambda }_1}{n{b}^2}-\frac{2}{{\lambda }_1}\right),}& \multicolumn{1}{c}{}\end{array}$$

(40)

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{P_2=\frac{\partial \mathrm{\Psi }}{\partial {\lambda }_2}=G\left(\frac{3l_2^2{\lambda }_2}{n{b}^2}-\frac{1}{{\lambda }_2}-\frac{3a{l}_2}{n{b}^2}\coth \left[\frac{3a{l}_2{\lambda }_2}{n{b}^2}\right]\right).}& \multicolumn{1}{c}{}\end{array}$$

(41)

The distinct stress–stretch relationship in each principal direction demonstrates the intrinsic anisotropy of the four-node coil-rod structure shown in Figure 4(a). The above stress components can be set to zero for ${\lambda }_1={\lambda }_2=1$ to ensure that the initial state is stress-free. This results in the following conditions for the stress-free size of the four-node coil-rod structure,

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{l_1=\sqrt{\frac{2}{3}n}b,}& \multicolumn{1}{c}{}\end{array}$$

(42a)

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{\frac{3l_2^2}{n{b}^2}-\frac{3a{l}_2}{n{b}^2}\coth \left[\frac{3a{l}_2}{n{b}^2}\right]=1.}& \multicolumn{1}{c}{}\end{array}$$

(42b)

Equation (42b) is an implicit function of a and typically requires numerical computations. However, for $a>\sqrt{n{b}^2}$ , the approximation $\coth \left[\frac{3a{l}_2}{n{b}^2}\right]\simeq 1$ can be utilized, which allows for an analytical solution for $l_2$ :

$$\begin{array}{ccc}\multicolumn{1}{c}{l_2=\frac{a}{2}+\frac{\sqrt{3a^2+4n{b}^2}}{\sqrt{12}}.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(43)

In this case, as $n\to 0$ , $l_2\to a$ verifying that the pure rod with no coil forms a rectangle with dimension $l_1=0$ and $l_2=a$ . In addition, by applying Padé approximation about $a=0$ to the second order for $\coth (x)$ , one can show that

$$\begin{array}{ccc}\multicolumn{1}{c}{l_2=\sqrt{\frac{n{b}^2}{6a^2}}\sqrt{7a^2-5n{b}^2+\sqrt{25n^2{b}^4-30a^{2}n{b}^2+49a^4}}}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(44)

holds for $a<0.5\sqrt{n{b}^2}$ . For $a\ll \sqrt{n{b}^2}$ , $l_2\to l_1$ , in this case the stress-free state of the rectangular cell becomes a square for the four-node structure with a cross-link rather than a junction zone. For a four-node coil-rod structure with $n=25$ and $a=5b$ , the contour of $\left(\mathrm{\Psi }+{\mathrm{\Psi }}_0\right)∕G$ is depicted in Figure 4(b) as a function of ${\lambda }_1{l}_1∕b$ and ${\lambda }_2{l}_2∕b$ . The minimum of Helmholtz free energy is indicated by the red point on the contour. By setting ${\lambda }_1={\lambda }_2=1$ , the red point corresponds to $l_1=0.82\sqrt{n{b}^2}$ and $l_2=1.27\sqrt{n{b}^2}$ , consistent with equations (42a) and (44), respectively, thereby signifying the stress-free state.

As it was alluded to, the applicability of Gaussian coils in the freely jointed chain model is constrained to small extensions of the coils and a large number of Kuhn segments. To assess the validity of such approximation in determining $l_1$ and $l_2$ , the exact Helmholtz free energy of the system is minimized based on equation (11) combined with (4) and (5). To obtain the exact result, all chains 15, 52, 63, and 46 in equation (11) with the same number of Kuhn segments n, are assumed to obey the following non-Gaussian probability distribution [18]

$$\begin{array}{ccc}\multicolumn{1}{c}{W_{n_{15}}^{\mathrm{C}}(\left|\mathbf{r}\right|)=W_{n_{52}}^{\mathrm{C}}(\left|\mathbf{r}\right|)=W_{n_{63}}^{\mathrm{C}}(\left|\mathbf{r}\right|)=W_{n_{46}}^{\mathrm{C}}(\left|\mathbf{r}\right|)=\frac{A_0\beta }{\left|\mathbf{r}\right|}{\left(\frac{\sinh \beta }{\beta }\right)}^n\exp \left[-\frac{\beta \left|\mathbf{r}\right|}{b}\right],}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(45)

with

$$\begin{array}{ccc}\multicolumn{1}{c}{\beta ={ℒ}^{-1}\left(\frac{\left|\mathbf{r}\right|}{\mathrm{nb}}\right),}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(46)

and $A_0$ being a normalization factor. $ℒ(x)=\coth x-1∕x$ denotes the Langevin function and to calculate its inverse, the Padé approximation [19] is utilized. The five-tuple integration (11) is performed by employing the adaptive quadrature techniques [20]. To obtain $l_1$ and $l_2$ associated with the minimum free energy, the differential evolution algorithm [21] is adopted, with a convergence tolerance of 0.01. The initial length $l_1$ and $l_2$ are depicted vs. n for $a=5b$ in Figures 5(a) and 5(b), respectively, where the results are compared with the Gaussian approximation (42a) and (42b). The excellent agreement suggests that the Gaussian approximation can be satisfactorily applied to determine $l_1$ and $l_2$ .

Figure 5.

(a) $l_1$ vs. n, (b) $l_2$ vs. n for $a=5b$ . Comparison is made between the Gaussian approximations (42a) and (42b) and exact result based on equation (11).

Next, the accuracy of the constitutive relations (40) and (41) is validated against the exact result without the Gaussian approximation. As mentioned earlier, after calculation of the exact Helmholtz free energy via equations (11), (4), (5), and (45), the stress components can be determined through equation (6). Since the numerical burden of evaluating (11) is heavy, the biaxial test is preferred over the uniaxial test. In the biaxial test, $P_1$ and $P_2$ are directly obtained from known ${\lambda }_1$ and ${\lambda }_2$ , unlike the uniaxial test, which requires solving for unknown stretches to achieve zero stress in the direction perpendicular to the loading axis. Figure 6(a) shows $P_1∕G$ as a function of ${\lambda }_1$ for $n=25$ and $a=5b$ , comparing the result based on Gaussian distribution (38) (dash-dotted black) with those based on non-Gaussian distribution (11) for ${\lambda }_2=1$ (solid blue) and ${\lambda }_2=4$ (dashed yellow). Note that under the Gaussian approximation (equation (40)), $P_1$ is a function of ${\lambda }_1$ only and does not depend on ${\lambda }_2$ . Compared to the solid blue curve for ${\lambda }_2=1$ , the Gaussian approximation exhibits minimal error in both compression and tension up to ${\lambda }_1=3$ . Beyond this point, as the coils within the cell undergo significant elongation, prediction from the Gaussian approximation starts deviating from the exact result. The discrepancy is even more pronounced when compared to the dashed yellow curve ( ${\lambda }_2=4$ ), where the Gaussian model fails to capture the effect of increased extension along direction 2. For the non-Gaussian distribution, a larger extension in direction 2 reduces the material’s deformability in direction 1, leading to a more rapid increase in stress. In contrast, the Gaussian probability distribution assigns nonzero values even at large extensions, making it incapable of capturing the coupling effects between stretches in different principal directions. However, in compression ( ${\lambda }_1<1$ ), such dependency vanishes, and the Gaussian approximation renders good agreement with the exact treatment. Figure 6(b) demonstrates similar curves for $P_2∕G$ as a function of ${\lambda }_2$ . Unlike direction 1 where the rod length is inconsequential (see equation (40)), in this case (equation (41)), the rod length plays a role, leading to increased stiffness under both tension ( ${\lambda }_2>1$ ) and compression ( ${\lambda }_2<1$ ) compared to Figure 6(a). Similar to Figure 6(a), the Gaussian approximation which removes the dependence of $P_2$ on ${\lambda }_1$ (equation (41)) agrees with the exact results in compression and small tension ( ${\lambda }_2<2$ ). However, the results deviate from the non-Gaussian curves for ${\lambda }_2>2$ , and the discrepancy is higher for larger ${\lambda }_1$ .

Figure 6.

The normalized principal stress vs. the principal stretch for biaxial test and $a=5b$ , $l_2=6.33b$ , $l_1=4.09b$ with and without Gaussian approximation, (a) $P_1∕G$ vs. ${\lambda }_1$ with ${\lambda }_2=1$ and 4 (b) $P_2∕G$ vs. ${\lambda }_2$ with ${\lambda }_1=1$ and 4.

3.2. Extension to 3D

Based on the formulation in Section 3.1, the four-node coil-rod structure shown in Figure 4(a) can be incorporated into the three-dimensional network to account for more general deformation states. One example is shown in Figure 7, where identical four-node coil-rod structures are placed, with different orientations, on three faces of the rectangular prism: 1234, 1265, and 1485. Suppose that the principal stretches are aligned with the prism edges. The arrangement of the four-node coil-rod structures is such that the prism has a cube shape with edge length h when ${\lambda }_1={\lambda }_2={\lambda }_3=1$ in the stress-free state. In the Cartesian coordinate system fixed at the center of the cell, the positions of ${\mathbf{R}}_i$ are specified as follows:

$$\begin{array}{cccccccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{{\mathbf{R}}_1=\left(-\frac{h{\lambda }_1}{2},-\frac{h{\lambda }_2}{2},-\frac{h{\lambda }_3}{2}\right),}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{{\mathbf{R}}_2=\left(\frac{h{\lambda }_1}{2},-\frac{h{\lambda }_2}{2},-\frac{h{\lambda }_3}{2}\right),}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

$$\begin{array}{cccccccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{{\mathbf{R}}_3=\left(\frac{h{\lambda }_1}{2},\frac{h{\lambda }_2}{2},-\frac{h{\lambda }_3}{2}\right),}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{{\mathbf{R}}_4=\left(-\frac{h{\lambda }_1}{2},\frac{h{\lambda }_2}{2},-\frac{h{\lambda }_3}{2}\right),}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{{\mathbf{R}}_5=\left(-\frac{h{\lambda }_1}{2},-\frac{h{\lambda }_2}{2},\frac{h{\lambda }_3}{2}\right),}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{{\mathbf{R}}_6=\left(\frac{h{\lambda }_1}{2},-\frac{h{\lambda }_2}{2},\frac{h{\lambda }_3}{2}\right),}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{{\mathbf{R}}_7=\left(\frac{h{\lambda }_1}{2},\frac{h{\lambda }_2}{2},\frac{h{\lambda }_3}{2}\right),}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{{\mathbf{R}}_8=\left(-\frac{h{\lambda }_1}{2},\frac{h{\lambda }_2}{2},\frac{h{\lambda }_3}{2}\right).}& \multicolumn{1}{c}{}\end{array}$$

(47)

Figure 7.

Under the Gaussian approximation, the topology of 3 four-node coil-rod structures in the faces of the rectangular prism is equivalent to combination of three three-chain networks.

Similar to (4), the entropy of the deformed cell per unit reference volume can be written as

$$\begin{array}{cccc}\multicolumn{1}{c}{S=c}& \multicolumn{1}{c}{+\frac{k_B}{h^3}\ln \left[{\left(\mathrm{d}{\mathbf{R}}_{12}\right)}^2\mathrm{d}{\mathbf{R}}_{13}{\left(\mathrm{d}{\mathbf{R}}_{14}\right)}^2{\left(\mathrm{d}{\mathbf{R}}_{15}\right)}^2\mathrm{d}{\mathbf{R}}_{16}\mathrm{d}{\mathbf{R}}_{18}\right]}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{+\frac{k_B}{h^3}\ln \left[P({\mathbf{R}}_{12},{\mathbf{R}}_{15},{\mathbf{R}}_{16})P({\mathbf{R}}_{15},{\mathbf{R}}_{14},{\mathbf{R}}_{18})P({\mathbf{R}}_{14},{\mathbf{R}}_{12},{\mathbf{R}}_{13})\right].}& \multicolumn{1}{c}{}\end{array}$$

(48)

By defining $\varrho =3∕h^3$ (three chains in a reference volume of $h^3$ ), and using (5), it can be shown that

$$\begin{array}{cccc}\multicolumn{1}{c}{\mathrm{\Psi }=-{\mathrm{\Psi }}_0}& \multicolumn{1}{c}{-G\ln \left[{\lambda }_1{\lambda }_2{\lambda }_3\right]+\frac{G{h}^2}{n{b}^2}\left({\lambda }_1^2+{\lambda }_2^2+{\lambda }_3^2\right)}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{-\frac{G}{3}\ln \left[\sinh \left(\frac{3\mathrm{ah}{\lambda }_1}{n{b}^2}\right)\sinh \left(\frac{3\mathrm{ah}{\lambda }_2}{n{b}^2}\right)\sinh \left(\frac{3\mathrm{ah}{\lambda }_3}{n{b}^2}\right)\right],}& \multicolumn{1}{c}{}\end{array}$$

(49)

where $G=\varrho k_{B}T$ , and ${\mathrm{\Psi }}_0$ is the free energy per unit reference volume ensuring that the condition $\mathrm{\Psi }({\lambda }_1={\lambda }_2={\lambda }_3=1)=0$ is satisfied. Now using (6), the principal stresses are

$$\begin{array}{ccc}\multicolumn{1}{c}{P_i=G\left(\frac{-1}{{\lambda }_i}+\frac{2h^2{\lambda }_i}{n{b}^2}-\frac{\mathrm{ah}}{n{b}^2}\coth \left[\frac{3\mathrm{ah}{\lambda }_i}{n{b}^2}\right]\right),i=1,2,3.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(50)

Analogous to the treatment in Section 3.1, the stress-free state at ${\lambda }_1={\lambda }_2={\lambda }_3=1$ can provide the value of h obtained by solving the following equation:

$$\begin{array}{ccc}\multicolumn{1}{c}{\frac{2h^2}{n{b}^2}-\frac{\mathrm{ah}}{n{b}^2}\coth \left[\frac{3\mathrm{ah}}{n{b}^2}\right]=1.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(51)

By decomposing the four-node coil-rod structure positioned at each prism face into one two-node coil-rod and two coil structures, one can conclude that the three-dimensional network model is equivalent to the summation of three three-chain networks (Figure 7): one containing three two-node coil-rod structures on the prism edges, each having rod length a and n Kuhn segments in the coil; and two identical networks each containing three Gaussian coils with $2n$ Kuhn segments. The integration of these three networks yields the Helmholtz free energy (49).

In the literature, there is a correspondence between three-chain network models and a well-known type of material referred to as Valanis-Landel [22], where the strain energy is written as the sum of three scalar functions, each evaluated independently for three principal stretches. Due to the simplicity of this model, there are some generalizations for such network models, as shown by Ehret and Stracuzzi [23]. However, more involved network models, such as the eight-chain model, demonstrate superior performance in handling different types of experiments [7]. The comparison and relationship between the three-chain model and the eight-chain model are provided by Carroll [24] for pure coils without the presence of a junction zone. To examine the performance of the proposed network model (Figure 7), we formulate the eight-chain network model of the two-node coil-rod structure based on Moosavian and Tang [14], but with Gaussian approximation applied to the coil. Following the original model of Arruda and Boyce [7], eight two-node coil-rod structures are placed inside a diagonals of the cell with edges ${\lambda }_{1}h$ , ${\lambda }_{2}h$ , and ${\lambda }_{3}h$ . The probability of finding a two-node coil-rod with one end fixed at the center and the other inside the volume element $\mathrm{d}{\mathbf{R}}_1$ is given by $W_n^{\mathrm{CR}}\left(\left|{\mathbf{R}}_1\right|\right)\mathrm{d}{\mathbf{R}}_1$ . Likewise, the probability of the other seven chains is established. Therefore, the entropy of the deformed state per unit reference volume is given by

$$\begin{array}{ccc}\multicolumn{1}{c}{S=c+\frac{\varrho k_B}{8}\sum _{i=1}^8\ln \left[W_n^{\mathrm{CR}}\left(\left|{\mathbf{R}}_i\right|\right)\mathrm{d}{\mathbf{R}}_i\right],}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(52)

where $\varrho =8∕h^3$ . Now by applying (3), (5), and relation (47), one can conclude that

$$\begin{array}{cccc}\multicolumn{1}{c}{\mathrm{\Psi }+{\mathrm{\Psi }}_0=}& \multicolumn{1}{c}{-G\ln [{\lambda }_1{\lambda }_2{\lambda }_3]-G\ln \left[\frac{\sinh \left(\frac{3\mathrm{ah}}{2n{b}^2}\sqrt{{\lambda }_1^2+{\lambda }_2^2+{\lambda }_3^2}\right)}{\sqrt{{\lambda }_1^2+{\lambda }_2^2+{\lambda }_3^2}}\right]}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{+\frac{3G{h}^2}{8n{b}^2}\left({\lambda }_1^2+{\lambda }_2^2+{\lambda }_3^2\right),}& \multicolumn{1}{c}{}\end{array}$$

(53)

where ${\mathrm{\Psi }}_0$ is the Helmholtz free energy at ${\lambda }_1={\lambda }_2={\lambda }_3=1$ yielding $\mathrm{\Psi }({\lambda }_1={\lambda }_2={\lambda }_3=1)=0$ . By employing (6), the principal components of stress are written as

$$\begin{array}{cccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{P_i=G\left(-\frac{1}{{\lambda }_i}+\frac{3h^2{\lambda }_i}{4n{b}^2}+\frac{{\lambda }_i}{{\lambda }_1^2+{\lambda }_2^2+{\lambda }_3^2}-\frac{3\mathrm{ah}{\lambda }_i}{2n{b}^2}\frac{\coth \left[\frac{3\mathrm{ah}}{2n{b}^2}\sqrt{{\lambda }_1^2+{\lambda }_2^2+{\lambda }_3^2}\right]}{\sqrt{{\lambda }_1^2+{\lambda }_2^2+{\lambda }_3^2}}\right),}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{i=1,2,3.}& \multicolumn{1}{c}{}\end{array}$$

(54)

Setting $P_1=P_2=P_3=0$ for ${\lambda }_1={\lambda }_2={\lambda }_3=1$ gives rise to the following relation between h, a, and $\sqrt{n}b$ :

$$\begin{array}{ccc}\multicolumn{1}{c}{\frac{2}{3}=\frac{3h^2}{4n{b}^2}-\frac{\sqrt{3}\mathrm{ah}}{2n{b}^2}\coth \left(\frac{3\sqrt{3}\mathrm{ah}}{2n{b}^2}\right).}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(55)

Using the above relation, the initial stress-free state of the network with cube edge length h can be determined for given values of a and n.

Figure 8 illustrates the difference between $P_1$ vs. ${\lambda }_1$ in the uniaxial test for equations (50) and (54) with $a=5b$ and $n=25$ . To ensure a fair comparison, for equation (50), $G=0.73\mathrm{kPa}$ , while for (54), $G=1\mathrm{kPa}$ , such that the shear modulus (curve slope) of the biopolymer networks is identical at ${\lambda }_1=1$ . Furthermore, the values of h are obtained by solving equations (51) and (55), and the results are $h=5.00b$ and $h=8.40b$ , respectively. Although, in compression the results show indistinguishable behavior, in tension the eight-chain model containing two-node coil-rod structures exhibits larger stress for ${\lambda }_1>1.9$ . The results are also compared with the eight-chain model containing non-Gaussian two-node coil-rod structures proposed by Moosavian and Tang [15], as shown by the dashed curve. The curve is generated by applying equations (5) and (47) in (52), with the non-Gaussian distribution for the coil rather than equation (3) (refer to equations (4) and (17) in [15]). Non-Gaussian effects emerge when the sample is stretched beyond ${\lambda }_1=2.0$ (separation of the dashed and solid curves), while the results remain consistent below this stretch and under compression. Based on these results, it can be concluded that the overall behavior of the proposed network model is similar to that of the eight-chain network model when Gaussian statistics are assumed for the coils.

Figure 8.

$P_1$ vs. ${\lambda }_1$ in the uniaxial test with $a=5b$ and $n=25$ . For equation (50), $G=0.73\mathrm{kPa}$ and $h=5.00b$ , while for the eight-chain network (equation (54) and the result of Moosavian and Tang [15]), $G=1\mathrm{kPa}$ and $h=8.40b$ .

4. Discussion

This study aimed to provide explicit modeling of interactions within junction zones, and it can be concluded that under the Gaussian statistics assumption, the two-node coil-rod structure network can be effectively applied. Based on the current framework, it is also possible to construct other three-dimensional network models, such as the symmetric eight-node structure with a rod, as elaborated in Appendix 3. However, a more detailed comparison between different networks requires testing them on real disordered biopolymer gels, where various mechanical tests are conducted on the samples.

It should be recalled that the phantom network model assumes that components can easily pass through each other. Therefore, the interaction between coils and rods can be extended by considering long-range interaction effects, such as excluded volume effects or entanglement between components [25–27], which would render the governing formulation more realistic but also computationally more challenging. In addition, the current work is based on the presence of a single rod within the unit cell. For a more comprehensive generalization, the inclusion of multiple rods and their associated configurations can be considered.

The unzipping behavior in disordered biopolymer gels is another intriguing topic that warrants future investigation. In this phenomenon, the rod can exchange segments with coils, resulting in a new configuration with longer coils and shorter rods. The conversion of the four-node coil-rod structure to a two-node coil-rod structure and two coils allows for the direct application of the previously developed unzipping formulation to the current framework [14]. Moreover, the Gaussian approximation for the coils is a more reasonable assumption in the presence of unzipping, as unzipping allows the system to adopt a new configuration to avoid high extension or entering the non-Gaussian regime. Finally, the current framework can incorporate different unzipping along different principal directions—an anisotropic feature that likely exists in reality but absent in the model developed by Moosavian and Tang [15].

5. Conclusion

This work provides preliminary insight into modeling junction zones and their explicit interactions with the amorphous regions in disordered biopolymer gels. The previous model of a two-node coil-rod structure is generalized to a four-node coil-rod structure, where two chains intertwine to form a common rod as a junction zone. The probability distribution function in the context of a phantom network is developed, which is mathematically more involved than the two-node coil-rod structure. Nevertheless, within the Gaussian statistics of the coils, the analytical probability distribution function can be derived with the aid of Fourier transform. It is demonstrated that, under this assumption, the four-node coil-rod structure is equivalent to a collection of two coils and one two-node coil-rod structure situated inside an affine network. This argument can also be generalized to multiple Gaussian coils sharing a single rod. To examine the characteristics of this structure, the four nodes are placed on the same plane to form a rectangular cell. The stress–stretch relationship of this network is derived, and the applicability of the Gaussian distribution is evaluated through comparison with the exact distribution. Finally, a three-dimensional network model is developed by arranging the four-node coil-rod structures on the faces of a rectangular prism. Results from this model is compared with the eight-chain network model combined with two-node coil-rod structure developed by Moosavian and Tang [14], which confirms that the network model with the two-node coil-rod structure is suitable for investigating the effects of more explicit interactions between coils and rods, as long as the Gaussian approximation for the coils remains valid. The proposed model provides a systematic framework for developing more complex networks that involve collections of coils connected by multiple rods.