Asymptotically consistent and computationally efficient modeling of short-ranged molecular interactions between curved slender fibers undergoing large 3D deformations

This article proposes a novel computational modeling approach for short-ranged molecular interactions between curved slender fibers undergoing large 3D deformations, and gives a detailed overview how it fits into the framework of existing fiber or beam interaction models, either considering microscale molecular or macroscale contact effects. The direct evaluation of a molecular interaction potential between two general bodies in 3D space would require to integrate molecule densities over two 3D volumes, leading to a sixfold integral to be solved numerically. By exploiting the short-range nature of the considered class of interaction potentials as well as the fundamental kinematic assumption of undeformable fiber cross-sections, as typically applied in mechanical beam theories, a recently derived, closed-form analytical solution is applied for the interaction potential between a given section of the first fiber (slave beam) and the entire second fiber (master beam). This novel approach based on a pre-defined section-beam interaction potential (SBIP) requires only one single integration step along the slave beam length to be performed numerically. In terms of accuracy, the total beam-beam interaction potential resulting from this approach is shown to exhibit an asymptotically consistent angular and distance scaling behavior. In addition to elementary two-fiber systems, carefully chosen to verify accuracy and asymptotic consistence of the proposed SBIP approach, a potential practical application in form of adhesive nanofiber-grafted surfaces is studied. Involving a large number of helicoidal fibers undergoing large 3D deformations, arbitrary mutual fiber orientations as well as frequent local fiber pull-off and snap-into-contact events, this example demonstrates the robustness and computational efficiency of the new approach.


Introduction
This work is motivated by the abundance and manifoldness of biological, fiber-like structures on the nano-and microscale, including filamentous actin, collagen, and DNA, among others.These slender, deformable fibers form a variety of complex, hierarchical assemblies such as networks (e.g.cytoskeleton, extracellular matrix, mucus) or bundles (e.g.muscle, tendon, ligament), which are crucial for numerous essential processes in the human body.Mainly due to the involved length and time scales and the complex composition of these systems, the design and working principles on the mesoscale often remain poorly understood and their significance regarding human physiology and pathophysiology can only be estimated so far.In this context, computational modeling and simulation is expected to complement theoretical and experimental approaches and significantly contribute to the scientific progress in this field [1][2][3][4][5][6][7][8][9][10][11].Moreover, the application of the accurate, efficient and versatile computational models and methods developed in this context may well be extended towards (future) technologies using, e.g., synthetic polymer, glass or carbon nanofibers in form of meshes, webbings, bundles or fibers embedded in a matrix material [12][13][14][15][16][17][18][19][20].
The targeted class of problems is characterized by involving length and time scales that are inaccessible for simulations with atomistic resolution, but require a level of detail, which precludes homogenized continuum models.Focusing on the efficient and accurate modeling of (short-ranged) molecular interactions such as van der Waals adhesion and steric repulsion, this article tackles one of the major challenges of the simulation-based investigation of these complex biophysical systems: In many cases, these interactions between atoms or charges are the key to the functionality and behavior on the system level.Yet, the inherent high dimensionality of interactions within an ensemble of objects in 3D poses a great challenge in terms of computational efficiency.Our strategy hereby is to rigorously derive models from the first principles of molecular interactions, which are formulated as interaction potentials of atoms or unit charges, and at the same time exploit the dimensionally reduced, slender structure of the fibers, which is inspired by mechanical beam theories, as well as the short-range nature of the considered class of interaction potentials.This leads to accurate, efficient and versatile beam-beam interaction formulations that allow for the computational study of so far intractable problems in complex biophysical systems of slender fibers.
There is a rich body of literature dealing with the analytical and computational modeling of molecular interactions between arbitrarily shaped, solid bodies in 3D space [21][22][23][24][25][26][27].A direct evaluation of the interaction potential between two general bodies in 3D space would require to integrate molecule densities over their volumes, leading to a sixfold integral (two nested 3D integrals) that has to be solved numerically.Even though dimensionally reduced models have been derived that are tailored for short-ranged interactions of such general-shaped 3D bodies and, thus, only require a numerical integration across the interacting surfaces [23], the solution of the remaining fourfold integral in combination with the sharp gradients characterizing these short-ranged interaction forces is still too demanding from a computational point of view to simulate representative 3D systems of curved slender fibers, which necessitates the development of reduced-order models for slender fibers that consistently account for their molecular interactions.While there is a large number of articles [12,13,[28][29][30][31][32][33][34][35][36] focusing on macroscale contact interaction between slender fibers respectively beams, comparable formulations for microscale molecular interactions are still missing.Important steps into this direction have been made by the works [4,37,38], however limited to the interaction of fibers respectively beams with a rigid half-space.
Based on the fundamental kinematic assumption of undeformable fiber cross-sections, as typically applied in mechanical beam theories, the authors recently proposed a generalized formalism to postulate section-section interaction potentials (SSIP) integrated into the framework of Cosserat beam theories, which allows to consider general interaction potentials between curved fibers in 3D space characterized by large deformations, arbitrary mutual orientations, initial curvatures and cross-section shapes as well as inhomogeneous molecule/charge distributions within the cross-sections [39].In a previous contribution of the authors, where this SSIP approach has originally been proposed, exemplary closed-form analytical solutions for the required SSIP laws could be derived for different long-ranged (e.g., electrostatic) and short-ranged (e.g., van der Waals adhesion and steric repulsion) interactions, relying on the additional assumptions of circular cross-section shapes and homogeneous molecule distributions and considering the asymptotic limits of either large distances for long-ranged or small distances for short-ranged interactions [40].Due to the pre-calculated analytical representation of section-section interaction potentials, this SSIP approach only requires the twofold integration along the fiber length directions to be performed numerically.While this approach is general in the sense that tailored SSIP laws could be derived for either long-or short-ranged interactions, the latter class turned out to be critical in terms of accuracy and algorithmic complexity, i.e., the derived SSIP laws for short-ranged interactions did not exhibit a consistent asymptotic scaling behavior and the required numerical solution of a twofold integral was still dominating the overall computational costs when simulating fiber systems.Following up these conclusions from our previous contribution [40], we aim to develop an enhanced approach for the specifically important and challenging case of very short-ranged interactions.
Remarkably, the novel approach proposed in this article turns out to be superior to the previous SSIP approach not only in terms of efficiency, but also regarding the model accuracy.This represents a major step forward enabled by exploiting the short-range nature of the considered class of interaction potentials and the kinematic assumption of undeformable fiber cross-sections.A recently derived, closed-form analytical solution [41] is applied for the interaction potential between a given section of the first fiber (slave beam) and the entire second fiber (master beam), whose geometry is linearly expanded at the point with smallest distance to the given slave beam section.This novel approach based on a pre-defined section-beam interaction potential (SBIP) requires only one single integration step along the slave beam length to be performed numerically.
Based on elementary two-fiber test cases considering pairs of either straight&straight, straight&curved or curved&curved beams, the total beam-beam interaction potential resulting from this approach is shown to exhibit an asymptotically consistent angular and distance scaling behavior in the decisive regime of small separations.Thus, remarkably, the newly proposed SBIP approach turns out to be superior to the previously derived SSIP approach, when applied to short-ranged interactions, not only in terms of computational efficiency but also in terms of model accuracy.Based on conservative estimates for the algorithmic complexity and by considering practically relevant parameter settings, it is demonstrated that the SBIP approach has the potential to reduce the computational costs by at least one order of magnitude as compared to the SSIP approach, and by at least five orders of magnitude as compared to the direct approach of sixfold numerical integration.
To formulate a specific computational model for fiber systems, this interaction model is combined with the geometrically exact beam theory, representing the mechanics of individual fibers.We present all necessary steps of deriving the virtual work contribution, its discretization on basis of the finite element method, and the consistent linearization required for tangent-based solvers for the resulting nonlinear system of equations.A particularly important additional algorithmic aspect is the proposed regularization of the interaction potential in the limit of zero separation, which remedies the inherent singularity of molecular interaction potential laws and only enables a robust numerical solution scheme.
This in turn allows to study a set of numerical examples.An elementary test case to verify the model accuracy of the proposed SBIP approach is the peeling and pull-off behavior of a pair of fibers consisting either of a straight and a curved beam or of two curved beams.Several aspects such as the influence of the strenghth of adhesion on the force-displacement curve and the important difference to the results of the previous SSIP approach with inconsistent asymptotic scaling behavior are studied in detail.Critically, by switching the master-slave assignment, the variant with one straight and one curved beam allows to verify the fundamental modeling assumption of approximating the geometry of one beam (master) as straight cylinder following a first-order expansion.
To demonstrate the readiness of the new approach for real-world use cases, the numerical examples include an application to adhesive nanofiber-grafted surfaces.It experiments with a newly suggested design principle using helical fibers to break the symmetry between strong adhesion in one configuration, where the adhesive fibers are mostly aligned, and seamless pull-off in a 90 • twisted configuration.Beside testing this hypothesis, the simulation results obtained for this numerical example showcase the efficacy, computational efficiency and robustness of the novel computational approach in combination with implicit time integration for large-scale, complex 3D systems.The challenging nature of this setup arises from the large number of fibers with a broad range of different mutual orientations, frequent local pull-off and snap-into-contact events as well as large deformations of the interacting, strongly curved fibers.
Additionally, we propose a new classification of beam interaction models, which puts the newly proposed approach into perspective of previously existing approaches by means of a unifying theoretical framework.Most importantly, this thought concept shows that the new SBIP approach fills a gap in the range from atomistic to macroscopic level of model resolution.The availability of several alternative models for beam interactions in turn motivates a comparison of models for beam contact, where we briefly discuss the advantages and drawbacks of each approach and come to a conclusion of which approach to apply depending on the type of application.
Finally, motivated by an in-depth analysis of the studied numerical examples, we derive the dimensionless key parameters that govern the fundamental behavior of any mechanical system with elastic, adhesive fibers.On this basic theoretical level, we find that there is a dimensionless parameter characterizing the strength of adhesion relative to the structural rigidity of the fibers and we therefore call this parameter "adhesive compliance".This nicely agrees with the numerical results obtained for the studied examples of peeling deformable fibers, where this ratio is indeed identified as the crucial influence determining how much the fibers will deform due to the exposure to adhesive interactions with other fibers.
The remainder of this article is structured as follows.First, fundamentals required for the modeling of molecular interactions and of slender fibers will be briefly summarized in Sect. 2. In the following, the novel SBIP approach will be presented in Sect. 3 and the applied closed-form analytical interaction law as well as the verification of its consistent scaling behavior by means of elementary straight-fiber-pair test cases with analytical solutions will be presented in Sect. 4. These aspects will be combined to derive the resulting virtual work contribution in Sect.5, its discretization based on finite elements as well as the consistent linearization required for tangentbased solution schemes (see Appendix B).Eventually, Sect.6 gives a detailed overview on how the resulting computational models according to the novel SBIP approach and the previously proposed SSIP approach fit into the framework of existing fiber or beam interaction models, either considering microscale molecular or macroscale contact effects.Important numerical and algorithmic aspects, such as the regularization of the interaction potential in the (singular) zero separation limit, search schemes for relevant pairs of interaction partners and overall algorithmic complexity are presented in Sect.7.This discussion paves the way for the set of numerical examples in Sect.8. Finally, this article will be concluded by a summary and outlook in Sect.9.For readers with a special interest in theoretical analysis, we would like to note that Appendix A derives the dimensionless key parameters governing the behavior of the underlying mechanical system of adhesive, elastic fibers.

Fundamentals
This section briefly summarizes essential aspects from the fields of molecular interactions and beam theory that will be referred to in the remainder of this article and are thus crucial for the overall understanding.

Two-body interaction from a molecular perspective
This section briefly recapitulates how the interaction of two extended, macromolecular or macroscopic bodies can be described starting from the first principles of atom-atom, i.e., point-pair interaction.The following summary is reproduced from our previous work [40] for the reader's convenience.Refer e.g. to the textbook [42] for further details.Fig. 1 schematically visualizes the distribution of elementary interaction partners, i.e., atoms or molecules, within two macromolecular or macroscopic bodies.
Consider a point pair interaction potential (r) as a function of the mutual distance r.Popular examples include the inverse-sixth power law valid for van der Waals (vdW) interactions and the LJ interaction law, which extends the attractive vdW part by a repulsive steric contribution: r LJ,eq r 12 − 2 r LJ,eq r 6 . ( Both are typical examples for short-ranged interactions and the generalized form of an inverse power law with high exponent m > 3 will serve as the prime example to be used throughout this article.To ease the comparison with the literature, a few equivalent forms using the most popular definitions and notations of the constant parameters k m , C vdW , r LJ,eq , LJ,eq are stated above.
Assuming additivity, we apply pairwise summation to arrive at the two-body interaction potential Note that the assumption of additivity is known to not hold unconditionally for instance in the important case of vdW interaction.However, the distance dependency obtained from pairwise summation is still valid and only the prefactor called Hamaker "constant" needs to be obtained from advanced Lifshitz theory, which yields Hamaker-Lifshitz hybrid forms [42,43] and motivates us to still apply pairwise summation here.Further assuming a continuous atomic density ρ i , i = 1, 2, the total interaction potential can alternatively be rewritten as nested integrals over the volumes V 1 , V 2 of both bodies B 1 and B 2 : It can be shown that this continuum approach is the result of coarse-graining, i.e., smearing out the discrete positions of atoms in a system into a smooth atomic density function ρ(x). [22]

General strategy to account for molecular interactions
The general approach to incorporate the effect of molecular interactions is identical to the one suggested for solid bodies in previous work [21,22] and has been summarized also in our recent contribution [40], which is repeated here for convenience.For a classical conservative system, the total potential energy of the system can be stated taking into account the internal and external energy int and ext .The additional contribution from molecular interaction potentials ia is simply added to the total potential energy as follows.
Note that the standard parts int and ext remain unchanged from the additional contribution.One noteworthy difference is that internal and external energy are summed over all individual bodies in the system whereas the total interaction free energy is summed over all pairs of interacting bodies.In order to shed some light on the basic characteristics of systems with adhesive elastic fibers, the dimensionless parameters of the governing Eq. ( 6) are identified by means of nondimensionalization and discussed in Appendix A.
The weak form of the equilibrium equations, which serves as basis for a subsequent finite element discretization of the problem, can now be derived by applying the principle of minimum of total potential energy and reads Alternatively, the same result is obtained by applying the more general principle of virtual work, which also holds for non-conservative systems.Clearly, the evaluation of the interaction potential ia , or rather its variation δ ia , is the crucial step here.Recall Eq. ( 5) to realize that it generally requires the evaluation of two nested 3D integrals 1 .The direct approach of incorporating this interaction potential in a computational model using 6D numerical quadrature turns out to be extremely costly and in fact inhibits any application to (biologically) relevant multi-body systems.See Sect.7.3 for more details on the algorithmic complexity and the computational cost of this naive, direct approach as well as the novel SBIP approach to be proposed in Sect.3.

Van der Waals interaction potential between two straight cylinders
One of the main objectives of this work is that the mutual orientation of two fibers, i.e. the angle α enclosed by the local tangent vectors of the (potentially curved) fiber centerlines, shall be taken into account in order to improve the accuracy of the reduced interaction law and the overall computational model.In this section, we review the angle dependency for the simple case of two straight fibers, which will serve as an analytical reference solution in order to verify the specific reduced interaction law to be derived in Sect. 4. To begin with, recall the analytical solutions for the special cases of parallel and perpendicular cylinders of infinite length (in the regime of small separations).These expressions agree with the following, more general relationship valid for all mutual angles α ∈ ]0, π/2] stated e.g. in the textbook [43, p. 173]: Here, R 1 and R 2 denote the radius of the first and second cylinder, and g −1 bl denotes their (bilateral) smallest surface-to-surface separation, also known as gap.For the limiting case of perpendicular cylinders α = π/2, this coincides with vdW,cyl⊥cyl,ss listed in the quick reference table of analytical solutions in our previous article [40] (alongside other expressions mentioned here).For the case of parallel cylinders, however, note that the total interaction potential of infinitely long cylinders is infinite and we obtain πvdW,cyl cyl,ss instead.Remark.Interestingly, the 1/ sin α-scaling behavior also holds true for screened electrostatic interaction of two cylinders [44], [45, p. 23].Indeed, it is shown in [42, p. 218], that this relation results from fundamental geometric considerations related to the socalled Derjaguin approximation, and is thus independent of the type of interaction, i.e., the specific form of the point interaction potential law (r).
In addition to these analytical expressions obtained by means of simple pairwise summation, further theoretical work that relax certain assumptions and consider more advanced aspects like interaction across inhomogeneous or anisotropic media, differences in optical material properties or retardation can be found in the literature.To give but one example, [46] studies cylinders with anisotropic optical properties considering the example of carbon nanotubes.A review of recent research activities on this topic is given in [47].This work however focuses on the extension towards curved slender fibers with arbitrary mutual separations/orientations due to their possibly large elastic deformations in 3D, and therefore uses the basic pairwise summation approach as mentioned and motivated in Sect.2.1.

