Multidimensional coupling: A variationally consistent approach to fiber-reinforced materials
Ustim Khristenko, Stefan Schu?, Melanie Krüger, Felix Schmidt, Barbara Wohlmuth, Christian Hesch
MMultidimensional coupling: A variationally consistentapproach to fiber-reinforced materials.
Ustim Khristenko a , Stefan Schuß b , Melanie Krüger b , Felix Schmidt b ,Barbara Wohlmuth a , Christian Hesch b ∗ a Faculty of Mathematics, Technical University of Munich, Garching, Germany b Chair of Computational Mechanics, University of Siegen, Siegen, Germany
Abstract
A novel coupling algorithm for the coupling of 1-dimensional fibers considered as beammodels to the surrounding 3-dimensional matrix material is presented. An overlappingdomain decomposition is used to embed the fibers within the matrix material. The coup-ling conditions acknowledge both, the forces and the moments of the beam model andtransfer them to the background material. A suitable static condensation procedure isapplied to remove the beam balance equations in the continuum setting and afterwardsin the coupling constraints in the discrete setting, reducing drastically the computationalcost. Eventually, a series of benchmark tests demonstrate the flexibility and robustnessof the proposed methodology.
Keywords : 1D-3D coupling, overlapping domain decomposition, immersed techniques,nonlinear beam model, condensation, second gradient material
Fiber reinforced materials are subjected to various physical mechanisms on differentscales, depending on the size, orientation and distribution within a suitable matrixmaterial. Applications contain, e.g. steel reinforced ultra-high performance concreteor fiber reinforced polymers. Beside such technical materials, many biological tissuesare reinforced by certain types of fibers. A model containing fibers fully resolved asa discretized continua, subsequently referred to as Cauchy continuum approach, is farout of the range of today’s computational capabilities. Therefore, we develop a modelcontaining the fibers as continuum degenerated 1-dimensional beams within the matrixmaterial in the sense of an overlapping domain decomposition method.This intermediate scale model can be considered as a material with microstructures.We refer here to dell’Isola et al. [11], where panthographic mechanisms as prototypicalmaterial with dedicated microstructure have been investigated. In Giorgio [17] a detailedanalysis demonstrates that a microscale model using a Cauchy continuum model can behomogenized using a macro second gradient model, see also dell’Isola et al. [10] fora suitable Piola homogenization of rod-like structures. More generally, we consider a ∗ Corresponding author. E-mail address: [email protected] a r X i v : . [ c s . C E ] F e b igher-gradient framework as proposed by Mindlin [30, 31], see also the work of Germain[15], Toupin [45, 46] as well as Eringen [13]. We refer also to Asmanoglo & Menzel [2, 3]for higher-order formulations used in the context of composites based on the early workof Spencer & Soldatos [43] and Soldatos [42]. For the application to Kirchhoff-Love shellelements see Schulte et al. [38] and Dittmann et al. [12] for the application of straingradient formulation on porous-ductile fracture.For the embedded beams, several models can be taken into account. Geometrically ex-act beam formulations using finite elements are presented foremost in the Simo-Reissnerbeam model, see Simo [39], Simo & Vu-Quoc [41] and Reissner [34]. For the parametriza-tion of the rotation, director interpolations can be employed to the Simo-Reissner beamtheory as shown in Romero & Armero [35] and Betsch & Steinmann [6, 5] and furtheradvanced in Eugster et al. [14]. Quaternions are used e.g. in McRobie & Lasenby [26],see also in Weeger et al. [47] for the application of quaternions in the context of iso-geometric collocation methods. In Meier et al. [28], geometrically exact Kirchhoff rodsare used and in Meier et al. [27] applied for the modeling of fiber based materials usingcontact algorithms.In contrast to the surface coupling of dedicated surfaces between beam and solid usingeither Dirichlet or Neumann conditions (e.g. within staggered algorithms or via Lagrangemultipliers using Mortar methods), overlapping domain decomposition allows for a flex-ible and efficient embedding of the fibers. This idea has been applied in the context offluid-structure interaction problems (FSI) using the terminus immersed technologies, see[33, 24, 23, 16, 19, 20]. From a methodological point of view, immersed methods for FSIare part of the so-called Fictitious Domain (FD) philosophy, introduced by Glowinskiet al. in [18] for the resolution of boundary value problems in complex geometrical set-tings. For the application on two elastic domains, we refer to Sanders & Puso [37]. Anoverlapping domain decomposition for beam-solid interaction problems using positionconstraints is presented in Steinbrecher et al. [44].To be more specific, the model proposed in this contribution derives the coupling condi-tions between a 3-dimensional matrix material and a dimensional reduced beam model.The model is derived from a full surface-to-volume beam-matrix coupling via Lagrangemultiplier, representing the coupling force. This full model is reduced via its Fourierexpansion as function of the angular coordinate in the beam cross-section. Taking intoaccount only the first term (average over cross-section) leads to the common displace-ment coupling constraint (see, e.g., [44]), which neglects the connection of the beamdirectors to the matrix material. Consideration of the second term allows to transferdilatation, shear and coupled stresses on the beam-matrix interface. It thus providesan additional constraint, which couples the beam directors with the deformation gra-dient, such that bending and torsion are transferred to the matrix as well. This alsoyields an additional tensor-valued Lagrange multiplier enforcing the matrix deformationto take into account the incompressibility of the beam cross-section. The constraintstypically involve the circular mean over the beam mantle (see, e.g., [9, 22]). Assumingenough regularity, we truncate the Taylor expansion of the matrix displacement field at2he beam centerline up to the linear term with respect to the beam radius r . Hence,the constraints are collapsed to the beam centerline, yielding a modeling error of or-der O ( r ). A subsequent condensation of the beam balance equations leads to a secondgradient formulation along the beam, which allows us to extract the so-called coupledforces.The paper is structured as follows. In Section 2 we present the fundamental formulationsfor the first and second gradient continuum as well as the beam. In Section 3, we derivethe coupling terms for the overlapping domain decomposition method and the conden-sation of the beam balance equations. The spatial discretization is given in Section 4,followed by a series of representative examples in Section 5. Eventually, conclusions aresummarised in Section 6. This section gives a brief summary on the used notation. A single contraction of twovectors will be understood as [ aaa · bbb ] = a i b i , where the Einstein summation convention onrepeated indices is used. For two second order tensors it holds [ AAA BBB ] ij = A ik B kj and thedouble contraction reads [ AAA : BBB ] = A ij B ij . Next, we define the gradient with respect tothe reference ∇ ( • ) of a vector field aaa and of a second-order tensor field AAA as[ ∇ aaa ] iJ = ∂ [ aaa ] i ∂ [ XXX ] J and [ ∇ AAA ] iJK = ∂ [ AAA ] iJ ∂ [ XXX ] K . (1)For the divergence operator it follows[ ∇ · AAA ] i = ∂ [ AAA ] iJ ∂ [ XXX ] J and [ ∇ · A ] iJ = ∂ [ A ] iJK ∂ [ XXX ] K . (2)The triple contraction for two third order tensors is given via [ A ... B ] = A ijk B ijk . Forthe double contraction we define[ A : AAA ] i = A iJK A JK and [ A T : AAA ] K = A KiJ A iJ . (3)The axial vector of a skew symmetric 3 × AAA ] i = − (cid:15) ijk [ AAA ] jk and the spin of a 3-dimensional vector [[ aaa ] × ] ij = − (cid:15) ijk [ aaa ] k such that aaa × bbb = [ aaa ] × bbb , i.e.axl AAA = 12 A − A A − A A − A and [ aaa ] × = − a a a − a − a a . (4)Here, we have made use of the Levi-Civita permutation tensor (cid:15) ijk with (cid:15) ijk (cid:15) ijl = 2 δ kl and the Kronecker symbol δ kl = [ III ] kl , such that axl [ aaa ] × = aaa . Hence, [ aaa ] × belongs to thevector space of skew-symmetric matricesso(3) = { AAA ∈ R × | AAA + AAA T = 000 } , (5)3he image of the matrix exponential map of so(3) is the special orthogonal group SO(3),i.e. exp : so(3) → SO(3) andSO(3) = { RRR ∈ GL(3) | RRR T RRR = RRR RRR T = III, det(
RRR ) = 1 } , (6)where GL(3) is the general linear group of 3 × ××× , defined as [ AAA ×××
BBB ] iJ = (cid:15) imn (cid:15) JP Q [ AAA ] mP [ BBB ] nQ forthe two-point second order tensors AAA and
BBB .Moreover, for any function f ( s ), s ∈ { , L } , we introduce the following notation:[[ f ]] | L := f ( L ) − f (0) . (7) We start with a short summary of non-linear continuum mechanics. Therefore, considera Lipschitz bounded domain Ω ⊂ R in its reference configuration with boundary ∂ Ω and outward unit normal NNN . The actual configuration Ω ⊂ R with boundary ∂ Ωand outward unit normal nnn is related to the reference configuration by a deformationmapping ϕϕϕ : Ω → R , such that Ω = ϕϕϕ (Ω ). Material points are labelled by XXX ∈ Ω with corresponding actual position xxx = ϕϕϕ ( XXX ). Kinematics.
For the matrix material, we introduce the deformation gradient as secondorder tensor field
FFF : Ω → R × , such that FFF = ∇ ϕϕϕ ( XXX ) , (8)which maps the infinitesimal vector d XXX at XXX ∈ Ω to the infinitesimal vector d xxx in theactual configuration. Next, the second order tensor field HHH : Ω → R × is introduced,where HHH = cof
FFF denotes the cofactor of
FFF , defined as follows
HHH = 12
FFF ×××
FFF , (9)which maps the infinitesimal oriented area element d
AAA = NNN ( XXX ) d A to the infinitesimaloriented area element d aaa = nnn ( XXX ) d a in the actual configuration. Eventually, the scalarfield J : Ω → R is introduced, where J = det FFF denotes the determinant, defined as J = 16 ( FFF ×××
FFF ) :
FFF , (10)which relates the infinitesimal volume element d V to the corresponding infinitesimalvolume element d v in the actual configuration. This last equation completes the set ofkinematic relations, to be used to define the strain energy of a non-linear material. Notethat we omit in the following the time dependency of the equations and postulate a zerodensity, as we focus here on static problems.4 arge strain elasticity. We assume the existence of a strain energy function of theform Ψ := Ψ(
FFF , cof
FFF , det
FFF ) , (11)where Ψ : R × × R × × R → R . This allows us to formulate the total stored energy inthe material configuration along with the external contributionsΠ int = Z Ω Ψ( FFF , cof
FFF , det
FFF ) d V, Π ext = − Z Ω BBB ext · ϕϕϕ d V − Z Γ σ TTT ext · ϕϕϕ d A. (12)Here, BBB ext are given body forces, and
TTT ext surface stresses on the Neumann boundary,where the boundary ∂ Ω is decomposed in Dirichlet boundary Γ ϕ with prescribed dis-placements ¯ ϕϕϕ and Neumann boundary Γ σ with the standard relationship ∂ Ω = Γ ϕ ∪ Γ σ and Γ ϕ ∩ Γ σ = ∅ . For ease of presentation, only dead loads with associated potentialenergy are taken into account, other choices are possible and do not restrict our ap-proach. Moreover, to ensure frame invariance, i.e., objectivity, we require that Ψ mustbe independent of rotational components of the deformation gradient and the cofactor.The principle of virtual work now readsdd (cid:15) (Π int ( ϕϕϕ (cid:15) ) + Π ext ( ϕϕϕ (cid:15) )) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 = 0 , (13)where ϕϕϕ (cid:15) ( XXX ) = ϕϕϕ ( XXX ) + (cid:15) δϕϕϕ ( XXX ). Here, we have made use of the configuration space ϕϕϕ ∈ C and the space of admissible variations δϕϕϕ ∈ V , defined as C = { ϕϕϕ ∈ (cid:16) H (Ω ) (cid:17) | J ( XXX ) > ∀ XXX ∈ Ω , ϕϕϕ = ¯ ϕϕϕ on XXX ∈ Γ ϕ } , (14)and V = { δϕϕϕ ∈ (cid:16) H (Ω ) (cid:17) | δϕϕϕ = 000 on XXX ∈ Γ ϕ } . (15)The virtual work of the internal contributions reads δ Π int ( ϕϕϕ ) = Z Ω PPP : ∇ δϕϕϕ d V, (16)where the first Piola-Kirchhoff stress tensor is given by PPP = ∂ Ψ ∂FFF + ∂ Ψ ∂HHH ××× FFF + ∂ Ψ ∂J HHH. (17)Now, we can reformulate the principle of virtual work (13) as Z Ω PPP : ∇ δϕϕϕ d V − Z Ω BBB ext · δϕϕϕ d V − Z Γ σ TTT ext · δϕϕϕ d A = 0 ∀ δϕϕϕ. (18) Remark 1.
The well-known Mooney-Rivlin constitutive model Ψ MR ( FFF , HHH, J ) = α FFF : FFF + β HHH : HHH + f ( J ) , (19) where f ( J ) = − α ln( J ) − β J + λ/ J − − α + β ) and α, β, λ are non-negativematerial parameters, is a suitable polyconvex constitutive model. In particular, Ψ MR isconvex with respect to its 19 variables, i.e. the components of FFF , HHH and J , which can beseen by using the Hessian of Ψ MR , see [7]. econd gradient material. Next we assume, that the second gradient of the deforma-tion ∇ FFF with respect to the reference configuration can be taken into account as well,see, among many others, [21, 31]. Thus, the strain energy function can be written asΨ := Ψ( ∇ FFF , FFF , cof
FFF , det
FFF ) (20)such that the internal virtual work reads δ Π int ( ϕϕϕ ) = Z Ω PPP : ∇ δϕϕϕ + P ... ∇ δϕϕϕ d V, (21)where P is a third order stress tensor, conjugated to ∇ δϕϕϕ . To identify and collectthe corresponding boundary conditions, we apply twice integration by parts and thedivergence theorem (see Appendix A.1 for details) yielding δ Π int ( ϕϕϕ ) = Z Ω ∇ · ( ∇ · P − PPP ) · δϕϕϕ d V + Z ∂ Ω δϕϕϕ · ( PPP − ∇ · P ) NNN + ∇ δϕϕϕ : ( P · NNN ) d A. (22)Using the orthogonal decomposition ∇ ⊥ · ( • ) = ∇ ( • ) : ( NNN ⊗ NNN ) and ∇ k · ( • ) = ∇ ( • ) :( III − NNN ⊗ NNN ), we obtain after some further technical steps (c.f. Appendix A.1) δ Π int ( ϕϕϕ ) = Z Ω ∇ · ( ∇ · P − PPP ) · δϕϕϕ d V + Z ∂ Ω δϕϕϕ · ( PPP − ∇ · P ) NNN d A − Z ∂ Ω h δϕϕϕ · ( K ( P NNN ) NNN + ∇ k · ( P NNN )) − ∇ ⊥ δϕϕϕ : ( P NNN ) i d A + Z ∂ Ω δϕϕϕ · ( P : ( ˆ NNN ⊗ NNN )) d S. (23)where ˆ NNN is the normal to ∂ Ω and the tangent to ∂ Ω . and ∂ Ω is defined by theunion of the boundary curves of the boundary surface patches. Moreover, K = −∇ k · NNN is the curvature of the surface. Since the fibers are embedded within the polymermatrix, the application of external boundaries for standard industrial components withgradient contributions as they arise in the second and third line of (23) is not feasible.Equilibrating this with the external contributions yields the principle of virtual work Z Ω PPP : ∇ δϕϕϕ + P ... ∇ δϕϕϕ d V − Z Ω BBB ext · δϕϕϕ d V − Z Γ σ TTT ext · δϕϕϕ d A = 0 ∀ δϕϕϕ, (24)The configuration space ϕϕϕ ∈ C and the space of admissible variations δϕϕϕ ∈ V are nowdefined as C = { ϕϕϕ ∈ (cid:16) H (Ω ) (cid:17) | J ( XXX ) > ∀ XXX ∈ Ω , ϕϕϕ = ¯ ϕϕϕ on Γ ϕ } , (25)and V = { δϕϕϕ ∈ (cid:16) H (Ω ) (cid:17) | δϕϕϕ = 000 on Γ ϕ } , (26)where we omit higher order Dirichlet boundary conditions. There are only a few exam-ples for higher-order strain energy formulations known, we refer to dell’Isola et al. [10]where a Piola homogenization procedure is used to derive a suitable formulation.6 .2 Continuum degenerate beam formulation In this section, we degenerate the general continuum mechanical framework as intro-duced above to a beam formulation. As we intend to embed fibers with a length-to-diameter ratio of 20 as standard for e.g. fiber reinforced polymers, it is reasonable thatwe restrict the kinematics of the 3-dimensional continuum along the fiber direction toa beam-like kinematic. In particular, we use the theory of geometrically exact beams,also known as Cosserat beam, introduced in [8], see also [36].
Beam kinematic.
Let us consider a beam as a 3-dimensional body, occupying ˜Ω ⊂ R with the following kinematical ansatz for the position field in the reference configura-tion ˜ XXX ( θ α , s ) = ˜ ϕϕϕ ( s ) + θ α DDD α ( s ) , (27)parametrized in terms of s ∈ [0 , L ] along the center line of the beam with length L .Moreover, an orthonormal triad [ DDD , DDD , DDD ] with convective coordinates ( θ , θ , s ) ∈ ˜Ω , is introduced, where DDD α , α = 1 , xxx ( θ α , s ) = ˜ ϕϕϕ ( s ) + θ α ddd α ( s ) . (28)Here, the orthonormal triad ddd i is related to the reference triad via the rotation tensor˜ RRR ∈ SO (3), i.e., ˜ RRR = ddd i ⊗ DDD i . The deformation gradient now reads [32]˜ FFF = ˜ ϕϕϕ ⊗ DDD + ddd α ⊗ DDD α + θ α ddd α ⊗ DDD , (29)where ( • ) represents the derivative with respect to s . Following the arguments in Ap-pendix A.2, the last equation can be rewritten as follows˜ FFF = ˜
RRR (cid:16)
ΓΓΓ ⊗ DDD + [ KKK ] × θ α DDD α ⊗ DDD + III (cid:17) , (30)using the usual strain measure (c.f. [6])ΓΓΓ = ˜ RRR T ˜ ϕϕϕ − DDD , [ KKK ] × = ˜ RRR T ˜ RRR . (31)The strain measure ΓΓΓ is known as the axial-shear strain vector, whereas the curvaturerepresented by KKK is called the torsional-bending strain vector. We refer to [14] fora detailed analysis of the contravariant components of the effective curvature. For thestrain energy to be defined subsequently, it is useful to introduce the polar decompositionof the deformation gradient via˜
FFF = ˜
RRR ˜ UUU , ˜ UUU = III + BBB ⊗ DDD , BBB = ΓΓΓ + ( KKK × θ α DDD α ) . (32)7ith regard to (6), the rotation tensor can be defined in terms of a rotation vector φφφ ∈ R via the exponential map ˜ RRR ( φφφ ) = e [ φφφ ] × . A closed form expression is given byRodrigues formula ˜ RRR ( φφφ ) = III + sin k φφφ kk φφφ k [ φφφ ] × + 12 sin ( k φφφ k / k φφφ k / ! [ φφφ ] × . (33)According to [29], the relationship between the spatial variation of δφφφ and δ ˜ RRR is givenby δ ˜ RRR = [ δφφφ ] × ˜ RRR . Strain energy for large deformation beam.
With regard to (32), we can write ourstrain energy function as Ψ(
FFF , HHH, J ) ˆ= Ψ(
BBB (ΓΓΓ , KKK )). In [32], we have shown that theMooney-Rivlin material model presented in (19) can be rewritten after some technicalcalculations in terms of the beam strain measures as˜Ψ MR ( BBB (ΓΓΓ , KKK )) = Z A ( s ) [ α ( BBB · BBB + 2
BBB · DDD + 3)+ β (2 BBB · BBB − ( DDD × BBB ) · ( DDD × BBB ) + 4
BBB · DDD + 3)+ f ( BBB · DDD + 1)] d θ, (34)where A ( s ) denotes the local cross section of the beam, such that the internal energycan be written as a line integral˜Π int = Z C ˜Ψ MR ( BBB (ΓΓΓ , KKK )) d s. (35)The most simplest model is written directly in terms of ΓΓΓ and KKK :˜Ψ(ΓΓΓ , KKK ) = 12 ΓΓΓ · ( K ΓΓΓ) + 12
KKK · ( K KKK ) , (36)where K = Diag[ GA , GA , EA ] and K = Diag[ EI , EI , GJ ]. The correspondingconstitutive relations in the reference and actual configuration reads˜ NNN = ∂ ˜Ψ ∂ ΓΓΓ , ˜ MMM = ∂ ˜Ψ ∂KKK , ˜ nnn = ˜ RRR ˜ NNN , ˜ mmm = ˜ RRR ˜ MMM . (37)The principle of virtual work of the beam is given in Appendix A.3. Introducing externalforces ¯˜ nnn and the couples ¯˜ mmm , where we assume the existence of an external energy˜Π ext = − Z C (cid:16) ¯˜ nnn · ˜ ϕϕϕ + ¯˜ mmm · φφφ (cid:17) d s, (38)we obtain from integration by parts the local form − Z C δ ˜ ϕϕϕ · (cid:16) ˜ nnn + ¯˜ nnn (cid:17) d s = 0 , − Z C δφφφ · (cid:16) ˜ mmm + ˜ ϕϕϕ × ˜ nnn + ¯˜ mmm (cid:17) d s = 0 , (39)8hich are the classical balance equations in the actual configuration, c.f. Antman [1].Additionally, we obtain information on the boundaries of the beam, i.e., on the endpointsas [[ δ ˜ ϕϕϕ e · (˜ nnn − nnn eext )]] | L = 0 and [[ δφφφ e · ( ˜ mmm − mmm eext )]] | L = 0 , (40)where nnn eext and mmm eext are appropriate external forces and moments at the endpoints. The method of weighted residual.
To avoid shear locking, a mixed, Hu-Washizutype method is employed. Therefore, we introduce additional, independent fields for theresultant spatial contact force vector n and contact torque vector m as well as the spatialaxial-shear strain vector g and the torsional-bending strain vector k . Now we assume astrain energy function of the form˜Ψ HW = ˜Ψ( g , k ) + n · ( γγγ − g ) + m · ( kkk − k ) (41)where γγγ = ˜ RRR
ΓΓΓ and kkk = ˜
RRR KKK . The principle of virtual work yields now Z C " − δ ˜ ϕϕϕ · (cid:16) n + ¯˜ nnn (cid:17) − δφφφ · (cid:16) m + ˜ ϕϕϕ × n + ¯˜ mmm (cid:17) + δ g · ∂ ˜Ψ( g , k ) ∂ g − n ! + δ k · ∂ ˜Ψ( g , k ) ∂ k − m ! + δ n · ( γγγ − g ) + δ m · ( kkk − k ) d s = 0 , ∀ δ ˜ ϕϕϕ, δφφφ, δ g , δ k , δ n , δ m . (42)Assuming that the last two variations with respect to δ n and δ m in (42) are fulfilledlocally, we can rewrite the mixed formulation as follows Z C " − δ ˜ ϕϕϕ · (cid:16) n + ¯˜ nnn (cid:17) − δφφφ · (cid:16) m + ˜ ϕϕϕ × n + ¯˜ mmm (cid:17) + δ g · ˜ RRR ∂ ˜Ψ(ΓΓΓ , KKK ) ∂ ΓΓΓ − n ! + δ k · ˜ RRR ∂ ˜Ψ(ΓΓΓ , KKK ) ∂KKK − m ! d s = 0 , (43)for ˜ ϕϕϕ , φφφ , n and m from H ( C ), and for all test functions δ ˜ ϕϕϕ , δφφφ , δ g and δ k from L ( C ). In this work, starting from a surface-to-volume coupling formulation for the matrix/beamsystem, we derive a reduced surrogate model, where the new coupling constraints definedon the beam centerline involve the deformation gradient and transfer both the linearforces and the moments of the beam to the matrix. The reduced model is based on thefollowing assumptions: 91.
Overlapping domains . The matrix continuum is extended inside the beam domain.This yields a modeling error due to additional stiffness, which depends on the fiberradius r and the ratio of the matrix stiffness to the fiber stiffness. Such type ofmodeling error was studied for diffusion problem in [22].A2. Form of the coupling force . We make an assumption on the form of the couplingforce (Lagrange multiplier), such that, applied to a beam cross-section, it presentsthe mean force, couple stresses, shear and dilatation forces.A3.
Sufficient regularity . We assume additional regularity of the solution, such that thetraces of deformation and of its gradient are defined on the beam centerline. Wewill then truncate Taylor expansion of the deformation in the cross-section withrespect to the radius r . Considering all terms up to the linear part yields a modelingerror of order O ( r ). C Ω DDD DDD DDD ϕϕϕ ddd ddd ddd Figure 1.
Deformation of the matrix/beam coupled system.
Let us consider the matrix/beam system illustrated in Figure 1. The coupling conditionsare defined via Lagrange multipliers on the beam mantle Γ C and at the end faces A , A L (see Figure 2). The work associated to the coupling forces over the whole matrix/beaminterface reads Π Γ = Z ∂ ˜Ω µµµ · ( ϕϕϕ − ˜ xxx ) d A = Π C + Π A , (44)where Π C = Z C Z C ( s ) µµµ · ( ϕϕϕ − ˜ xxx ) d C d s, Π A = Z A s µµµ · ( ϕϕϕ − ˜ xxx ) d A (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) L . (45)Above, A s and C ( s ) are respectively the beam cross-section and its boundary corre-sponding to the arc length s . The area and the circumference of the cross-section arerespectively denoted as | A | := π r and | C | := 2 π r . Here and further, for the sake ofshortness, we use the following abuse of notations: for any function f ( XXX ) defined in Ω, f ( θ α , s ) := f ( ˜ XXX ( θ α , s )) means the restriction of f onto the beam domain ˜Ω in convectivecoordinates; and f ( s ) := f (˜ ϕϕϕ ( s )), s ∈ C , is understood in the sense of the trace f | C of f on the beam centerline. 10 A A L Γ C C ( s ) r DDD Figure 2.
3D beam reference configuration.
The Lagrange multiplier µµµ , physically interpreted as the interface load, is defined onthe beam mantle Γ C = { ( r cos θ, r sin θ, s ) | s ∈ [0 , L ] , θ ∈ [0 , π ] } , where r is the beamradius and θ = arctan( θ /θ ) denotes the angle along the beam cross-section circumfer-ence. Thus, µµµ = µµµ ( θ, s ) can be decomposed in Fourier series as periodic function of theangle θ ∈ [0 , π ]. Here, we base our reduced model on truncation of this series after thefirst two terms. That is, we assume the interface load in the form: µµµ ( θ, s ) = ¯ µµµ ( s ) + ˜ µµµ ( s ) cos θ + ˜ µµµ ( s ) sin θ, (46)where ¯ µµµ ( s ) = π r R π µµµ ( θ, s ) d θ is the mean cross-sectional load, and ˜ µµµ α , α = 1 ,
2, arerespectively the first cosine and sine Fourier coefficients of µµµ . These two terms composethe load fluctuation ˜ µµµ ( θ, s ) := ΣΣΣ( s ) NNN ( θ ) = ˜ µµµ ( s ) cos θ + ˜ µµµ ( s ) sin θ, (47)where ΣΣΣ = ˜ µµµ α ⊗ DDD α denotes the interface stress tensor, and NNN = DDD cos θ + DDD sin θ is the unit outer normal to the interface in the reference configuration. Note that R π ˜ µµµ ( θ, s ) d θ = 0.While ¯ µµµ corresponds to the resulting force in the cross-section and, therefore, couplesposition of the beam center-line to the matrix, ˜ µµµ contains the coupled stresses, shearsand hydrostatic pressure, and is thus responsible for transition of the beam bending andtorsion to the matrix. In the following remark, we discuss in detail the structure of ˜ µµµ and ΣΣΣ. Remark 2.
Due to the definition (47) , the stress tensor
ΣΣΣ has only six degrees offreedom. In particular, denoting
DDD := DDD j ⊗ eee j , we write ΣΣΣ in the following form,involving coupled stresses, shears and dilatation:
ΣΣΣ = ˜
RRR
DDD ˜ µ n, ˜ µ n, − ˜ µ τ, µ n, + ˜ µ τ, ˜ µ n, − µ τ, µ τ, DDD T , (48)11 here ˜ µ n, , ˜ µ n, and ˜ µ n, are associated to the dilatation and the shear of the cross-section, ˜ µ τ, and ˜ µ τ, to the bending, and ˜ µ τ, to the torsion. Let ˜ µµµ n = ˜ µ n,j ddd j and ˜ µµµ τ = ˜ µ τ,j ddd j = axl (cid:18) ΣΣΣ ˜
RRR T (cid:19) . Then, we have ΣΣΣ = (
PPP α ˜ µµµ τ + QQQ α ˜ µµµ n ) ⊗ DDD α , (49) where PPP = ddd ⊗ ddd − ddd ⊗ ddd , QQQ = ddd ⊗ ddd + ddd ⊗ ddd , PPP = − ddd ⊗ ddd + 2 ddd ⊗ ddd , QQQ = ddd ⊗ ddd + ddd ⊗ ddd . (50) Moreover, the tensor [˜ µµµ n ] s := ˜ RRR T QQQ α ˜ µµµ n ⊗ DDD α = ˜ µ n, DDD ⊗ DDD + ˜ µ n, DDD ⊗ DDD + ˜ µ n, ( DDD ⊗ DDD + DDD ⊗ DDD ) (51) is symmetric. The following lemma gives the representation of the energy contribution due to thecoupling on the beam mantle Γ C . Lemma 1. If µµµ is in form (46) , the coupling energy term (45) writes Π C = Z C µµµ · ( ϕϕϕ − ˜ xxx ) | C | d s = Z C ¯ µµµ · (¯ ϕϕϕ − ˜ ϕϕϕ ) | C | d s + 12 Z C ΣΣΣ : (cid:16)
ΦΦΦ − ˜ RRR (cid:17) | C | d s = Z C ¯ µµµ · (¯ ϕϕϕ − ˜ ϕϕϕ ) | C | d s + 12 Z C ˜ µµµ α · (ˆ ϕϕϕ α − ddd α ) | C | d s, (52) with ¯ ϕϕϕ = 1 | C | Z C ( s ) ϕϕϕ d C, ΦΦΦ = 2 | C | Z C ( s ) ϕϕϕ ⊗ NNN d C, (53) and Fourier cosine and sine coefficients respectively ˆ ϕϕϕ = 2 | C | Z C ( s ) ϕϕϕ cos θ d C, ˆ ϕϕϕ = 2 | C | Z C ( s ) ϕϕϕ sin θ d C. (54) Proof.
Substitution of (46) and (47) to (45) yieldsΠ C = Z C Z C ( s ) (¯ µµµ + ΣΣΣ NNN ) | {z } µµµ · ϕϕϕ − (cid:16) ˜ ϕϕϕ + ˜ RRR NNN r (cid:17)| {z } ˜ xxx ! d C d s. (55)12iven NNN = DDD cos θ + DDD sin θ , we have R C ( s ) NNN d C = 000 and R C ( s ) NNN ⊗ NNN d C = π DDD α ⊗ DDD α .In particular, π Z NNN ⊗ NNN d θ = DDD ⊗ DDD π Z cos θ d θ | {z } = π + DDD ⊗ DDD π Z sin θ d θ | {z } = π + [ DDD ⊗ DDD + DDD ⊗ DDD ] π Z sin θ cos θ d θ | {z } =0 = π DDD α ⊗ DDD α . (56)Thus, taking into account that ΣΣΣ DDD α ⊗ DDD α = ΣΣΣ, (55) writesΠ C = Z C ¯ µµµ · Z C ( s ) ( ϕϕϕ − ˜ ϕϕϕ ) d C d s + Z C ΣΣΣ : Z C ( s ) (cid:16) ϕϕϕ − ˜ RRR NNN r (cid:17) ⊗ NNN d C d s = Z C ¯ µµµ · Z C ( s ) ϕϕϕ d C − ˜ ϕϕϕ | C | d s + Z C ΣΣΣ : Z C ( s ) ϕϕϕ ⊗ NNN d C − π r ˜ RRR d s. (57)Hence if follows (52) . Finally, (52) is obtained by substitution of ΣΣΣ = ˜ µµµ α ⊗ DDD α into(52) . Remark 3.
Let us note that if ϕϕϕ ∈ H (Ω) , then according to the trace theorem, ¯ ϕϕϕ ∈ H ( C ) and thus is in L ( C ) . Likewise, ΦΦΦ ∈ L ( C ) × . Thus, for ˜ ϕϕϕ ∈ H ( C ) , ˜ RRR ∈ H ( C ) × , ¯ µµµ ∈ L ( C ) and ΣΣΣ ∈ L ( C ) × , the integral in (52) is defined. Thus, we obtained the energy term (52) coupling the 1D Cosserat beam to the 3Dmaterial matrix to the circular integrals of the matrix deformation ϕϕϕ on the lateralarea Γ C of the beam. Such approach is common for 1D-3D coupling models (see, e.g.,[9, 22]), since the trace of ϕϕϕ ∈ H (Ω) on the beam centerline C is not well-defined.However, in practice, computation of the circular integrals can be subtle. Moreover,Galerkin approximation using smooth basis functions allows us to deal with the tracesof the matrix on the beam centerline. In this paper, we are not addressing mathematicalaspects of the variational problem, but rather focus on the computational model. So,considering Galerkin projection, we assume necessary regularity and make use of thefollowing lemma. Lemma 2.
Let ϕϕϕ be regular enough, such that the traces ϕϕϕ c ( s ) := ϕϕϕ (˜ ϕϕϕ ( s )) and FFF c ( s ) := FFF (˜ ϕϕϕ ( s )) are square integrable. Then, if µµµ is in form (46) , the coupling energy term (45) writes Π C = Z C (cid:20) ¯ µµµ · ( ϕϕϕ c − ˜ ϕϕϕ ) + 12 ΣΣΣ : (cid:16) FFF c − ˜ RRR (cid:17)(cid:21) | C | d s + O ( r ) . (58)13 roof. Using Fourier representation of the trace of ϕϕϕ onto the plane containing C ( s ) andthe formal power series of the exponential, we can write ϕϕϕ (cid:12)(cid:12)(cid:12) C ( s ) = Z R ˆ ϕϕϕ ( ωωω, s ) e i ωωω · NNN ( θ ) r d ωωω = ∞ X k =0 Z R ˆ ϕϕϕ ( ωωω, s ) (i ωωω · NNN ( θ ) r ) k k ! d ωωω = ϕϕϕ c ( s ) + FFF c ( s ) · NNN ( θ ) r + O ( r ) . (59)Then, substitution of (59) to the equation (52) of Lemma 1 yields the statement.Under the assumption that the fiber is thin, we neglect in (58) the terms with powersof r greater than one, keeping only the linear term.Eventually, one can consider the end faces A s = { ( ρ cos θ, r sin θ, s ) | ρ ∈ [0 , r ] , θ ∈ [0 , π ] } , s ∈ { , L } . Similar to (46)-(47) with replacing a circular mean by a cross-sectional one, we treat the Lagrange multipliers µµµ on the end faces A and A L : µµµ ( ρ, θ, s ) = ¯ µµµ ( ρ, s ) + ΣΣΣ( ρ, s ) NNN ( θ ) , (60)where ¯ µµµ ( ρ, s ) = π ρ R π µµµ ( ρ, θ, s ) d θ , and ΣΣΣ has the same form as in (48). Then, usingthe same arguments as in Lemma 1 and Lemma 2, the coupling energy term at the endfaces writesΠ A = (cid:20)(cid:20)Z r (cid:20) ¯ µµµ · ( ϕϕϕ c − ˜ ϕϕϕ ) + 12 ΣΣΣ : (cid:16) FFF c − ˜ RRR (cid:17)(cid:21) | C ( ρ ) | d ρ (cid:21)(cid:21)(cid:12)(cid:12)(cid:12)(cid:12) L + O ( r ) . (61)However, since this term scales as r , we also neglect it. Remark 4.
Remark that further consideration of higher order terms in the expan-sion (59) will lead to higher gradients models, reducing the modeling error.
Thus, the virtual work associated to the kinematic constraints is obtained as variationof (58) in the following form: δ Π Γ ≈ δ Π C ≈ Z C h δ ¯ µµµ · ( ϕϕϕ c − ˜ ϕϕϕ ) + ¯ µµµ · ( δϕϕϕ − δ ˜ ϕϕϕ ) i | C | d s + Z C (cid:20)
12 ΣΣΣ : (cid:16) ∇ δϕϕϕ − [ δφφφ ] × ˜ RRR (cid:17) + 12 δ ΣΣΣ : (cid:16)
FFF c − ˜ RRR (cid:17)(cid:21) | C | d s, (62)where we note that ΣΣΣ : (cid:16) [ δφφφ ] × ˜ RRR (cid:17) = ˜ µµµ τ · δφφφ .14 .2 Total coupled system Finally, we obtain the principle of virtual work for the coupled matrix/beam system δ Π int + δ Π ext + δ ˜Π int + δ ˜Π ext + δ Π Γ = 0 , (63)in terms of (18) and (43) along with the virtual work of the coupling forces in (62).Regarding the matrix material, i.e., the variation with respect to ϕϕϕ , we obtain Z Ω " ∂ Ψ ∂FFF : ∇ δϕϕϕ − BBB ext · δϕϕϕ d V − Z Γ σ TTT ext · δϕϕϕ d A + Z C (cid:18) ¯ µ ¯ µ ¯ µ · δϕϕϕ + 12 ΣΣΣ : ∇ δϕϕϕ (cid:19) | C | d s = 0 , (64)and the beam contributions yield − Z C (cid:16) n + ¯˜ nnn + ¯ µµµ | C | (cid:17) · δ ˜ ϕϕϕ d s = 0 , (65) − Z C (cid:16) m + ˜ ϕϕϕ × n + ¯˜ mmm + ˜ µµµ τ | C | (cid:17) · δφφφ d s = 0 , (66)with boundary conditions [[ δ ˜ ϕϕϕ · ( n − nnn eext )]] | L = 0 and [[ δφφφ · ( m − mmm eext )]] | L = 0, and sup-plemented by the weak constitutive equations Z C δ g · ˜ RRR ∂ ˜Ψ(ΓΓΓ , KKK ) ∂ ΓΓΓ − n ! d s = 0 , (67) Z C δ k · ˜ RRR ∂ ˜Ψ(ΓΓΓ , KKK ) ∂KKK − m ! d s = 0 , (68)as well as the coupling constraints Z C δ ¯ µ ¯ µ ¯ µ · ( ϕϕϕ c − ˜ ϕϕϕ ) | C | d s = 0 , (69)12 Z C δ ΣΣΣ : (cid:16)
FFF c − ˜ RRR (cid:17) | C | d s = 0 . (70) The above system is large, as we have to deal with the degrees of freedom of the ma-trix material (three per node in the 3-dimensional continuum), and the 21 unknownsincluding 9 Lagrange multipliers per node along the beam center line. Thus, we aim at atwo step static condensation procedure: First, we condense (65)-(66) and eliminate thecorresponding Lagrange multipliers in the continuous system. The remaining equationsfor the constraints and the constitutive laws for the beam are condensed in the discretesetting, such that finally only the matrix degrees of freedom remain and the beam isfully condensed. 15 heorem 1 ( Condensed system ) . Let V = { v ∈ H (Ω ) , s.t. DDD α · ∇ v | C ∈ L ( C ) } .The system (64) - (70) can be reduced to the following system of the unknowns ϕϕϕ ∈ V , ˜ ϕϕϕ, φφφ ∈ (cid:16) H ( C ) (cid:17) , n , m , ˜ µµµ n ∈ (cid:16) L ( C ) (cid:17) : Z Ω (cid:16) PPP : ∇ δϕϕϕ − BBB ext · δϕϕϕ (cid:17) d V − Z Γ σ TTT ext · δϕϕϕ d A + Z C P ... ∇ δϕϕϕ + (cid:16) PPP n + PPP m + PPP g + 12 FFF c [˜ µµµ n ] s | C | (cid:17) : ∇ δϕϕϕ − ¯˜ nnn · δϕϕϕ ! d s − (cid:20)(cid:20) nnn eext · δϕϕϕ + 12 ( PPP α mmm eext ⊗ DDD α ) : ∇ δϕϕϕ (cid:21)(cid:21)(cid:12)(cid:12)(cid:12)(cid:12) L = 0 , Z C δ g · ˜ RRR ∂ ˜Ψ(ΓΓΓ , KKK ) ∂ ΓΓΓ − n ! d s + Z C δ k · ˜ RRR ∂ ˜Ψ(ΓΓΓ , KKK ) ∂KKK − m ! d s = 0 , Z C δ ¯ µ ¯ µ ¯ µ · ( ϕϕϕ c − ˜ ϕϕϕ ) | C | d s + 12 Z C PPP T α ( FFF c DDD α − ddd α ) · δ ˜ µµµ τ | C | d s + 12 Z C (cid:16) FFF T c FFF c − III (cid:17) : [ δ ˜ µµµ n ] s | C | d s = 0 , (71) for all δϕϕϕ ∈ { v ∈ H (Ω ) , s.t. DDD α · ∇ v | C ∈ H ( C ) } , δ g , δ k , δ ¯ µ ¯ µ ¯ µ, δ ˜ µµµ τ , δ ˜ µµµ n ∈ (cid:16) L ( C ) (cid:17) ,where the corresponding stresses are defined as follows: PPP := ∂ Ψ ∂FFF , (72) PPP n := n ⊗ DDD , (73) PPP m := 12 h PPP α m − PPP α (cid:16) ˜ ϕϕϕ × n + ¯˜ mmm (cid:17)i ⊗ DDD α , (74) PPP g := 12 PPP α m ⊗ DDD α , (75) P := 12 PPP α m ⊗ DDD α ⊗ DDD . (76) Recall that the tensors
PPP α and [ δ ˜ µµµ n ] s are defined in Remark 2. For the proof of Theorem 1, let us first consider the following remark along with tworelated lemmas.
Remark 5.
Note that the variation of
ΣΣΣ = µµµ α ⊗ DDD α has the form δ ΣΣΣ = δµµµ α ⊗ DDD α . Thus,the constraints (70) can be rewritten as follows: Z C (cid:16) FFF c DDD − ddd (cid:17) · δ ˜ µµµ | C | d s = 0 , Z C (cid:16) FFF c DDD − ddd (cid:17) · δ ˜ µµµ | C | d s = 0 . (77)16 lternatively, in line with the representation (49) , let δ ΣΣΣ = (
PPP α δ ˜ µµµ τ + QQQ α δ ˜ µµµ n ) ⊗ DDD α .Then, the constraints (70) rewrite as Z C (cid:16) PPP T α [ FFF c DDD α − ddd α ] (cid:17) · δ ˜ µµµ τ | C | d s = 0 , (78)12 Z C (cid:16) QQQ T α [ FFF c DDD α − ddd α ] (cid:17) · δ ˜ µµµ n | C | d s = 0 , (79) such that the hydrostatic pressure and shear are separated from the bending and torsionalparts. In fact, the last condition does not involve the beam at all, restraining only thestretches and the shears of the matrix, which is shown in the following lemma. Lemma 3.
Condition (79) can be rewritten as Z C (cid:16) FFF T c FFF c − III (cid:17) : [ δ ˜ µµµ n ] s | C | d s = 0 , (80) with [ δ ˜ µµµ n ] s = δ ˜ µ n, DDD ⊗ DDD + δ ˜ µ n, DDD ⊗ DDD + δ ˜ µ n, ( DDD ⊗ DDD + DDD ⊗ DDD ) .Proof. Integrating in (70) with δ ΣΣΣ = ˜
RRR [ δ ˜ µµµ n ] s and δ ΣΣΣ =
FFF c [ δ ˜ µµµ n ] s , and using symmetryof [ δ ˜ µµµ n ] s , we obtain12 Z C RRR T ( FFF c − ˜ RRR ) : [ δ ˜ µµµ n ] s | C | d s = 0 , Z C ( FFF c − ˜ RRR ) T FFF c : [ δ ˜ µµµ n ] s | C | d s = 0 , (81)respectively. Given FFF T c FFF c − III = ˜
RRR T ( FFF c − ˜ RRR ) + (
FFF c − ˜ RRR ) T FFF c , summing up the integralsin (81) yields the statement. Lemma 4.
The following equality holds: Z C QQQ α ˜ µµµ n ⊗ DDD α : ∇ δϕϕϕ | C | d s = Z C FFF c [˜ µµµ n ] s : ∇ δϕϕϕ | C | d s (82) Proof.
Integrating in (70) with δ ΣΣΣ = ∇ δϕϕϕ [˜ µµµ n ] s yields the statement. Proof of Theorem 1.
Using (65) with test functions δ ˜ ϕϕϕ = δϕϕϕ | C , we obtain (64). Hence,after integration by parts and noting that ∂ s ≡ DDD · ∇ and the endpoint conditions n = nnn eext at s = 0 , L , we obtain Z Ω ∂ Ψ ∂FFF : ∇ δϕϕϕ − BBB ext · δϕϕϕ ! d V − Z Γ σ TTT ext · δϕϕϕ d A + Z C (cid:16) [ n ⊗ DDD ] : ∇ δϕϕϕ − ¯˜ nnn · δϕϕϕ (cid:17) d s − [[ nnn eext · δϕϕϕ ]] | L + Z C
12 ΣΣΣ : ∇ δϕϕϕ | C | d s = 0 . (83)17hus, we have eliminated the set of Lagrange multipliers ¯ µµµ . For the second set ofLagrange multipliers ΣΣΣ, we integrate in (66) with the test functions δφφφ = PPP T α ∇ δϕϕϕ DDD α ∈ L ( C ). Taking the decomposition (49) and Lemma (4) into account, yields − Z C h PPP α (cid:16) m + ˜ ϕϕϕ × n + ¯˜ mmm (cid:17) ⊗ DDD α + (ΣΣΣ − FFF c [˜ µµµ n ] s ) | C | i : ∇ δϕϕϕ d s = 0 . (84)Then, multiplying the above equation by , adding it to (83) and integrating it by partswith the endpoint conditions m = mmm eext at s = 0 , L , we obtain Z Ω ∂ Ψ ∂FFF |{z} PPP : ∇ δϕϕϕ − BBB ext · δϕϕϕ ! d V − Z Γ σ TTT ext · δϕϕϕ d A + Z C n ⊗ DDD | {z } PPP n : ∇ δϕϕϕ − ¯˜ nnn · δϕϕϕ ! d s − (cid:20)(cid:20) nnn eext · δϕϕϕ + 12 ( PPP α mmm eext ⊗ DDD α ) : ∇ δϕϕϕ (cid:21)(cid:21)(cid:12)(cid:12)(cid:12)(cid:12) L + Z C
12 [
PPP α m ⊗ DDD α ⊗ DDD ] | {z } P ... ∇ δϕϕϕ d s + Z C
12 [
PPP α m ⊗ DDD α ] | {z } PPP g : ∇ δϕϕϕ d s + Z C h PPP α m − PPP α (cid:16) ˜ ϕϕϕ × n + ¯˜ mmm (cid:17)i ⊗ DDD α | {z } PPP m + 12 FFF c [˜ µµµ n ] s | C | ! : ∇ δϕϕϕ d s = 0 , (85)supplemented with constitutive equations (67)-(68) and coupling constraints (69), (78)and (80).The existence of solution to the system (71) is an open question and beyond the scopeof this work. For the remainder of this section, a series of remarks with further de-tails on technical issues as well as the interpretation of the arising coupled stresses areprovided. Remark 6.
Note that the term
PPP g in (75) is due to the initial curvature of the beam.In the case of a straight initial configuration, i.e., when DDD α ≡ , this term vanishes.Moreover, according to Appendix A.4, we have in this case PPP α = [ κκκ ] × PPP α − PPP α [ κκκ ] × . Remark 7.
We now obtain in (71) an explicit representation for the third order stresstensor introduced in (24) , which emanates from the condensation of the moment equationof the beams. Note that this term directly depends on the resultant torque m , derivedfrom, e.g., (36) , which is written in terms of the second and the polar areal moment,respectively. This third order stress tensor is restricted to the line C at the scale of thefiber, which is, in general, in terms of [µm] for short fiber reinforced polymers. For amore general strain gradient framework at a macroscopic scale, suitable homogenizationprocedures have to be applied. emark 8. For the physical interpretation of the third order stress tensor P given in (76) , let us consider equation (84) , such that the stress contributions of the condensedbeam equations along the line integral are written as ΣΣΣ = (cid:16)
PPP α h − ( m + ˜ ϕϕϕ × n + ¯˜ mmm ) | C | − i| {z } ˜ µµµ τ + QQQ α ˜ µµµ n (cid:17) ⊗ DDD α . (86) Using the component ˜ µ τ,j and ˜ µ n,j with regard to the basis ddd j in the actual configurationand (30) we find ΣΣΣ ˜
FFF T = ˜ RRR
DDD ˜ µ n, ˜ µ n, − ˜ µ τ, µ n, + ˜ µ τ, ˜ µ n, − µ τ, µ τ, | {z } σσσ b DDD T ˜ RRR T . (87) The decomposition of σσσ b in symmetric and skew symmetric parts reads σσσ bsym = ˜ µ n, ˜ µ n, − ˜ µ τ, ˜ µ n, ˜ µ n, ˜ µ τ, − ˜ µ τ, ˜ µ τ, , σσσ bskw = − ˜ µ τ, ˜ µ τ, ˜ µ τ, − ˜ µ τ, − ˜ µ τ, ˜ µ τ, . (88) Note, that n ⊗ DDD induces additional contributions to both parts if n is not aligned with ddd . Summarized, we obtain skew symmetric stress contributions, which can be consideredas the result of the line divergence of the coupled forces. Concerning the spatial approximation, suitable isogeometric discretization schemes areemployed for the different fields. In particular, a standard IGA approach is applied tothe matrix as well as to the beam.
Matrix material.
A standard displacement-based finite element approach employs ap-proximations of the deformation field ϕϕϕ and its variation of the form ϕϕϕ h = X A ∈I R A qqq A and δϕϕϕ h = X A ∈I R A δqqq A , (89)respectively, where qqq A ∈ R and δqqq A ∈ R . Here, R A : B h0 → R are NURBS basedshape functions associated of order p with control points A ∈ I = { , . . . , M } , where M denote the overall number of control points. The discrete form of (18) reads now δqqq A · Z Ω PPP h ∇ R A d V − Z Ω R A BBB h ext d V − Z Γ σ R A TTT h ext d A = 0 , ∀ δqqq A , (90)19here PPP h , BBB h ext and TTT h ext are the appropriately approximated Piola-Kirchhoff stress andthe external contributions. Beam.
For the beam, we require in a first step the approximations of the centerlineand the rotations, where we have make use of quaternions q for the parametrization, c.f.Appendix B. The approximation along the centerline is given by˜ ϕϕϕ h = X A ∈J ˜ R A ˜ qqq A and q h = X A ∈J ˜ R A q A . (91)Here, ˜ R A : C h0 → R are 1-dimensional, NURBS based shape functions of order p withassociated control points A ∈ J = { , . . . , O } . For the mixed formulation as introducedbefore, additional approximations are required for the resultants n h = X A ∈L ˜ M A n A and m h = X A ∈L ˜ M A m A , (92)where ˜ M A : C h0 → R are 1-dimensional, NURBS based shape functions of order p − A ∈ L = { , . . . , O − } . The test functions δ ˜ ϕϕϕ h and δφφφ h and in the context of the applied mixed method, δ g h and δ k h are discretized via aclassical Bubnov-Galerkin approach. Thus, we obtain δ ˜ ϕϕϕ h = X A ∈J ˜ R A δ ˜ qqq A , δφφφ h = X A ∈J ˜ R A δφφφ A , (93)and δ g h = X A ∈L ˜ M A δ g A , δ k h = X A ∈L ˜ M A δ k A . (94)As the balance equations are condensed within the matrix equation, we will show thediscrete contributions of the mixed beam formulation subsequently within the coupledsystem. Coupling conditions.
Next, we introduce suitable approximations for the variations ofthe Lagrange multipliers as δ ¯ µµµ h = X A ∈J ˜ N A δ ¯ µµµ A , (95)along with the 6 variations δ ˜ µµµ h τ = X A ∈J ˜ N A δ ˜ µµµ τ,A , δ ˜ µµµ h n = X A ∈L ˜ N A δ ˜ µµµ n,A , (96)using the discretization of the beam, along which we define the Lagrange multipliers.Here, ˜ N A : C h0 → R are 1-dimensional, NURBS based shape functions, this time oforder p − A ∈ L = { , . . . , O − } . Note that ¯ µµµ and˜ µµµ τ are already condensed and ˜ µµµ h n = P A ∈L ˜ N A ˜ µµµ n,A follows immediately from (96) for aBubnov-Galerkin approach. 20 emark 9. We have to note that the coupling condition in the form (79) explicitlyincludes the rotation tensor ˜ RRR . Due to its orthogonality, this condition can be writtenonly in terms of
FFF (see Lemma 3). However, the quaternion approach, which is freeof singularities, produces an error in the orthogonality of ˜ RRR h , as the unity constraintsprovide their own approximation error. Thus, (79) couples q h and FFF h providing anadditional numerical error. These issues are avoided when using the condition in formof (80) . Discrete system.
The discretized condensed system as proposed in (71) reads Z Ω PPP h ∇ R A d V − Z Ω R A BBB h ext d V − Z Γ σ R A TTT h ext d A + Z C (cid:18) P h : ∇ R A + (cid:18) PPP n , h + PPP m , h + PPP g, h + 12 FFF h c [˜ µµµ h n ] s (cid:12)(cid:12)(cid:12) C h (cid:12)(cid:12)(cid:12)(cid:19) ∇ R A − R A ¯˜ nnn h (cid:19) d s − (cid:20) R A nnn e, h ext + 12 (cid:16) PPP h α mmm e, h ext ⊗ DDD h α (cid:17) ∇ R A (cid:21) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) L = 000 , (97)evaluated at every node A . Additionally, we require000 = Z C ˜ M A " ˜ RRR h ∂ ˜Ψ ∂ ΓΓΓ (cid:16)
ΓΓΓ h , KKK h (cid:17) − n h d s,
000 = Z C ˜ M A " ˜ RRR h ∂ ˜Ψ ∂KKK (cid:16) ΓΓΓ h , KKK h (cid:17) − m h d s, Z C ˜ R A h q h · q h − i d s (98)along with 000 = Z C ˜ N A h ϕϕϕ h c − ˜ ϕϕϕ h i (cid:12)(cid:12)(cid:12) C h (cid:12)(cid:12)(cid:12) d s,
000 = 12 Z C ˜ N A PPP T , h α (cid:20) FFF h c DDD h α − ˜ RRR h DDD h α (cid:21) (cid:12)(cid:12)(cid:12) C h (cid:12)(cid:12)(cid:12) d s,
000 = 12 Z C ˜ N A QQQ T , h α ˜ RRR h h FFF T , h c FFF h c − III i DDD h α (cid:12)(cid:12)(cid:12) C h (cid:12)(cid:12)(cid:12) d s, (99)Note that the discrete values of P h , PPP n , h and PPP m , h , are given as the discrete counterpartsof (85). Moreover, the unity constraint of the quaternions is evaluated in an integralsense, i.e., this is not necessarily fulfilled at each Gauss point where we use them toconstruct the rotation matrix. 21 iscrete reduction scheme. Within the Newton-Raphson iteration, we have to solveat every step
KKK ϕϕϕϕϕϕ
KKK ϕϕϕ ˜ ϕϕϕ KKK ϕϕϕ n KKK ϕϕϕ m KKK ϕϕϕ q KKK ϕϕϕ ˜ µµµ n KKK ˜ ϕϕϕϕϕϕ MMM ˜ ϕϕϕ ˜ ϕϕϕ
000 000 000 000000
KKK n ˜ ϕϕϕ MMM nn KKK nq MMM mm KKK mq KKK q ϕϕϕ
000 000 000
KKK qq KKK ˜ µµµ n ϕϕϕ
000 000 000
KKK ˜ µµµ n q ∆ ϕϕϕ ∆˜ ϕϕϕ ∆ n ∆ m ∆ q ∆˜ µµµ n = RRR ϕϕϕ
RRR ˜ ϕϕϕ RRR n RRR m RRR q RRR ˜ µµµ n , (100)where we have made use of (36). Here, RRR ϕϕϕ correlates to (97),
RRR ˜ ϕϕϕ to (99) , RRR n to (98) , RRR m to (98) , RRR q and RRR p to (98) , (99) and (99) , respectively. The matrices on thediagonal with components[ MMM ˜ ϕϕϕ ˜ ϕϕϕ ] AB = − Z C ˜ N A ˜ R B d s, and [ MMM nn ] AB = [ MMM mm ] AB = − Z C ˜ M A ˜ M B d s (101)are invertible, as long as MMM ˜ ϕϕϕ ˜ ϕϕϕ is quadratic and invertible, depending on the order of δ ¯ µµµ .This is obvious for ˜ N A = ˜ R A , for all other cases the invertibility has to be shown forthe specific case. Different choices are compared in the subsequently following numericalexamples. With this, we can solve the second, third and fourth line in (100) via∆˜ ϕϕϕ = MMM − ϕϕϕ ˜ ϕϕϕ ( RRR ˜ ϕϕϕ − KKK ˜ ϕϕϕϕϕϕ ∆ ϕϕϕ ) , ∆ n = MMM − nn ( RRR n − KKK n ˜ ϕϕϕ ∆˜ ϕϕϕ − KKK nq ∆ q ) , ∆ m = MMM − mm ( RRR m − KKK mq ∆ q ) , (102)and obtain ¯ KKK ϕϕϕϕϕϕ ¯ KKK ϕϕϕ q KKK ϕϕϕ ˜ µµµ n KKK q ϕϕϕ KKK qq KKK ˜ µµµ n ϕϕϕ KKK ˜ µµµ n q ∆ ϕϕϕ ∆ q ∆˜ µµµ n = ¯ RRR ϕϕϕ
RRR q RRR ˜ µµµ n , (103)where ¯ KKK ϕϕϕϕϕϕ = KKK ϕϕϕϕϕϕ − ( KKK ϕϕϕ ˜ ϕϕϕ − KKK ϕϕϕ n MMM − nn KKK n ˜ ϕϕϕ ) MMM − ϕϕϕ ˜ ϕϕϕ KKK ˜ ϕϕϕϕϕϕ , ¯ KKK ϕϕϕ q = KKK ϕϕϕ q − KKK ϕϕϕ n MMM − nn KKK nq − KKK ϕϕϕ m MMM − mm KKK mq , (104)and the modified residual vector¯ RRR ϕϕϕ = RRR ϕϕϕ − ( KKK ϕϕϕ ˜ ϕϕϕ − KKK ϕϕϕ n MMM − nn KKK n ˜ ϕϕϕ ) MMM − ϕϕϕ ˜ ϕϕϕ RRR ˜ ϕϕϕ − KKK ϕϕϕ n MMM − nn RRR n − KKK ϕϕϕ m MMM − mm RRR m . (105) Remark 10.
Without further proof we remark here, that the quaternions can be removedfrom the system using (99) with regard to (99) . Therefore, we evaluate ddd α = FFF c DDD α xz M y xz M Figure 3.
Reference configuration of the bending test (left) and the torsion test (right). and calculate ddd = ddd × ddd at the respective Gauss point along the beam center line, suchthat the rotation tensor is given by ˜ RRR = ddd i ⊗ DDD i . Eventually, we need ˜ RRR = ddd i ⊗ DDD i assuming again straight initial fibers, and we obtain immediately ddd α = ∇ FFF : DDD α ⊗ DDD , and ddd = ( ∇ FFF : DDD ⊗ DDD ) × ddd + ddd × ( ∇ FFF : DDD ⊗ DDD ) . (106) This leads finally to a system to be solved with respect to the matrix degrees of freedom ϕϕϕ and the Lagrange multipliers ˜ µµµ n . If required, suitable methods like augmented Lagrangecan be considered to solve solely with respect to the matrix unknowns. For the numericalexamples in the next section, quaternions are used. In this section, we investigate the accuracy of the proposed formulation. In particular, weconsider a benchmark test with results for a surface coupling between beam and matrixinstead of a multidimensional coupling from [44]. Afterwards we investigate a torsionaltest, where it is obvious that pure position constraints are not suitable. Eventually, wedemonstrate the applicability towards larger representative volume elements (RVE) forthe analysis of multiple embedded fibers.
In this first numerical example, we consider a model problem from [44], Section 4.2: Abeam of length 5 m and radius r = 0 .
125 m with Young modulus 4346 N / m is embeddedinto the 1 m × × / m . Poisson ratio is zero for both materials. The geometry and materialparameters for both, the matrix and the beam are chosen such that both systems havesimilar properties, i.e., if the matrix material would be considered as beam, we wouldobtain the same moment of inertia.The matrix and the beam are both fixed at x = 0 m, and we apply a moment mmm Lext =[0 , , . x = 5 m as a dead load, see Figure 3, left, for details.23 igure 4. Bending of a beam: Von Mises stress distribution.
This benchmark test allows us to illustrate the stability of the proposed approach, itsconvergence and the model error, since the solution for a surface coupling is providedin [44]. The matrix is discretized with B-Splines of order p = [ p x , p y , p z ]. We considerthe same order in y - and z -direction, p y = p z , and the order in x -direction is the sameas for the beam. The number n of elements in y - and z -direction is the same and by afactor of five smaller than in x -direction. For the beam, the same number of elementsin x -direction is used.The deformation and the von Mises stress distribution for the bending problem is pre-sented in Figure 4. As expected, the stresses concentrate at the right hand side of thebeam where the external moment is applied to. In Figure 5, the displacement of the tipof the beam is presented for a second and fourth order example. The numbers in thelegend denote the order of the B-splines for [ ϕϕϕ, ˜ ϕϕϕ, q , n , m , ¯ µµµ, ˜ µµµ τ , ˜ µµµ n ]. As can be seen,both combinations of different orders converge to the same result.In Figure 6, the results for n , m and ˜ µµµ n along the beam center line are displayed. Thedifferent order chosen for [ ϕϕϕ, ˜ ϕϕϕ, q , n , m , ¯ µµµ, ˜ µµµ τ , ˜ µµµ n ] are displayed as well. Note thatwe fixed ˜ µµµ n at the first node equal to zero, to avoid spurious oscillations leading toa divergence within the Newton iteration. Moreover, the components [ n ] , [ m ] , [ m ] and [˜ µµµ n ] are close to zero, as expected for this example. Finally, we have plotted theresults for an explicit enforcement of (70) at both endpoints of the beam additionallyto the integral enforcement of the constraints along C . Note that no differences can beobserved. 24 igure 5. Bending of a beam: Convergence of the tip displacement for different combinations ofshape functions.
Routine expl. EP expl. EP expl. EP expl. EP cond condq-unity weak direct direct direct direct directorders 44433222 44433222 44444333 44444222 44433442 44433443expl. EP weak 44433222 1.13E-4 % 1.32E-4 % 1.13E-4% 1.17E-4 % 1.15E-4 %expl. EP direct 44433222 1.13E-4 % 1.90E-5 % 1.00E-7 % 4.00E-6% 2.00E-6 %expl. EP direct 44444333 1.32E-4 % 1.90E-5 % 1.90E-5 % 1.50E-5 % 1.70E-5 %expl. EP direct 44444222 1.13E-4 % 1.00E-7 % 1.90E-5 % 4.00E-6% 2.00E-6 %cond direct 44433442 1.17E-4 % 4.00E-6 % 1.50E-5 % 4.00E-6 % 2.00E-6 %cond direct 44433443 1.15E-4 % 2.00E-6 % 1.70E-5 % 2.00E-6 % 2.00E-6 %
Table 1: Bending of a beam:
Deviation of the beam tip displacement in percentage for differentorders and approaches for the quaternion unity constraint. igure 6. Bending of a beam: From top to bottom, n , m and ˜ µµµ n , from left to right 1-, 2- and3-components are displayed. Eventually, we compare the deviation of the beam tip displacement in percentage for dif-ferent approaches. To be more specific, we evaluate the absolute value of | k (˜ ϕϕϕ ( L ) − ˜ ϕϕϕ ( L ) k a k (˜ ϕϕϕ ( L ) − ˜ ϕϕϕ ( L ) k b − | in percentage for approach a compared to b , where the absolute tip displacement ofthe first approach displaced in Table 1 is 0 . ϕϕϕ, ˜ ϕϕϕ, q , n , m , ¯ µµµ, ˜ µµµ τ , ˜ µµµ n ], see Table 1. We remark here, thatthe relative deviation between [4 , , , , , , ,
2] compared to [4 , , , , , , , − %, i.e. negligible and already includes the better approxima-tion of the higher order shape functions. However, the implementation simplifies for[4 , , , , , , ,
3] significant, as only two different order of shape functions have tobe evaluated at each Gauss point. 26 igure 7. Torsion Test:
Von Mises stress distribution.
Next, we use the same setup as in Section 5.1 using a Mooney-Rivlin material modelΨ(
J, I , I ) = c ( J − − d ln( J ) + c ( I −
3) + c ( I − , (107)as we expect to obtain large deformations. Here, J = det( FFF ), I = tr( FFF T FFF ) and I = cof( FFF T FFF ). Moreover, c = 2 / c + c ), d = 2 ( c + 2 c ), c = 2 and c = 1. Weapply an external moment mmm Lext = [0 . , ,
0] Nm, such that we obtain a pure torsion(see Figure 3, right). Note that using only position constraints is not possible for thisexample, as the immersed beam can rotate without influencing the matrix material. InFigure 7, the von Mises stress distribution is shown. For visualization, we have plottedthe (virtual) surface of the beam inside the matrix material. As can be observed, thematrix rotates within the virtual area of the beam geometry due to the incorporationof geometrical data in terms of the second moment of inertia within the beam materialmodel in (36). In fact, the gradient around the 1-dimensional beam is restricted; a fullreconstruction of the beam geometry would require all higher-order terms.In Figure 8, the results for n , m and ˜ µµµ n along the beam center line are displayed.Note that the n , n , m and m are again near the numerical limit, which corre-sponds to the physics of this example. The chosen order for [ ϕϕϕ, ˜ ϕϕϕ, q , n , m , ¯ µµµ, ˜ µµµ τ , ˜ µµµ n ] is[4 , , , , , , , igure 8. Torsion test: From top to bottom, n , m and p , from left to right 1-, 2- and 3-componentsare displayed. ration equal zero, we obtain a slight elongation in x − direction, resulting in the displayedforces obtained in n ( L ). A comparative study using the Saint-Venant Kirchhoff modelof the sub-section before results in an exact zero elongation and zero forces in n ( L ). Asthe results are not physically meaningful due to the obtained large deformation, onlythe results for the Mooney-Rivlin material are plotted. The last example using the above described geometry can be denoted as shear test.Therefore, we fix the lower surface of the matrix material and prescribe the displacementof the top surface. In Figure 9 the resulting von Mises stress distribution with (left) andwithout (right) constraints on hydrostatic pressure and shear are presented., see Remark5 and Lemma 3.The shear test does not bend the beam, nor does the beam elongate. Additionally,28 igure 9. Shear Test:
Von Mises stress distribution with (left) and without (right) constraints onhydrostatic pressure and shear.
Figure 10. Application of multiple beams:
Reference configuration of the RVE, isotropic (left)and anisotropic (right) case with 21 beams. the rotation in the center line is nearly constant, such that without the constraints thematrix is not affected by the beam. This can be observed in Figure 9, right. In contrast,using the constraints in (80), restricting the stretches in terms of the right Cauchy-Greentensor on the cross section of the beam, leads to a redistribution of the stresses aroundthe virtual beam geometry, see left plot in Figure 9.
In this final example, we investigate a representative volume element (RVE) for fiberreinforced ultra high performance concrete (UHPFRC). The concrete has a Young mod-ulus of E = 50 .
000 N / mm , a Poisson ratio of ν = 0 . d = 0 . l = 14 mm, a Young modulus of E = 200 .
000 N / mm and a Poisson ratio of ν = 0 . ×
20 mm ×
20 mm and we made use of 20 × × ϕϕϕ, ˜ ϕϕϕ, q , n , m , ¯ µµµ, ˜ µµµ τ , ˜ µµµ n ] is set to [4 , , , , , , , FFF is applied as Dirichlet boundary on the six surfaces, resulting ina volumetric average first Piola-Kirchhoff ¯
PPP stress tensor via¯
PPP = 1 V Z Ω PPP d V. (108)With this, we obtain ¯ FFF = . . − . − . . − . − . . . , ¯ PPP iso = − . . − . . − . . − . . − . , ¯ PPP ani = − . . − . . − . . − . . − . . (109) Figure 11 displays the corresponding stress distribution, where we obtain peak values ofabout 50 N / mm for the isotropic and 18 N / mm for the anisotropic case compared toan average value of ≈ / mm . The framework provided in this contribution utilizes an overlapping domain decompo-sition technique similar to immersed techniques in fluid and solid mechanics to embed1-dimensional fibers into a suitable 3-dimensional matrix material. For the first time,the coupling terms not only consider the forces, but also the bending and torsional mo-ments of the beam model, leading to an enhanced formulation, which has shown to besuperior to previous approaches. The chosen benchmark tests for bending and torsiondemonstrate the flexibility as well as the accuracy of the proposed model. Especially thelatter example regarding a torsional moment cannot be obtained from a simple couplingof forces. 30 igure 11. Application of multiple beams:
Von Mises stress distribution, isotropic (left) andanisotropic (right) case, plotted in the reference configuration.
The proposed static condensation procedures in the continuum as well as in the dis-crete setting reduces the computational effort on the calculation of the matrix materialwith regard to the enforcement of the shear and pressure terms of the beam, whichkinematically constrains the cross sectional area. This condensation procedures for themultidimensional coupling of the immersed fibers give rise to a third order stress tensor,which is related to the concept of coupled stresses. This can be considered as first stepof a homogenization, as we can now avoid the highly inefficient resolution of the beamas a classical Cauchy continuum.
Acknowledgements
Support for this research was provided by the Deutsche Forschungsgemeinschaft (DFG)under grant DI2306/1-1. The author C. Hesch gratefully acknowledge support by theDFG. U. Khristenko and B. Wohlmuth gratefully acknowledge support by the EuropeanUnion’s Horizon 2020 research and innovation programme under grant agreement No800898, the German Research Foundation by grants WO671/11-1 and WO671/15-2.We would like to thank the group of Alexander Popp at the Universität der BundeswehrMünchen and his Co-Worker Ivo Steinbrecher for providing us additional data for thebenchmark test. We also thank the group of Torsten Leutbecher at the University ofSiegen for providing all necessary information on fiber reinforced ultra high performanceconcrete in the final example. We also thank Patrick Le Tallec (École PolytechniqueParis) for helpful discussions. 31
Proofs
A.1 Surface terms
The internal virtual work reads δ Π int ( ϕϕϕ ) = Z Ω (cid:18) PPP : ∇ δϕϕϕ + P ... ∇ δϕϕϕ (cid:19) d V. (110)The first integration by parts gives then δ Π int ( ϕϕϕ ) = Z Ω (cid:16) −∇ · PPP · δϕϕϕ + ∇ · ( PPP T δϕϕϕ ) − ∇ · P : ∇ δϕϕϕ + ∇ · ( P T : ∇ δϕϕϕ ) (cid:17) d V, (111)and in the second step δ Π int ( ϕϕϕ ) = Z Ω (cid:16) ∇ · ( ∇ · P − PPP ) · δϕϕϕ + ∇ · [( PPP T − ∇ · P T ) δϕϕϕ ] + ∇ · ( P T : ∇ δϕϕϕ ) (cid:17) d V. (112)Application of the divergence theorem on the second and third term yields δ Π int ( ϕϕϕ ) = Z Ω ∇ · ( ∇ · P − PPP ) · δϕϕϕ d V + Z ∂ Ω (cid:16) δϕϕϕ · [( PPP − ∇ · P ) NNN ] + ∇ δϕϕϕ : ( P NNN ) (cid:17) d A. (113)Next, we decompose the last term as follows ∇ δϕϕϕ : ( P NNN ) = ∇ · ( δϕϕϕ P NNN ) − δϕϕϕ · ∇ · ( P NNN )= ∇ ⊥ · ( δϕϕϕ P NNN ) + ∇ k · ( δϕϕϕ P NNN ) − δϕϕϕ · ∇ ⊥ · ( P NNN ) − δϕϕϕ · ∇ k · ( P NNN ) . (114)Using the identity ∇ ⊥ · ( δϕϕϕ P NNN ) = δϕϕϕ · ∇ ⊥ · ( P NNN ) + ∇ ⊥ δϕϕϕ : ( P NNN ) and the divergencetheorem: Z ∂ Ω ∇ k · ( δϕϕϕ P NNN ) d A = Z ∂ Ω ( δϕϕϕ P NNN ) · ˆ NNN d S − Z ∂ Ω K ( δϕϕϕ P NNN ) · NNN d A (115)we obtain Z ∂ Ω ∇ δϕϕϕ : ( P NNN ) d A = − Z ∂ Ω δϕϕϕ · ( K ( P NNN ) NNN + ∇ k · ( P NNN )) − ∇ ⊥ δϕϕϕ : ( P NNN ) d A + Z ∂ Ω δϕϕϕ · ( P : [ ˆ NNN ⊗ NNN ]) d S. (116)32 .2 Deformation gradient of the beam The deformation gradient of the beam (29) is given here once again:˜
FFF = ˜ ϕϕϕ ⊗ DDD + ddd α ⊗ DDD α + θ α ddd α ⊗ DDD . (117)The next step is to multiply the deformation gradient with the identity matrix and tosubstitute the product ddd α ⊗ DDD α :˜ RRR ˜ RRR T = ˜ RRR T ˜ RRR = III, ddd α ⊗ DDD α = ˜ RRR − ddd ⊗ DDD . (118)Furthermore, we use ˜ RRR T ddd α = ˜ RRR T ˜ RRR DDD α and DDD = ˜ RRR T ddd such that (117) ends up in:˜ FFF = ˜
RRR (cid:18) ˜ RRR T ˜ ϕϕϕ ⊗ DDD + ˜ RRR T (cid:16) ˜ RRR − ddd ⊗ DDD (cid:17) + θ α ˜ RRR T ˜ RRR DDD α ⊗ DDD (cid:19) , = ˜ RRR (cid:18)(cid:20) ˜ RRR T ˜ ϕϕϕ − DDD (cid:21)| {z } ΓΓΓ ⊗ DDD + θ α ˜ RRR T ˜ RRR | {z } [ KKK ] × DDD α ⊗ DDD + III (cid:19) . (119) A.3 Virtual work of the beam
To consider again the principle of virtual work, we introduce the internal virtual workof the beam δ ˜Π int = Z C ˜ NNN · δ ΓΓΓ + ˜
MMM · δKKK d s (120)with the variations of the beam strain measures: δ ΓΓΓ = δ (cid:18) ˜ RRR T ˜ ϕϕϕ − DDD (cid:19) = δ ˜ RRR T ˜ ϕϕϕ + ˜ RRR T δ ˜ ϕϕϕ = − ˜ RRR T [ δφφφ ] × ˜ ϕϕϕ + ˜ RRR T δ ˜ ϕϕϕ (121)and [ δKKK ] × = δ (cid:18) ˜ RRR T ˜ RRR (cid:19) = δ ˜ RRR T ˜ RRR + ˜ RRR T δ ˜ RRR = − ˜ RRR T [ δφφφ ] × ˜ RRR + ˜ RRR T (cid:16) [ δφφφ ] × ˜ RRR (cid:17) = ˜ RRR T [ δφφφ ] × ˜ RRR = − (cid:15) ijk δφ k DDD i ⊗ DDD j , (122)where we make use of the skew-symmetric matrix [ δφφφ ] × , which is given by:[ δφφφ ] × = − (cid:15) ijk δφ k ddd i ⊗ ddd j . (123)33he axial vector of [ δKKK ] × is then given by:[ δKKK ] = − (cid:15) pij ( − (cid:15) ijk δφ k ) DDD p = δ pk δφ k DDD p = ˜ RRR T δφφφ (124)Inserting (121) and (124) into the internal virtual work (120), we get: δ ˜Π int = Z C (cid:16) ˜ RRR ˜ NNN (cid:17) · δ ˜ ϕϕϕ − (cid:16) ˜ ϕϕϕ × ˜ RRR ˜ NNN (cid:17) · δφφφ + (cid:16) ˜ RRR ˜ MMM (cid:17) · δφφφ d s, (125)where we used once again the relation h [ δφφφ ] × i jk = − (cid:15) jkl δφφφ l . Here, C = [0 , L ] ⊂ R and δφφφ ∈ R are the virtual rotations. A.4 Derivative of the terms in (50)
With regard to (50) , we write PPP = 12 ( ddd ⊗ ddd + ddd ⊗ ddd ) − ( ddd ⊗ ddd + ddd ⊗ ddd ) . (126)Noting that ddd i = ˜ RRR DDD i = ˜ RRR ˜ RRR T ddd i , where ˜ RRR ˜ RRR T is skew symmetric, such that we candefine κκκ = axl( ˜ RRR ˜ RRR T ). Thus, we obtain PPP = 12 (cid:16) [ κκκ ] × ddd ⊗ ddd + ddd ⊗ [ κκκ ] × ddd (cid:17) − (cid:16) [ κκκ ] × ddd ⊗ ddd + ddd ⊗ [ κκκ ] × ddd (cid:17) , = 12 (cid:16) [ κκκ ] × ddd ⊗ ddd − ddd ⊗ ddd [ κκκ ] × (cid:17) − (cid:16) [ κκκ ] × ddd ⊗ ddd − ddd ⊗ ddd [ κκκ ] × (cid:17) , = [ κκκ ] × PPP − PPP [ κκκ ] × . (127)All further terms presented in (50) follow analoguesly. B Parametrization of the rotation tensor in terms ofquaternions
Quaternions q = ( q , qqq ) can be written in terms of the scalar part q ∈ R and the vectorpart qqq ∈ R . Restricting the quaternions to unit length, i.e., | q | = 1 such that weconsider a unit sphere in R , we form a group as S = { q ∈ R | | q | = 1 } , (128)see [4] for details. This allows us to introduce a mapping ˜ RRR : S → SO (3)˜ RRR ( q ) = ( q − qqq · qqq ) III + 2 qqq ⊗ qqq + 2 q ˆ qqq, (129)called Euler-Rodrigues parametrization, see Marsden and Ratiu [25].34 Rigid bodies
If a Cosserat point instead of a Cosserat beam is taken into account, i.e. a 0D-3Dcoupling, we have to constraint the whole sphere instead of the cross section. Withregard to (47), we decompose ˜ µµµ into torque and pressure˜ µµµ ( θ i ) = µµµ τ × h ˜ RRR NNN ( θ i ) i + µ n ˜ RRR NNN ( θ i ) . (130)Next we insert (130) in (55) and note that (cid:16) µµµ τ × h ˜ RRR NNN i(cid:17) · (cid:16)h FFF c − ˜ RRR i θ i DDD i (cid:17) = (cid:18) ˜ RRR T µµµ τ (cid:19) · (cid:20) NNN × (cid:18) ˜ RRR T h FFF c − ˜ RRR i NNN (cid:19)(cid:21) k θ i DDD i k , = (cid:18) ˜ RRR T µµµ τ (cid:19) · (cid:20) NNN × (cid:18) ˜ RRR T FFF c NNN (cid:19)(cid:21) k θ i DDD i k . (131)Making use of a spherical parametrization NNN = DDD cos θ sin φ + DDD sin θ sin φ + DDD cos φ, (132)we obtain for an arbitrary matrix AAA π Z π Z − π NNN × (cid:16) AAA NNN (cid:17) d θ d φ = π Z (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) DDD DDD DDD cos θ sin φ sin θ sin φ cos φDDD T1 AAA NNN DDD T2 AAA NNN DDD T3 AAA NNN (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) d θ = [ DDD , DDD , DDD ] · π Z DDD T3 AAA NNN sin θ sin φ − DDD T2 AAA NNN cos φ − DDD T3 AAA NNN cos θ sin φ + DDD T1 AAA NNN cos φDDD T2 AAA NNN cos θ sin φ − DDD T1 AAA NNN sin θ sin φ d θ = π [ DDD , DDD , DDD ] · DDD T3 AAA DDD − DDD T2 AAA DDD − DDD T3 AAA DDD + DDD T1 AAA DDD DDD T2 AAA DDD − DDD T1 AAA DDD = 2 π [ DDD , DDD , DDD ] · axl (cid:16) [ DDD , DDD , DDD ] T AAA [ DDD , DDD , DDD ] (cid:17) , (133)and thus π Z π Z − π NNN × (cid:16) ˜ RRR T FFF c NNN (cid:17) d θ d φ = 2 π [ DDD , DDD , DDD ] · axl (cid:18) [ DDD , DDD , DDD ] T ˜ RRR T FFF c [ DDD , DDD , DDD ] (cid:19) . (134)Finally, we can stateskew (cid:16) [ DDD , DDD , DDD ] T ˜ RRR T FFF c [ DDD , DDD , DDD ] (cid:17) = [ DDD , DDD , DDD ] T skew( ˜ RRR T FFF c ) [ DDD , DDD , DDD ] , (135)i.e., we have to enforce that skew( ˜ RRR T FFF c ) is zero. Note that this is equivalent to theformulation considered in [40]. 35 eferenceseferences