A Computational Model for Molecular Interactions Between Curved Slender Fibers Undergoing Large 3D Deformations With a Focus on Electrostatic, van der Waals and Repulsive Steric Forces
AA Computational Model for Molecular Interactions Between CurvedSlender Fibers Undergoing Large 3D Deformations With a Focus onElectrostatic, van der Waals and Repulsive Steric Forces
Maximilian J. Grill ∗ , Wolfgang A. Wall, Christoph Meier Technical University of Munich, Institute for Computational Mechanics, Boltzmannstr. 15, 85748 Garching b. M¨unchen,Germany
Abstract
This contribution proposes the first 3D beam-to-beam interaction model for molecular interactions betweencurved slender fibers undergoing large deformations. While the general model is not restricted to a specificbeam formulation, in the present work it is combined with the geometrically exact beam theory and dis-cretized via the finite element method. A direct evaluation of the total interaction potential for general 3Dbodies requires the integration of contributions from molecule or charge distributions over the volumes ofthe interaction partners, leading to a 6D integral (two nested 3D integrals) that has to be solved numer-ically. Here, we propose a novel strategy to formulate reduced section-to-section interaction laws for theresultant interaction potential between a pair of cross-sections of two slender fibers such that only two 1Dintegrals along the fibers’ length directions have to be solved numerically. This section-to-section interactionpotential (SSIP) approach yields a significant gain in efficiency, which is essential to enable the simulationof relevant time and length scales for many practical applications. In a first step, the generic structureof SSIP laws, which is suitable for the most general interaction scenario (e. g. fibers with arbitrary cross-section shape and inhomogeneous atomic/charge density within the cross-section) is presented. Assumingcircular, homogeneous cross-sections, in a next step, specific analytical expressions for SSIP laws describingshort-range volume interactions (e. g. van der Waals or steric interactions) and long-range surface interac-tions (e. g. Coulomb interactions) are proposed. Besides ready-to-use expressions for the total interactionpotential, also the resulting virtual work contributions, its finite element discretizations as well as a suitablenumerical regularization for the limit of zero separation are derived. The validity of the SSIP laws as wellas the accuracy and robustness of the general SSIP approach to beam-to-beam interactions is thoroughlyverified by means of a set of numerical examples considering steric repulsion, electrostatic or van der Waalsadhesion.
Keywords: slender continua, molecular interactions, geometrically exact beam theory, finite elementmethod, intermolecular potentials, van der Waals interaction, electrostatic interaction, steric exclusion
1. Introduction
Biopolymer fibers such as actin, collagen, cellulose and DNA, but also glass fibers or carbon nanotubes areubiquitous examples for slender, deformable structures to be found on the scale of nano- to micrometers. Onthese length scales, molecular interactions such as electrostatic or van der Waals (vdW) forces are of utmostimportance for the formation and functionality of the complex fibrous systems they constitute [1, 2, 3].Biopolymer networks such as the cytoskeleton or the extracellular matrix, muscle fibers, Gecko spatulae or ∗ Corresponding author
Email address: [email protected] (Maximilian J. Grill)
Preprint submitted to International Journal for Numerical Methods in Engineering August 1, 2019 a r X i v : . [ c s . C E ] J u l hromosomes are some of the most popular examples. To foster the understanding of such systems, whichin turn allows for innovations in several fields from medical treatment to novel synthetic materials, there isan urgent need for powerful simulation tools. Finite element formulations based on the geometrically exactbeam theory [4, 5, 6] are known to model the transient (elastic) deformation of these slender structuresin an accurate and efficient manner. However, no corresponding numerical methods for above mentionedmolecular interactions between deformable fibers have been published yet. We thus aim to develop methodsthat both accurately as well as efficiently describe these molecular phenomena based on the geometricallyexact beam theory in order to ultimately solve relevant practical problems on the scale of complex systemsconsisting of a large number of fibers in arbitrary arrangement.A comprehensive review of the origin, characteristics and mathematical description of intermolecularforces can nowadays be found in (bio)physical textbooks [2, 3]. The critical point is to transfer the firstprinciples formulated for the interaction between atoms or single molecules to the interaction betweenmacromolecules such as slender fibers. Here, the analytical approaches to be found in textbooks and alsoin recent contributions [7, 8, 9, 10] from the field of theoretical biophysics are (naturally) restricted toundeformable, rigid bodies with primitive geometries such as spheres, half spaces or, most relevant in ourcase, cylinders. Some computational approaches can be found in the literature, but rather aim at includingmore complex phenomena such as retardation and solvent effects in vdW interactions [11], still limited torigid bodies. All-atom simulation methods like molecular dynamics do not suffer from this restriction, buttheir computational cost is by orders of magnitude too high to be applied to the relevant, complex biologicalsystems mentioned in the beginning and thus currently limited to time scales of nano- to microseconds [2].On the other hand, studying the deformation of elastic, slender bodies has a long history in mechanics andtoday’s geometrically exact finite element formulations for shear-deformable (Simo-Reissner) as well as shear-rigid (Kirchhoff-Love) beams have proven to be both highly accurate and efficient [4, 5, 6]. Moreover, contactinteraction between beams has been considered in a number of publications, e. g. [12, 13, 14, 15, 16, 17, 18,19, 20, 21]. However, all these methods are motivated by the macroscopic perspective of non-penetratingsolid bodies rather than the microscopic view considering first principles of intermolecular repulsive forces.The combination of elastic deformation of general 3D bodies and intermolecular forces has first been con-sidered by Argento et. al. [22] for small deformations, by Sauer and Li [23] for large deformations and finallyby Sauer and Wriggers [24] also for three-dimensional problems. In order to reduce the high computationalcost associated with the required high-dimensional numerical integrals, a possible model reduction frombody to surface interaction in case of sufficiently short-ranged interactions as e. g. predominant in (adhesive)contact scenarios has already been addressed in these first publications and has been the focus of subsequentpublications [25, 26]. However, since these formulations aim to describe the interaction between 3D bodiesof arbitrary shape, the inherent complexity of the problem still requires a four-dimensional integral overboth surfaces in case of surface interactions and a six-dimensional integral over both volumes for volumeinteractions, respectively. In contrast, beam theory describes a slender body as a 1D Cosserat continuum,such that a further reduction in the dimensionality and thus computational cost can be achieved. So far,this idea has been applied to describing the interaction between a beam and an infinite half-space in 2D asa model for the adhesion of a Gecko spatula on a rigid surface [24, 27] and later also for the interaction of acarbon nanotube with a Lennard-Jones wall in 3D [28]. In both cases, the influence of the rigid half spacecan be evaluated analytically and formulated as a distributed load on the beam.To the best of the authors’ knowledge, no approach for describing molecular interactions between curved3D beams for arbitrary configurations and large deformations has been proposed yet. Notable previousapproaches to similar problems have made simplifying assumptions. Ahmadi and Menon [29] study theclumping of fibers due to vdW adhesion by means of an analytical 2D beam method, yet only includevdW interaction between the hemispherical tips based on an analytical expression for the interaction of twospheres. A numerical study of the influence of inter-fiber adhesion on the mechanical behavior of 2D fibernetworks assumes an effective adhesion energy per unit length of perfectly parallel fiber segments and solvesfor the unknown contact length in a second, nested minimization algorithm [30].In this article, we propose the first model specifically for molecular interactions between arbitrarily curvedand oriented slender fibers undergoing large deformations in 3D. While the general model is not restrictedto a specific beam formulation, in the present work it is combined with the geometrically exact beam theory2nd discretized via the finite element method. This novel approach is based on reduced section-to-sectioninteraction potential (SSIP) laws that describe the resulting interaction potential between a pair of cross-sections as a closed-form analytical expression. Thus, the two-body interaction potential follows from twonested 1D integrals over this SSIP law along both fibers’ axes, which are evaluated numerically. In this way,the proposed, so-called SSIP approach significantly reduces the dimensionality of the required numericalintegration from six to two, and hence the associated computational cost. As compared to methods for 3Dsolid bodies or even to all-atom methods, this gain in efficiency opens up new fields of applications, e. g., thecomplex biological systems mentioned above. Regarding the practicability of the SSIP approach, it is alsoimportant to emphasize that it can be seamlessly integrated into an existing finite element framework for solidmechanics. In particular, it does neither depend on any specific beam formulation nor the applied spatialdiscretization scheme and in the context of the present work, it has exemplarily been used with geometricallyexact Kirchhoff-Love as well as Simo-Reissner type beam finite elements. Likewise, it is independent of thetemporal discretization and we have used it along with static and (Lie group) Generalized-Alpha timestepping schemes as well as inside a Brownian dynamics framework.For the proposed SSIP laws, which can either be derived analytically or postulated and fitted e. g. toexperimental data, we first present the most general form describing the interaction between arbitrarilyshaped cross-sections with inhomogeneous distribution of the elementary interacting points (e. g. atoms orcharges). Subsequently, we focus on homogeneous, circular cross-sections and propose specific, ready-to-use SSIP laws for vdW adhesion, steric repulsion and electrostatic interaction. Based on the fundamentaldistinction into either short-range or long-range interactions, we present the required steps and theoreticalconsiderations underlying the analytical derivation of the SSIP laws in a general manner, starting fromfirst principles in form of a point pair interaction potential that is described by a power law with generalexponent. Besides the expressions for the total interaction potential, also the corresponding virtual workcontributions, its finite element discretization and the consistent linearization are presented. Due to thecharacteristic singularity of molecular interactions in the limit of zero separation, also a suitable numericalregularization of the SSIP laws will be proposed.The remainder of this article is structured as follows. Section 2 briefly summarizes the fundamentalconcepts and theory of molecular interactions. Along with the fundamentals of the geometrically exact beamtheory to be introduced in Section 3, this forms the basis for the novel SSIP approach to be proposed inSection 4. This general approach will then be applied to specific types of physical interactions, namely vdW,steric and electrostatic interactions in Section 5. In Section 6, we turn to the finite element discretizationof the newly developed numerical methods and discuss some important algorithmic aspects such as theregularization of the reduced interaction laws and the algorithmic complexity of the reduced approach.Finally, the accuracy of the proposed SSIP laws as well as the general SSIP approach to beam-to-beaminteraction is validated by means of analytical as well as numerical reference solutions for academic testcases in Section 7.1. In a series of numerical examples including steric repulsion, electrostatic or vdWadhesion, the effectiveness and robustness of the novel approach is verified in the remainder of Section 7.We conclude the article in Section 8 and provide an outlook to promising future enhancements of the novelapproach.
2. Fundamentals of Intermolecular Forces and Potentials
Interactions between molecules may result from various physical origins and are a complex and highlyactive field of research within the community of theoretical as well as experimental physics. The methods tobe derived in this work make use of the most essential and well established findings as summarized e. g. inthe textbooks [2] and [3]. This section briefly presents a selection of aspects relevant for this work.
To begin with, a number of universal aspects characterizing molecular interactions, especially with regardto the numerical methods to be developed in this work, shall be presented. A few simple facts about theexamples mostly considered throughout this work, namely electrostatic and vdW interaction, are presented3traight away, whereas the details on these and further types of molecular interactions are to be discussedin the subsequent Section 2.2.
A collection of characteristics of molecular interactions with high relevance for this work. • Type of elementary interaction partners
Interaction may originate from unit charges as in thecase of electrostatics. Another popular example are vdW effects that are caused by fluctuatingdipole interactions occurring in every molecule and hence are related to the molecular density ofthe material. • Spatial distribution of elementary interaction partners
Thinking of the resulting interactionbetween two bodies as accumulation of all molecular interactions, the question for the locationsof all elementary interaction partners arises. Charges can often be found on the bodies’ surfaceswhereas molecules relevant for vdW interactions spread over the entire volume of the bodies.This work focuses on solid bodies (i. e. condensed matter) that are non-conducting such thatinteraction partners will not redistribute, i. e., change their position within a body. • Distance-dependency of the fundamental potential law
Generally, the strength of molecularinteractions decays with increasing distance. Most frequently, inverse power laws with differentexponents or exponential decay can be identified. • Range of interactions
As a result of the previous aspects, a range of significant strength of aninteraction can be defined. Rather than an inherent property, the classification of long- versusshort-ranged interactions is a theoretical concept to judge the perceptible impact in specificscenarios. Moreover, it is a decisive factor in the derivation of well-suited numerical methods. • Additivity and higher-order contributions
Many approaches including the one presented in thiswork make use of superposition, i. e., accumulating all the individual contributions from ele-mentary interaction partners to obtain the total effect of interaction. This assumes that theinteractions behave additively, i. e., that the sum of all pair-wise interactions describes the overallinteraction sufficiently well. More specifically, the presence of other elementary interaction part-ners in the surrounding must not have a pronounced effect as compared to an isolated system ofan interacting pair. Otherwise, the sum of all pair-wise interactions would need to be extendedby contributions from sets of three, four and more elementary interaction partners.As a matter of course, this list is not exhaustive but represents a selection of the most relevant aspectsconsidered in the development of our methods throughout Section 4 and 5.
Interaction potential and corresponding force.
An interaction potential Φ( r ), also known as (Gibbs) free energy of the interaction, is defined as the amountof energy required to approach the interaction partners starting from a reference configuration with zeroenergy at infinite separation. Hence, the following relations between the interaction potential Φ( r ) and themagnitude of the force f ( r ) acting upon each of the partners, each in terms of the distance between bothinteracting partners r , hold true:Φ( r ) = r (cid:90) ∞ − f ( r ) d r ↔ f ( r ) = − d Φ( r )d r (1)Although the final quantities of interest are often the resulting, vectorial forces on slender bodies, it is yetconvenient and sensible to consider the scalar interaction potential throughout large parts of this work. Thisis underlined by the fact that nonlinear finite element methods in the context of structural dynamics canbe formulated on the basis of energy and work expressions. Equation (1) expresses the direct and inherentrelation between force and potential. Note that the forces emanating from such interaction potentials areconservative and the integral value in (1) is path-independent. Furthermore, the interaction is symmetric4n a sense that the force acting upon the first interacting partner f ( r ) has the same magnitude yet op-posite direction as compared to the force acting on the second partner f ( r ). Using the partners’ positionvectors x , x ∈ R , we can formulate the vectorial equivalent of the formula above: f ( r ) = − d Φ( r )d r x − x (cid:107) x − x (cid:107) and f ( r ) = d Φ( r )d r x − x (cid:107) x − x (cid:107) with r = (cid:107) x − x (cid:107) (2) Disambiguation.
In order to particularize the very general term molecular interactions , we may note that we solely considerinteractions between distinct, solid (macro-)molecules, i. e., no covalent or other chemical bonds, but ratherwhat is sometimes referred to as physical bonds . Thus, we restrict ourselves to intermolecular forces asopposed to intramolecular ones.
First principles describing molecular interactions are formulated for a pair of atoms, molecules or pointcharges. In the following, all types of interactions to be considered in this work are thus first presented fora minimal system consisting of one pair of these elementary interaction partners.
Coulomb’s law is one of the most fundamental laws in physics and describes the interaction of a pair ofpoint charges under static conditions byΦ elstat ( r ) = Q Q πε ε r , (cid:107) f elstat ( r ) (cid:107) = Q Q πε ε r , f elstat , ( r ) = Q Q πε ε x − x (cid:107) x − x (cid:107) (3)where ε and ε are the vacuum and dielectric permittivity, respectively. For the sake of brevity in anylater usage, let us define the abbreviation C elstat := (4 πε ε ) − . Depending on the signs of Q and Q ,electrostatic forces may either be repulsive or attractive. Besides a pair of point charges, Coulomb’s lawlikewise holds for a pair of spherically symmetric charge distributions with resulting charges Q and Q ,respectively. This is an important insight, since ultimately we are interested in interactions between twobodies with finite extension rather than points. Furthermore, interactions between rigid spheres and rigidbodies are of interest for applications such as particle diffusion in hydrogels. Throughout the entire work, noelectrodynamic effects shall be considered. This is a valid assumption as long as bodies are non-conductiveand the motion of bodies carrying the attached charges happens on much larger time scales than relevanteigenfrequencies in electrodynamics.Due to the inverse-first power law, the electrostatic potential has quite a long range, meaning thattwo point charges at a large distance still experience a considerable interaction force as compared to smalldistances. This behavior is even more pronounced for the interaction of two extended bodies, where thewhole lot of all distant point pairs dominates the total interaction energy as compared to the few closestpoint pairs. This property is crucially different as compared e. g. to vdW interactions considered in the nextsection. We will account for and indeed make use of this important property in the development of themethods to be presented in this work. forces originate from charge fluctuations, thus being an electrodynamic effect caused byquantum-mechanical uncertainties in positions and orientations of charges. Depending on the interactionpartners, three subclasses can be distinguished as Keesom (two permanent dipoles), Debye (one permanentdipole, one induced dipole) and London dispersion interactions (two transient dipoles). The ubiquitousnature of van der Waals interactions is due to the fact that the latter contribution even arises in neutral,nonpolar, yet polarizable matter that means basically every atom or molecule. All three kinds of dipole5nteractions can be unified in that their interaction free energy follows an inverse-sixth power law in theseparation [3]: Φ vdW ( r ) = − C vdW r (4)This is a pleasantly simple expression, yet intricate when it comes to transferring it to two-body interactions,as we will discuss in Section 2.3.2. In general, van der Waals forces are always attractive for two identicalor similar molecules, yet may be repulsive for other material combinations. Two approaching atoms or molecules will at some very small separation suddenly experience a seeminglyinfinite repulsive force. This effect is attributed to the overlap of electron clouds and referred to as stericrepulsion , steric exclusion or hard core repulsion . Without thorough theoretical foundation, several (almost)infinitely steep repulsive potential laws are empirically used to model this phenomenon. The first option isa hard wall/core/sphere potential which has a singularity at zero separationlim r → Φ c,hs ( r ) = ∞ and Φ c,hs ( r ) ≡ r > . (5)Other common choices include a power-law potential with a large integer exponent n c,pow Φ c,pow ( r ) = C c,pow r − n c,pow (6)and finally an exponential potential Φ c,exp ( r ) = C c,exp e − r/r c,exp . (7)Note that the former two coincide in the limit n c,pow → ∞ and increase indefinitely for r → In many systems of interest, any two or more of the aforementioned effects may be relevant at the sametime such that a combination of the pair potentials is required. This is typically done by summation ofthe individual potential contributions and leads to a total intermolecular pair potential . Among the largenumber of possible combinations , the Lennard-Jones (LJ) potential is probably the most commonly usedvariant (see Figure 1).Φ LJ ( r ) = k r − + k r − = − Φ LJ,eq (cid:18)(cid:16) r LJ,eq r (cid:17) − (cid:16) r LJ,eq r (cid:17) (cid:19) (8)It is a special case of the Mie potential Φ Mie ( r ) = C Mie,m r − m − C Mie,n r − n with exponents being chosen tomodel the inverse-sixth van der Waals attraction on the one hand and a strong repulsion on the other hand.The parameters can be identified as the minimal value Φ LJ,eq < r LJ,eq >
0, i. e., at the separation where the resulting force is zero.Other important quantities characterizing the LJ force law f LJ ( r ) = −
12 Φ
LJ,eq r LJ,eq (cid:18)(cid:16) r LJ,eq r (cid:17) − (cid:16) r LJ,eq r (cid:17) (cid:19) (9)are the minimal force value f LJ,min and corresponding distance r LJ,f min f LJ,min ≈ . LJ,eq r LJ,eq and r LJ,f min = (cid:18) (cid:19) / r LJ,eq ≈ . r LJ,eq . (10) See [2][p. 138] for a comprehensive list of combinations used as total pair potential laws. igure 1: Lennard-Jones interaction potential for a pair of points, i. e. atoms. The minimal force, i. e., the maximal adhesive force, is commonly referred to as pull-off force . Israelachvili [2]also points out the chance of a fortunate cancellation of errors in total pair potentials, especially close tothe limit r →
0. In this regime, attractive forces tend to be underestimated by the simplified inverse-sixthterm but likewise the steric repulsion is probably stronger than estimated from the power law. Both errorscancel rather than accumulate and increase the model accuracy.
Remarks.
1. Many of the presented point-point interaction potentials decay rapidly with the distance as shownexemplarily for a law Φ( r ) ∝ r − in Figure 2. In anticipation of the numerical methods to be proposedin this work we can already state at this point that these extreme gradients are very challenging fornumerical quadrature schemes that will therefore be discussed in a dedicated Section 6.4.2. In molecular dynamics, a force field is typically used instead of the potential law to model the totalinteraction of a pair of atoms. Specific forms are being proposed for (coarse-grained) force fieldsmodeling the interaction of macromolecules such as DNA instead of atoms. Since these all-atomapproaches follow an entirely different approach as compared to the continuum model proposed here,we will not discuss force fields any further at this point. (a) (b)Figure 2: Example of a point-point interaction potential Φ( r ) ∝ r − plotted over a circular domain (blue circle) with (a)small and (b) large distance to the point-like interaction partner (red dot). Note the huge difference in scales. In this section, we take the important step from interacting point pairs to interactions between twobodies with defined spatial extension containing many of the fundamental point-like interaction partners7onsidered throughout the preceding Section 2.2. The obvious question for the spatial distribution of theinteraction partners leads to the important distinction between surface and volume interactions . As thename suggests, in the first case, the elementary interaction partners are distributed over the surface of thebodies but not in the interior. The most important example from this category are electrostatic interactionsbetween bodies where the charges sit on the surfaces and are not free to move around. This applies to a largenumber of charged, non-conductive biopolymer fibers such as actin or DNA. In the second case of volumeinteractions, the elementary interacting partners are distributed over the entire volume of the bodies. Themost important examples here are van der Waals interactions and steric exclusion. As compared to surfaceinteractions, this further increases the dimension of the problem making it more challenging to tackle, bothby analytical as well as numerical means. In terms of notation one may also find the expressions body forces or bulk interaction referring to this category of interactions.Let us briefly look at volume and surface interactions as an abstract concept, leaving aside the specificsof the underlying physical effects that are to be discussed in the subsequent Section 2.3.1 and Section 2.3.2.Likewise, we assume additivity here and discuss the applicability later with the physical type of interaction.Since volume interactions are the more general and challenging case, we will discuss most aspects and ap-proaches first for volume interactions and later only point out the differences considering surface interactionsthroughout this article. Figure 3 schematically visualizes the distribution of elementary interaction partnerswithin two macromolecular or macroscopic bodies. B : ρ , V B : ρ , V Figure 3: Two arbitrarily shaped, deformable bodies B and B with volumes V , V and continuous particle densities ρ , ρ . Assuming additivity, we apply pairwise summation to arrive at the two-body interaction potentialΠ ia = (cid:88) i ∈B (cid:88) j ∈B Φ( r ij ) . (11)Further assuming a continuous atomic density ρ i , i = 1 ,
2, the total interaction potential can alternativelybe rewritten as integral over the volumes V , V of both bodies B and B :Π ia = (cid:90) (cid:90) V ,V ρ ( x ) ρ ( x )Φ( r ) d V d V with r = (cid:107) x − x (cid:107) (12)It can be shown that this continuum approach is the result of coarse-graining , i. e., smearing-out the discretepositions of atoms in a system into a smooth atomic density function ρ ( x ) [23].In the case of surface interactions, the dimensionality of the problem reduces and summation or integra-tion is carried out over both bodies’ surfaces ∂ B , ∂ B :Π ia = (cid:90) (cid:90) S ,S σ ( x ) σ ( x )Φ( r ) d S d S with r = (cid:107) x − x (cid:107) (13)Accordingly, surface densities σ i ( x i ), i = 1 ,
2, replace the volume densities in this case.8 he range of two-body interaction forces originating from point pair potentials. Let us assume a general inverse power law Φ( r ) = kr − m for the point pair interaction potential. It is obvious,that the potential becomes infinitely large if the separation of the two individual points r approaches zeroand, on the other hand, that the potential rapidly decays with increasing distance. Turning to two bodiesof finite size, i. e., two clouds of points, things are more involved as the following theoretical considerationsdemonstrate. In short, it can be shown that there is a fundamental difference between potentials with anexponent m ≤ m > m >
3, e. g.,vdW interactions, the two-body interaction potential goes to infinity if the bodies approach until theirsurfaces touch each other. This can be illustrated by the simple example of two spheres of radius R , wherethe vdW interaction potential scales with Π vdW ∝ g − (cf. [2, p.255]) with surface-to-surface separationor gap g = d − R and the distance between the spheres’ centers d . This singularity of the two-bodyinteraction potential in the limit of zero separation g → m > m > m ≤
3, e. g., Coulombic interactions. Here,the total contribution of all distant point pairs dominates over the few closest point pairs and the totalinteraction potential remains finite even if both bodies are in contact. Looking once again at the simpleexample of two spheres, Coulomb’s law (cf. eq. (3)) directly shows Π elstat ∝ d − and thus no singularityoccurs for (nearly) contacting surfaces g →
0, i. e., d → R . Also, in contrast to the case of vdW mentionedabove, the Coulomb interaction potential of two perpendicular cylinders would increase if their length isincreased. The underlying theoretical derivations revealing also the turning point, i. e., the exponent m = 3were first noted by Newton and can be found e. g. in [2, p. 11]. Due to this crucial difference, potentialswith m > short-range interactions (e. g. repulsive as well as attractive part of LJ) andpotentials with m ≤ long-range interactions (e. g. Coulomb) throughout this work. The Coulomb interaction is additive such that the net force acting on an individual point charge in asystem of point charges can be calculated from superposition of all pair-wise computed force contributions[2]. Equivalently, the net interaction potential results from summation of all pair potentials. A large body ofliterature deals with the problem of electrostatic multi-body interaction. One concept of high relevance forthe present work is a well-known strategy called multipole expansion, which aims to express the resultantelectrostatic potential of a (continuous) charge distribution as an (infinite) series (see e. g. [31] for details).The individual terms of the series expansion generally are inverse power laws in the distance with increasingexponent and referred to as mono-, di-, quadru-, up to n -pole moments. At points far from the location of thecharge cloud, the series converges quickly and can thus be truncated in good approximation. Regarding thetotal interaction potential of two charged bodies as formulated in (12) or (13), this already outlines the wayhow to determine Π elstat for trivial geometries of the interacting bodies, where the integrals can potentiallybe solved analytically. We will return to this concept in the context of deformable slender fibers whenproposing the general SSIP approach in the beginning of Section 4 and make use of a (truncated) multi-poleexpansion of the charged cross-sections for the (simplified) SSIP law for long-range surface interactions tobe proposed in Section 5.3. Here, we want to discuss vdW interactions as one example of physically relevant volume interactionsthat is based on the inverse-sixth power law (4). However, very similar considerations and formulae applyto steric interactions as well as LJ interactions. 9oday, we know that vdW interactions are generally non-additive. The latest and most accurate modelsfor two-body vdW interactions are based on Lifshitz theory and, among other effects, include retardation,anisotropy and differences in polarizability. Nevertheless, a “happy convergence” of old, i. e., Hamaker’spairwise summation, and new, i. e. Lifshitz, theory allows to determine the distance dependency from pairwisesummation and then estimate the prefactor, i. e., Hamaker “constant” A Ham , from more advanced moderntheory. This approach yielding a so-called Hamaker-Lifshitz hybrid form [3], [2, p. 257] is what motivates usto use pairwise summation in the derivation of the numerical methods to be proposed in the present work.Also, there are some special scenarios (negligible retardation, negligible difference in optical properties of thebodies, interaction in vacuum, ...), where additivity can be assumed as a good approximation even withoutadaption of the Hamaker constant.Generally, even the simple approach of pairwise summation requires two nested 3D integrals over bothbodies’ volumes, i. e., six-dimensional integration. Mainly due to this high dimensionality of the problem,unfortunately, (closed-form) analytical solutions can only be obtained for some simple special cases. Still,careful considerations and selection allow us to exploit some of these analytical expressions in order todevelop efficient, reduced methods in Section 4.2. To get a concise overview of all expressions relevant forthe remainder of this work, we will provide a collection of closed-form analytical solutions in the following.First, we want to look at two cylinders representing the simplest model for two interacting straight, rigidfibers with circular cross-section. A number of publications consider this scenario and due to the simplicity ofthe geometry they were able to derive analytical solutions for some special cases. The resulting expressionsare summarized in Table 1 and will be used for verification purposes in Section 7.1. A second, highlyrelevant scenario is the one considering two disks. These analytical expressions, summarized in Table 2,will be beneficial, and in fact provide the main ingredient, for the SSIP approach to describe molecularinteractions between deformable fibers modeled as 1D Cosserat continua.
Two Cylinders.
To begin with, we consider the cases of parallel and perpendicular cylinders. Generally, the cylinders areassumed to be infinitely long, such that the boundary effects at its ends may be neglected. As the inter-action potential for parallel cylinders would be infinite, one typically considers a length-specific interactionpotential ˜ π vdW,cyl (cid:107) cyl with dimensions of energy per unit length. This quantity thus describes the interactionof one infinitely long cylinder with a section of unit length of the other infinitely long cylinder. For perpen-dicular orientation (and all other mutual angles apart from α = 0) on the other hand, the total interactionpotential Π vdW,cyl ⊥ cyl remains finite.Even for this simple case of two cylinders, no closed-form analytical solution for the vdW interactionenergy can be found for all mutual angles and all separations. One thus resorts to the consideration ofthe limits of small and large surface-to-surface separations for which the general solution, an infinite series,converges to the expressions presented in the following Table 1.Despite the different dimensions of the quantities for parallel and perpendicular cylinders, we can stillcompare these expressions as becomes clear by the following thought experiment. Considering two “suffi-ciently long” cylinders of length L in parallel orientation, the total interaction potential is well describedby Π vdW,cyl (cid:107) cyl = ˜ π vdW,cyl (cid:107) cyl · L and thus shows the same scaling behavior in the separation as (14) and(15). In addition, (16) and (17) are also a good approximation for the perpendicular orientation of thesecylinders of finite length L since the difference in the distant point pairs is negligible.We would like to point out just a few interesting aspects of these equations. First, it is remarkable howthe expressions differ in the exponent of the power law describing the distance dependency of the potential.This relates to a diverse and highly nonlinear behavior already for this simplest model system of fiber-fiberinteractions composed of two cylinders. Second, the parallel orientation is a very special orientation thatgives rise to the strongest possible adhesive forces between two cylinders and at the same time is the onlystable equilibrium configuration. Third, the distance scaling behavior of two parallel cylinders at smallseparations ˜ π vdW,cyl (cid:107) cyl,ss ∝ g − lies between the fundamental solutions known for two infinite half spaces˜˜ π ∝ g − (double tilde indicates a potential per unit area) and two spheres Π ∝ g − . Note that againmultiplication of these laws by a length or area does not alter the scaling law in the distance. Lookingat the equations for large separations, we see similar relations, once again with a stronger distance scaling10imit of small separations Limit of large separations g (cid:28) R , R g, d (cid:29) R , R parallel ˜ π vdW,cyl (cid:107) cyl,ss = − A Ham (cid:113) R R R + R g − (14) ˜ π vdW,cyl (cid:107) cyl,ls = − π A Ham R R d − (15) (cid:104) energylength (cid:105) see [2][p. 255],[3][p. 172] see [3][p. 16, p. 172]perp. Π vdW,cyl ⊥ cyl,ss = − A Ham √ R R g − (16) Π vdW,cyl ⊥ cyl,ls = − π A Ham R R d − (17)[energy] see [2][p. 255] see [3][p. 16] Table 1: A collection of analytical solutions for the cylinder-cylinder interaction potential derived via pairwise summation.Here, R i denotes the cylinder radii, d denotes the closest distance between the cylinder axes, g denotes the surface-to-surfaceseparation, i. e., gap, A Ham := π ρ ρ C vdW is the commonly used abbreviation known as Hamaker constant where ρ i denotesthe particle densities and C vdW denotes the constant prefactor of the point-pair potential law (see eq. (4)). behavior in the parallel case. Generally, the solutions for large separations are expressed more naturally inthe inter-axis separation d rather than the surface-to-surface separation g . Two Disks.
This problem has been studied in literature on vdW interaction of straight, rigid cylinders of infinite length[32]. In analogy to the cylinder-cylinder scenario, it turns out that not even in the simplest case of paralleloriented disks, i. e., two disks with parallel normal vectors, a closed analytical solution can be found for allseparations. Instead, two expressions for the limit of small and large separations of the disks g as comparedto their radii R , R are presented in the following Table 2, respectively.Limit of small separations Limit of large separations g (cid:28) R , R g, d (cid:29) R , R parallel ˜˜ π vdW,disk (cid:107) disk,ss = − A Ham (cid:113) R R R + R g − (18) ˜˜ π vdW,disk (cid:107) disk,ls = − A Ham R R d − (19) (cid:104) energylength (cid:105) see [32] see [32] Table 2: A collection of analytical solutions for the disk-disk interaction potential derived via pairwise summation. Here, R i denotes the disk radii, d denotes the closest distance between the disk midpoints, g denotes the surface-to-surface separation,i. e., gap, A Ham := π ρ ρ C vdW is the commonly used abbreviation known as Hamaker constant where ρ i denotes the particledensities and C vdW denotes the constant prefactor of the point-pair potential law (see eq. (4)). To summarize, a closed-form expression for the two-body vdW interaction potential is only known forsome rare special cases and the ones relevant for fiber-fiber interactions have been identified in the voluminousliterature on this topic and presented here in a brief and concise manner.We would like to conclude this section on two-body vdW interactions with a note on the analogy tosteric exclusion, i. e., contact interactions, as already discussed for point pairs in Section 2.2.3. This class ofphysical interactions shares the two central properties of being extremely short in range and being volumeinteractions. Starting from an inverse-twelve power law as in the repulsive part of the LJ interaction law, onemay apply very similar solution strategies and finally obtain very similar expressions as the ones presented11n this section. For the sake of brevity, we refer to the derivations in Appendix A and the analysis of theresulting total LJ interaction that will also be used for the regularization of the reduced potential laws inSection 6.3.
3. Fundamentals of Geometrically Exact 3D Beam Theory
This section aims to provide a brief and concise introduction to well-known concepts of beam theoryto be used in the remainder of this article. As commonly used in engineering mechanics, we refer tothe term beam as a mathematical model for a three-dimensional, slender, deformable body for which thefollowing assumption can be made. The much larger extent of the body in its axial direction as comparedto all transverse directions often justifies the Bernoulli hypothesis of rigid and therefore undeformable cross-sections. This in turn allows for a reduced dimensional description as a 1D Cosserat continuum embeddedin the 3D Euclidian space.The so-called Simo-Reissner beam theory dates back to the works of Reissner [33], Simo [34], and Simoand Vu-Quoc [35], who generalized the linear Timoshenko beam theory [36] to the geometrically nonlinearregime. Since the Simo-Reissner model, which accounts for the deformation modes of axial tension, bending,torsion and shear deformation, is the most general representative of geometrically exact beam theories, wechoose it as the one to be used exemplarily throughout this work. Nevertheless, the novel approach to beproposed is not restricted to a specific beam formulation. We have likewise applied it to formulations ofKirchhoff-Love type, which are known to be advantageous in the regime of high slenderness ratios where theunderlying assumption of negligible shear deformation is met [6, 37]. Refer to Section 6.1 for more details.
Geometry representation.
A certain configuration of the 1D Cosserat continuum is uniquely defined by the centroid position and theorientation of the cross-section at every point of the continuum. The set of all centroid positions is referredto as centerline or neutral axis and expressed by the curve s, t (cid:55)→ r ( s, t ) ∈ R (20)in space and time t ∈ R . Each material point along the centerline is represented by a corresponding valueof the arc-length parameter s ∈ [0 , l ] =: Ω l ⊂ R . Note that this arc-length parameter s is defined inthe stress-free, initial configuration of the centerline curve r ( s ) = r ( s, t = 0). Thus, the norm of the initialcenterline tangent vector yields (cid:107) r (cid:48) ( s ) (cid:107) := (cid:107) ∂ r ( s ) ∂ s (cid:107) ≡ , (21)but generally, in presence of axial tension, (cid:107) r (cid:48) ( s, t ) (cid:107) = (cid:107) ∂ r ( s, t ) /∂ s (cid:107) (cid:54) = 1.Furthermore, the cross-section orientation at each of these material points is expressed by a right-handedorthonormal frame often denoted as material triad : s, t (cid:55)→ Λ ( s, t ) := [ g ( s, t ) , g ( s, t ) , g ( s, t )] ∈ SO (22)The second and third base vector follow those material fibers representing the principal axes of inertia ofarea. Such a triad can equivalently be interpreted as rotation tensor transforming the base vectors of aglobal Cartesian frame E i ∈ R , i ∈ , , g i ∈ R , i ∈ , , g i ( s, t ) = Λ ( s, t ) E i ( s, t ) . (23)In summary, a beam’s configuration may be uniquely described by a field of centroid positions r ( s, t ) and afield of associated material triads Λ ( s, t ), altogether constituting a 1D Cosserat continuum (see Figure 4).According to this concept of geometry representation, the position x of an arbitrary material point P of theslender body is obtained from x P ( s, s , s , t ) = r ( s, t ) + s g ( s, t ) + s g ( s, t ) . (24)12 E E r (cid:48) ( s ) r ( s ) → r ( s ) Λ ( s ) → Λ ( s ) r ( s ) r ( s ) s s g g g g g g Figure 4: Geometry description and kinematics of the Cosserat continuum formulation of a beam: Initial, i. e., stress-free (blue)and deformed (black) configuration. Straight configuration in initial state is chosen exemplarily here without loss of generality.
Here, the additional convective coordinates s and s specify the location of P within the cross-section, i. e.,as linear combination of the unit direction vectors g and g . For a minimal parameterization of the triad,e. g. the three-component rotation pseudo-vector ψ may be used such that we end up with six independentdegrees of freedom ( r , ψ ) to define the position of each material point in the body by means of (24). Remark on notation.
Unless otherwise specified, all vector and matrix quantities are expressed in the globalCartesian basis E i . Differing bases as e. g. the material frame are indicated by a subscript [ . ] g i . Quantitiesevaluated at time t = 0, i. e., the initial stress-free configuration, are indicated by a subscript 0 as e. g. in r ( s ).Differentiation with respect to the arc-length coordinate s is indicated by a prime, e. g., for the centerlinetangent vector r (cid:48) ( s, t ) = d r ( s, t ) / d s . Differentiation with respect to time t is indicated by a dot, e. g., forthe centerline velocity vector ˙ r ( s, t ) = d r ( s, t ) / d t . For the sake of brevity, the arguments s, t will often beomitted in the following. Remark on finite 3D rotations.
To a large extent, the challenges and complexity in the numerical treatmentof the geometrically exact beam theory can be traced back to the presence of large rotations. In contrastto common vector spaces , the rotation group SO (3) is a nonlinear manifold (with Lie group structure) andlacks essential properties such as additivity and commutativity, which makes standard procedures such asinterpolation or update of configurations much more involved. While Section 4 introduces the concept ofsection-to-section interaction laws in the most general manner, in Section 5, some additional (practicallyrelevant) assumptions are made that allow to formulate the interaction laws as pure function of the beamcenterline configuration. In turn, this strategy will allow to avoid the handling of finite rotations and toachieve simpler and more compact numerical formulations. Kinematics, deformation measures and potential energy of the internal, elastic forces.
Figure 4 summarizes the kinematics of geometrically exact beam theory. Based on these kinematic quantities,deformation measures as well as constitutive laws can be defined. Finally, the potential 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. [4, 5, 6] for a detailed presentation of these steps.13 . The Section-to-Section Interaction Potential (SSIP) Approach Based on the fundamentals of molecular interactions (Section 2) as well as geometrically exact beamtheory (Section 3), this section will propose the novel SSIP approach to model various types of molecularinteractions between deformable fibers undergoing large deflections in 3D.
For a classical conservative system, the total potential energy of the system can be stated taking intoaccount the internal and external energy Π int and Π ext . The additional contribution from molecular inter-action potentials Π ia is simply added to the total potential energy as follows.Π TPE = Π int − Π ext + Π ia ! = min. (25)Note that the existing parts remain unchanged from the additional contribution. One noteworthy differenceis that internal and external energy are summed over all bodies in the system whereas the total interactionfree energy is summed over all pairs of interacting bodies.According to the principle of minimum of total potential energy , the weak form of the equilibriumequations is derived by means of variational calculus. The very same equation may alternatively be derivedby means of the principle of virtual work which also holds for non-conservative systems: δ Π int − δ Π ext + δ Π ia = 0 (26)Clearly, the evaluation of the interaction potential Π ia , or rather its variation δ Π ia , is the crucial stephere. Recall (12) to realize that it generally requires the evaluation of two nested 3D integrals . The directapproach using 6D numerical quadrature turns out to be extremely costly and in fact inhibits any applicationto (biologically) relevant multi-body systems. See Section 6.4 for more details on the complexity and thecost of this naive, direct approach as well as the novel SSIP approach to be proposed in the following. We propose a split of the integral in the length dimensions l , l on the one hand and the cross-sectionaldimensions A , A on the other hand:Π ia = (cid:90) (cid:90) l ,l (cid:90) (cid:90) A ,A ρ ( x ) ρ ( x )Φ( r ) d A d A (cid:124) (cid:123)(cid:122) (cid:125) =: ˜˜ π ( r − , ψ − ) d s d s with r = (cid:107) x − x (cid:107) . (27)Exploiting the characteristic slenderness of beams, the 4D integration over both undeformable cross-sectionsshall be tackled analytically and only the remaining two nested 1D integrals along the centerline curves shallbe evaluated numerically to allow for arbitrarily deformed configurations. Generally speaking, we followthe key idea of reduced dimensionality from beam theory and thus aim to express the relevant informationabout the cross-sectional dimensions by the point-wise six degrees of freedom ( r i , ψ i ) of the 1D Cosseratcontinua without loss of significant information. To this end we need to consider the resulting interac-tion between all the elementary interaction partners within two cross-sections expressed by an interactionpotential ˜˜ π ( r − , ψ − ) that depends on the separation r − of the centroid positions and the relative rota-tion ψ − between both material frames attached to the cross-sections. For this reason, the novel approachis referred to as the section-to-section interaction potential (SSIP) approach. The SSIP ˜˜ π is a double length-specific quantity in the way that it measures an energy per unit length of beam 1 per unit length of beam 2,which is indicated here by the double tilde. The sought-after total interaction potential Π ia of two slender It is important to mention that, assuming additivity of the involved potentials, systems with more than two bodies can behandled by superposition of all pair-wise two-body interaction potentials. It is thus sufficient to consider one pair of beams inthe following. The same reasoning applies to more than one type of physical interaction, i. e., potential contribution. ia = l (cid:90) l (cid:90) ˜˜ π ( r − , ψ − ) d s d s (28)This relation suggests another, alternative interpretation of the SSIP ˜˜ π . In analogy to the term inter-surfacepotential , introduced by [22], ˜˜ π can be understood as an inter-axis potential , i. e., describing the interactionof two spatial curves (with attached material frames).To further illustrate this novel concept, a simple, demonstrative example is shown in Figure 5. In this E E E ξ ξ ψ − r − r ( ξ ) r ( ξ ) ξ ξ Figure 5: Illustration of the novel SSIP approach: Two cross-sections at integration points ξ / , GP of beam 1 and 2, respectively,their separation r − and relative rotation ψ − . scenario of two beams with circular cross-section, the SSIP ˜˜ π ( r − , ψ − ) describes the interaction of twocircular disks at arbitrary mutual distance and orientation. To evaluate the two nested 1D integrals along thebeam axes numerically, the SSIP needs to be evaluated for all combinations of integration points (denotedhere as Gaussian quadrature points (GP), without loss of generality). For one of these pairs ( ξ , GP , ξ , GP ),the geometrical quantities are shown exemplarily.While analytical integration of the inner 4D integral of (27) has already been suggested above as oneway to find a closed-form expression for the SSIP ˜˜ π , we would like to stress the generality of the SSIPapproach at this point. The question of how to find ˜˜ π is independent of the strategy to determine theinteraction energy Π ia of two slender bodies via numerical double integration as proposed in this section.This is important to understand because the SSIP ˜˜ π will obviously depend on the type of interaction as wellas the cross-section shape and a number of other factors and there might also be cases where no analyticalsolution can be obtained and one wants to resort to relations fitted to experimental data. In the scope ofthis work, several specific expressions of ˜˜ π , e. g., for vdW as well as electrostatic interactions, will be derivedanalytically in the following Section 5.In its most general form ˜˜ π will be a function of the relative displacement r − and the relative rota-tion ψ − between both cross-sections, i. e., three translational and three rotational degrees of freedom.This becomes clear if one recalls that the position x P of every material point in a slender body can beuniquely described by the six degrees of freedom of a 1D Cosserat continuum (cf. eq. (24)). Thus, keepingone cross-section fixed, the position x P1 − P2 of every material point in the second cross-section relative to the(centroid position and material frame of the) first cross-section is again uniquely described by six degrees offreedom ( r − , ψ − ). This insight naturally leads to the interesting question under which conditions the15SIP ˜˜ π can be described by a smaller set of degrees of freedom, thus simplifying the expressions. Rotationalsymmetry of the interacting cross-sections is one common example where the SSIP would be invariant underrotations around the cross-section’s normal axis. We will return to this topic in Section 5.1 as a preparationfor the following derivation of specific expressions for the SSIP. Remark on the included special case of surface interactions.
It is very convenient that the practically highlyrelevant case of surface potentials is already included as a simpler, special case in the proposed SSIP approachto model molecular interactions between the entire volume of flexible fibers. In simple words, it is sufficientto omit one spatial dimension of analytical integration on each interacting body in the analytical derivationof the required SSIP ˜˜ π . More specifically, this means that ˜˜ π may be obtained from solving analytically twonested 1D integrals along both, e. g. ring-shaped, contour lines of the fiber cross-sections.
5. Application of the General SSIP Approach to Specific Types of Interactions
At this point we would like to return to the fact that the SSIP approach proposed in Section 4 is generalin the sense that it does not depend on the specific type of physical interaction. This section provides thenecessary information and formulae to apply the newly proposed, generally valid approach from Section 4.2to certain types of real-world, physical interactions such as electrostatics or vdW. As mentioned above, theapproach requires a closed-form expression for the SSIP ˜˜ π . We basically see two alternative promising waysto arrive at such a reduced interaction law ˜˜ π :1. analytical integration, e. g., as presented in Section 5.22. postulate ˜˜ π as a general function of separation r − and mutual orientation ψ − and determine thefree parameters via fitting to(a) experimental data for specific section-to-section configurations, i. e., discrete values of separa-tion r − and mutual orientation ψ − (b) data from (one-time) numerical 4D integration for specific section-to-section configurations, i. e.,discrete values of r − and ψ − (c) experimental data for the global system response, e. g., of the entire fiber pair or a fiber networkAs a starting point we will restrict ourselves to the first option based on analytical integration throughoutthe remainder of this work. See Section 5.2 for an example of the further steps required to derive the final,ready-to-use expressions in case of vdW interactions. To give but one example for a recent experimentalwork, which could serve as the basis for the second option listed above, we refer to [38] measuring cohesiveinteractions between a single pair of microtubules. Postulating an SSIP and studying the global systemresponse in numerical simulations could also be used as a verification of theoretical predictions for thesystem behavior. To give an example we refer to the work of theoretical biophysicists studying the structuralpolymorphism of the cytoskeleton resulting from molecular rod-rod interactions [39], which is based on apostulated model potential “that captures the main features of any realistic potential”. In summary, we seea large number of promising future use cases for the proposed SSIP approach. Recall that the most general form of the SSIP is uniquely described by a set of six degrees of freedom,three for the relative displacement and three for the relative orientation of the two interacting cross-sections,as presented in the preceding Section 4.2. The following assumptions turn out to significantly simplify thismost general form of the SSIP law by reducing the number of relevant degrees of freedom from six to four,two or even one. This in turn eases the desirable derivation of analytical closed-form solutions of the SSIP ˜˜ π based on the point pair potentials Φ( r ) presented in Section 2.2. Specifically, these assumptions are:1. undeformable cross-sections2. circular cross-section shapes3. homogeneous (or, more generally, rotationally symmetric) particle densities ρ , ρ in the cross-sectionsor surface charge densities σ , σ over the circumference16he first assumption is typical for geometrically exact beam theory and the second and third assumption arereasonable regarding our applications to biopolymer fibers such as actin or DNA that can often be modeledas homogeneous fibers with circular cross-sections. Based on these three assumptions, we can conclude thatthe interaction between two cross-sections is geometrically equivalent to the interaction of two homogeneous,circular disks (or rings in case of surface interactions). The rotational symmetry of the circular disks thenimplies that the interaction potential is invariant under rotations around their own axes and thus reducesthe number of degrees of freedom to four. The relative importance of the remaining degrees of freedom, i. e.,modes, will be the crucial point in the following discussion, where we turn to the interaction of two slenderbodies, i. e., consider the entirety of all cross-section pairs. At this point, recall the fundamental distinctionbetween either short-range or long-range interactions as outlined in Section 2.3.2. β β γ γ α r (cid:48) ( s ) r (cid:48) ( s ) r − (a) Ω (b)Figure 6: Sketches to illustrate the simplifications resulting for (a) short- and (b) long-range interactions. In the case of short-range interactions, the cross-section pairs in the immediate vicinity of the mutualclosest points of the slender bodies dominate the total interaction. As is known from macroscopic beamcontact formulations [12, 20, 21], the criterion for the closest point is that the distance vector r − isperpendicular to both centerline tangent vectors r (cid:48) i , i. e., (assuming small shear angles) the normal vectorsof the disks (see Figure 6(a)). Since only cross-section pairs in the direct vicinity of the closest pointsare relevant, arbitrary relative configurations (i. e. separations and relative rotations) between those cross-sections shall in the following be discussed on the basis of six alternative degrees of freedom as illustratedin Figure 6(a). By considering the cross-sections A and A at the closest points as reference, the relativeconfiguration between cross-sections in the direct vicinity of A and A can be described via (small)rotations of A around the axis r (cid:48) (angle γ ) and r − × r (cid:48) (angle β ), (small) rotations of A around theaxes r (cid:48) (angle γ ) and r − × r (cid:48) (angle β ), (small) relative rotations between A and A around the axis r − (angle α ), and (small) changes in the (scalar) distance d = (cid:107) r − (cid:107) . As a consequence of assumptions1–3 discussed above, the considered interaction potentials are invariant under rotations γ and γ . Fromthe remaining four degrees of freedom, the scalar distance d clearly has the most significant influence onthe interaction potential because changes in d directly affect the mutual distance r of all point pairs inthe body and, most importantly, the smallest surface separation g between both bodies. The second mostsignificant influence is expected for the scalar relative rotation angle α between the cross-section normalvectors, i. e., cos( α ) = n · n ≈ r (cid:48) · r (cid:48) / ( (cid:107) r (cid:48) (cid:107) (cid:107) r (cid:48) (cid:107) ). A change in α does not alter the gap g , but influencesthe distance of all next nearest point pairs in the immediate vicinity of the closest surface point pair. For theremaining two relative rotations β and β , arguments for both sides, either significant or rather irrelevantinfluence on the total interaction potential, can be found at this point. On the one hand, the orthogonalityconditions r (cid:48) i · r − = 0, i = 1 , β / could be considered negligible.On the other hand, even small rotations β / change the smallest separation of any two point pairs in theimmediate neighborhood of the closest point pair as soon as the centroid distance vector r − rotates outof the two cross-section planes. Therefore, it seems hard to draw a final conclusion with respect to theinfluence and thus importance of β / based on the qualitative theoretical considerations of this section.17o summarize, the scalar distance d between the cross-section centroids, the scalar relative rotationangle α between the cross-section normal vectors and possibly also the relative rotation components β / are supposed to have a perceptible influence on the short-range interaction between slender beams fulfillingassumptions 1–3, with a relative importance which decreases in this order. In favor of the simplest possiblemodel, we will therefore assume at this point that the effect of this scalar relative rotations α, β and β is negligible as compared to the effect of the scalar separation d . This allows us to directly use theanalytical, closed-form expression for the disk-disk interaction potential as presented in Section 2.3.2. Theerror for arbitrary configurations associated with this model assumption will be thoroughly analyzed inSection 7.1.1. In this context, it is a noteworthy fact, that the first published method for 2D beam-to-rigidhalf space LJ interaction [24] likewise neglects the effect of cross-section orientation. In the subsequentpublication [27], the effect of cross-section rotation, i. e., interaction moments, has finally been included anda quantitative analysis considering a peeling experiment of a Gecko spatula revealed that the differences inthe resulting maximum peeling force and bending moment are below 8% and 2%, respectively. However, itis unclear whether this assessment also holds for beam-to-beam interactions modeled via the proposed SSIPapproach. Including the orientation of the cross-sections thus is a work in progress and will be addressed ina subsequent publication. Finally, it is emphasized that by the discussed assumptions the SSIP law ˜˜ π as wellas the total two-body interaction potential Π ia can be formulated as pure function of the beam centerlines r and r without the necessity to consider cross-section orientations via rotational degrees of freedom. Thisis considered a significant simplification of the most general case of the SSIP approach and thus facilitatesboth the remaining derivations in the present work as well as potential future applications. Remark on configurations with non-unique closest points.
It is well-known from the literature on macro-scopic beam contact that the location of the closest points is non-unique for certain configurations of twointeracting beams, e. g., the trivial case of two straight beams, where an infinite number of closest pointpairs exists (see e. g. [20]). Note however that the reasoning presented above also holds in these cases, sincethe cross-section pairs in either one or several of these regions will dominate the total interaction potential.In the case of long-range interactions, the situation is fundamentally different. Recall that here the largenumber of cross-section pairs with large separation d (cid:29) R outweighs the contributions from those few pairsin the vicinity of the closest point and dominates the total interaction. Thus, the regime of large separationsis decisive in this case and it has already been shown in the literature considering disk-disk interaction(see the brief summary in Appendix A.1.1) that in this regime the exact orientation of the disks can beneglected as compared to the centroid separation d . In simple terms, this holds because the distance x P1 − P2 between any point in disk 1 and any point in disk 2 may be approximated by the centroid separation d , if d is much larger than the disk radii R i , which - again - holds for the large majority of all possible cross-sectionpairs. The validity of this assumption will be thoroughly verified by means of numerical reference solutionsin Section 7.1.2. Remark.
The following, similar reasoning from the perspective of slender continua comes to the same con-clusion. As visualized in Figure 6(b), even pure (rigid body) rotations of slender bodies always entail largedisplacements of the centerline in the region far away from the center of rotation. The displacement of anymaterial point due to cross-section rotation will be in the order of Ω R , where Ω is the angle of rotation and R denotes the cross-section radius, whereas the displacement due to centerline displacement will be in the orderof Ω L , where L is the distance from the center of rotation and thus in the order of the beam length l . Dueto the high slenderness l/R (cid:29) L (cid:29) R , which is the decisive one here, because it includes the largemajority of all possible cross-section pairs, as outlined above. The original, analogous reasoning has been ap-plied to the relative importance of translational versus rotational contributions to the mass inertia of beams. Disregarding rotations around its own axis, which are irrelevant here due to rotational symmetry, as mentioned above.
18o conclude, we have discussed the possibility of defining and using SSIP laws ˜˜ π as a function of thescalar separation of the centroids d instead of the six degrees of freedom in the most general form. Thissignificantly simplifies the theory because the analytical solutions for the planar disk-disk interaction fromliterature can directly be used and the complex treatment of large rotations is avoided. Having consideredthe additional assumptions above in the context of short-range interactions, the relative importance of cross-section rotations still needs to be verified in the subsequent quantitative analysis of Section 7.1.1. In thecase of long-range interactions between slender bodies, we have argued that the application of such simpleSSIP laws ˜˜ π ( d ) is expected to be a good approximation which will be confirmed by the quantitative analysisof Section 7.1.2. In the following, a generic short-range volume interaction described by the point-pair potential lawΦ m ( r ) = k m r − m , m > m = 6 (cf. eq. (4)) as well as stericrepulsion as modeled by LJ for exponent m = 12 (cf. eq. (8)). As outlined already in the preceding section,only the regime of small separations is practically relevant in this case of short-range interactions and weneglect the effect of cross-section rotations throughout this article. At this point, we can thus return to theresults for the disk-disk scenario obtained in literature on vdW interactions [32] and summarized in Table 2.In particular, we make use of expression (18) or rather the more general form (A.12). The latter is valid forall power-law point pair interaction potentials with a general exponent m > /
2, i. e., all interactions wherethe strength decays “fast enough”.First, let us introduce the following abbreviation containing all constants in the lengthy expression: c m,ss := k m ρ ρ π ( m − (cid:114) R R R + R Γ( m − )Γ( m − )Γ( m − m −
1) (30)Using eq. (A.12) in combination with the general SSIP approach (28) from Section 4.2, we directly obtain anexpression for the total interaction potential of two deformable fibers in the case of short-range interactions:Π m,ss = (cid:90) l (cid:90) l c m,ss g − m + 72 d s d s for m > (31)with g ( s , s ) = (cid:107) r ( s ) − r ( s ) (cid:107) − R − R (32)Here, the so-called gap g is the (scalar) surface-to-surface separation, i. e., the beams’ centerline curves r ( s ) and r ( s ) minus the two radii R i , as visualized in Figure 5. In general, the particle densities ρ / may depend on the curve parameters s / , i. e., vary along the fiber, without introducing any additionalcomplexity at this point. For the sake of brevity, these arguments s / will be omitted in the remainder ofthis section.The variation of the interaction potential required to solve eq. (26) finally reads δ Π m,ss = ( − m + ) (cid:90) l (cid:90) l c m,ss (cid:0) δ r T − δ r T (cid:1) r − r d g − m + 52 d s d s for m > . (33)Here, we used the variation of the gap δg , which is a well-known expression from the literature on macroscopicbeam contact [12] and is identical to the variation of the separation of the beams’ centerlines δd to be usedin (36), because the cross-sections are assumed to be undeformable: δg = δd = (cid:0) δ r T − δ r T (cid:1) r − r d (34)Solving eq. (26) generally requires two further steps of discretization and subsequent linearization of thisadditional contribution δ Π m,ss to the total virtual work. The resulting expressions will be presented inSection 6.1.1 and Appendix B.1, respectively. As discussed along with the general SSIP approach inSection 4.2, the remaining two nested 1D integrals are evaluated numerically, e. g., by means of Gaussianquadrature. See Section 6.4 for details on this algorithmic aspect.19 emark on the regularization of the integrand. The inverse power law in the integrand of eq. (33) has asingularity for the limit of zero surface-to-surface separation g →
0. Consequently, a so-called regularization of the potential law is needed to numerically handle (the integration of) this term robustly as well as suf-ficiently accurate. This approach is well-known e. g. from (beam) contact mechanics (see e. g. [40, 25, 20])and will be further discussed and elaborated in Section 6.3.At the end of this section, we can conclude that we found specific, ready-to-use expressions for theinteraction free energy as well as virtual work of generic short-range interactions described via the SSIPapproach. Thus, vdW interaction or steric exclusion of slender, deformable continua can now be modeledin an efficient manner, reducing the numerical integral to be evaluated from six to two dimensions. Adetailed quantitative study of the approximation quality with regard to the assumptions discussed in thepreceding Section 5.1 is content of Section 7.1.1.
Having discussed short-range volume interactions, we now want to consider one example of long-rangesurface potentials . Since electrostatic interaction is the prime example of surface potential interaction and atthe same time of high interest for the application to biopolymers we have in mind, we will focus on this casethroughout the following section and mostly speak of point charges as the elementary interaction partners.However, the required steps and formulae will be presented as general as possible in order to allow for asmooth future transfer to other applications.Especially in this context, it is important to stress again that within this model the elementary interactionpartners, i. e., charges must not redistribute within the bodies. Hence, only non-conducting materials can bemodeled with the SSIP approach. This however covers our main purpose to model electrostatic interactionsbetween bio-macromolecules such as protein filaments and DNA because charges are not free to move therein.According to the SSIP approach proposed in Section 4.2, we aim to use analytical expressions for the twoinner integrals over the cross-section circumferences, while the integration along the two beam centerlineswill be evaluated numerically (cf. eq. (27) in combination with the remark on surface interactions at the endof the corresponding section). As discussed in Section 5.1, the regime of large separations is the decisiveone for beam-to-beam interactions in this case of long-range interactions and the SSIP law can be simplifiedin good approximation to depend only on the centroid separation d , which will be confirmed numericallyin Section 7.1.2.At this point, we again return to the expressions for the disk-disk interaction based on a generic point pairpotential Φ m ( r ) = k r − m , as derived in the literature on vdW interactions [32] and summarized in AppendixA.1.1. In particular, the relation (A.6) will be used, which is the same approximation used to derive eq. (19)that describes the practically rather irrelevant scenario of short-range vdW interactions in the regime of largeseparations. Note that in the context of electrostatics, this result is well-known as the first term, i. e., zerothpole or monopole of the multipole expansion of the ring-shaped charge distribution on each of the disks’circumference, which represents the effect of the net charge of a (continuous) charge distribution and has noangular dependence (see Section 2.3.1). In simple terms, this monopole-monopole interaction means thatthe point pair interaction potential Φ( r ) is evaluated only once for the distance between the centers of thedistributions r = d = (cid:107) r − r (cid:107) and weighted with the number of all point charges on the two circumferencesof the circular cross-sections. The expression for the SSIP law to be used throughout this work would thusbe exact for the scenario of the net charge of each cross-section concentrated at the centroid position (ordistributed spherically symmetric around the centroid position). If the accuracy of the SSIP approach needsto be improved beyond the level resulting from this simplified SSIP law (see Section 7.1.2 for the analysis),one could simply include more terms from the multipole expansion of the (ring-shaped) charge distributionsto the SSIP law, which would take the relative rotation of the cross-sections into account. Throughout thiswork and for the applications we have in mind, the simplified SSIP law, which is based on the monopole-monopole interaction of cross-sections, turns out to be an excellent approximation for the true electrostaticinteraction law and we thus restrict ourselves to this variant.20wo nested 1D integrals over the beams’ length dimensions then yield the two-body interaction potentialfor two fibers with arbitrary centerline shapesΠ ia,ls = (cid:90) l (cid:90) l πR σ πR σ Φ( r = d ) d s d s with d ( s , s ) = (cid:107) r ( s ) − r ( s ) (cid:107) . (35)The surface (charge) densities σ j , j = 1 , λ j = 2 πR j σ j , j = 1 ,
2, and iscommonly referred to as linear charge density . Note however that eq. (35) holds for all long-range point pairpotential laws Φ( r ), e. g., all power laws Φ m ( r ) = k r − m with m ≤
3. In order to obtain the weak form ofthe continuous problem, the variation of this total interaction energy needs to be derived. This variationalform can immediately be stated as δ Π ia,ls = (cid:90) l (cid:90) l λ λ ∂ Φ( r = d ) ∂ d δd d s d s with δd = (cid:0) δ r T − δ r T (cid:1) r − r d (36)as the consistent variation of the separation of the beams’ centerlines d , which is well-known from macroscopicbeam contact formulations [12]. By inserting the generic (long-range) power lawΦ m ( r ) = k r − m , m ≤ δ Π m,ls = (cid:90) l (cid:90) l km λ λ (cid:124) (cid:123)(cid:122) (cid:125) =: c m,ls (cid:0) − δ r T + δ r T (cid:1) r − r d m +2 d s d s . (38)The specific case of Coulombic surface interactions follows directly for m = 1 and k = C elstat (cf. eq. (3)). Atthis point, we have once again arrived at the sought-after contribution to the weak form (26) of the space-continuous problem. The steps of finite element discretization and linearization will again be presented later,in Section 6.1.2 and Appendix B.2, respectively. Remark on volume interactions.
Note that there is no conceptual difference if long-range volume interactionswere considered instead of the long-range surface interactions presented exemplarily in this section. Theonly difference lies in the constant prefactor c m,ls , which would read c m,ls = kmA A ρ ρ instead. Ratherthan the spatial distribution of the elementary interaction points in the volume or on the surface, it is thelong-ranged nature of the interactions, which is important for the derivations in this section and allows theuse of approximations for large separations (refer to the extensive discussion in Section 5.1). Remark on intra- versus inter-body interactions.
The electrostatic interaction of point charges on the sameslender body may cause unexpected effects. Assuming equal charges along the beam length leads to repul-sive forces which in turn cause tensile axial forces in the beam. At the start of a dynamic simulation, asimply supported beam will undergo axial strain oscillations before eventually an equilibrium state is found.Alternatively, these interactions of charges within the same body may be included in the constitutive modelused for the continuous body and to this end be modeled by an increased effective stiffness as has e. g. beensuggested by [23, 41]. The latter approach has been applied in the numerical examples of Section 7.
6. Finite Element Discretization and Selected Algorithmic Aspects
Having discussed the space-continuous theory in Section 4 and 5, we now turn to the step of spatialdiscretization by means of finite elements. Subsequently, the most important aspects of the required algo-rithmic framework will be presented briefly and discussed specifically in the light of the novel SSIP approach.This includes the applied regularization technique, multi-dimensional numerical integration, an analysis ofthe algorithmic complexity as well as the topics of search for interaction partners and parallel computing.21 .1. Spatial discretization based on beam finite elements
As presented in Section 3, the centerline position r and the triad Λ arise as the two primary fields ofunknowns. Within Simo-Reissner beam theory, both fields are uncorrelated and their discretization canhence be considered independently as follows. The Simo-Reissner finite beam element used throughout thiswork originates from [4, 5], although we apply a different centerline interpolation scheme here. We employcubic Hermite polynomials based on nodal position vectors ˆ d , ˆ d and tangent vectors ˆ t , ˆ t as the primaryvariables. See [42] for a detailed discussion of Hermite centerline interpolation in the context of geometricallyexact (Kirchhoff) beams and [37] for the details on the Hermitian Simo-Reissner beam element that is usedwithin this article. Applying this interpolation scheme results in the following discretized centerline geometryand variation: r ( ξ ) ≈ r h ( ξ ) = (cid:88) i =1 H id ( ξ ) ˆ d i + l (cid:88) i =1 H it ( ξ ) ˆ t i =: H ˆ d ,δ r ( ξ ) ≈ δ r h ( ξ ) = (cid:88) i =1 H id ( ξ ) δ ˆ d i + l (cid:88) i =1 H it ( ξ ) δ ˆ t i =: H δ ˆ d (39)Here, all the degrees of freedom of one element relevant for the centerline interpolation, i. e., nodal posi-tions ˆ d i and tangents ˆ t i , i = 1 , d and H is the accordingly assembled matrixof shape functions, i. e., Hermite polynomials H id and H it . The newly introduced element-local parame-ter ξ ∈ [ −
1; 1] is biuniquely related to the arc-length parameter s ∈ [ s ele,min ; s ele,max ] describing the verysame physical domain of the beam as follows and the scalar factor defining this mapping between bothinfinite length measures is called the element Jacobian J ( ξ ):d s = d s d ξ d ξ =: J ( ξ ) d ξ with J ( ξ ) := (cid:107) d r , h ( ξ )d ξ (cid:107) . (40)Our motivation to use Hermite interpolation is that it ensures C -continuity, i. e., a smooth geometryrepresentation even across element boundaries. This property turned out to be crucial for the robustness ofsimulations in the context of macroscopic beam contact methods [37], and is just as important if we includemolecular interactions as proposed in this article. See [43] for a comprehensive discussion of (non-)smoothgeometries and adhesive, molecular interactions using 2D solid elements. Note however that neither the SSIPapproach proposed in Section 4 nor the specific expressions for the interaction free energy and the virtualwork are limited to this Hermite interpolation scheme. In fact, all of the following discrete expressions willbe equally valid for a large number of other beam formulations, where the discrete centerline geometry isdefined by polynomial interpolation, which can generally be expressed in terms of the generic shape functionmatrix H introduced above.Recall also, that the proposed SSIP laws from Section 5 solely depend on the centerline curve description,i. e., the rotation field does not appear in the additional contributions and hence its discretization is notrelevant in the context of this work. It is therefore sufficient to apply the discretization scheme for thecenterline field stated in eq. (39) to the expressions for the virtual work contributions δ Π ia presented inSection 5 and finally end up with the discrete element residual vectors r ia , / . The latter need to beassembled into the global residual vector R as it is standard in the (nonlinear) finite element method. Notethat the linearization of all the expressions presented in this Section 6.1 is provided in Appendix B. Discretization of the centerline curves according to (39), i. e., r j ≈ r h ,j = H j ˆ d j and δ r T j ≈ δ r Th ,j = δ ˆ d T j H T j , for both beam elements j = 1 , δ Π m,ss into its discrete counterpart δ Π m,ss,h = − ( m − ) (cid:90) l (cid:90) l c m,ss (cid:16) δ ˆ d T1 H T1 − δ ˆ d T2 H T2 (cid:17) r h, − r h, d h g − m + 52 h d s d s . (41)22efer to (30) for the definition of the constant c m,ss . Note that (41) only contributes to those scalar residuaassociated with the centerline, i. e., translational degrees of freedom ˆ d . This is a logical consequence of thefact that the SSIP law solely depends on the beams’ centerline curves, as discussed in detail in Section 5.For the sake of brevity, the index ’h’, indicating all discrete quantities, will be omitted from here on sinceall following quantities are considered discrete. In eq. (41), the discrete element residual vectors of the twointeracting elements j = 1 , r m,ss , = − ( m − ) (cid:90) l (cid:90) l c m,ss H T1 ( r − r ) d g − m + 52 d s d s and (42) r m,ss , = ( m − ) (cid:90) l (cid:90) l c m,ss H T2 ( r − r ) d g − m + 52 d s d s . (43)See Section 6.4 for details on the numerical quadrature required to evaluate these expressions. In analogy to the previous section, we discretize (38) and obtain the discrete element residual vectors r m,ls , = − (cid:90) l (cid:90) l c m,ls H T1 ( r − r ) d m +2 d s d s and r m,ls , = (cid:90) l (cid:90) l c m,ls H T2 ( r − r ) d m +2 d s d s . (44)As mentioned already in Section 5.3, the discrete element residual vectors in the specific case of Coulombicinteractions follow directly for m = 1 and c m,ls = C elstat λ λ . See Section 2.2.1 for the definition of C elstat and Section 5.3 for the definition of the linear charge densities λ i . Again, as mentioned in Section 5.3, thecase of long-range volume interactions only requires to adapt the constant prefactor via c m,ls = kmA A ρ ρ . It can be shown that the proposed SSIP approach from Section 4 in combination with the SSIP lawsfrom Section 5 fulfills the essential mechanical properties of objectivity, global conservation of linear andangular momentum as well as global conservation of energy. Due to the equivalent structure of the resultingspace-discrete contributions, e. g., equation (41), as compared to the terms obtained in macroscopic beamcontact formulations, we refer to the proof and detailed discussion of these important aspects in [21, Ap-pendix B]. The fulfillment of conservation properties will furthermore be verified by means of the numericalexamples in Section 7.2 and Section 7.4.
The singularity of inverse power laws for zero separation is a well-known pitfall when dealing with thiskind of interaction laws. See e. g. [2, p.137] for a discussion of this topic in the context of point-point LJinteraction as compared to a hard-sphere model. In numerical methods, one therefore typically appliesa regularization that cures the singularity and ensures the robustness of the method. Sauer [43] gives anexample for a regularized LJ force law between two half-spaces, where the force is linearly extrapolated belowa certain separation, which is chosen as 1 .
05 times the equilibrium spacing of the two half spaces. Also,existing, macroscopic beam contact formulations rely on the regularization of the seemingly instantaneousand infinite jump in the contact force when two macroscopic beams come into contact (see e. g. [21, 44]).However, the SSIP laws derived for disk-disk vdW or LJ interaction from Section 5.2 have not yet beenconsidered in literature. Note that LJ is the most general and challenging case considered in this work, sincestrong adhesive forces compete with even stronger repulsive forces whenever two fibers are about to comeinto contact. To be more precise, it is not only the strength of these competing forces, but also the highgradients in the force-distance relation that lead to a very stiff behavior of the governing partial differentialequations. This alone places high demands on the nonlinear solver, which in combination with the alreadymentioned singularity at zero separation g = 0, and the fact that LJ interaction laws are not defined forconfigurations g < g ≤ g reg,LJ in a manner very similar to [43] with the only difference that it is applied to the length-specific disk-diskforce law instead of the force law between two half spaces. Linear extrapolation means that the originalexpression ( m − / c m,ss g − m + in (42) and (43) is replaced by a linear equation a g + b in the gap g for all g < g reg,LJ . The two constants a and b are determined from the requirements that the force valueas well as the first derivative of the original and the linear expression are identical for the regularizationseparation g = g reg,LJ . Figure 7 shows both the original (blue) and the regularized (red) LJ disk-disk forcelaw as a function of the smallest surface separation g . Figure 7: Comparison of regularized (red) and full (blue) LJ disk-disk force law. Here, g reg,LJ = g LJ,eq,disk (cid:107) disk is shownexemplarily . The numerical experiment of adhesive fibers studied in [45] reveals that this regularization yields thealready mentioned great enhancement in terms of robustness as well as efficiency without any change in thesystem response. As shown in the comparison of the force-displacement curves therein, the results obtainedwith the full LJ and with the regularized LJ force law do indeed coincide down to machine precision. This isreasonable and expected, because we chose a regularization parameter g reg,LJ ≤ g LJ,eq,cyl (cid:107) cyl that is smallerthan any separation value g occurring anywhere in the system in any converged equilibrium state. Thus, thesolution never “sees” the modification to the vdW force law in the interval g < g reg,LJ and the results areidentical. However, since during the nonlinear iterations also non-equilibrium configurations with g < g reg,LJ occur, the nonlinear solution procedure is influenced in an extremely positive way, leading to an overall savingof a factor of five in the number of nonlinear iterations as compared to the full LJ interaction without anyregularization. For more details on the comparison including all parameter values we kindly refer the readerto [45]. See eq. (A.20) for an analytical expression of g LJ,eq,disk (cid:107) disk . .4. Numerical evaluation of n-dimensional integrals of intermolecular potential laws Generally, we use n nested loops of a 1D Gauss-Legendre quadrature scheme which is the well-establishedand de-facto standard method in nonlinear finite element frameworks and has been used also in previouspublications in the context of molecular interactions [22, 25]. Due to the strong nonlinearity, i. e., highgradients of the power laws, a large number of quadrature points is required in each dimension to achievesufficient accuracy. This effect is most critical for high exponents of the potential law, i. e., vdW andsteric interactions, and small separations of the interacting bodies. We thus implemented the possibility tosubdivide the domain of a finite element into n IS integration segments and apply an n GP -point Gauss ruleon each of them in order to achieve sufficient density of quadrature points in every case. Multi-dimensional numerical integration of the intermolecular potential laws as discussed above turnsout to be the crucial factor in terms of efficiency. For the following analysis of efficiency, we consider theassociated algorithmic complexity. Generally, all possible pairs of elements need to be evaluated, whichhas O (cid:0) n (cid:1) complexity. Let us assume we apply a total of n GP,tot,ele-length integration points along theelement length and n GP,tot,transverse integration points in the transversal, i. e., cross-sectional in-plane direc-tions. Thus, the complexity of an approach based on full 6D numerical integration over the 3D volumes ofthe two interacting bodies (cf. (12)) can be stated as O (cid:0) n · n · n (cid:1) . (45)In contrast to that, the novel SSIP approach proposed in Section 4.2 reduces the dimensionality of numericalintegration from six to two (cf. (28)) and thus yields O (cid:0) n · n (cid:1) (46)complexity. The resulting difference between both clearly depends on the problem size, type of interactionand other factors. To get an impression, typical numbers for the total number of quadrature points intransverse dimensions based on the numerical examples of Section 7 are given as n GP,tot,transverse = 10 . . . and can be as large as a factor of 10 . In additionto this tremendous saving from the inherent algorithmic complexity, the power law integrand has a smallerexponent due to the preliminary analytical integration in case of the SSIP approach. This in turn allowsfor a smaller number of integration points n ele · n GP,tot,ele-length for each of the two nested 1D integrationsalong the centerline, given the same level of accuracy. To give an example, the vdW interaction forcescales with an exponent of − − / m = 6). This makes another significant difference, especially ifvery small separations - as typically observed for contacting bodies - are considered. The combination ofhigh dimensionality and strong nonlinearity of the integrand renders the direct approach of six-dimensionalnumerical quadrature to evaluate eq. (12) unfeasible for basically any problem of practical relevance. Infact, even a single evaluation of the vdW potential of two straight cylinders to serve as a reference solutionturned out to be too computationally costly below some critical, small separation. See Section 7.1.1 fordetails on this numerical example. Note that although there might be more elaborate numerical quadratureschemes for these challenging integrands consisting of rational functions (see e. g. [46]), the basic problemand the conclusions drawn from this comparison of algorithmic complexities remain the same.These cost estimates based on theoretical algorithm complexity and the experience from rather smallacademic examples considered in Section 7 show that the SSIP approach indeed makes the difference betweenfeasible and intractable computational problems. This directly translates to the applicability to complexbiopolymer as well as synthetic fibrous systems that we have in mind and thus significantly extends therange of (research) questions that are accessible by means of numerical simulation.25 .6. Search algorithm and parallel computing In order to find the relevant pairs of interaction partners, the same search algorithms as in the case ofmacroscopic contact (between beams or 3D solids) may be applied, however, the obvious difference lies inthe search radius. For contact algorithms, a very small search radius covering the immediate surroundingof a considered body is sufficient, whereas for molecular interactions the search radius depends on the typeof interaction and must be at least as large as the so-called cut-off radius. Only at separations beyond thecut-off radius, the energy contributions from a particular interaction are assumed to be small enough toneglect them. Depending on the interaction potential and partners, the range and thus cut-off radius canbe considerably large which underlines the importance of an efficient search algorithm. In the scope of thiswork, a so-called bucket search strategy has been used, that divides the simulation domain uniformly intoa number of cells or buckets and assigns all nodes and elements to these cells to later determine spatiallyproximate pairs of elements based on the content of neighboring cells. This leads to an algorithmic complexityof O ( n ele ) and the search thus turned out to be insignificant in terms of computational cost as compared tothe evaluation of pair interactions as discussed in the preceding section. See [47] for an overview of searchalgorithms in the context of computational contact mechanics.To speed up simulations of large systems, parallel computing is a well-established strategy of ever in-creasing importance. Key to this concept is the ability to partition the problem such that an independentand thus simultaneous computation on several processors is enabled. In our framework, this partitioning isbased on the same bucket strategy that handles the search for interaction partners. Regarding the evalu-ation of interaction forces, a pair (or set) of interacting beam elements is assigned to the processor whichowns and thus already evaluates the internal and external force contribution of the involved elements. Atprocessor boundaries, i. e., if the two interacting elements are owned by different processors, one processoris chosen to evaluate the interaction forces and the required information such as the element state vectorof the element owned by the other processor is communicated beforehand. Upon successful evaluation ofthe element pair interaction, the resulting contribution to the element residual vector and stiffness matrix isagain communicated for the element whose owning processor was not responsible for the pair evaluation.
7. Numerical Examples
The set of numerical examples studied in this section aims to verify the effectiveness, accuracy androbustness of the proposed SSIP approach and the corresponding SSIP laws as a computational modelfor either steric repulsion, electrostatic or vdW adhesion and also a combination of those. Supplementaryinformation on the code framework and the algorithms used for the simulations is provided in Appendix C.
As a follow-up to the general discussion of using simplified SSIP laws in Section 5.1 and the proposalof specific closed-form analytic expressions in Section 5.2 and Section 5.3, this section aims to analyze theaccuracy in a quantitative manner. The minimal examples of two disks or two cylinders are considered inorder to allow for a clear and sound analysis of either the isolated SSIP laws or its use within the generalSSIP approach to modeling beam-to-beam interactions, respectively.
Throughout this section, we consider the example of vdW interaction, but analogous results are expectedfor steric interaction or any other short-range volume interaction. Specifically, we will focus on the approxi-mation quality of the proposed SSIP law from Section 5.2, which is based on the assumptions and resultingsimplifications discussed in-depth in Section 5.1. Recall that, beside the obviously most important surface-to-surface separation g , the rotation of the cross-sections around the closest point α (quantified by the angleenclosed by their tangent vectors) and potentially also the rotation components β , β (see Figure 6(a))have been identified as relevant degrees of freedom, yet are neglected in the simplified SSIP law proposedin Section 5.2. The influence of these factors, separation and rotation, on the approximation quality willthus be analyzed numerically in the following. 26ecall also from the discussions in Section 5.1 and 5.2 that only the regime of small separations will be ofpractical relevance in the case of short-range interactions considered here. However, we include the regimeof large separations in the following analyses, mainly because it will be interesting to see the transition fromsmall to large separations and confirm that potential values indeed drop by several orders of magnitude ascompared to the regime of small separations. Moreover, it is a question of theoretical interest and has beenconsidered in literature on vdW interactions [32]. This regime of large separations can be treated withoutany additional effort as described for the case of long-range interactions in Section 5.3 (where this regime isthe decisive one) if we take into account the corresponding remark on volume interactions.As presented in Section 2.3.2, analytical solutions for the special cases of parallel and perpendicularcylinders, for the regime of small and large separations, respectively, can be found in the literature [2, 3]and thus serve as reference solutions in this section. To the best of the authors’ knowledge, no analyticalreference solution has yet been reported for the intermediate regime in between the limits of large and smallseparations. Another source for reference solutions is the full numerical integration of the point pair potentialover the volume of the interacting bodies, however it is limited due to the tremendous computational cost.Only a combination of both analytical and numerical reference solutions thus allows for a sound verificationof the novel SSIP approach and the proposed SSIP laws.In the following analyses, either the SSIP, i. e., the interaction potential per unit length squared ˜˜ π ofa pair of circular cross-sections, the interaction potential per unit length ˜ π of a pair of parallel cylindersor the interaction potential Π of a pair of perpendicular cylinders will be plotted as a function of thedimensionless surface-to-surface separation g/R , respectively. For simplicity, the radii of the beams are setto R = R =: R = 1. Parallel disks and cylinders.
Figure 8(a) shows the SSIP ˜˜ π of two disks in parallel orientation, i. e., their normal vectors are parallelwith mutual angle α = 0, as a function of the normalized separation g/R . This is the simplest geometricalconfiguration and forms the basis of the proposed SSIP laws from Section 5. We thus begin our analysis withthe verification of the used analytical solutions in the limit of small (green line, cf. eq. (18)) and large (redline, cf. eq. (19)) separations by means of a numerical reference solution (black dashed line with diamonds)obtained from 4D numerical integration of the point pair potential law (4).
1 -6 1 -2.5 (a) (b)Figure 8: (a) VdW interaction potential per unit length squared ˜˜ π of two disks in parallel orientation over normalized surfaceseparation g/R . The analytical expressions (18) (green line) and (19) (red line) used as SSIP laws throughout this work areverified by means of a numerical reference solution (black dashed line with diamonds). (b) Subdivision of circular cross-sectionsinto integration sectors used to compute numerical reference solution. For rapidly decaying potentials, only the areas highlightedin dark and light gray considerably contribute to the total interaction potential. Figure 8(a) confirms that both analytical expressions match the numerical reference solution perfectlywell in the limit of large and small separations, respectively. As predicted, the interaction potential of27wo circular disks follows a power law with (negative) exponent 2 . .Note that all plots in this section are normalized with respect to the length scale R and the energy scale ρ ρ C vdW . It is remarkable that the obtained values span several orders of magnitude which illustrates thenumerical challenges associated with power laws, especially in the context of numerical integration schemesas discussed already in Section 6.4. Moreover, it underlines that the regime of large separations is practicallyirrelevant in the case of short-range interactions, because the potential values are basically zero as comparedto those obtained in the small separation regime.Regarding the full range of separations, one may ask where either of the two expressions may be usedgiven a maximal tolerable error threshold. As can be concluded from Figure 8(a), the resulting error is smallfor separations g/R < . g/R >
10 with a relative error below 7%. Inthe region of intermediate separations, the analytical solution for small separations seems to yield an upperbound, whereas the one for large separations seems to yield a lower bound for the interaction potential.Let us have a look at the efficiency gain from using the analytical solutions. The numerical referencesolution requires the evaluation of a 4D integral over both cross-sectional areas for a given separation g andhas been carried out in polar coordinates. Assuming Gaussian quadrature with the same number of Gausspoints n GP,tot,transverse in radial and circumferential dimension and for both cross-sections, this requires atotal of ( n GP,tot,transverse ) function evaluations. In contrast, the analytical expressions for the large andsmall separation limit, respectively, require only one function evaluation. This significant gain in efficiency ismost pronounced for small separations, where the number of required Gauss points increases drastically dueto the high gradient of the power law that needs to be resolved (see Section 6.4 for details). If the numberof Gauss points is not sufficient, this leads to so-called underintegration and we observed that the obtainedcurve of the numerical reference solution erroneously flattens (because the contribution of the closest-pointpair is not captured) or becomes steeper (because the contribution of the closest-point pair is overrated).For these reasons, the computation of an accurate numerical reference solution shown in Figure 8(a)requires quite some effort. The integration domains were subdivided into integration sectors (see Figure 8(b))in order to further increase the Gauss point density. But even in this planar disk-disk scenario requiringonly 4D integration, we reached a minimal separation of g/R ≈ −
3, below which the affordable numberof Gauss points was not sufficient to correctly evaluate the SSIP ˜˜ π by means of full numerical integration .For these very small separations, only the exact analytical dimensional reduction from 4D to 2D accordingto Langbein (cf. [32] and eq. (A.7)) allowed to compute an accurate numerical reference solution. Theanalytical solutions for the disk-disk interaction potential (18) (and (19)), used as SSIP law in eq. (33)(and (36)), thus realize a significant increase in efficiency and indeed only enable the accurate evaluationof the interaction potential in the regime of very small separations. Note that such small separations arehighly relevant if we consider fibers in contact since surface separations are expected to lie on atomic lengthscale in this case. For instance, the work of Argento et al. [22] mentions g = 0 . R varies from several nm for DNA to mm for synthetic polymerfibers, resulting in a potentially very small normalized separation g/R . An example for the simulation ofadhesive fibers in contact can be found in the authors’ recent contribution [45], which studies the peelingand pull-off behavior of two fibers attracting each other either via vdW or electrostatic forces.As a next step, the interaction potential per unit length ˜ π of two parallel straight beams is considered.The length of the beams is chosen sufficiently high such that it has no perceptible influence on the resultsand meets the assumption of infinitely long cylinders made to derive the analytical reference solution from[32]. Accordingly, a slenderness ratio ζ = l/R = 50 is used in the regime of small separations, whereas ζ = l/R = 1000 is used for large separations. Based on the experience from the disk-disk scenario, it is notsurprising that the full 6D numerical integration in this case exceeds the affordable computational resourcesby orders of magnitude and thus can not serve as a reliable reference solution. In fact, we were not able to Note that in the double logarithmic plot, a power law with exponent m is a linear function with slope m . The maximum number of n GP,tot,transverse = 8 ×
32 = 256 considered in the scope of this work led to several hours ofcomputation time on a desktop PC for the evaluation of ˜˜ π as a numerical reference solution for Figure 8(a). -1.5 1 -2 1-5 Figure 9: VdW interaction potential per unit length ˜ π of two parallel cylinders over normalized surface separation g/R . reproduce the theoretically predicted power law scaling in the regime of small separations despite using anumber of Gauss points that led to computation times of several days. However, instead of the numericalreference solution, the analytical solution for infinitely long cylinders in the limit of very small (black dashedline, cf. eq. (14)) and very large separations (blue dashed line, cf. eq. (15)) serves as a reference in Figure 9.Note that as compared to the case of two circular disks the exponent of the power laws and thus the slopeof the curves drops by one due to the integration over both cylinders’ length dimension.Interestingly, the SSIP approach using the simplified SSIP law from Section 5.2 (green line with crosses)does not yield the correct scaling behavior even in this case of parallel cylinders. This confirms the con-cerns from Section 5.1 that the simplified SSIP law neglecting any relative rotations of the cross-sectionsdeteriorates the accuracy of the approach in the case of short-ranged interactions in the regime of smallseparations. Due to this specific scenario of parallel cylinders, this deterioration can be attributed solely tothe rotation components β / (see Figure 6(a)) since the included angle of the cross-section normal vectors α is zero for every of the infinitely many pairs of cross-sections. Despite the correct trend of the resultinginteraction potential as an inverse power law of the surface separation, it must thus be stated that thesimplified SSIP law overestimates the strength of interaction and that the error increases with decreasingseparation . In the regime of large separations, however, the results for the SSIP approach (red line withcircles) perfectly match the analytical reference solution (blue dashed line). This confirms the hypothesisfrom Section 5.1 that the relative rotation of cross-sections is negligible in this regime and a high accuracycan be achieved with the simplified SSIP law. Although being of little practical importance here due to thenegligible absolute values, this is a first numerical evidence for the validity of the SSIP approach in generaland its high accuracy even in combination with simplified SSIP laws in the particular case of long-rangeinteractions to be considered in the following Section 7.1.2. Perpendicular disks and cylinders.
Up to now, we have only discussed the situation of parallel orientation of disks and cylinders. In thefollowing, the accuracy of the simplified SSIP laws as well as the SSIP approach for twisted configurationswill be analyzed by considering the most extreme configuration of perpendicular disks and cylinders. Again,computing a reference solution by means of full numerical integration was only affordable for the 4D caseof two disks. The results for perpendicular disks shown in Figure 10(a) confirm that there is no differencebetween perpendicular and parallel orientation for large separations and the scaling behavior of the numericalreference solution (black dashed line with diamonds) with exponent 6 is met by the simplified SSIP law (redline, cf. eq. (19)). On the other hand, there is a remarkable difference in the scaling behavior for small Note that the numerical integration error has been ruled out as cause for this behavior by choosing a high number ofGauss points n GP,tot,ele-length = 2 ×
50 = 100 for each of the 64 elements used to discretize each cylinder. A further increaseof n GP by a factor of five does not change the results using double precision. -61 -2 1 -2.5 (a) (b)Figure 10: (a) VdW interaction potential per unit length squared ˜˜ π of two perpendicular disks and (b) vdW interactionpotential Π of two perpendicular cylinders, plotted over the normalized surface separation g/R , respectively. separations. While the interaction potential of two parallel disks, which is the underlying assumption ofthe proposed SSIP law (green line, cf. eq. (18)), follows an inverse 2 . α , i. e.,the angle included by the cross-section normal vectors and again the error of the proposed simplified SSIPlaw increases with decreasing separation.Finally, the scenario of perpendicular cylinders is considered and Figure 10(b) shows the total interactionpotential Π as a function of the normalized smallest surface separation g/R . As discussed before, thecomputational cost of the full 6D numerical integration is too high to compute a reliable reference solutionin the case of two 3D bodies and we resort to the analytical solutions for the limits of very small andvery large separations, respectively. Note that in contrast to the case of infinitely long parallel cylinders(cf. Figure 9) the total interaction potential of infinitely long perpendicular cylinders is finite and the resultthus has dimensions of energy instead of energy per length. Perpendicular cylinders are worth to considerbecause they trigger both of the sources of errors that have been analyzed individually so far - neglectingrelative rotations α as well as β / in the simplified SSIP law. In short, the resulting accuracy is similaras for either perpendicular disks or parallel cylinders. In the decisive regime of small separations, the SSIPapproach based on the simplified SSIP law (green line with crosses) fails to reproduce the correct scalingbehavior of the analytical reference solution (black dashed line, cf. eq. (16)), whereas in the regime of largeseparations, the SSIP approach based on the simplified SSIP law (red line with circles) perfectly matchesthe analytical reference solution (blue dashed line, cf. eq. (17)). Conclusions.
First, this section reveals that full 6D numerical integration to compute the total interaction potential ofslender continua is by orders of magnitude too expensive and can not reasonably be used as a numericalreference solution even in minimal examples of one pair of cylinders. At most, 4D numerical integrationrequired for disk-disk interactions allows to compute numerical reference solutions for the intermediateregime of separations where no analytical solutions are known. This underlines the importance of reducingthe dimensionality of numerical integration to 2D as achieved by the proposed SSIP approach in order toenable the simulation of large systems as well as a large number of time steps.Second, the thorough analysis of the accuracy resulting from using the proposed simplified SSIP law,neglecting the cross-section rotations, reveals that one has to distinguish between the regime of small andlarge separations. In the decisive regime of small separations, we find that the scaling behavior deviates fromthe analytical prediction for perpendicular disks and parallel as well as perpendicular cylinders and that30he resulting error increases with decreasing separation. A remedy of this limitation could be a calibration,i.e. a scaling of the prefactor k in the simple SSIP law, to fit a given reference solution within a smallrange of separations (e.g. around the equilibrium distance of the LJ potential). In the authors’ recentcontribution [45], this pragmatic procedure is shown to reproduce the global system response very well.Still, it would be valuable to include the relative rotations of the cross-sections in the applied SSIP lawto obtain the correct asymptotic scaling behavior. To the best of the authors’ knowledge, no analyticalclosed-form expression has been published yet and the like is far from trivial to derive. Therefore, we leavethis as a promising enhancement of the novel approach, which we are currently working on and will addressin a future publication. As mentioned before, the regime of large separations is of little practical relevance inthe case of short-range volume interactions, however it is of some theoretical interest and the correspondingfindings and conclusions of this regime will hold true also for long-range interactions such as electrostaticsto be considered in the following section. Here, the results are in excellent agreement with the theoreticallypredicted power laws for parallel as well as perpendicular disks and cylinders.
Turning to long-range interactions, again parallel and perpendicular disks and cylinders will be consideredin order to analyze the accuracy of the simplified SSIP law from Section 5.3 both individually as well asapplied within the general SSIP approach proposed in Section 4.2. As before, Coulombic surface interactionsare chosen as specific example, however the conclusions are expected to hold true also for other types oflong-range interactions. As compared to the preceding section, the computation of a numerical referencesolution simplifies mainly due to the reduction from volume to surface interactions but also due to thesmaller gradient values, which need to be resolved in the regime of small separations thus requiring lessintegration points. This allows for a verification by means of a numerical reference solution also in the caseof cylinder-cylinder interaction. (a) (b)Figure 11: Electrostatic interaction potential per length squared ˜˜ π of (a) two parallel disks and (b) two perpendicular disks,plotted over the normalized surface separation g/R , respectively. Figure 11 shows the results for the simplified SSIP law obtained from the monopole-monopole interactionof two disk-shaped cross-sections in Section 5.3 (red line) and a numerical reference solution (black dashedline with diamonds). As expected, the proposed SSIP law excellently matches the reference solution in theregime of large separations, both for the parallel as well as perpendicular configuration. In both cases, therelative error is below 7% already for g/R = 1. The most important and remarkable result of this sectionhowever is the following. The inevitable error of the simplified SSIP law in the regime of small separationsdoes not carry over to beam-to-beam interactions as shown in Figure 12. For both parallel as well asperpendicular cylinders, the results from the SSIP approach using this simplified SSIP law from Section 5.3(red line with crosses) agree very well with the numerical reference solution (black dashed line with diamonds)31ver the entire range of separations. This confirms the theoretical considerations from Section 5.1 arguingthat the beam-beam interaction will be dominated by the large number of section pairs with large separations,which outweigh the contributions of the few section pairs with smallest separations.A closer look reveals that the relative error for the parallel cylinders is below 0 .
3% even for the smallestseparation g/R = 10 − considered here. For the presumably worst case of perpendicular cylinders, thisdeviation is even smaller with a relative error of 0 . (a) (b)Figure 12: Electrostatic interaction potential Π of (a) two parallel cylinders and (b) two perpendicular cylinders, plotted overthe normalized surface separation g/R , respectively. The slenderness ratio of the cylinders is ζ = L/R = 50.
Note that unlike in the case of short-range interactions, here the total interaction potential is consideredalso for the parallel cylinders. Due to the long range of interactions, the interaction potential per lengthdepends on the length of the cylinders and is thus no representative quantity. For Figure 12, a slendernessratio of ζ = L/R = 50 is chosen exemplarily. The two nested 1D integrals along the cylinder lengthdimensions are evaluated using n GP,ele-length = 5 Gauss points for each of the 64 elements used to discretizeeach cylinder. Additionally, n GP,circ = 8 ×
32 = 256 Gauss points over the circumference of each disk are usedto compute the numerical reference solution. In all cases, it has been verified that the numerical integrationerror does not influence the results noticeably.At this point, recall from Section 5.3 that the accuracy of the applied SSIP law can still be increasedwhenever deemed necessary by including more terms of the multipole expansion of the cross-sections. How-ever, because the results of this section show a high level of accuracy and the resulting simplification issignificant, the simplified SSIP law seems to be the best compromise for our purposes. To conclude thissection it can thus be stated that the novel SSIP approach as proposed in Section 4.2 in combination with thesimplified SSIP law from Section 5.3 is a simple, efficient, and accurate computational model for long-rangeinteractions of slender fibers. In the following, it will be applied to first numerical examples of deformable slender fibers in Section 7.3 and 7.4.
This numerical example aims to demonstrate the general ability of our proposed method to preclude pen-etration of two slender bodies that come into contact under arbitrary mutual orientation in 3D. No adhesiveforces are considered in this example. The setup is inspired by Example 1 in [21] where the macroscopic,so-called all-angle beam contact (ABC) formulation is used to account for the non-penetrability constraint.Here, we model the contact interaction based on the repulsive part of the LJ interaction potential (8). More32pecifically, we apply the novel SSIP approach as proposed in Section 4.2 in combination with the SSIP lawproposed in Section 5.2. The parameter specifying the strength of repulsion is set to be kρ ρ = 10 − .To be consistent throughout this article, we apply Hermitian Simo-Reissner beam elements instead of thetorsion-free Kirchhoff elements used in [21]. As compared to the original example, this requires us to replacethe hinged support of the upper beam by clamped end Dirichlet boundary conditions in order to eliminateall rigid body modes in this quasi-static example. The same number of three finite elements for the upper,deformable beam and one element for the lower, rigid beam is used. (a) initial configuration (b) time t = 1 . t = 1 . t = 2 . t = 2 . t = 3 . A sequence of the resulting simulation snapshots is shown in Figure 13. As expected, the two beams donot penetrate each other in any of the various mutual orientations throughout the simulation. Figure 14 (a) time t = 1 .
95 (b) time t = 1 .
97 (c) time t = 1 .
99 (d) time t = 2 . visualizes the contact force distributions in the most interesting time span before the beams reach the More precisely, the vectorial line load with dimensions of force per unit length is visualized as an arrow at each integrationpoint. The force resultant therefore equals the integral over the contour curve defined by the arrows’ tips (i. e. the area underthis curve), and not the vector sum of all arrows shown. This is important to understand because the number of visible arrowsper unit length depends on the discretization and is thus higher for the upper, deformable beam. t = 2 .
0. The force distribution quickly changes from a point-like force forlarge mutual angles to a broad distributed load for parallel beam axes. Note also that the line load hasa three dimensional shape where the out-of-plane component decreases with decreasing mutual angle untilboth beam axes and thus also the line loads lie in one plane at t = 2 .
0. Another remarkable result is thesymmetry between the line loads on both fibers. It nicely confirms that the novel approach indeed fulfills theexpected local equilibrium of interaction forces in good approximation. In contrast to existing, macroscopicformulations for beam contact, this is not postulated a priori in our approach and hence is a valuableverification at this point. See [25] for a comprehensive discussion of this important topic in the context ofcontact between 3D solids described by inter-surface potentials. The global equilibrium of contact forces onthe other hand is fulfilled exactly, as can be concluded from the global conservation of linear momentumthat can be shown analytically as outlined in Section 6.2. In this numerical example, we found that the sumof all reaction forces in either of the spatial dimensions is indeed zero with a maximal residuum of 10 − throughout all simulations considered here, which confirms the statement numerically. (a) (b)Figure 15: (a) Reaction force and (b) interaction potential over time. Figure 15 shows the resulting vertical reaction force as well as the interaction potential over time. Dueto the inverse-twelve power law and the extremely small separations of the interacting bodies, the numericalintegration of the disk-disk interaction forces is very challenging and we studied the influence of the numberof Gauss points. For this purpose, the number of integration segments per element with five Gauss pointseach is set to 30, 50, or 64, respectively. Interestingly, the interaction potential shown in Figure 15(b)seems to be more sensitive with respect to the integration error than the vertical reaction force shown inFigure 15(a) despite the fact that the latter has a higher inverse power law exponent. Presumably, this isdue to the fact that the reaction force is dominated by the bending deformation of the beams. For reference,the reaction force obtained by using the macroscopic ABC formulation is shown as well and is in excellentagreement with the one resulting from the repulsive part of the LJ interaction potential.A more comprehensive comparison of this novel SSIP approach to model contact between beams basedon (the repulsive part of) the molecular LJ interaction and existing, macroscopic formulations based onheuristic penalty force laws is a highly interesting subject that is worth to investigate in the future.
The following example consists of two initially straight and parallel, deformable fibers that attract eachother due to their surface charge of opposite sign. Its setup is kept as simple as possible to allow for anisolated and clear analysis of the physical effects as well as the main characteristics of the proposed SSIPapproach. In a first step presented here, the interplay of elasticity and electrostatic attraction in the regimeof large separations is studied. Additionally, the authors’ recent contribution [45] considers the scenario ofseparating these adhesive fibers starting from initial contact and studies a variety of physical effects andinfluences in depth, which would go beyond the scope of this work.34n this numerical example, we are interested in the static equilibrium configurations for varying attractivestrength. As shown in Figure 16(a), two straight beams of length l = 5 are aligned with the global y -axisat an inter-axis separation d = 5. Both are simply supported and restricted to move only within the xy -plane and rotate only around the global z -axis. The beams have a circular cross-section with radius R =0 .
02 which results in a slenderness ratio of ζ = 250. Cross-section area, area moments of inertia andshear correction factor are computed using standard formula for a circle. A hyperelastic material law withYoung’s modulus E = 10 and Poisson’s ratio ν = 0 . dl y x (a) Problem setup: undeformed configuration. (b) Static equilibrium configurations for varying attractivestrength. Solution for beam centerlines and corresponding valueof the potential law prefactor k shown in the same color.Figure 16: Two parallel beams with constant surface charge density (left beam positive, right beam negative). long-range Coulomb interactions in Section 5.3. Both beams are nonconducting with a constant surfacecharge density of σ = 1 . σ = − .
0, respectively. For simplicity, we vary the prefactor k of theunderlying Coulomb law Φ( r ) = k r − to vary the strength of attraction. However, as becomes clearfrom (38), this is equivalent to a variation of surface charge densities because in our case the product of thesequantities is a constant prefactor in all relevant equations. In order to evaluate the electrostatic force andstiffness contributions, Gauss quadrature with two integration segments per element and ten Gauss pointsper integration segment is applied. This turns out to be fine enough to not change the presented resultsperceptibly. More precisely, the difference in the displacement of the beam midpoint for n GP = (2 × ascompared to (2 × is below 10 − . No cut-off radius is applied here, i. e., the contributions of all Gausspoint pairs are evaluated and included.Figure 16(b) finally shows the resulting static equilibrium configurations for different levels of attractivestrength. As expected, the beams are increasingly deflected and pulled towards each other if the prefactorof the applied Coulomb law k and thus the attractive strength is increased. Like the problem definition, alsoall the solutions are perfectly symmetric with respect to the vertical axis of symmetry located at x = d/ k = 1 . (a) Resulting electrostatic forces evaluated at the Gauss points.(b) Individual electrostatic force contributions of all element pairs.(c) Electrostatic force contributions of the closest element pair. (d) Electrostatic force contributions of the next-nearest pair.Figure 17: Electrostatic forces acting on the beams for k = 1 .
0. Color indicates force magnitude.
To conclude this example of two charged, attractive beams, we briefly look at the nonlinear solver. New-ton’s method without any adaptations is used here to allow for a clear and meaningful analysis of nonlinear36onvergence behavior. The solutions for k ≤ . k = 0 .
1. As convergence criteria, we enforced both for the Euclidean norm of theresidual vector (cid:107) R (cid:107) < − and for the norm of the iterative displacement update vector (cid:107) ∆ X (cid:107) < − .In fact, this combination leads to (cid:107) R (cid:107) < − in almost all equilibrium configurations shown here. Due to the high gradients in the inverse power laws, molecular interactions give rise to highly dynamicsystems. This is a first, simple example for a dynamic system consisting of two oppositely charged fiberswith a hinged support at one end each, that will snap into contact. In the initial configuration shownin Figure 18(a), the straight fibers include an angle of 45 ◦ and their axes are separated by 5 R in the out-of-plane direction z . With a cross-section radius R = 0 .
02 and length l set to l = 5, they have a highslenderness ratio of ζ = 250 and 354. Each of the fibers is discretized by 10 Hermitian Simo-Reissner beamelements and the material parameters are chosen to be E = 10 , ν = 0 .
3, and ρ = 10 − . The fibers carry aconstant, opposite surface charge σ / = ± . C elstat = 0 .
4. In order to start from a stress-free initial configuration, the chargeof one of the fibers is ramped up linearly within the first 100 time steps. We apply the SSIP approachas proposed in Section 4.2 and applied to Coulomb interactions in Section 5.3. A total of 5 integrationsegments per element with 10 Gauss points each is used to evaluate these electrostatic contributions. Thecontact interaction between the fibers is modeled by the line contact formulation proposed in [20], using apenalty parameter ε = 10 and 20 integration segments per element with 5 Gauss points each for numericalintegration. An undetected crossing of the fiber axes is prevented by applying the modified Newton methodlimiting the maximal displacement increment per nonlinear iteration to R/ ( l , , R )( l , l, R ) ( l, l, y xz (a) Problem setup. (b) Energy over time.Figure 18: Two oppositely charged, crossed beams dynamically snap into contact. In terms of temporal discretization, we apply the Generalized-Alpha scheme for Lie groups as proposedin [48] and set the spectral radius at infinite frequencies to ρ ∞ = 0 . t = 10 − is applied to account for the highly dynamic behavior of this system. Figure 19shows a sequence of simulation snapshots where the electrostatic forces on both fibers are visualized asgreen arrows. We observe a large variety of mutual orientations of the two fibers and a strong couplingof adhesive, repulsive and elastic forces that demonstrate the effectiveness and robustness of the proposedSSIP approach. Most importantly, we see that the total system energy is preserved with very little deviationof ±
2% as shown in Figure 18(b). Note that the negative energy values result from defining the zero level ofthe interaction potential at infinite separation as described in Section 2.1. Based on this numerical example,37e can thus conclude that the novel SSIP approach proves to be effective as well as robust in a highlydynamic example with arbitrary mutual orientations of the fibers in three dimensions. (a) initial configuration (b) time t = 1 × − , ramp-up of chargecompleted (c) time t = 5 × − (d) time t = 1 × − (e) time t = 2 × − (f) time t = 4 × − (g) time t = 6 × − (h) time t = 8 × − (i) time t = 1 × − Figure 19: Sequence of simulation snapshots. Electrostatic forces acting on both fibers shown in green.
8. Conclusions and Outlook
This contribution proposes the first 3D beam-to-beam interaction model for molecular interactions suchas electrostatic, van der Waals (vdW) or repulsive steric forces between curved slender fibers undergoinglarge deformations. While the general model is not restricted to a specific beam formulation, in the presentwork it is combined with the geometrically exact beam theory and discretized via the finite element method.A direct evaluation of the total interaction potential for general 3D bodies requires the integration of con-tributions from molecule or charge distributions over the volumes of the interaction partners, leading toa 6D integral (two nested 3D integrals) that has to be solved numerically. The central idea of our novelapproach is to formulate reduced interaction laws for the resultant interaction potential between a pair ofcross-sections of two slender fibers such that only the two 1D integrals along the fibers’ length directions haveto be solved numerically. This section-to-section interaction potential (SSIP) approach therefore reduces thedimensionality of the required numerical integration from 6D to 2D and yields a significant gain in efficiency,which only enables the simulation of relevant time and length scales for many practical applications. Beingthe key to this SSIP approach, the analytical derivation of the specific SSIP laws is based on careful consid-eration of the characteristics of the different types of molecular interactions, most importantly their pointpair potential law and the range of the interaction. In a first step, the most generic form of the SSIP law,which is valid for arbitrary shapes of cross-sections and inhomogeneous distribution of interacting points(e. g. atoms or charges) within the cross-sections has been presented before the assumptions and resultingsimplifications for the specific SSIP laws have been discussed in detail. For the practically relevant case of38omogeneous, disk-shaped cross-sections, specific, ready-to-use SSIP laws for short-range volume interac-tions such as vdW or steric interactions and for long-range surface interactions such as Coulomb interactionshave been proposed. We would like to stress that postulating the general structure of the SSIP law andfitting the free parameters to e. g. experimental data is one of the promising alternatives to the strategyof analytical derivation of the SSIP law as applied in this article. It is also important to emphasize thatthe general SSIP approach can be seamlessly integrated into an existing finite element framework for solidmechanics. In particular, it does neither depend on any specific beam formulation nor the applied spatialdiscretization scheme and in the context of the present work, we have exemplarily used it with geometricallyexact Kirchhoff-Love as well as Simo-Reissner type beam finite elements. Likewise, it is independent ofthe temporal discretization and we have used it along with static and (Lie group) Generalized-Alpha timestepping schemes as well as inside a Brownian dynamics framework.The accuracy of the proposed SSIP laws as well as the general SSIP approach has been studied in athorough quantitative analysis using analytical as well as numerical reference solutions for the case of vdWas well as electrostatic interactions. We find that a very high level of accuracy is achieved for long-rangeinteractions such as electrostatics both for the entire range of separations as well as all mutual angles ofthe fibers from parallel to perpendicular. In the case of short-range interactions, however, the derived SSIPlaw without cross-section orientation information slightly overestimates the asymptotic power-law exponentof the interaction potential over separation. As a pragmatic solution, a calibration of the simple SSIP lawhas been proposed to fit a given reference solution in the small yet decisive range of separations aroundthe equilibrium distance of the Lennard-Jones (LJ) potential. In the authors’ recent contribution [45],this strategy led to very good agreement in the force response on the system level. While this accuracymight already be sufficient for certain real world applications, our future research work will focus on thederivation of enhanced interaction laws including information about the cross-section orientation with theaim to achieve higher accuracy and the exact asymptotic scaling behavior.The presented set of numerical examples finally demonstrates the effectiveness and robustness of theSSIP approach to model steric repulsion, electrostatic or vdW adhesion. Several important aspects such asthe influence of the Gauss integration error and the spatial discretization error as well as local and globalequilibrium of forces and conservation of energy are studied in these simulations, including quasi-staticand dynamic scenarios as well as arbitrary mutual orientations and separations of the interacting fibers. Inorder to remedy the characteristic singularity of inverse power interaction laws in the limit of zero separation,we have proposed a numerical regularization of the LJ SSIP law, which leads to a significant increase inrobustness and efficiency, saving a factor of five in the number of nonlinear iterations while yielding identicalresults.
Appendix A. Examples for the Derivation and Analysis of the Two-Body Interaction Poten-tial and Force Laws for Parallel Disks and Cylinders
The aim of this appendix is to present the mathematical background of analytical solutions for two-bodyinteraction potential as well as force laws. Generally, the strategy of pairwise summation, i. e., integrationof a point pair potential, is applied. See Section 2.3.2 for a discussion of the applicability of this approach.Exemplarily, we consider the interaction between two parallel disks and two parallel cylinders since thesescenarios proved to be most important throughout the derivation of SSIP laws as well as their verificationin Section 5 and 7.1, respectively. In addition, we are interested in the total LJ interaction potential andforce law in the limit of small separations, because the regularization proposed in Section 6.3 is based onthese theoretical considerations. Finally, also the equilibrium spacing g LJ,eq,cyl (cid:107) cyl of two infinitely longcylinders interacting via the LJ potential will be derived and has proven helpful in order to choose an almoststress-free initial configuration of two deformable, straight fibers e. g. in the authors’ recent contribution [45]studying the peeling and pull-off behavior.
Appendix A.1. A generic interaction potential described by an inverse power law
Instead of Φ vdW ( r ) from (4) or any other particular interaction type, here, we rather use the more generalpower law Φ m ( r ) = k m r − m for the point pair potential. As noted already in [32], this does not introduce any39dditional complexity in the derivations and the solutions can directly be used for other exponents m . Wewill make use of this fact when considering LJ interaction between two disks and two cylinders analyticallyin Appendix A.2.1 and Appendix A.2.2, respectively. These findings are to be used in the context ofderiving a proper regularization of the potential laws in Section 6.3. Appendix A.1.1. Disk-disk interaction
The following refers to the analytical solutions for the disk-disk vdW interaction potential from literaturethat is summarized in Section 2.3.2. Let us first state the underlying mathematical problem. We would liketo find an analytical solution for the required 4D integral C m over the circular area of each disk C m := (cid:90) (cid:90) A ,A Φ m ( r ) d A d A with Φ m ( r ) = k m r − m (A.1)in order to arrive at the disk-disk interaction potential˜˜ π m = ρ ρ C m . (A.2) Details on 2(a) in Section 2.3.2: The regime of large separations.
For the limit of large separations g (cid:29) R , R , the solution is quite straightforward and can be explained insimple words as follows. The distance of any point in a disk to its center is of order O ( R ) and thus muchsmaller than the disks’ surface-to-surface separation g :˜ r / = O ( R / ) (cid:28) O ( g ) (A.3)Figure A.20 illustrates the introduced geometrical quantities. The distance r between any two points x in b b g R R rA A d ˜ r ˜ r b b Figure A.20: Two circular cross-sections, i. e. disks in parallel alignment disk 1 and x in disk 2 may therefore be approximated by the inter-axis distance d = g + R + R : r = (cid:107) x − x (cid:107) = (cid:107) r + ˜ r − r − ˜ r (cid:107) ≈ (cid:107) r − r (cid:107) = d (A.4)Double integration over both disks is hence equivalent to a multiplication with the disks’ areas A , A C m ,ls ≈ A A Φ m ( r = d ) (A.5)and finally we end up with the sought-after expression for the general disk-disk interaction potential in thelimit of large separations ˜˜ π m, disk (cid:107) disk,ls ≈ ρ ρ A A Φ m ( r = d ) . (A.6)Note that this approximation is valid for arbitrary pair interaction functions Φ( r ). Moreover, this solutiondoes not even depend on the parallel orientation of the disks. It is valid for all mutual angles of the diskswhich is important because we will apply it to arbitrary configurations of deflected beams. For the specialcase of parallel disks, this result can alternatively be obtained by the sound mathematical derivation of [32,eq. (10)]. The leading term of his hypergeometric series is identical to the right hand side of equation (A.6).40 b g R R p tϕψ Figure A.21: Integration over the cross-sections at small separations, figure taken from [32] with adapted notation.
Details on 1(a) in Section 2.3.2: The regime of small separations.
Now, we consider the limit of small separations g (cid:28) R , R . The problem has been studied by Langbein [32]in the context of vdW attraction of rigid cylinders, rods or fibers. In the following, we will briefly presentthe central mathematical concept of his derivations. The basic idea is to choose a favorable set of integrationvariables p, t, ϕ, ψ as shown in Figure A.21. In this way, the four dimensional integral (A.1) can be reducedto a double integral because the integrand Φ m ( p ) does not depend on the angles ϕ and ψ : C m = (cid:90) A (cid:90) A Φ m ( p ) d A d A = (cid:90) p (cid:90) t (cid:90) ϕ (cid:90) ψ Φ m ( p ) d ψ d ϕ d t d p = (cid:90) p (cid:90) t Φ m ( p ) 2 pϕ ( p, t ) 2 tψ ( p, t ) d t d p (A.7)where cos( ϕ ) = p + t − R pt , cos( ψ ) = t + d − R td , d = g + R + R For a general potential law Φ m ( r = p ) = k m p − m , this reads C m = 4 k m (cid:90) p (cid:90) t p − m +1 ϕtψ d t d p (A.8)Making use of g (cid:28) R , R and introducing reduced variables ¯ p = p/g and ¯ t = t/g leads to C m ,ss = 4 k m (cid:114) R R R + R (cid:90) g +2 R +2 R g p − m +1 (cid:90) pg √ t − g arccos (cid:18) tp (cid:19) d t d p (A.9)= 4 k m (cid:114) R R R + R g − m +7 / (cid:90) ∞ ¯ p − m +1 (cid:90) ¯ p (cid:112) ¯ t − (cid:18) ¯ t ¯ p (cid:19) d¯ t d ¯ p (A.10)Another substitution of variables x = ¯ t/ ¯ p and interchanging the order of integration finally yields the solu-tion C m ,ss = g − m + 72 k m π ( m − (cid:114) R R R + R Γ( m − ) Γ( m − )Γ( m −
2) Γ( m −
1) for m > (A.11)Here, Γ denotes the gamma function which is defined by Γ( z ) = (cid:82) ∞ w z − e − w d w .Multiplication with the particle densities finally results in the sought-after general disk-disk interactionpotential for the regime of small separations˜˜ π m, disk (cid:107) disk,ss = ρ ρ C m ,ss for m > (A.12)that can be further specified by means of m = 6 and k = − C vdW to end up with ˜˜ π vdW,disk (cid:107) disk,ss as in (18). Note that in the original article [32], the final form of C m (eq. (15) on p. 65) seems to be incorrect. A comparison with [3,p. 172] for the case of vdW potential with m = 6 confirms the solution presented here. Additionally, this solution is verifiedby means of numerical quadrature in Section 7.1.1 (cf. Figure 8(a)). emarks.
1. Note that this solution is valid for exponents m > / r ).2. Note however the conceptual similarity of this expression to the one valid for the limit of large sepa-rations (A.6). Here, we also find a power law, however in the surface-to-surface distance g instead ofthe inter-axis distance d and with a different exponent. Appendix A.1.2. Cylinder-cylinder interaction
Considering the case of two parallel cylinders, we are interested in the length-specific interaction potential˜ π m, cyl (cid:107) cyl = lim l →∞ l (cid:90) l / − l / (cid:90) ∞−∞ (cid:90) (cid:90) A ,A ρ ρ Φ m ( r ) d A d A d s d s with Φ m ( r ) = k m r − m . (A.13)The integral over s = − l / . . . l / l since the integrand is constant along s and thusimmediately cancels with the normalization factor 1 /l .Exemplarily, we want to discuss the more interesting and challenging regime of small separations here.Following [32, p.63], one can interchange the order of integration, solve the integral over the infinitely longcylinder length analytically in a first step, and then make use of the generic solution for C m, ss from (A.11),but this time with reduced exponent m −
1, to end up with˜ π m, cyl (cid:107) cyl ,ss = (cid:90) (cid:90) A ,A (cid:90) ∞−∞ ρ ρ Φ m ( r ) d s d s d A d A (A.14)= 3 π ρ ρ k m C m − , ss k m − . (A.15)Plugging in m = 6 for vdW interaction directly yields the two-body interaction potential per unit lengthfor two parallel cylinders in the regime of small separations ˜ π vdW,cyl (cid:107) cyl,ss as stated in eq. (14). This genericexpression (A.15) will be exploited when deriving the total LJ interaction law in Appendix A.2.2. Appendix A.2. Lennard-Jones force laws in the regime of small separations
As compared to the preceding sections, we now want to turn to the LJ interaction consisting of twopower law contributions, one adhesive and one repulsive, respectively. Our motivation is to study thecharacteristics of the resulting, superposed force laws for disk-disk as well as cylinder-cylinder interactionsby means of theoretical analysis of the analytical expressions. These findings shall prove valuable whenderiving an effective yet accurate regularization of the LJ potential law for the limit of zero separation inSection 6.3. We therefore focus on the regime of small separations throughout this section.Coming from the expressions for the two-body interaction potential ˜˜ π m, disk (cid:107) disk,ss and ˜ π m, cyl (cid:107) cyl,ss derivedfor a generic point pair potential Φ m in Appendix A.1, we will now sum the adhesive contribution m = 6and the repulsive contribution m = 12 and differentiate once to arrive at the desired LJ force laws. Appendix A.2.1. Disk-disk interaction
As outlined above, we make use of (A.12) for both parts of the LJ interaction and immediately obtain˜˜ π LJ,disk (cid:107) disk,ss = ˜ k g − + ˜ k g − (A.16)where the following abbreviations for the constant prefactors have been introduced:˜ k := π k ρ ρ (cid:114) R R R + R Γ (cid:18) (cid:19) and ˜ k := k ρ ρ (cid:114) R R R + R . − k , k specifying the point pair LJ potential to the other commonly used set Φ LJ,eq , r
LJ,eq according to (8): k = 2Φ LJ,eq r and k = − Φ LJ,eq r (A.18)42ifferentiation with respect to the separation yields the disk-disk LJ force law˜˜ f LJ,disk (cid:107) disk,ss = − d ˜˜ π LJ,disk (cid:107) disk,ss d g = 52 ˜ k g − + 172 ˜ k g − . (A.19)See Section 6.3 for a plot of the function. This expression allows us to determine some very interesting,characteristic quantities like the equilibrium spacing g LJ,eq,disk (cid:107) disk , i. e., the distance where the force van-ishes: g LJ,eq,disk (cid:107) disk = (cid:32) −
175 ˜ k ˜ k (cid:33) ≈ .
653 513 r LJ,eq . (A.20)Due to the fact, that repulsive contributions from proximate point pairs decay faster than the adhesivecontributions, we obtain a smaller equilibrium spacing as compared to the scenario of a point pair. Anotherdifferentiation allows us to determine the value of the force minimum, i. e., the maximal adhesive force, andthe corresponding separation˜˜ f LJ,disk (cid:107) disk,min ≈ .
904 115 ρ ρ (cid:114) R R R + R r LJ,eq Φ LJ,eq (A.21) g ˜˜ f LJ,disk (cid:107) disk,min = (cid:32) − k ˜ k (cid:33) ≈ .
771 844 8 r LJ,eq ≈ .
181 07 g LJ,eq,disk (cid:107) disk . (A.22)These quantities turn out to be decisive for the choice of a regularized, i. e., altered force law that is to beused instead of the original one in order to cure the numerical problems that come with the singularity atzero separation g = 0.In summary, we have found an analytical, closed-form expression for the disk-disk LJ force law (A.19),valid in the regime of small separations and for parallel disks. By means of elementary algebra, we werethus able to determine analytical expressions for the characteristic equilibrium spacing as well as value andspacing of the force minimum. Appendix A.2.2. Cylinder-cylinder interaction
As in the previous section (Appendix A.1.2), we want to restrict ourselves to parallel, infinite cylindersand consider the length-specific interaction potential as well as force law. Again, starting from the expressionfor a generic interaction potential (A.15), superposition yields˜ π LJ,cyl (cid:107) cyl,ss = ˜ k cyl , g − + ˜ k cyl , g − (A.23)where the following abbreviations have been introduced:˜ k cyl , := π k ρ ρ (cid:114) R R R + R and ˜ k cyl , := 5 .
818 68e − k π ρ ρ (cid:114) R R R + R (A.24)Differentiation with respect to the separation yields the cylinder-cylinder LJ force law˜ f LJ,cyl (cid:107) cyl,ss = − d ˜ π LJ,cyl (cid:107) cyl,ss d g = 32 ˜ k cyl , g − + 152 ˜ k cyl , g − . (A.25)that shall be further analyzed in the following. To begin with, the equilibrium spacing for two parallelcylinders interacting via a LJ potential can be derived as g LJ,eq,cyl (cid:107) cyl = (cid:32) − k cyl , ˜ k cyl , (cid:33) ≈ .
571 69 r LJ,eq . (A.26)43his is an extremely interesting and important result, since it leads the way to the non-trivial stress-freeconfiguration of two flexible, initially straight fibers. We make use of this knowledge e. g. in [45]. Again,since the repulsive contribution of proximate point pairs decays faster than the adhesive contribution, thisequilibrium spacing is smaller than g LJ,eq,disk (cid:107) disk for the disks, that in turn is smaller than r LJ,eq in thefundamental case of a point pair. The very same value of 57% of the point pair equilibrium spacing hasalready been mentioned as a side note by Langbein [32, p. 62], however, without presenting the detailed,comprehensive derivation. In addition to the equilibrium spacing, we can again determine and look at thevalue and location of the force minimum˜ f LJ,cyl (cid:107) cyl,min ≈ .
116 34 ρ ρ (cid:114) R R R + R r LJ,eq Φ LJ,eq (A.27) g ˜ f LJ,cyl (cid:107) cyl,min = (cid:32) − k cyl , ˜ k cyl , (cid:33) ≈ .
701 04 r LJ,eq ≈ .
226 25 g LJ,eq,cyl (cid:107) cyl . (A.28)Here, we find that the force minimum, i. e., the maximal adhesive force is slightly shifted towards a smallerseparation as compared to the disk-disk interaction. However, expressed in terms of the respective equilib-rium spacing g LJ,eq,cyl (cid:107) cyl , the value is slightly larger as compared to 1 . g LJ,eq,disk (cid:107) disk from (A.22). Withthese results we conclude the derivation and analysis of LJ force laws in the regime of small separations andsummarize the most important results in the following table to serve as a quick access reference.
Appendix A.2.3. Summary
The following table gives an overview of some important quantities characterizing the LJ force laws forpoint-point, parallel disk-disk, and parallel cylinder-cylinder interaction.equilibrium spacing location of force min. min. force value r LJ,eq / g LJ,eq r f LJ,min / g f LJ,min f LJ,min / ˜˜ f LJ,min / ˜ f LJ,min point-point 1 [ r LJ,eq ] 1 .
11 [ r LJ,eq ] 2 . (cid:104) Φ LJ,eq r LJ,eq (cid:105) disk (cid:107) disk 0 .
65 [ r LJ,eq ] 0 .
77 [ r LJ,eq ] 0 . (cid:104) ρ ρ (cid:113) R R R + R r LJ,eq Φ LJ,eq (cid:105) cylinder (cid:107) cylinder 0 .
57 [ r LJ,eq ] 0 .
70 [ r LJ,eq ] 2 . (cid:104) ρ ρ (cid:113) R R R + R r LJ,eq Φ LJ,eq (cid:105)
Table A.3: Comparison of characteristic quantities of LJ force laws for a pair of points, parallel disks and parallel cylinders.
Appendix B. Linearization of the Virtual Work Contributions from Molecular Interactions
Generally, the discrete residual vectors r ia ,j from molecular interactions between two beam elements j =1 , x k of both beam elements k = 1 ,
2. Consistent linearization thus yieldsthe following four sub-matrices k jk to be considered and assembled into the global stiffness matrix, i. e.,system Jacobian K : k := d r ia , d ˆ x , k := d r ia , d ˆ x , k := d r ia , d ˆ x , k := d r ia , d ˆ x (B.1)Note that the linearization with respect to the primary variables ˆ x k of both interacting beam elementssimplifies due to the fact that the residuals r ia ,j do not depend on the cross-section rotations as discussedalong with the derivation of the specific SSIP laws in Section 5. Thus, only the linearization with respectto the centerline degrees of freedom ˆ d k yields non-zero entries and are therefore presented in the remainderof this section. 44 ppendix B.1. Short-range volume interactions such as van der Waals and steric repulsion The linearization of the residual contributions with respect to the primary variables ˆ x of both interactingbeam elements is directly obtained from differentiation of eq. (42) and (43): k m,ss , = ( m − ) (cid:90) l (cid:90) l c m,ss (cid:18) − d − g − m + 52 H T1 H + (cid:18) d − g − m + 52 + ( m − ) d − g − m + 32 (cid:19) H T1 ( r − r ) ⊗ ( r − r ) T H (cid:19) d s d s (B.2) k m,ss , = ( m − ) (cid:90) l (cid:90) l c m,ss (cid:18) d − g − m + 52 H T1 H − (cid:18) d − g − m + 52 + ( m − ) d − g − m + 32 (cid:19) H T1 ( r − r ) ⊗ ( r − r ) T H (cid:19) d s d s (B.3) k m,ss , = ( m − ) (cid:90) l (cid:90) l c m,ss (cid:18) d − g − m + 52 H T2 H − (cid:18) d − g − m + 52 + ( m − ) d − g − m + 32 (cid:19) H T2 ( r − r ) ⊗ ( r − r ) T H (cid:19) d s d s (B.4) k m,ss , = ( m − ) (cid:90) l (cid:90) l c m,ss (cid:18) − d − g − m + 52 H T2 H + (cid:18) d − g − m + 52 + ( m − ) d − g − m + 32 (cid:19) H T2 ( r − r ) ⊗ ( r − r ) T H (cid:19) d s d s . (B.5)See eq. (30) for the definition of the constant c m,ss and eq. (39) for the definition of the shape functionmatrices H j . Note that the ‘mixed’ matrix products H T1 ( . . . ) H and H T2 ( . . . ) H lead to off-diagonal entriesin the tangent stiffness matrix of the system which couple the corresponding degrees of freedom. This isreasonable and necessary because these couplings represent the interaction between the respective bodies. Appendix B.2. Long-range surface interactions such as electrostatics
In analogy to the previous section, differentiation of eq. (44) yields k m,ls , = − (cid:90) l (cid:90) l c m,ls H T1 d m +4 (cid:16) H d m +2 − ( m + 2) d m ( r − r ) ⊗ ( r − r ) T H (cid:17) d s d s = (cid:90) l (cid:90) l c m,ls (cid:18) − d m +2 H T1 H + ( m + 2) d m +4 H T1 ( r − r ) ⊗ ( r − r ) T H (cid:19) d s d s (B.6) k m,ls , = (cid:90) l (cid:90) l c m,ls (cid:18) d m +2 H T1 H − ( m + 2) d m +4 H T1 ( r − r ) ⊗ ( r − r ) T H (cid:19) d s d s (B.7) k m,ls , = (cid:90) l (cid:90) l c m,ls (cid:18) d m +2 H T2 H − ( m + 2) d m +4 H T2 ( r − r ) ⊗ ( r − r ) T H (cid:19) d s d s (B.8) k m,ls , = (cid:90) l (cid:90) l c m,ls (cid:18) − d m +2 H T2 H + ( m + 2) d m +4 H T2 ( r − r ) ⊗ ( r − r ) T H (cid:19) d s d s (B.9)See again eq. (39) for the definition of the shape function matrices H j . As mentioned before, the discreteelement residual vectors in the specific case of Coulombic interactions directly follow for m = 1 and c m,ls = C elstat λ λ . See Section 2.2.1 for the definition of C elstat and Section 5.3 for the definition of the linearcharge densities λ i . Again, as mentioned already in Section 5.3, the case of long-range volume interactionsonly requires to adapt the constant prefactor via c m,ls = kmA A ρ ρ .45 ppendix C. Supplementary information on algorithms and code framework used for thesimulations Implementation
All novel methods have been implemented in C++ within the framework of the multi-purpose and multi-physics in-house research code BACI [49].
Integration into existing code framework
The novel SSIP approach can be integrated very well in an existingnonlinear finite element solver for solid mechanics. In particular, it does not depend on a specific beam(finite element) formulation and has been used with geometrically exact Kirchhoff-Love as well as Simo-Reissner beam elements. Also, it is independent of the temporal discretization and has been used alongwith statics, Lie group Generalized-Alpha as well as Brownian dynamics.
Load/time stepping
We either applied a fixed step size or an automatic step size adaption that is outlinedin the following. Starting from a given initial step size, a step is repeated with half of the previous stepsize if and only if the nonlinear solver did not converge within a prescribed number of iterations. Thisprocedure may be repeated until convergence is achieved (or until a given finest step size is reachedwhich aborts the algorithm). After four subsequent converging steps with a new, small step size, thestep size is doubled. Again, this is repeated until the initial step size is reached.
Nonlinear solver
The Newton-Raphson algorithm used throughout this work is based on the package NOXwhich is part of the Trilinos project [50]. Unless otherwise stated, the Euclidean norms of the dis-placement increment vector and of the residual vector are used as convergence criteria. Typically, thecorresponding tolerances were chosen as 10 − and 10 − , respectively. In some of the numerical exam-ples, an additional Newton step size control is applied. It restricts the step size such that a specifiedupper bound of the displacement increment per nonlinear iteration is not exceeded. In simple terms,it is meant to prevent any two points on two beams from moving too far and eventually crossing eachother without being detected from one iteration to the other. For this reason, the value for this upperbound is typically chosen as half of the beam radius. Linear solver
We use the algorithm UMFPACK [51] which is a direct solver for sparse linear systems ofequations based on LU-factorization and included in the package Amesos which is part of the Trilinosproject [50].
Parallel computing
The implementation of the novel methods supports parallel computing and is basedon the package Epetra which is part of the Trilinos project [50]. See Section 6.6 for details on thepartitioning of the problem in the context of the search algorithm applied to identify spatially proximateinteraction partners.
Post-processing and visualization
The computer program MATLAB [52] was used to post-process and plotsimulation data. All visualizations of the simulation results were generated using Paraview [53].
References [1] R. H. French, V. A. Parsegian, R. Podgornik, R. F. Rajter, A. Jagota, J. Luo, D. Asthagiri, M. K. Chaudhury, Y. M.Chiang, S. Granick, S. Kalinin, M. Kardar, R. Kjellander, D. C. Langreth, J. Lewis, S. Lustig, D. Wesolowski, J. S.Wettlaufer, W. Y. Ching, M. Finnis, F. Houlihan, O. A. Von Lilienfeld, C. J. Van Oss, T. Zemb, Long range interactionsin nanoscale science, Reviews of Modern Physics 82 (2) (2010) 1887–1944.[2] J. N. Israelachvili, Intermolecular and surface forces, 3rd Edition, Academic press, 2011.[3] V. A. Parsegian, Van der Waals forces: a handbook for biologists, chemists, engineers, and physicists, Cambridge UniversityPress, 2005.[4] G. Jelenic, M. A. Crisfield, Geometrically exact 3D beam theory: implementation of a strain-invariant finite element forstatics and dynamics, Computer Methods in Applied Mechanics and Engineering 171 (1–2) (1999) 141–171.[5] M. A. Crisfield, G. Jelenic, Objectivity of strain measures in the geometrically exact three-dimensional beam theory andits finite-element implementation, Proceedings of the Royal Society of London. Series A: Mathematical, Physical andEngineering Sciences 455 (1999) 1125–1147.[6] C. Meier, A. Popp, W. A. Wall, Geometrically Exact Finite Element Formulations for Slender Beams: Kirchhoff–LoveTheory Versus Simo–Reissner Theory, Archives of Computational Methods in Engineering 26 (1) (2019) 163–243.
7] H. Ohshima, A. Hyono, Electrostatic interaction between two cylindrical soft particles, Journal of Colloid and InterfaceScience 333 (1) (2009) 202–208.[8] R. P. Jaiswal, S. P. Beaudoin, Approximate Scheme for Calculating van der Waals Interactions between Finite CylindricalVolume Elements, Langmuir 28 (22) (2012) 8359–8370.[9] T. Stedman, D. Drosdoff, L. M. Woods, Van der Waals interactions between nanostructures: Some analytic results fromseries expansions, Physical Review A 89 (1) (2014) 012509.[10] H. Maeda, Y. Maeda, Orientation-Dependent London-van der Waals Interaction Energy between Macroscopic Bodies,Langmuir 31 (26) (2015) 7251–7263.[11] D. M. Dryden, J. C. Hopkins, L. K. Denoyer, L. Poudel, N. F. Steinmetz, W.-Y. Ching, R. Podgornik, A. Parsegian,R. H. French, van der Waals Interactions on the Mesoscale: Open-Science Implementation, Anisotropy, Retardation, andSolvent Effects, Langmuir 31 (37) (2015) 10145–10153.[12] P. Wriggers, G. Zavarise, On contact between three-dimensional beams undergoing large deflections, Communications inNumerical Methods in Engineering 13 (6) (1997) 429–438.[13] P. Litewka, The penalty and Lagrange multiplier methods in the frictional 3d beam-to-beam contact problem, Civil andEnvironmental Engineering Reports 1 (2005) 189–207.[14] D. Durville, Simulation of the mechanical behaviour of woven fabrics at the scale of fibers, International Journal of MaterialForming 3 (2) (2010) 1241–1251.[15] A. Kulachenko, T. Uesaka, Direct simulations of fiber network deformation and failure, Mechanics of Materials 51 (2012)1–14.[16] M. Chamekh, S. Mani-Aouadi, M. Moakher, Stability of elastic rods with self-contact, Computer Methods in AppliedMechanics and Engineering 279 (2014) 227–246.[17] A. Gay Neto, P. M. Pimenta, P. Wriggers, A master-surface to master-surface formulation for beam to beam contact. PartI: Frictionless interaction, Computer Methods in Applied Mechanics and Engineering 303 (2016) 400–429.[18] A. Konyukhov, O. Mrenes, K. Schweizerhof, Consistent Development of a Beam-To-Beam Contact Algorithm via theCurve-to-Solid Beam Contact - Analysis for the Nonfrictional Case, International Journal for Numerical Methods inEngineering 113 (7) (2018) 1108–1144.[19] O. Weeger, B. Narayanan, L. De Lorenzis, J. Kiendl, M. L. Dunn, An isogeometric collocation method for frictionlesscontact of Cosserat rods, Computer Methods in Applied Mechanics and Engineering 321 (2017) 361–382.[20] C. Meier, A. Popp, W. A. Wall, A finite element approach for the line-to-line contact interaction of thin beams witharbitrary orientation, Computer Methods in Applied Mechanics and Engineering 308 (2016) 377–413.[21] C. Meier, W. A. Wall, A. Popp, A unified approach for beam-to-beam contact, Computer Methods in Applied Mechanicsand Engineering 315 (2017) 972–1010.[22] C. Argento, A. Jagota, W. C. Carter, Surface formulation for molecular interactions of macroscopic bodies, Journal of theMechanics and Physics of Solids 45 (7) (1997) 1161–1183.[23] R. A. Sauer, S. Li, A contact mechanics model for quasi-continua, International Journal for Numerical Methods inEngineering 71 (8) (2007) 931–962.[24] R. A. Sauer, Multiscale modelling and simulation of the deformation and adhesion of a single gecko seta, ComputerMethods in Biomechanics and Biomedical Engineering 12 (6) (2009) 627–640.[25] R. A. Sauer, L. De Lorenzis, A computational contact formulation based on surface potentials, Computer Methods inApplied Mechanics and Engineering 253 (2013) 369–395.[26] H. Fan, S. Li, A three-dimensional surface stress tensor formulation for simulation of adhesive contact in finite deformation,International Journal for Numerical Methods in Engineering 107 (3) (2016) 252–270.[27] R. A. Sauer, J. C. Mergel, A geometrically exact finite beam element formulation for thin film adhesion and debonding,Finite Elements in Analysis and Design 86 (2014) 120–135.[28] M. G. Schmidt, A. E. Ismail, R. A. Sauer, A continuum mechanical surrogate model for atomic beam structures, Interna-tional Journal for Multiscale Computational Engineering 13 (5) (2015) 413–442.[29] S. Ahmadi, C. Menon, A new model for predicting fiber clumping phenomenon in bio-inspired dry adhesives, Journal ofAdhesion 90 (2) (2014) 135–155.[30] V. Negi, R. C. Picu, Mechanical behavior of cross-linked random fiber networks with inter-fiber adhesion, Journal of theMechanics and Physics of Solids 122 (2018) 418–434.[31] K. Prytz, Electrodynamics: The Field-Free Approach, Springer, 2015.[32] D. Langbein, Van der Waals attraction between cylinders, rods or fibers, Physik der kondensierten Materie 15 (1) (1972)61–86.[33] E. Reissner, On finite deformations of space-curved beams, Zeitschrift f¨ur Angewandte Mathematik und Physik (ZAMP)32 (6) (1981) 734–744.[34] J. C. Simo, A Finite Strain Beam Formulation. The Three-Dimensional Dynamic Problem. Part I, Computer Methods inApplied Mechanics and Engineering 49 (1985) 55–70.[35] J. C. Simo, L. Vu Quoc, A Three Dimensional Finite Strain Rod Model Part II: Computational Aspects, ComputerMethods in Applied Mechanics and Engineering 58 (1986) 79–116.[36] S. Timoshenko, LXVI. On the correction for shear of the differential equation for transverse vibrations of prismatic bars,The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 41 (245) (1921) 744–746.[37] C. Meier, M. J. Grill, W. A. Wall, A. Popp, Geometrically exact beam elements and smooth contact schemes for themodeling of fiber-based materials and structures, International Journal of Solids and Structures 154 (2018) 124–146.[38] F. Hilitski, A. R. Ward, L. Cajamarca, M. F. Hagan, G. M. Grason, Z. Dogic, Measuring cohesion between macromolecularfilaments one pair at a time: Depletion-induced microtubule bundling, Physical Review Letters 114 (13) (2015) 138102.
39] I. Borukhov, R. F. Bruinsma, W. M. Gelbart, A. J. Liu, Structural polymorphism of the cytoskeleton: A model of linker-assisted filament aggregation, Proceedings of the National Academy of Sciences of the United States of America 102 (10)(2005) 3673–3678.[40] D. Durville, Finite Element Simulation of Textile Materials at Mesoscopic Scale, in: Finite element modelling of textilesand textile composites, Saint-Petersbourg, Russian Federation, 2007, pp. 1–14.[41] C. J. Cyron, K. W. M¨uller, A. R. Bausch, W. A. Wall, Micromechanical simulations of biopolymer networks with finiteelements, Journal of Computational Physics 244 (2013) 236–251.[42] C. Meier, A. Popp, W. A. Wall, An objective 3D large deformation finite element formulation for geometrically exactcurved Kirchhoff rods, Computer Methods in Applied Mechanics and Engineering 278 (2014) 445–478.[43] R. A. Sauer, Enriched contact finite elements for stable peeling computations, International Journal for Numerical Methodsin Engineering 87 (6) (2011) 593–616.[44] D. Durville, Contact-friction modeling within elastic beam assemblies: an application to knot tightening, ComputationalMechanics 49 (6) (2012) 687–707.[45] M. J. Grill, C. Meier, W. A. Wall, Investigation of the Peeling and Pull-off Behavior of Adhesive Elastic Fibers via a NovelComputational Beam Interaction Model, arXiv preprint, arXiv:1907.13021 .[46] W. Gautschi, The use of rational functions in numerical quadrature, Journal of Computational and Applied Mathematics133 (1) (2001) 111–126.[47] P. Wriggers, Computational Contact Mechanics, 2nd Edition, Springer, 2006.[48] O. Br¨uls, A. Cardona, On the Use of Lie Group Time Integrators in Multibody Dynamics, Journal of Computational andNonlinear Dynamics 5 (2010) 1415–1555.[49] W. A. Wall, A. Popp, M. Kronbichler, M. Mayr, C. Meier, A.-T. Vuong, C. Ager, F. Br¨au, M. J. Grill, BACI: A multiphysicssimulation environment (2018).[50] M. A. Heroux, J. M. Willenbring, A New Overview of the Trilinos Project, Scientific Programming 20 (2) (2012) 83–88.[51] T. A. Davis, Algorithm 832: UMFPACK V4.3 – an unsymmetric-pattern multifrontal method, ACM Transactions onMathematical Software 30 (2) (2004) 196–199.[52] MATLAB R2017b, The MathWorks, Inc., Natick, Massachusetts, United States.[53] Paraview 5.4.1, Kitware, Inc., Clifton Park, New York, United States..[46] W. Gautschi, The use of rational functions in numerical quadrature, Journal of Computational and Applied Mathematics133 (1) (2001) 111–126.[47] P. Wriggers, Computational Contact Mechanics, 2nd Edition, Springer, 2006.[48] O. Br¨uls, A. Cardona, On the Use of Lie Group Time Integrators in Multibody Dynamics, Journal of Computational andNonlinear Dynamics 5 (2010) 1415–1555.[49] W. A. Wall, A. Popp, M. Kronbichler, M. Mayr, C. Meier, A.-T. Vuong, C. Ager, F. Br¨au, M. J. Grill, BACI: A multiphysicssimulation environment (2018).[50] M. A. Heroux, J. M. Willenbring, A New Overview of the Trilinos Project, Scientific Programming 20 (2) (2012) 83–88.[51] T. A. Davis, Algorithm 832: UMFPACK V4.3 – an unsymmetric-pattern multifrontal method, ACM Transactions onMathematical Software 30 (2) (2004) 196–199.[52] MATLAB R2017b, The MathWorks, Inc., Natick, Massachusetts, United States.[53] Paraview 5.4.1, Kitware, Inc., Clifton Park, New York, United States.