Steric repulsion-mechanical contact
The prevailing notion of contact between bodies in biophysics is commonly described as excluded volume effect, which means that bodies may approach each other without any influence on each other and only as soon as their surfaces touch, the repulsive contact forces that inhibit any overlap of the bodies' volumes may rise to infinite strength.
Throughout this work, repulsive contact forces will be modeled based on the repulsive part of the LJ potential law (Eq.2), which is an inverse-twelve power law in the separation of the point-like atoms.A number of alternative force-distance laws can be found in literature, however this approach seems to be most consistent with the modeling of vdW interactions as inverse-six point-pair potential.In particular, in Sect.6, the approach of using a repulsive interaction potential will be compared to macroscopic formulations well-established for the mechanical contact interaction of slender fibers respectively beams.

Geometrically exact 3D beam theory and corresponding finite element formulations
We use geometrically exact 3D beam theory to model large elastic deformations of slender fibers in 3D space.See e.g.Ref. [48] for a recent review of both space-continuous beam theories as well as suitable temporal and spatial discretization schemes.The interaction Fig. 2 Example configurations and kinematics of the Cosserat continuum formulation of a beam to illustrate the field of centroid positions r(s, t) and material triad field (s, t): Initial, i.e., stress-free (blue) and deformed (black) configuration.Straight configuration in the initial state is chosen here as an example, although the applied beam theory is more general approach to be proposed in this article is both independent of the specific beam formulation and the discretization schemes used to describe the mechanics of individual fibers.In this work the proposed interaction approach is exemplarily applied in combination with geometrically exact beam elements of both shear-deformable (Simo-Reissner) and shear-free (Kirchhoff-Love) type, which have been previously described in Ref. [16] and the aforementioned review article [48].The following two subsections briefly summarize the fundamentals required for the remainder of this article.

Space-continuous beam theory
The configuration of a beam at time t is uniquely defined by the controid position vector r(s, t) ∈ R 3 and the orthonormal frame (s, t) = [g 1 , g 2 , g 3 ] ∈ SO(3) describing the cross-section orientation at each point s along the 1D Cosserat continuum.Note that the arc-length parameter s is hereby defined in the stress-free, initial configuration of the centerline curve r 0 (s) := r(s, t = 0).See Fig. 2 for an illustration of these geometrical quantities and the resulting kinematics.
According to this concept of geometry representation, the position x of an arbitrary material point P of the slender body is obtained from Here, the additional convective coordinates s 2 and s 3 specify the location of P within the cross-section, i.e., as linear combination of the orthonormal directors g 2 and g 3 .For a minimal parameterization of the triad, e.g. the three-component rotation pseudo-vector ψ may be used, i.e. (s, t) = (ψ(s, t)), such that we end up with six independent degrees of freedom (r, ψ) at every centerline location s to define the position of each material point in the body by means of Eq. ( 9).Again refer to Fig. 2 for a sketch of the kinematics of geometrically exact beam models.Based on these kinematic quantities, deformation measures as well as constitutive laws can be defined.Finally, the potential energy of the internal (elastic) forces and moments int is expressed uniquely by means of the set of six degrees of freedom (r, ψ) at each point of the 1D Cosserat continuum.See e.g.[48][49][50] for a detailed presentation of these steps.Remark on notation.Unless otherwise specified, all vector and matrix quantities are expressed in the global Cartesian basis E i .Differing bases as e.g. the material frame are indicated by a subscript [.] g i .Quantities evaluated at time t = 0, i.e., the initial stress-free configuration, are indicated by a subscript 0 as e.g. in r 0 (s).Differentiation with respect to the arc-length coordinate s is indicated by a prime, e.g., for the centerline tangent vector r (s, t) = d r(s, t)/d s .Differentiation with respect to time t is indicated by a dot, e.g., for the centerline velocity vector ṙ(s, t) = d r(s, t)/d t .For the sake of brevity, the arguments s, t will often be omitted in the following.Remark on finite 3D rotations.To a large extent, the challenges and complexity in the theoretical as well as numerical treatment of the geometrically exact beam theory can be traced back to the presence of large rotations.In contrast to the much more common vector spaces, the rotation group SO( 3) is a Lie group with associated Lie algebra so(3) consisting of all skew-symmetric 3D tensors.The fact that this Lie group is non-commutative and non-additive renders standard procedures such as the interpolation or the update of configurations quite intricate (see e.g.Ref. [51] for details).We thus follow a two-part strategy.First, we aim to develop and formulate the novel approach in the most general form in Sect.3, including also cases such as arbitrary cross-section shapes or inhomogeneous atomic densities, where the involvement of large 3D rotations is inevitable.In a second step, however, we aim to abstain from the handling of finite 3D rotations wherever possible when proposing specific reduced interaction laws for instance for the case of homogeneous, circular cross-sections considered in Sect. 4. As a result, this will allow to avoid the handling of finite rotations in the interaction potentials and to achieve simpler and more compact numerical formulations whenever possible.

Spatial discretization based on beam finite elements
A smooth enough, i.e., at least C 1 -continuous, discrete representation of the centerline curve is inevitable in the context of molecular interaction laws with its typical high gradients in order to ensure robustness of the numerical method also for reasonably coarse discretizations.This is a well-known general issue and has been discussed in the context of macroscopic beam contact interaction [16] and (surface enrichment of) 2D and 3D solid contact elements based on the LJ interaction potential [52] before.In the scope of this work, it is addressed by applying a third order Hermite interpolation scheme to discretize the centerline curve r(s).The corresponding beam finite element formulations of both Simo-Reissner type and Kirchhoff-Love type have been presented in [16] and [48], respectively, and only the most essential aspects required later in this work will be briefly summarized in the following.
The spatial centerline curve r is approximated by means of the discrete set of the i-th node's position vector di ∈ R 3 and tangent vector ti ∈ R 3 (i = {1, 2}) as primary degrees of freedom and the four scalar cubic Hermite polynomials H i d/t used for interpolation as follows.
Here, l ele denotes the initial length of the element.The newly introduced elementlocal parameter ξ ∈ [−1; 1] is biuniquely related to the arc-length parameter s ∈ [s ele,min ; s ele,max ] describing the very same physical domain of the beam and the scalar factor defining this mapping between both length measures in differential form is called the element Jacobian J (ξ ) with ds =: J (ξ ) dξ .Note that on the right hand side of Eq. ( 10), all the centerline degrees of freedom of one beam element, i.e., the nodal positions di and tangents ti are collected in one vector d for a more compact notation.Accordingly, H is the assembled matrix of shape functions, i.e., Hermite polynomials H i d and H i t .Following a Bubnov-Galerkin scheme, this very same interpolation is applied to the test functions, i.e., the variation of the centerline curve As mentioned already in the last section, the specific reduced interaction law to be proposed in Sect. 4 will be described by the centerline curve (and tangent field) only and avoid the use of the rotation field in favor of a simple and efficient formulation.Thus, the spatial discretization of the rotation field will not be required and omitted here for the sake of brevity.
Remark on notation.Note that a derivative with respect to the element-local parameter ξ will be indicated by a short upright prime, e.g., r (ξ , t) = d r(ξ , t)/d ξ in order to differentiate it from the derivative r (s, t) with respect to s.
At the end of this section, we would like to point out that both the general SBIP approach and the reduced interaction law to be proposed in this article are not limited to a specific interpolation scheme and that the presented Hermite interpolation is just one example that we employ throughout this work.

