## 2 Development of Beam Element-Based Simulation
### Model Description
#### 2.1.2 Multipoint Constraints  
Because the beads are modeled using beams, which are one-dimensional elements, additional strategies must be employed to transfer forces between the beads. Vertical forces are transferred via vertical springs connecting a beam's end nodes to those of the beams immediately above and below it. Interface shear force transfer between beads is tackled with multipoint constraints (MPCs). For example, MPCs have been used for shear transfer in a composite beam [71], shear transfer between a bridge deck and girder [72] and during the simulation of warping of a composite beam's cross-section [73]. In the current model, multipoint constraints are employed to ensure strain compatibility at every bead-bead interface. A pair of constrained elements, with notation used throughout this paper, is shown in Figure 4.  
Full composite action is assumed between layers, meaning beads do not slip relative to each other. This implies that from the time of activation of the upper element, the change in longitudinal strain at the top of the layer below (the lower element in Figure 4) is equal to the change in longitudinal strain at the bottom of the layer above (the upper element in Figure 4). The total longitudinal strain at any point within a beam element \(\epsilon(x,y)\) is given in Equation (6),  
\[\epsilon(x,y)=\epsilon_{a}-v^{\prime\prime}(x)y, \tag{6}\]  
where \(e_{a}\) is the axial strain, \(v^{\prime\prime}\) is the element bending curvature and \(y\) is the signed distance from the neutral axis of the beam (centerline of the bead). Using the conventional Euler-Bernoulli beam element, axial strain is constant along the length of the beam and is given by Equation (7),  
\[e_{a}=\frac{u_{j}-u_{i}}{L}, \tag{7}\]  
where \(u_{i}\) and \(u_{j}\) are the axial displacements of element starting and ending nodes \(i\) and \(j\), respectively, and \(L\) is the element length. The curvature of a beam element is given in Equation (8),  
\[v^{\prime\prime}(x)=N_{b}^{\prime\prime}(x)U^{el}, \tag{8}\]  
where \(N_{b}^{\prime\prime}\) as defined previously is the second derivative of the element bending shape functions with respect to \(x\) and \(U^{el}=\left[ \begin{array}{cccc}u_{j}&v_{i}&\theta_{i}&u_{j}&v_{j}&\theta_{j }\end{array} \right]^{T}\) is the element displacement vector of \(x\)-direction nodal displacements \(u\), \(y\)-direction nodal displacements \(v\) and rotations \(\theta\) at nodes \(i\) and \(j\) of an Euler beam element. Since axial strain is a function of \(x\)-direction nodal displacements, the corresponding coefficients can be added to \(yN_{b}^{\prime\prime}\) in the correct locations to obtain \(B\), the matrix that relates strains to \(U^{el}\). Total longitudinal strain at any point in a beam element can then be expressed as shown in Equation (9),  
\[\epsilon(x,y)=B(x,y)U^{el}. \tag{9}\]  
Figure 4: Constrained elements. To enforce strain compatibility between two beads in the model at a location \(x\) at any time \(t\) after the activation of the upper element at \(t_{u}^{act}\), the strain at the bottom of the upper element \(e_{u}\) must equal the strain accumulated in the top of the lower element \(e_{l}\) following activation of the upper element. Equation (10) offers the full constraint expression at any point \(x\) within a pair of upper and lower elements that must be enforced to ensure inter-bead shear transfer and \(\epsilon_{u}=\epsilon_{l}\).