The section-beam interaction potential (SBIP) approach
When considering two slender, fiber-like bodies with lengths l i and cross-sections A i (i = 1, 2), it is reasonable to split the two-fold volume integral of the interaction potential Eq. ( 5) in integrals across the fiber lengths and cross-sections as follows: In our previous work [40], the double length-specific interaction potential π(r 1−2 , ψ 1−2 ) between two cross-sections characterized by a distance vector r 1−2 and mutual orientation vector ψ 1−2 [39] has been approximated analytically by exploiting the shortrange nature of the considered class of interaction potentials and the fundamental kinematic assumption Eq. ( 9) characterizing the fiber deformation.Thus, in the proposed computational model only the two integrals across the fiber lengths l 1 and l 2 Fig. 3 Illustration of the section-beam interaction potential (SBIP) approach.The actual, deformed volume of the second interaction partner ("master") is approximated by a surrogate body (blue) located at the closest point ξ 2c to a given integration point ξ 1,GP of the first interaction partner ("slave").Distance vector r 1−2c and relative rotation vector ψ 1−2c uniquely describe the mutual separation and orientation had to be solved numerically, a procedure that was denoted as section-section interaction potential (SSIP) approach referring to the physical meaning of π(r 1−2 , ψ 1−2 ).The present works aims to follow this path one (essential) step further, by approximating the second of the interacting fibers/beams as a (cylinder-shaped) surrogate body constructed at the position of smallest distance with respect to a given point on the first fiber such that the single length-specific interaction potential π (r 1−2c , ψ 1−2c ) between a cross-section of the first fiber and the entire second fiber can be approximated analytically.This procedure will be denoted as section-beam interaction potential (SBIP) approach, again referring to the physical meaning of π(r 1−2c , ψ 1−2c ), and will only require a 1D integral, i.e. integration across the first fiber's length l 1 , to be solved numerically.This general approach will be motivated in the following, before a specific expression for the section-beam interaction potential π (r 1−2c , ψ 1−2c ) will be presented in Sect. 4. Consider a point-pair interaction potential with a very steep gradient as for example the inverse power laws with exponent six or twelve from the popular LJ interaction law (Eq.2).The rapid decay of the potential with increasing distance implies that among all possible point pairs between both bodies only those with the smallest separation contribute significantly to the total interaction potential of both bodies.When looking at the interaction of two deformable slender bodies such as fibers, this consideration gives rise to an approach where the geometry of the second body is approximated by a surrogate body with simplified geometry located at the point of closest distance from a given point on the first body.See Fig. 3 for an illustration of the approach using the example of circular cross-sections and therefore a cylinder-shaped surrogate body.
In the region around the closest point, this straight cylinder is expected to be a good approximation for the actual, possibly deformed, beam geometry.Note, however, that the general SBIP approach to short-ranged beam-beam interactions is not limited to the circular cross-section shape shown in this example.
In accordance with formulations for macroscopic beam contact, the body which is projected onto, i.e., here the one with approximated geometry is referred to as master beam (indicated by subscript m) whereas the first body is called slave beam (indicated by subscript s).Without loss of generality, the beam with index 1 is assumed to be the slave beam whereas index 2 is used as a synonym for master.
From a mathematical point of view, the geometrical approximation used in this context is equivalent to a Taylor series expansion of the centerline curve r 2 (ξ 2c ) of the master beam at the closest point ξ 2c truncated after the second, i.e., linear term: Here, the linear term represents the orientation of the surrogate body in the direction of the master beam's tangent vector at the closest point r 2 (ξ 2c ).Recall from the previous section that the short prime denotes a differentiation with respect to the element parameter coordinate, i.e., r i ( As stated above, we assume an interaction potential π(r 1−2c , ψ 1−2c ) for the interaction between all the points within one cross-section of the slave beam and the entire master beam surrogate, i.e., the tangential straight cylinder in the example above.Following Eq. ( 12), the total interaction potential is evaluated as an integral along the slave beam centerline curve as follows: Generally, such a section-beam interaction potential (SBIP) is a length-specific quantity with dimensions of energy per unit length (of the slave beam).It is an analytical expression uniquely defined by the mutual configuration, i.e., the distance vector r 1−2c and relative rotation vector ψ 1−2c as illustrated in Fig. 3.
To evaluate the remaining 1D integral in Eq. ( 14), Gaussian quadrature is applied throughout this work.The reduction starting from 6D integration (cf.Eq. ( 5)) to 1D integration already indicates the overall gain in efficiency that will be further analyzed and discussed in Sec.7.3.Recall also the 2D integral (two nested 1D integrals) resulting for the general SSIP approach to realize the superior efficiency of this novel SBIP approach that will also be verified in the numerical experiments of Sec. 8.
In analogy to the previously presented section-section interaction potential (SSIP) approach from Ref. [40], the question of how to find an analytical, closed-form expression for the reduced interaction law π(r 1−2c , ψ 1−2c ) can be considered separately from the generally valid SBIP approach proposed in this section.Just as for the SSIP laws, such an effective SBIP law π will depend on the considered type of interaction, the cross-section shape(s) and dimensions, the atom density distributions, and possibly other interaction-specific factors.An example for how to determine this (single) lengthspecific potential π(r 1−2c , ψ 1−2c ) analytically by means of 5D integration starting from the point-point interaction potential law (Eq.2) in case of LJ interaction is presented in Sect. 4.However, other strategies such as postulating the general form of the SBIP law and fitting of the parameters to experimental results e.g. for the force response in two-fiber systems are considered to be promising alternative ways that could enable a broad variety of future applications of this SBIP approach.
Discussion of the choice of master and slave side.Starting from the problem of twobody interaction that is symmetric with respect to the two interaction partners, the SBIP approach introduces the notion of master and slave, which causes a bias in the formulation and asks for a criterion how to assign these roles.Both from the mathematical description as a linear Taylor approximation and the illustration in Fig. 3, it becomes clear that the resulting model error and bias will depend on the magnitude of curvature of the master beam's centerline.This would give rise to a criterion that chooses the beam with the smaller (maximum or average) curvature as the master beam.However, such criteria might lead to sudden changes of master and slave over time, which is numerically unfavorable.Alternatively, one could consider to evaluate the pair of beams in two half passes, where the roles of master and slave switch and only the contributions on the slave side are evaluated in each of the passes (see e.g.[24,53]), such that the bias in the formulation is avoided.However, it can be argued that the model error introduced by the use of the surrogate beam on the master side is negligible, because first, the curvature anyway is assumed to be limited in the underlying beam theory (typically compared to the inverse radius 1/R in case of circular cross-sections) and second, the very short range of the considered interactions naturally limits the impact of the master beam's shape deviation from the surrogate shape, because only the immediate surrounding of the expansion point will contribute noticeably to the total interaction.Following this assumption that the corresponding model error will be negligible, we apply the simple heuristic that the beam (element) with the smaller (global) identification number (ID) will generally be the slave beam throughout this work and validate this assumption in the numerical example of Sect.8.1.The resulting maximal relative difference in the force response on system level turns out to be below 1.5% even for relatively large curvatures, which is considered to be a reasonably small model error.In addition, this simple criterion based on element IDs ensures a unique decision that does not change in the course of the simulation.Remark on self-interactions.As already discussed for the SSIP approach, self-interactions, i.e., the interaction of distinct parts of the same beam, can be treated naturally also within the SBIP approach.Leaving everything else unchanged, the search for and evaluation of (non-neighboring) beam element pairs from one and the same physical beam directly allows to incorporate the effect of self-interaction.This is considered to be important for long, flexible fibers showing the tendency to large deformations.

Closed-form expression for the disk-cylinder interaction potential
There are different ways to arrive at a closed-form expression for the required SBIP law π .One of them is the analytical integration of a point-pair potential over all point pairs in the section-beam (surrogate) system.This strategy has been applied in our recent contribution [41] considering a generic inverse power law m (r) = k m r −m with exponent m ≥ 6. Due to the generality, the resulting reduced interaction law π can be used to model both the adhesive vdW part (m = 6) and the repulsive part (m = 12) of the LJ potential.Moreover, it considers the practically relevant case of circular, undeformable cross-sections and homogeneous densities of the fundamental interacting points in both fibers.From a geometrical point of view, this leads to the scenario of a disk and a cylinder (cf.Fig. 3) with arbitrary mutual configuration, i.e., separation and orientation.It is of great importance for the accuracy of the SBIP approach that the applied disk-cylinder potential law is accurate for all mutual configurations.
Note that instead of the six degrees of freedom (r 1−2c , ψ 1−2c ) in the most general scenario [39] considered in the preceding section, this system consisting of a homogeneous circular disk and cylinder (radii R 1 and R 2 ) can be described by only three degrees of freedom.Out of different analytical approximations for (r 1−2c , ψ 1−2c ) as derived in Ref. [41], one variant has been identified as particularly appealing and will also be considered throughout the present work.This variant allows for sufficient model accuracy and at the same time for a pleasantly simple and compact closed-form representation of (r 1−2c , ψ 1−2c ) that only requires two degrees of freedom to describe the fiber interaction: the angle α := arccos (r 1 (ξ 1,GP ) • r 2 (ξ 2c ))/(||r 1 (ξ 1,GP )|| ||r 2 (ξ 2c )||) enclosed by the tangent vectors r 1 (ξ 1,GP ) and r 2 (ξ 2c ) as well as the surface gap function )||, between the position vectors r 1 (ξ 1,GP ) and r 2 (ξ 2c ) defined at a given location ξ 1,GP on the centerline of the slave beam and the associated closest projection point ξ 2c on the centerline of the (surrogate) master beam.
The derivation of this closed-form expression for the disk-cylinder potential law πm,disk-cyl requires a 5D analytical integration of the point-pair potential m , which is a challenging theoretical problem for itself.Here, we only briefly repeat the problem statement and the final expression found in [41] for the reader's convenience.

Problem statement
In the derivation, first an interaction potential m,pt-cyl between a single point on the disk (i.e., the slave beam cross-section) and the cylinder (i.e., the surrogate master beam) is derived before solving a double integral across the disk area.πm,disk-cyl := Here, x 1 ∈ A disk denotes any point on the disk.On the master side, x 2 ∈ V cyl denotes any point on the surrogate body assumed as an infinitely long auxiliary cylinder oriented along the (normalized) tangent vector t 2 = r 2 / r 2 .

General strategy
The general strategy follows the one generally known as point-pairwise summation (see e.g.[42,43] for details and a discussion) and e.g.applied in [54] for the analytical calculation of vdW forces for certain geometric configurations, e.g., a cylinder and a perpendicular disk.Since already for such specific scenarios, no exact analytical solution can be found for the integrals, also the proposed derivation made use of the common approach of series expansions in order to find an analytical, closed-form expression for the integral of the leading term(s) of the series.Due to the rapid decay of the inverse power laws, this is known to yield good approximations for the true solution of the integral.

Approximative solution
The final form of the disk-cylinder interaction potential to be used as reduced interaction law in the context of this article reads: For convenience in later reference, we explicitly state the most common prefactors for the vdW part m = 6 and the repulsive part m = 12 of the LJ potential as follows:

Verification
An immediate verification of these expressions for the special case α = 0 confirms that both π6,disk-cyl and π12,disk-cyl are identical to the independently derived analytical solutions for the interaction potential per unit length π6,cyl cyl and π12,cyl cyl of two infinitely long, parallel cylinders (cf.Eqs.(A23) and (A24) in our previous contribution [40]).This is an important finding, as it shows the consistency of the more general expression (Eq.17) valid for all mutual angles α with previously derived expressions for the important special case α = 0.A comprehensive analysis of the accuracy of Eq. ( 17) as well as a comparison to alternative, increasingly complex expressions can be found in Ref. [41].Here, we only briefly summarize the conclusions and finally present the important comparison to the accuracy of the reduced disk-disk interaction potential laws π used together with the SSIP approach in Ref. [40].Most importantly, the pleasantly simple expression from Eq. (17) shows the correct asymptotic scaling behavior in the decisive regime of small separations and small angles.It is thus considered the optimal compromise between accuracy and complexity of the expression for the purposes of this work.In particular, also the theoretically predicted 1/ sin α angle dependence (for non-parallel cylinders α = 0) has been successfully verified.
Figure 4 demonstrates that the novel SBIP law πm,disk-cyl from Eq. ( 17) in combination with the general SBIP approach from Sect. 3 is significantly more accurate than the previous method using the SSIP law πm,disk disk in combination with the general SSIP approach as proposed in our previous contribution [40].
Recall the important result of the SSIP verification therein that the simple and readily available section-section interaction law πm,disk disk from Eq. ( 18) in Ref. [40] does in general not yield the correct asymptotic scaling behavior in the limit of small separations, which is decisive in case of short-ranged interactions, and that the orientation of the (disk-shaped) cross-sections would need to be included in the reduced interaction law π to improve this crucial characteristic.According to Fig. 4, the alternative SBIP approach specialized on short-range interactions in combination with the SBIP law πm,disk-cyl derived in [41] has finally accomplished the goal of reproducing the correct asymptotic scaling behavior as demonstrated for the special cases of two parallel cylinders (expected asymptotic scaling of order 1.5; see Fig. 4a) and two perpendicular cylinders (expected asymptotic scaling of order 1; see Fig. 4b).
The plots in Fig. 5 complement the analysis above by showing the dimensionless interaction potential as a function of the mutual angle α.  18) in Ref. [40] (used together with the SSIP approach proposed in the same article; dark green line with crosses) with the analytical expression for the disk-cylinder potential π6,disk-cyl Eq. ( 17) (used together with the SBIP approach from Sect.3; red line with triangles).The numerical reference solution obtained via 3D Gaussian quadrature of the point-half space potential 6,pt-hs from Ref. [41] (green line with diamonds) and the analytical reference solutions summarized in Sect.for the disk-cylinder potential π6,disk-cyl from Eq. ( 17) (used together with the SBIP approach from Sect.3; red line with triangles) by means of a numerical reference solution obtained via 3D Gaussian quadrature of the point-half space potential 6,pt-hs from Ref. [41] (green line with diamonds) and by means of the analytical reference solution summarized in Sect.2.3 (black dashed line) The considered scenario of two cylinders and the three different solutions for the twocylinder interaction potential are identical to the previous Fig. 4. Most importantly, the theoretically predicted scaling behavior is confirmed by the numerical reference solution and reproduced by the disk-cylinder potential law π6,disk-cyl (option C) from Eq. ( 17).

Virtual work contribution
Recall from Eq. ( 7) that it is the variation of the two-body interaction energy δ ia , which needs to be evaluated to incorporate the effects of molecular interactions into both the theoretical and computational framework of nonlinear continuum mechanics.According to the general SBIP approach proposed in Sect. 3 combined with the generic disk-cylinder interaction potential law πm,disk-cyl from Sect. 4, the variation of the two-fiber interaction potential reads Note that δ(ds 1 ) vanishes due to the fact that ds 1 = r 01 (ξ 1 ) dξ 1 only depends on the element parameter coordinate ξ 1 and the initial ("0") configuration of the slave beam, but not on the current configuration, i.e., current values of the primary degrees of freedom.
For the sake of brevity, the subscripts "m" and "disk-cyl" as well as the function arguments (g ul , α) of π will be omitted throughout this section.The variation of π consists of the summands, each defined by two factors which will be determined subsequently in the next steps.First, the derivatives of π can be expressed in recursive manner: Note that the remaining two factors δg ul and δ(cos α) are known from the macroscopic line contact formulation proposed in [34] and its combination with a point contact formulation presented in a unified ABC beam contact formulation [35], respectively.Both are reproduced here for the sake of completeness and a unified notation: In the previous equations, we have introduced the unit tangent vectors t i := r i / r i , i = 1, 2, the unilateral unit normal vector n ul := d ul /d ul with unilateral distance vector d ul := r 1 (ξ 1 ) − r 2 (ξ 2c ) as well as the auxiliary vectors Note the difference between the notations δ r 2 (ξ 2c ) and δr 2 (ξ 2c ) (see Eq. 27), which originates from the fact that ξ 2c is the result of a (closest) point-to-curve projection, i.e., it depends on the primary variables of our problem.Thus, δ r 2 (ξ 2c ) represents a total variation, and δr 2 (ξ 2c ) a partial variation at fixed ξ 2c .In contrast to δg ul in Eq. ( 23) 2 , the additional contribution from δξ 2c must actually be computed and included to ensure a variationally consistent formulation in Eq. (24).Also for the later reference, all the expressions required in this respect are given here as using the following expression for the variation of the slave beam parameter coordinate where the derivative of the scalar orthogonality condition At this point, we have gathered all the required pieces that allow us to evaluate the virtual work contribution from molecular interactions δ ia according to Eq. ( 19).As discussed along with the general SBIP approach in Sect.3, the 1D integral along the slave beam is evaluated numerically, e.g., by means of Gaussian quadrature.Note that the correctness of the presented and implemented expressions of this section has been verified to be consistent with the corresponding interaction energy ia (see Eq. ( 14) with Eq. ( 17)) by means of an automatic differentiation tool [55].
In a next step, the contribution δ ia of the interaction potential to the weak form Eq. ( 7) needs to be discretized in space.The discrete counterpart δ ia,h of the space-continuous form δ ia is obtained via substitution of the centerline interpolation scheme from Eqs. Eq. (10) and Eq. ( 11) into Eq.(20).This allows to identify the discrete residual vectors r ia,i of the interacting beam elements i = {1, 2} that finally result from the molecular interactions.For quasi-static problems, this step is sufficient to transfer our mechanical problem into a discrete set of nonlinear algebraic equations that need to be solved numerically for the discrete (nodal) primary variables d.If tangent-based solution schemes, e.g.Newton-Raphson, shall be applied for this purpose, the required linearization of these residual vectors [δ ia ] with respect to the vector of primary degrees of freedom d is provided in Appendix B. Discussion of the special treatment required for master beam endpoints Recall from Sect. 3 and Sect. 4 that the cylinder used as the surrogate for the master beam has been assumed to have an infinite length.Due to the very short range of the interactions considered here, this is an excellent approximation in almost all cases.In the special case that the result of the closest-point projection is a master beam endpoint, however, this approximation overestimates the true contribution to the interaction potential approximately by a factor of two.Again, given the short range of interactions considered here, the resulting model error can be interpreted as if the master beam was slightly longer 3  than it actually should be.
Due to this short range of interactions and the rarity of this event involving the endpoints of slender fibers among all those cases involving the points between the endpoints, the influence of this model error on the total two-body interaction potential is expected to be negligible in almost all applications.Nevertheless, we can think of the worst case scenario, where two straight, parallel, adhesive fibers of finite length (with equilibrium inter-axis separation) slide along each other in axial direction and the only effective force would be the one at the fiber endpoints, where the influence of the second beam on an exemplarily considered cross-section of the first beam rapidly decreases, because the second beam ends.Whereas the unmodified SBIP approach using πm,disk-cyl would yield zero force, it could be augmented by a special treatment of master beam endpoints that subtracts half of the interaction potential contribution at any integration point where the result of the closest-point projection is a master beam endpoint.This procedure probably needs to be smooth such that the transition from the full disk-cylinder interaction potential contribution to half that value at the master beam endpoint needs to be smeared out over a small, yet finite length of the beam.Due to the expected negligible effect in almost all applications, this augmentation of the SBIP approach is left for future work, but this discussion as well as the described worst-case scenario should prove useful when implementing, calibrating and verifying this model enhancement.

Beam interaction formulations from a meta-level perspective
This section aims to take a step back and look at beam interaction formulations from a meta-level perspective in order to get an overview of the previously existing approaches and the new one proposed in this article.

Classification and comparison of approaches for beam-beam interactions
A classification and comparison of beam-beam interaction formulations is provided in Fig. 6.
The leftmost column depicts the approach of point-pairwise summation -or corresponding nested volume integration-of the fundamental interaction potentials for pairs of molecules or charges.The second and third column show the formulations based on section-section interaction potentials (SSIP) [40], and based on section-beam interaction potentials (SBIP) proposed in Sect.3, respectively.Finally, the rightmost column depicts a possible further dimensional reduction of the problem to the interaction of two beam surrogates, which would allow to evaluate the two-body interaction potential by a single function evaluation. 4From left to right, the resolution of details decreases and likewise the algorithmic complexity determined by the dimensionality of the underlying problem decreases.This overview of four distinct, logical categories of beam-beam interaction formulations therefore illustrates the tradeoff between resolution of details ranging from atomistic view to trivial meso/macroscopic bodies on the one hand and the dimensionality and thus main driver for the computational cost on the other hand.The ultimate goal for the derivation and choice of (a class of) formulations however is to outsmart this natural conflict of objectives by a consistent dimensional reduction of the fully resolved problem (on the left) to the minimal possible description that is yet able to reproduce the most important, characteristic features.Based on the detailed theoretical considerations in Ref. [40] and Sect.3, this optimal choice is given as the SSIP approach for long-range and the SBIP approach for short-range molecular interactions of beams.
This new overall picture of beam interaction formulations also allows to classify previous approaches to macroscopic contact of beams and interpret them in the context of molecular interactions, which are also the origin of the macroscopic contact forces and resulting non-penetrability of objects that we observe in everyday life.Interestingly, the very first numerical formulation of (macroscopic) beam contact is based on the concept of determining the one bilateral closest-point pair between both beams and evaluating an heuristic penalty force law as a function of the closest point-pair separation in order to preclude any (noticeable) penetrations [28].Given this new overall picture of beam interaction formulations, such an approach can be interpreted as the consistent, most extreme dimensional reduction of the problem motivated by the short range of interactions and the resulting possibility to evaluate the total interaction potential for a pair of surrogate bodies approximating the shape of the actual beams.However, this new perspective likewise reveals the well-known limitations of this kind of approach with respect to describing arbitrary mutual configurations as the illustrative examples of two parallel straight beams or one straight beam and a surrounding helical beam demonstrate (see e.g. the Discussion in Ref. [34]).The non-uniqueness of the bilateral closest-point pair in such situations is a confirmation of the oversimplification of the general beam-beam interaction problem by such an approach.Nevertheless, this category of surrogate-surrogate interaction formulations is the most efficient theoretically possible class of formulations and due to its validity for a certain range of mutual orientations, this efficiency can be exploited in combined approaches such as the ABC formulation [35] that handle the problematic mutual configurations differently.This recognition of the superior efficiency of the existing, combined macroscopic contact formulation asks for a more detailed discussion of the applicability of such an heuristic approach to preclude penetration also in micro-and nanoscale problem settings, which will thus be given in the following section.

Brief comparison of micro-and macroscopic approaches to beam contact
Modeling the steric repulsion that prevents a penetration of distinct fibers has a long history in the field of computational contact mechanics and has first been addressed in Ref. [28].The paradigm of these macroscopic continuum models is that the smallest surface separation or gap must be equal to or greater than zero which is formulated as an inequality constraint.With the development of the novel SSIP and SBIP approaches to molecular interactions of fibers, an alternative modeling strategy has arisen, which is motivated by the rather microscale perspective of LJ interactions between all material points in the slender continua.This asks for a brief assessment and comparison of the modeling approaches.
On the one hand, penalty-based models for beam contact are well-established formulations with optimized efficiency as well as robustness.On the other hand, the SSIP and SBIP approaches are based on first principles in form of the LJ law and are thus expected to better resolve the actual contact force distributions, especially in the case of nano-to micro-scale applications.This becomes obvious if we recall the purely heuristic nature of the penalty force law and the resulting (small) negative gap values, i.e., tolerated penetrations.It will most likely depend on the specific application whether the associated model error has a significant or rather negligible influence on the results.In order to answer this question in the context of the authors' recent work e.g. on biopolymer filament mechanics, it seems most important to look at the adhesive force laws to be applied in combination with the models for steric repulsion.The SSIP law modeling long-ranged electrostatic attraction [40,Eq. 35] is an inverse power law in the inter-axis separation d rather than the smallest surface separation g = d − R 1 − R 2 and thus expected not to be very sensitive to small changes in the gap values in case of contacting fibers g ≈ 0. On the contrary, the SSIP law [40,Eq. 31] as well as the SBIP law Eq.( 17) modeling the short-range vdW adhesion are inverse power laws in the gap itself and therefore highly sensitive with respect to the gap g.Indeed, the thorough validation of both adhesion models using the example of two straight slender fibers in Ref. [56] as well as an unsuccessful attempt to use penalty beam contact in combination with vdW adhesion to model the peeling process of two slender fibers5 confirm these considerations.Moreover, refer to Sect.7.1 and Ref. [40] for a detailed discussion of the importance to correctly resolve the regime of small gap values by means of a suitable regularization strategy to remedy the inherent singularity of the (reduced) vdW interaction laws for zero separation g → 0.
To conclude, the choice of a proper computational model for steric repulsion between contacting fibers is closely linked to the type of adhesion and will most likely also depend on the specific application.For the reasons outlined above, the authors decided to apply the penalty-based line contact formulation together with the rather long-ranged SSIP law for electrostatic attraction, whereas the SSIP/SBIP approach based on the repulsive part of the LJ law will be applied in combination with the short-ranged SSIP/SBIP law for vdW adhesion, respectively.Nevertheless, a more detailed analysis of the similarities and differences of existing, macroscopic beam contact formulations and the novel approaches based on molecular steric repulsion is considered an interesting aspect of future work in this field.

Regularization and selected further aspects
This section discusses the numerical regularization scheme as well as further (algorithmic) aspects that are of special importance for the application of the novel SBIP approach and the proposed interaction law π.

Regularization of the reduced disk-cylinder interaction law in the limit of zero separation
Due to the inherent singularity of molecular interaction potentials in the limit of zero separation, a numerical regularization is required in order to solve the governing, nonlinear system of equations resulting from Eq. ( 7) in a robust manner.Generally, such a numerical regularization is a standard approach in (beam) contact formulations (see e.g.[35,57]) and in the specific context of the LJ potential considered here, it has first been applied in [52].In analogy to the regularization of the section-section interaction potential law in our previous contribution [40], a quadratic/linear extrapolation of the section-beam interaction potential/force law will be applied here in the range of very small gap values g ul < g ul,reg below a certain regularization parameter g ul,reg ∈ R + .All the additionally required expressions to compute the regularized section-beam interaction potential law πreg are provided in Appendix C.
If this regularization parameter is chosen sufficiently small, which means smaller than any gap value occurring in any converged configuration of any time/load step throughout the entire simulation, such a regularization scheme will not influence the results at all and can thus be considered a mere auxiliary means to enable and improve the iterative process of solving the nonlinear system of equations.Altogether, the necessity of a suitable regularization scheme due to its significantly positive effect on robustness and efficiency will be confirmed also by the numerical examples to be presented in Sect.8.

Objectivity and conservation properties
Recall from Sect. 5 that the final contributions to the discrete element residual vectors r ia resulting from the general SBIP approach in combination with the reduced interaction law from Sect. 4 have the same abstract form as in the case of macroscopic beam contact formulations [35].Most importantly, they are functions of the unilateral gap g ul and the mutual angle of the tangent vectors α.Due to this fact, the proofs presented in [35, Appendix B] likewise hold in this case and it is thus straightforward to conclude that also the SBIP approach in combination with the here proposed reduced interaction law fulfills objectivity, global conservation of linear and angular momentum as well as global conservation of energy.

Algorithmic complexity
The following discussion focuses on the part of evaluating the total interaction potential δ ia (and likewise its linearization [δ ia ]) as this is the one determined by the computational approach to molecular interactions of fibers.Depending on many other factors, this part may or may not be the dominating one in the entire algorithmic framework required for a nonlinear finite element solver in structural dynamics.Based on the experiences with the novel SBIP approach, the previously proposed SSIP approach, and the attempt of directly evaluating δ ia via 6D numerical integration (see Eq. ( 5)), it can however be stated that in a best case scenario the computational cost of evaluating this part is still a considerable one and in the worst case scenario-using direct 6D numerical integration-it becomes so costly that it is actually unfeasible for any system of practical relevance.This initial general assessment shows both the urgent need for an efficient approach and also the high leverage of any potential improvement in this respect, which has actually been the main motivation for the development of the novel SBIP approach.
To narrow down the broad topic of algorithm efficiency, this analysis can be restricted to the evaluation of one element pair, because the number of interacting element pairs can be considered a fixed number for now.This number mainly depends on the spatial distribution of the fibers the range of interaction, i.e., the radius, which means that it will be limited due to both the short range of interactions considered in this article and the non-penetrability constraint restricting the closest packing of fibers.Note that the associated important question of an efficient search for element pairs and the selection of element pairs to finally evaluated will be discussed in the following Sect.7.4 and Sect.7.5.
At this point, recall from the analysis in our previous contribution [40] the algorithmic complexity of the evaluation of one element pair in case of direct 6D numerical integration of Eq. ( 5) will be O(n 2  GP,tot,ele-length • n 4 GP,tot,transverse ).Here, n GP,tot,ele-length and n GP,tot,transverse denote the number of integration points in axial and transverse direction, respectively.The SSIP approach already reduces this complexity significantly to O(n 2 GP,tot,ele-length ) due to the of the 4D numerical integration over the cross-section areas by an effective section-section potential law.By the novel SBIP approach, this is finally reduced even further and for the remaining 1D integral along the slave beam (cf.Eq. ( 14)), we obtain O n GP,tot,ele-length (30) complexity.Bearing in mind the typically large number of integration points required to integrate the (short-ranged) molecular interaction laws with its high gradients with sufficient accuracy, this linear complexity makes a significant difference as compared to the quadratic complexity of the SSIP approach.Based on the experience of the numerical examples studied throughout this work, typical values are given as n GP,tot,ele-length = 10 . . .100.This is thus the factor we can expect to save from the reduced dimensionality of numerical integration when comparing the proposed SBIP approach to the previously derived SSIP approach, which itself offers potential savings by a factor of 10 4 . . . 10 8 as compared to direct 6D numerical integration [40].Of course, this comes at the cost of the closest point-to-curve projection required in case of the novel SBIP approach.This projection consists of solving the scalar nonlinear orthogonality condition (cf.Eq. ( 29)) e.g. by means of Newton's method, which however turns out to be rather insignificant as compared to evaluating the terms of the integrand.The net saving will thus be slightly smaller than the number of (axial) integration points per element, but still significant.
In addition to that, there is another positive effect to be considered.Due to the additional analytical integration step in the derivation of the reduced SBIP law from Sect. 4 as compared to the SSIP laws, the (inverse) exponent of the effective interaction law and thus integrand is reduced by one: cf.−m+9/2 in Eq. ( 17) vs. −m+7/2 in Eq. (A12) of Ref. [40], for m ≥ 6.In turn, this makes the integrand smoother and less integration points are required to achieve the same accuracy of the numerical integration.Especially for the short-ranged interactions e.g. from the LJ interaction with m = 6 and m = 12, this makes a significant difference in the decisive regime of small separations and contributes to the superior efficiency of the SBIP approach as compared to the SSIP approach or even the direct 6D numerical integration.
In this it seems noteworthy to point out the clear separation of general SSIP/SBIP approach and the applied reduced interaction law.Generally, the complexity of the specific expression does not necessarily depend on whether it is a SSIP or SBIP law.However, some conclusions like the one just made for the resulting exponent of the power law-if consistently derived from an inverse-power point pair potential law-are possible.Likewise, assuming homogeneous, circular cross-sections we can state that the mutual configuration of the disk-disk system has four degrees of freedom whereas the disk-cylinder system can be described by three degrees of freedom as discussed in Ref. [41].However, this does not allow to estimate the complexity of the specific expressions even in the hypothetical case of exact analytical interaction laws.Given the concrete examples of expressions for short-range interactions presented in Sect. 4 and our previous contribution [40], respectively, it is important to underline that they are based on different simplifying assumptions and thus naturally have a different accuracy [41].The fact that the SSIP law expressions are simpler than their SBIP law counterparts must thus be seen in the light of this tradeoff between simplicity and accuracy.Nevertheless, when comparing the performance in the numerical example to be presented in Sect.8.1, one will notice this effect of less operations being required to evaluate the simpler yet less accurate specific SSIP law as compared to the SBIP law. 6This contrary effect diminishes the observable net speed-up resulting from the superior efficiency of the general SBIP approach over the general SSIP approach described above.
To conclude this important assessment of the algorithm's efficiency, it can be stated that the general SBIP approach is significantly more efficient than the SSIP approach (which in turn is still significantly more efficient than a direct 6D numerical integration).This holds even despite the superior accuracy of the applied SBIP law (Eq.17), which is therefore also slightly more complex as compared to the simple SSIP law from Eq. (A12) of our previous contribution [40].In the numerical example of peeling two adhesive fibers to be presented in Sect.8.1, the combination of the novel SBIP approach and the specific SBIP law turns out to be approximately a factor of 4 faster than its SSIP counterpart.

Search for interacting pairs and partitioning for parallel computing
The search for interacting beam element pairs follows an efficient standard algorithm commonly referred to as bucket search (see e.g.Ref. [58] for details), which has already been used in combination with the SSIP approach [40].Due to the very short range of the interactions such as vdW and repulsive steric forces considered in this article, the requirements for the search algorithm are very similar to those from macroscopic (beam) contact formulations.The resulting small cut-off radius is beneficial with respect to both minimizing the number of interaction pairs to be evaluated and an effective partitioning of the problem to parallelize the evaluation on multiple processors without excessive cost for communication between the processors.Hereby, the partitioning of the problem is based on the spatial arrangement of the beam elements and uses the same subdivision of the computational domain into buckets already obtained from the search algorithm.A repartitioning and thus redistribution of the interaction pairs to the processors is done only if the spatial distribution of the beam elements has changed so much that-considering the cut-off radius -there is a chance that new interaction pairs need to be identified and evaluated.Generally, the computational cost of these steps of search and partitioning turned out to be insignificant as compared to the evaluation of the interaction pairs.Therefore, the parallelization of the pair evaluation on multiple processors indeed reduces the overall computation time significantly.

A criterion to sort out element pairs separated further than the cut-off radius before the actual evaluation
The following is applied as an additional step after the search for spatially proximate and thus potentially interacting element pairs.Motivated by the critical influence of the pair evaluation on the overall computational cost, this additional filtering step aims to sort out element pairs that are identified by the rather rough and conservative bucket search algorithm, but are further separated than the cut-off radius and will thus not contribute to the total interaction energy.The key idea is thus to skip the entire evaluation of the element pair, i.e., the loop over the integration points on the slave side, which otherwise would only after the closest point-to-curve projection check the cut-off radius and skip the evaluation of terms for this specific integration point.
To achieve high net savings, the applied criterion must be cheap to evaluate and is thus only based on the nodal positions and an estimate of the actual, deformed centerline geometry of the elements by means of so-called spherical bounding boxes.A very similar filter criterion has been applied in the context of macroscopic beam contact [35], where further details can be found.The only difference lies in the distance threshold value that is used.Here, we skip the pair evaluation if the minimal distance between the spherical bounding boxes is more than twice the cut-off radius, which should be on the safe side to not miss any contributions also in the case of strongly deformed elements.Still, the resulting decrease in the overall runtime observed for the numerical examples from Sect. 8 was up to 30%, which is quite remarkable and underlines the effectiveness of this additional filtering step.As an additional benefit, it has been observed that the number of non-unique or ill-posed closest point-to-curve projections significantly decreased and in fact vanished for the numerical examples considered in the context of this work.

Numerical examples
All presented formulations and algorithms have been implemented in C++ and integrated into the existing computational framework of the in-house research code BACI [59].Note that the implementation of the SBIP approach as well as the disk-cylinder potential has been verified by means of a second, independent implementation in MATLAB [60].Furthermore, the correct implementation of the discrete element residual vectors as well as tangent stiffness matrices has been verified by means of an automatic differentiation tool provided via the package Sacado, which is part of the Trilinos project [55].At this point, it is also important to emphasize that the novel SBIP approach seamlessly integrates into an existing nonlinear finite element solver for solid and structural mechanics.It does neither depend on any specific beam (finite element) formulation nor time discretization scheme, which underlines the versatility of this novel approach.In the following numerical examples, it has been used in combination with geometrically exact, Hermitian Simo-Reissner beam elements [16] and both in statics as well as an (implicit) time integration framework.The resulting nonlinear system of equations is highly challenging to solve, mainly due to the competition of the strongly nonlinear, deformationdependent adhesive and repulsive forces, which also leads to physical instabilities such as snapping free and snapping into contact.We used a Newton-Raphson algorithm with an additional step size control, which is described in more detail in Appendix C of our previous contribution [40].

Peeling and pull-off behavior of two adhesive elastic fibers
This numerical example has first been studied in the authors' previous contribution [56] where the SSIP approach [40] has been used to model vdW adhesion and steric repulsive forces.Fig. 7a shows the problem setup consisting of two parallel, straight fibers interacting via a LJ potential.
The idea is to study the entire process of separating these adhesive fibers starting from contact along the entire length up to the point where they would suddenly snap free.Therefore, a displacement u x in x-direction is prescribed at both ends of the right fiber and the sum of the reaction forces F x := F tr x + F br x in this direction is measured in order to obtain the characteristic, quasi-static force-displacement curve.
Most importantly, this simple setup with two initially straight beams allows to verify the accuracy of the SBIP approach by means of analytical reference solutions.the theoretical work for the scenario of infinitely long and parallel rigid cylinders presented in Ref. [40,Appendix A.2.2] is able to predict both the equilibrium spacing g LJ,eq,cyl cyl and the maximal magnitude of adhesive forces per unit length fmin,LJ,cyl cyl , i.e. the effective local pull-off force per unit length.It will be shown that these local quantities are excellently reproduced by the SBIP approach whereas the previously used simple SSIP law fails to do so (without additional calibration).While the results from both approaches agree very well in a qualitative manner, only the SBIP approach is thus able to yield the quantitatively correct pull-off force and other important global quantities characterizing this numerical example.
Moreover, it will be shown that the bias introduced by the choice of master and slave is negligible even for the considerably large fiber deformations in this example.First, this is confirmed by the symmetries of the example, which are excellently preserved in the numerical solutions.And second, using a modified setup with one rigid, straight beam and one deformable beam, we can quantify the error introduced by approximating the master beam geometry as a straight cylinder, because one of the two possible choices for master and slave will be the exact solution of the problem.For the other choice, we obtain a maximum relative error of 1.4% in the global pulling force even for the rather extreme scenario of maximal adhesive force and thus strongly bent fibers to be considered here.Altogether, this is a very important confirmation that the fundamental modeling approach of approximating one of the fibers as cylinder when calculating the two-fiber interaction potential leads to a very high accuracy in the case of very shortranged interactions considered here.

Setup and parameters
At this point, only the differences in the setup compared to the original numerical peeling experiment mentioned above will be presented.Refer to [56,Sec. 4] for a complete presentation of the setup, numerical methods and parameter values.Most importantly, the LJ interaction between the fibers is modeled by means of the novel SBIP approach from Sect. 3 in combination with the proposed disk-cylinder interaction potential law πm,disk-cyl from Eq. ( 17), which is used for both the attractive m = 6 as well as the repulsive m = 12 part of the LJ interaction.The two prefactors k 6 and k 12 specifying the LJ point-pair potential law will be varied to study their influence on the system response.Instead of using these prefactors, however, it seems more meaningful and intuitive in this context to use an equivalent set of parameters, which describe the equilibrium spacing g LJ,eq,cyl cyl and the maximal magnitude of adhesive forces per unit length fmin,LJ,cyl cyl , i.e. the effective local pull-off force per unit length, of straight, parallel fibers.According to the theoretical work for the scenario of infinitely long and parallel rigid cylinders presented in [40, Appendix A.2.2], these two alternative sets of two LJ parameters are bi-uniquely related 7 as follows: Again refer to [40, Appendix A.2.2] for the rather lengthy exact expressions of the scalar prefactors that are given here in their approximate floating point representation.For the sake of brevity, the subscript "cyl cyl" will be omitted in the remainder of this section, however, keep in mind that these quantities describe the academic case of infinitely long and parallel, rigid cylinders.Since the meaning of these two parameters is much more intuitive in the context of fiber adhesion, we can now argue that the equilibrium surface spacing g LJ,eq is a fundamental property of the type of physical interaction (and maybe also the material combination of both fibers and surrounding media) and keep it fixed for now with an exemplarily chosen value of g LJ,eq = 10 −3 = 0.05 • R = 2 × 10 −4 • l, which corresponds to five percent of the fiber radius.The minimal LJ force per unit length fmin,LJ on the other hand is a viable measure for the strength of adhesion and will be varied to study this important influence on the system behavior.Unless otherwise stated, two integration segments with ten Gauss points each is used for numerical integration of the LJ contributions in each of the 64 beam elements per fiber.It has been verified by means of refinement studies that the influence of the spatial discretization error and numerical integration error is negligible for all results to be presented in the following.The regularization strategy proposed in Sect.7.1 has been applied with a regularization separation g ul,reg = 8 × 10 −4 = 0.04 • R, which is smaller than the equilibrium spacing g LJ,eq given above and thus-as has been argued in Sect.7.1-led to identical results as compared to the non-regularized interaction law, yet ensures robustness and a significant reduction of the number of required nonlinear iterations by almost a factor of five.In order to further reduce the computational cost, the very short range of the LJ interaction has been exploited by applying a cut-off radius of r c = 0.1 = 5R, which again had no influence on the results.As mentioned above, all other parameters including the geometrical and material properties of the fibers remain unchanged as compared to [56,Sec. 4].

Qualitative quasi-static system behavior
Let us first look at the case of fmin,LJ = −1, which is the strongest attractive strength considered in this numerical experiment.The resulting equilibrium configurations for exemplarily chosen displacement values u x are shown in Fig. 8.Note that the configurations are symmetric with respect to both a vertical and a horizontal middle axis, which also implies for the individual reaction force components on the left/right ("l"/"r") side and top/bottom ("t"/"b") that F tl x = F bl x = −F tr x = −F br x and F y ≡ 0, as expected from the symmetric setup of the experiment and already observed in [56].It seems worth mentioning that this symmetry is not broken by the choice of master and slave beam as required by the SBIP approach-more on this important aspect in Sect.8.1.4.The corresponding force-displacement curve is shown in Fig. 7b (blue line) and a magnified view of the small displacement value range is provided in Fig. 7c.Overall, the shape of the curve is very similar to the one obtained using the previous SSIP approach in [56] and also the three characteristic phases of initiation of fiber deformation, peeling, and pull-off as obtained and analyzed for the case of electrostatic attraction in [56] are recognizable here.Moreover, the sharp force maximum for a very small displacement value in the initiation phase is confirmed by this study using the novel SBIP approach and Fig. 9 Detail study of the resulting line force distributions for two displacement values obtained for the case fmin,LJ = −1.For clarity, the fibers are depicted as their centerlines and forces are shown for the left fiber only.Note that both pictures show details of the entire system and are not to scale more accurate disk-cylinder interaction law πm,disk-cyl .It can thus be concluded that the known limitation regarding the accuracy of the previously applied, simple SSIP law π does not affect the qualitative analysis and conclusions drawn in Ref. [56]. 8A detailed comparison, including a quantitative analysis of the differences in the global system response for the identical set of parameter values will follow in a dedicated section below.Following up on the qualitative results, note that the very short range of the LJ interactions leads to the fact that fibers interact almost exclusively if or where they are in contact and that there is no perceptible far range effect as observed and measured in form of a second branch ("separated fibers") in the force-displacement plot e.g. for the case of electrostatic attraction in [56].This observation is also in agreement with the common notion that short-ranged interactions such as vdW adhesion only have an influence on the state of being/remaining in contact, and not the process of coming into contact from an initially separated state.
A quantitative analysis of the results will be given in the following section discussing the influence of the strength of adhesion, however, two brief aspects are stated here as an immediate verification of the results.First, the specified parameter value for the equilibrium gap of infinitely long and parallel, rigid cylinders g LJ,eq = 10 −3 is indeed recovered as a simulation result in the middle parts of the fibers wherever they are approximately parallel.And second, as can be seen from the visualization of the resulting LJ force distributions in Fig. 9, the maximum magnitude of the observed interaction forces per unit length agrees very well with the specified parameter value fmin,LJ = −1.
The interesting shape of the interaction force distribution shown in Fig. 9a reveals the competing nature of the adhesive and repulsive forces.At the point, where the smallest surface separation of the fibers transitions from below to beyond the equilibrium distance, there is a sharp transition from strong repulsive to strong adhesive forces.Due to the bending resistance of the fibers, there must be a repulsive region when separating the two fibers.Only in the middle part, the fibers maintain the net force-free equilibrium surface separation g LJ,eq .Fig. 9b finally shows a state where the fibers have been pulled further apart than g LJ,eq even at theirs closest points such that solely adhesive forces prevail.

Influence of the strength of adhesion
A more general, theoretical study of parameter influences in systems of adhesive fibers by means of nondimensionalization of the governing partial differential equations is provided in Appendix A. Here, it is complemented by a numerical study considering the influence of the adhesive strength for this fundamental two-fiber peeling and pull-off example.Thus, keeping all other parameters unchanged, the minimal LJ force per unit length (i.e. the maximal adhesive force magnitude) of infinitely long and parallel cylinders fmin,LJ is varied over two orders of magnitude such that both limits of "strong adhesion/low fiber stiffness" and "weak adhesion/high fiber stiffness" can be observed.The resulting forcedisplacement curves for fmin,LJ = {−1, −0.1, −0.01} (blue solid line, red solid line, black dashed line) are shown in Fig. 7b, c and the corresponding final equilibrium configurations ultimately before snapping free are displayed in Fig. 10.
As expected, an increase of the adhesive strength generally leads to increased reaction force values and higher displacement values before the fibers would finally snap free.Specifically, we obtain normalized force peak values of Fx,max ≈ 0.32 at u x /l ≈ 3.4 × 10 −4 for fmin,LJ = −0.01,Fx,max ≈ 1.8 at u x /l ≈ 3.6 × 10 −4 for fmin,LJ = −0.1, and Fx,max ≈ 9.9 at u x /l ≈ 3.6 × 10 −4 for fmin,LJ = −1.The peak force values thus increase by a factor of approximately 5.6 and 31 if the minimal LJ force per unit length is increased by a factor 10 and 100, respectively.However, the location of the force peak is not influenced by the varying adhesive strength, which meets our expectations from keeping the equilibrium surface separation g LJ,eq fixed.Finally, the maximum normalized displacement values u x,max /l before the fibers would snap free are observed to be approximately 0.13, 0.48, and 0.82 for fmin,LJ = {−0.01,−0.1, −1}, respectively.Note however that these maximum displacement values are the result of a quasi-static analysis and that the exact point of snapping free will depend on the dynamics of the system.A final interesting observation based on the force-displacement curves in Fig. 7b is that both pairs of curves for fmin,LJ = {−0.01,−0.1} and fmin,LJ = {−0.1,−1} approximately correspond to each other in a certain, limited displacement range u x /l ≈ [0.08, 0.13] and u x /l ≈ [0.25, 0.48], respectively, before the system with weaker adhesive strength would snap free.

Influence of the choice of master and slave
Following up on the discussion of how to assign the roles of master and slave in Sect.6, this numerical example is used to study the influence of this choice and thus bias in the SBIP approach.For this sake, all three cases shown in Fig. 7b have been repeated with switched roles of master and slave and the resulting differences in the reaction force values along the entire force-displacement curve have been evaluated.The maximal relative differences along the entire curve are approximately {4 × 10 −6 , 5 × 10 −6 , 3 × 10 −5 } for fmin,LJ = {−0.01,−0.1, −1}, respectively, and thus are of the same order of or even smaller than the spatial discretization error that has been verified by doubling the number of elements.This is a very important confirmation that the introduction of master and slave and the corresponding approximation of the master beam's geometry as cylinder in the calculation of the two-fiber interaction potential does not break the inherent symmetries of the problem setup and is thus a reasonable modeling assumption.
To quantify the model error associated with replacing the master beam by a surrogate cylinder, the setup of the example is now slightly modified by constraining the left fiber's deformation to zero and thus preserving the initial straight cylinder shape.Using the left fiber as the master beam will thus eliminate the corresponding model error and serves as a reference solution for the subsequent simulation run with switched roles of master and slave.This results in maximal relative differences in the pulling force values of approximately {0.13%, 0.43%, 1.42%} for fmin,LJ = {−0.01,−0.1, −1}, once again measured along the entire force-displacement curve.Given the considerable magnitude of curvature for the highest adhesion force value, which is even slightly more than in the original example visualized in Fig. 10, this is considered to be a reasonably small error level, which will presumably be negligible in most practical applications of the novel SBIP approach.

Comparison of SSIP and SBIP approach and corresponding reduced interaction laws
At this point, recall from the prior analysis of the model accuracy using the example of two cylinders in Sect. 4 that the previously used simple SSIP law πm,disk disk from Ref. [40] overestimates both the total interaction potential and force in the decisive regime of small separations and that the novel SBIP law πm,disk-cyl from Eq. ( 17) on the contrary ensures the correct scaling behavior and is thus significantly more accurate.The influence of this model error on the global system behavior shall be briefly investigated by means of the numerical peeling and pull-off experiment considered in this section.Fig. 11 shows the corresponding force-displacement curves for the original set of LJ parameters k 6 = −10 −7 , k 12 = 5 × 10 −25 from [56], using either the previous SSIP law πm,disk-cyl (red solid line) or the proposed SBIP law πm,disk-cyl from Eq. ( 17) and corresponding approach from Sect. 3 (blue solid line).
The results confirm the assessment summarized above because the peak force value Fx,max is overestimated by a factor of approximately 2.6 and the maximum displacement value before snapping free u x,max /l is overestimated by a factor of approximately 1.3 if using the simple SSIP law from the authors' previous contribution [40].On the other hand, the location of the force peak at u x /l ≈ 3.0 × 10 −4 , which is determined by the ratio of adhesive and repulsive contributions, is identical for both approaches.Most importantly, however, is that the qualitative shape of the force-displacement curve is very similar for both approaches and indeed almost the identical force-displacement curve  Fig. 11 Comparison of the results obtained in [56] via the previous SSIP law πm,disk disk from [40] with the results for the proposed SBIP law πm,disk-cyl from Eq. ( 17) and corresponding approach from Sect. 3 can be obtained from the consistent SBIP law πm,disk-cyl if a different parameter value set k 6 = −4 × 10 −7 , k 12 = 2 × 10 −24 is used (black dashed line) 9 .This confirms the argument given in [40] that the simple SSIP law πm,disk disk can be calibrated to compensate for the overestimation of the interaction potential in the small yet decisive range of separations around the equilibrium separation of the LJ potential.In this manner, the simple SSIP law πm,disk disk yields results, which -in good approximation -agree with the consistent SBIP law on the system level, e.g. the reaction force-displacement curve studied in this numerical example.
After the accuracy, let us now briefly compare the efficiency of both approaches, which has been discussed on the theoretical level of algorithmic complexity in Sect.7.3, by looking at the computational cost for one simulation run of the numerical example considered here.Keeping in mind the general limitations of such a simple performance comparison and especially the unoptimized implementation of both approaches, we found that the SBIP approach is approximately a factor of 3.8 faster than the SSIP approach 10 .The advantage of the SBIP approach that less integration points will be required for the same level of numerical integration error due to the smaller exponent of the reduced interaction law as argued in Sect.7.3 has not been exploited in this comparison to isolate the net effect of the decreased algorithmic complexity of the approaches on the one hand and increased complexity of the reduced interaction law on the other hand.It is thus confirmed by this numerical experiment that the novel SBIP approach in combination with the proposed disk-cylinder interaction law πm,disk-cyl is significantly more efficient than the previous SSIP approach for the case of short-ranged interactions.

Adhesive nanofiber-grafted surfaces
This set of numerical examples deals with bioinspired adhesive surfaces, which have recently gained a lot of attention (see e.g.[61] for a review).This topic is motivated by Fig. 12 Initial, undeformed configuration of the numerical experiments studying adhesive nanofiber-grafted surfaces the fascinating skills of geckos [62], spiders [63], mussels [64], and other animals to stick to surfaces and defy significant detachment forces such as for instance their body weight when sitting on the ceiling or steeply inclined surfaces.The origin of this extraordinary adhesion is known to lie on the molecular scale and vdW forces between hierarchical arrays of nano-hairs and the surface have been shown to play a major role [65,66].
In this context, the following computational study suggests and analyzes a possible design for artificial adhesive surfaces based on grafted helical nanofibers.Hereby, nature's ubiquitous pattern of helical fibers is used to achieve the desired large ratio between strong adhesion under load on the one hand and easy removal of surfaces on the other hand.This numerical example has been studied in full detail in [67] and only a selection of results will be presented in the following.These results aim to demonstrate the capability of the novel SBIP approach to handle arbitrary mutual configurations of fibers in 3D by means of an initially curved spatial geometry as well as large, non-trivial deformations of fibers.Also the efficiency of the novel approach will be showcased by applying it to a large-scale system regarding both the number of fibers and time steps.

Setup and parameters
Figure 12 shows the setup of this numerical example from different view angles.As mentioned above, the system is mainly composed of helical fibers with a helix diameter of D h = 2R h = 1 and a pitch, i.e., height of one complete turn, of P h = 1, respectively.The grafting onto surfaces is modeled by respective Dirichlet boundary conditions applied to the fibers and the dark colored surfaces shown in Fig. 12 are only used for visualization and not considered in the simulation.In this first minimal setup, 2x2x2 fiber loops are considered, which refers to 2 surfaces with 2 helices each with 2 complete turns, i.e., loops each.Two different scenarios will be studied.In the normal loading, also referred to as "Pull" scenario, the surfaces will be moved together until the surfaces of the fibers touch and subsequently pulled apart without any rotation.In the second, the "Twist & Pull" scenario, the surfaces will be twisted by 75 • before they are being pulled apart.The full specification of this numerical example including the time curves used for the prescribed boundary conditions can be found in Ref. [67].

Results and discussion
Figure 13 shows selected simulation snapshots for the normal loading referred to as "Pull" scenario.The configuration in Fig. 13a corresponds to the end of the approach (and equilibration) phase at t = 2, where the helical fibers of both surfaces are perfectly aligned and their surfaces touch.In the subsequent pull-off phase shown in Fig. 13b-d, the fibers are continuously peeled and the contact length decreases while the fibers are strongly deformed.Only after approximately t = 4.09, the last contact points of the fibers detach and the separation process is completed.
Figure 14 likewise shows selected simulation snapshots for the "Twist & Pull" scenario, which differs from the "Pull" scenario in the phase t > 2. Figure 14a-d cover the twisting phase and show that twisting the surfaces by 75 • indeed separates the fibers almost entirely, which agrees with the idea of a subsequent, relatively easy separation of surfaces in normal direction.A closer look reveals that the strength of adhesion suffices to keep the fibers in contact almost along the entire length up to approximately 37.5 • , which leads to quite large and complex deformations of the fibers.Note also that Fig. 14d shows a limitation of our current setup in a way that the fibers can penetrate the surfaces, because the surfaces are not explicitly included in the simulation and only accounted for by means of Dirichlet boundary conditions on the grafted fiber endpoints as mentioned above.This is however not expected to have a major influence on the pull-off forces to be discussed next, as it occurs rather rarely and from Fig. 14d it can be estimated that a preclusion of the observed penetration would not change the overall deformed shape of the fibers in a significant manner.
Finally, the force-displacement curves for both scenarios are shown in Fig. 15.Specifically, the reaction force component F z in the surface normal, i.e., z-direction, is summed over all nodes grafted onto the top surface and normalized by means of the maximum force F z,max,2×2×2 ≈ 0.356 that is observed for this 2 × 2 × 2 loop system.As desired, the  Fig. 15 Force-displacement curves for the decisive phase, i.e., 2 ≤ t < 5 in the "Pull" scenario and 2 ≤ t < 7 in the "Twist & Pull" scenario.Force values are normalized with respect to the maximal pull-off force F z,max adhesive connection of surfaces is able to withstand approximately two times larger pulloff forces in the normal loading scenario (red line) as compared to the scenario meant to be used for the removal of surfaces (blue dashed line).Generally, a rich and highly nonlinear system behavior can be observed from these force-displacement curves, which reflects the complex interplay of strongly adhesive LJ interactions and large and complex structural deformations in 3D.Particularly for the "Pull" scenario, some characteristics of the fundamental problem of peeling two initially straight adhesive fibers (cf.our previous article [56]) such as the sharp force peak, the extended pull-off phase and the sudden snap-free are however still recognizable.Note also the considerable compressive forces that occur in the "Twist & Pull" for ūz = 0, i.e., during the twisting phase 2 ≤ t < 7.In this context, it seems worth mentioning that the reaction forces in x-and y-direction are below 0.1 and thus more than ten times smaller than the maximal pull-off force observed here.

Upscaling the computational experiment
Increasing the size of the system is important to judge the performance of the novel SBIP approach.Therefore, a 2 × 8 × 8 and a 2 × 16 × 16 loop system will be considered now.With respect to the original 2 × 2 × 2 loop system, this corresponds to a scaling factor of 16 and 64, respectively, which leads to a total of 2048 elements/4112 nodes and 8192 elements/16416 nodes, respectively.Except for the system size, the setup and parameters stated in Sect.8.2.1 remain unchanged.Figure 16 shows selected simulation snapshots of the "Pull" scenario for the 2 × 16 × 16 loop system and confirms that the periodicity in the problem setup is also observable in the resulting deformed configurations.Once again, refer to Ref. [67] for a detailed investigation of the force-displacement curves obtained for these scaled system sizes.

Entanglement of fibers
Turning to the second considered scenario including the twisting of surfaces, an-at first glance unintended-entanglement of fibers has been observed, which also has a noticeable influence on the force-displacement curve.This can be traced back to the (unphysical) behavior that during the twisting phase the outermost loop from the top surface dives under the fiber endpoint grafted onto the bottom surface and thus leaves the fibers entangled for the remaining course of the simulation, as can be seen from the simulation snapshots in Fig. 17.While this is obviously neither physically correct nor desirable in the initial context of separating the adhesive fibers from each other, it is, however, the correct numerical solution given that no beam-solid contact has been considered between the fibers and the surfaces.And, most importantly, it nicely demonstrates the effectiveness and robustness of the proposed SBIP approach even for highly complex, large 3D deformations in large-scale systems that include a broad variety of mutual orientations and separations of fibers.

Summary, conclusions and outlook
Following up on the development of the first 3D beam-beam interaction formulation for molecular interactions-the so-called section-section interaction potential (SSIP) approach [40]-this article proposes an enhanced approach, which is specialized on the most challenging case of short-ranged interactions.Examples for such short-ranged interactions include the van der Waals (vdW) adhesion and the steric repulsion of the Lennard-Jones (LJ) interaction.Exploiting the characteristic rapid decay of the interaction potential with increasing distance allows to further reduce the dimensionality of the beam-beam interaction problem and thus achieve a significantly more efficient formulation than the more general SSIP approach.The key idea is to approximate the second, arbitrarily deformed beam by a surrogate body with trivial geometry, which is located at the closest point from a given point on the first beam and oriented according to the centerline tangent vector of the second beam.Mathematically, this surrogate body is the linear series expansion of the second, so-called master beam around the closest point.In this manner, the interaction of the entire second beam with a given section of the first, so-called slave beam can be described by an analytical, closed-form interaction potential law, which replaces the numerical integration along the master beam by a single func-tion evaluation.This novel, specialized approach is therefore called the section-beam interaction potential (SBIP) approach.
In addition to being significantly more efficient, the novel formulation developed in this article is also significantly more accurate than the one used before.This is due to the fact that the underlying reduced interaction law is superior to the previously used simple SSIP law, because it considers the relative rotation, i.e., the orientation of the interacting bodies in addition to their scalar separation.The specific SBIP law used in the present work has been derived in the authors' recent contribution [41], considering the case of circular, homogeneous cross-sections and a generic inverse power law with exponent m ≥ 6 as fundamental point-pair interaction potential.Most importantly for this work, the SBIP law ensures the correct asymptotic scaling behavior in the decisive regime of small separations, which has been identified as the most critical aspect when using the simple SSIP law from Eq. (A12) of our previous contribution [40], which neglects the mutual orientation of cross-sections.This important property has been verified numerically by an exemplarily chosen test case considering two straight fibers at different relative orientations.Moreover, using the numerical example of peeling two adhesive elastic fibers, it has been shown that only the SBIP law reproduces the quantitatively correct local adhesive force as well as equilibrium surface separation and therefore also the decisive global pull-off force.
These two novelties-the general SBIP approach from Sect. 3 and the specific analytical interaction law from Ref. [41]-have been combined in Sect. 5 to obtain the resulting virtual work contribution.The subsequent discretization by means of finite elements and its consistent linearization finally give rise to a novel beam-beam interaction formulation that can be seamlessly integrated in an existing nonlinear finite element framework for structural mechanics.In this respect, it is important to point out that the approach does neither depend on a specific beam (element) formulation nor time integration scheme and is thus considered to be highly versatile.Further methodological aspects such as a suitable numerical regularization of the characteristic singularity at zero separation in the reduced interaction law and a criterion to sort out beam element pairs with a larger separation than the cutoff radius before the actual evaluation to save computational resources have been presented in Sect.7. By means of a numerical example studying the peeling and pull-off behavior of elastic adhesive fibers, the higher efficiency of the novel, specialized SBIP approach as compared to the previous, more general SSIP approach has been quantified.In this example and even without using all savings the proposed formulation offers, the novel approach is approx.4 times faster than the previous approach.This is slightly less than the factor of 10, which would be expected from the comparison of algorithmic complexities, because the applied reduced SBIP interaction law (Eq.17) is more complex due to the significantly enhanced accuracy if compared to the simple SSIP law (Eq.(A12) of our previous contribution [40]).It can be concluded that the combination of a better accuracy and higher efficiency leads to a drastically improved applicability-for instance with respect to the maximum feasible system size and time scales-and therefore opens the door to tackle a multitude of yet intractable problems in science and technology.This is demonstrated by the comprehensive numerical example studying adhesive nanofiber-grafted surfaces, which confirms the effectiveness, efficiency and robustness of the novel formulation also for large-scale, complex systems with arbitrary mutual configurations and large deformations of the interacting fibers.In addition, this example outlines the working principle of the proposed design of the nanofiber-grafted surfaces based on arrays of helical fibers by showing the desired large ratio of strong adhesive connection under load and easy removal of the surfaces if desired.This showcases how suitable simulation tools such as the one developed in this article could contribute to the design and manufacturing of bioinspired artificial adhesives in the long run.
A future derivation of other SBIP laws for instance for the case of screened electrostatics or hydrophobic/-philic interactions would further extend the range of applications of the novel SBIP approach.Given the importance of the mentioned interaction types in many biological systems on the nano-and microscale, this is considered a highly promising extension of the work presented here.Moreover, the identification and calibration of the parameters of reduced interaction laws based on experimental results (such as the forcedisplacement curves for simple characteristic setups with two interacting fibers) would be an interesting avenue of future research that is expected to be of great benefit to the prediction quality of computational experiments.

Appendix A: Dimensionless key parameters of a system with adhesive elastic fibers
This section aims to identify the dimensionless parameters that govern the fundamental behavior of a mechanical system with elastic, adhesive fibers.Considering the spacecontinuous, static problem without external loads, the governing equation for the total potential energy (Eq.6) simplifies to the internal energy contributions from elastic forces and moments int , which is well-known from beam theory, and the energy contribution from the interaction ia .For the latter, we use the most general form (Eq. 5) based on two nested volume integrals of the point-pair potential (r) and consider the LJ potential with its adhesive and repulsive component as an example.The entire set of parameters carrying units is thus given by the length L and radius R of the fibers, its Young's modulus E, the atom density ρ, and finally the prefactors of the adhesive and repulsive part of the LJ law k 6 and k 12 , respectively.Note that instead of these two prefactors, an equivalent set of parameters r LJ,eq and LJ,eq specifying the equilibrium distance and the corresponding minimal potential value of the LJ law is commonly used.Nondimensionalization of this equation by means of normalization of the primary variables with suitable length measures allows to identify the following (non-unique) set of dimensionless parameters: The slenderness ratio ζ is known to be the only dimensionless parameter of the static elastic problem of beams and e.g.determines the amount and thus relevance of shear deformation.The second and third parameter are specific for the interactions between beams and given for the two alternative sets of LJ parameters mentioned above.Using r LJ,eq and LJ,eq , these parameters nicely disclose their meaning as equilibrium distance measure and measure for the strength of adhesion relative to the fibers' structural rigidity.The latter is thus named "adhesive compliance", because it will determine how much the fibers where g ul = n T ul ( r 1 − r 2 ), (B7) and are analogous to the variations of these quantities given in Eq. ( 23) and Eq. ( 24).The required second derivatives of the disk-cylinder interaction potential π from Eq. ( 17) can-like the first derivatives in Eq. ( 21) and Eq. ( 22)-be conveniently expressed in recursive manner as (B11) and the missing second derivative of π with respect to cos α: The linearization of the variation of the unilateral gap reads Note the analogy to the corresponding variation δξ 2c from Eq. ( 26) and the definition of p 2,ξ 2 given in Eq. (29).
Remark.Up to this point, the number and complexity of the expressions is comparable to the ones required in the macroscopic line contact formulation proposed in [34] and the combination, i.e., blending 11 thereof with a point contact formulation in [35].Indeed, only the expressions that depend on the applied disk-cylinder interaction potential law deviate from the line contact case, where most commonly a quadratic penalty potential π = 0.5 ε g 2 ul is used.Here, the additional expression for the linearization of the variation of the cosine of the mutual angle Note that the correctness of the presented and implemented analytical linearization has been verified by means of automatic differentiation [55] of the virtual work expression from Eq. (20).Due to the high complexity of [δξ 2c ], it would be interesting to investigate how a neglect of this contribution would influence the computational cost on the one hand and the performance of the nonlinear solver on the other hand.In the context of this work, however, always the fully consistent linearization has been used.

Appendix C: Additionally required expressions for the regularization of the reduced disk-cylinder interaction law
The equation describing the proposed quadratic extrapolation below g ul,reg for a general section-beam interaction potential law π is given as follows:  Note that this regularization needs to be applied to both the adhesive m = 6 and the repulsive m = 12 part of the LJ potential law.

Fig. 4
Fig.4 Interaction potential of two cylinders as a function of the dimensionless minimal surface separation g bl /R at different mutual angles α.Comparison of the previously used SSIP law π6,disk disk from Eq. (18) in Ref.[40] (used together with the SSIP approach proposed in the same article; dark green line with crosses) with the analytical expression for the disk-cylinder potential π6,disk-cyl Eq. (17) (used together with the SBIP approach from Sect.3; red line with triangles).The numerical reference solution obtained via 3D Gaussian quadrature of the point-half space potential 6,pt-hs from Ref.[41] (green line with diamonds) and the analytical reference solutions summarized in Sect.2.3 (black dashed line) are plotted as reference

Fig. 5
Fig.4 Interaction potential of two cylinders as a function of the dimensionless minimal surface separation g bl /R at different mutual angles α.Comparison of the previously used SSIP law π6,disk disk from Eq. (18) in Ref.[40] (used together with the SSIP approach proposed in the same article; dark green line with crosses) with the analytical expression for the disk-cylinder potential π6,disk-cyl Eq. (17) (used together with the SBIP approach from Sect.3; red line with triangles).The numerical reference solution obtained via 3D Gaussian quadrature of the point-half space potential 6,pt-hs from Ref.[41] (green line with diamonds) and the analytical reference solutions summarized in Sect.2.3 (black dashed line) are plotted as reference

Fig. 6
Fig. 6 Classification of formulations for beam-beam interactions Quasi-static force-displacement curve.Force values to be interpreted as multiple of a reference point load that causes a deflection of ∕4 if applied at the fiber midpoint.
Detail view for small displacement values.

Fig. 7
Fig. 7 Numerical peeling experiment with two adhesive elastic fibers interacting via the LJ potential

Fig. 8 A
Fig. 8 A selection of characteristic equilibrium configurations obtained for the case fmin,LJ = −1

Fig. 10
Fig. 10 Comparison of the final configurations before snapping free for three values of the min.LJ force per unit length fmin,LJ Quasi-static force-displacement curve.Force values to be interpreted as multiple of a reference point load that causes a deflection of ∕4 if applied at the fiber midpoint.
Detail view for small separations.

Fig. 13
Fig.13 Sequence of simulation snapshots of the "Pull" scenario.Top surface is hidden for better visibility of the fibers

Fig. 14
Fig.14 Sequence of simulation snapshots of the "Twist & Pull" scenario.Top surface is hidden for better visibility of the fibers.Note that until time t = 2, this scenario is identical to the "Pull" scenario shown in Fig.13

Fig. 16
Fig. 16 Selected simulation snapshots of the "Pull" scenario with 2 × 16 × 16 loops.Top surface is hidden for better visibility of the fibers

Fig. 17
Fig. 17 Sequence of simulation snapshots of the "Twist & Pull" scenario with 2 × 8 × 8 loops.Top surface is hidden for better visibility of the fibers (and outlined by the black rectangle in the views from the top)

1 2 d 2 π d g 2
ul,reg ) + d π d g ul g ul,reg (g ul − g ul,reg ) + ul g ul,reg (g ul − g ul,reg ) 2 , g ul < g ul,reg π (g ul ), g ul ≥ g ul,reg .(C26) Note that the required derivatives of this regularized potential law πreg accordingly change to g ul,reg + d 2 π d g 2 ul g ul,reg (g ul − g ul,reg ), g ul < g ul,g ul,reg , g ul < g ul,reg d π d g ul,reg (g ul ), g ul ≥ g ul,reg α) d g ul g ul,reg (g ul − g ul,reg ) . . .α) d g 2 ul g ul,reg (g ul − g ul,reg ) 2 , g ul < g ul,reg d π d(cos α) (g ul ), g ul ≥ g ul,reg α) d g ul g ul,reg + d 3 π d(cos α) d g 2 ul g ul,reg (g ul − g ul,reg ), g ul < g ul,reg d 2 π d(cos α) d g ul (g ul ), g ul ≥ g ul,reg α) 2 d g ul g ul,reg (g ul − g ul,reg ) . . .α) 2 d g 2 ul g ul,reg (g ul − g ul,reg ) 2 , g ul < g ul,regd 2 π d(cos α) 2 (g ul ), g ul ≥ g ul,reg .(C31)These expressions for the regularized general interaction potential π and its first and second derivatives replace the original expressions presented in the context of Sect. 4 and Sect. 5 and Appendix B, respectively.For the specific disk-cylinder interaction potential from Eq. (17), most of the derivatives required on the right hand side of these equations have already been presented in Sect. 5 and Appendix B, and the additionally required third and fourth derivative follow (recursively) as d3 πm,disk-cyl d(cos α) 2 d g ul = r 2 (ξ 2c )]+ v T α2 [δr 2 ] + (r T 2 v α2 ) [δξ 2c ] sgn(t T 1 t 2 ),(B19)is required for a consistent linearization and in fact most of the complexity comes from this contribution or, to be more precise, from the linearization of the variation of the element parameter of the closest-point on the master side