Extension of a shell‐like RVE model to represent transverse shear effects in textile reinforced composites

Textile reinforcements have long been used in shell‐like components like tires, belts, hoses and, in particular, air spring bellows. These bellows consist of a certain number of cords embedded in soft rubber, resulting in high membrane stiffness to absorb tensile forces and low resistance to changes in curvature. The embedded cords are produced by twisting multifilament yarns. Modelling and simulation of bellows is a challenging task. Thereby, the complex internal geometry, large deformations as well as the strongly anisotropic and nonlinear material behavior have to be taken into account in a computationally efficient manner. One approach is the application of multiscale analyses in conjunction with representative volume elements (RVE). Within this study, a through‐thickness RVE based on in‐plane periodic microstructures is applied. The local behavior of the cords is represented by an anisotropic constitutive model that regards directions and mechanical properties of filaments, which result from previous cord twisting simulations. In this contribution, an extension of the classical periodic boundary conditions (PBC) is presented, which considers curvature, pressure load, and, in particular, transverse shear. The RVE is used to simulate the stresses and strains in cord‐rubber composites under realistic load conditions. Moreover, it is employed for homogenizing the mechanical behavior and, thus, for developing a suitable constitutive model for shells.


INTRODUCTION
Air spring bellows feature a highly complex internal geometry and are characterized by strong anisotropic, nonlinear material behavior, generally large deformations, and contact.The main objective is to provide novel high-resolution models that account for these characteristics for the structural mechanical simulation of air springs.Consequently, the models have to be both efficient and realistic.Figure 1 illustrates a cross-ply reinforced bellows.The rubber matrix contains two layers of textile cords, which are oriented at the angle  relative to the axial direction and are significantly stiffer than the cord matrix.However, the cord owns a complex geometry and filament structure.The filament structure plays an important role in determining the stiffness properties of the cord and, consequently, the resulting stiffness properties of the composite.This composite structure allows for the absorption of significant membrane forces, which are crucial for F I G U R E 1 Representation of a rolling lobe air spring (left) as well as microscopic images of bellows cross-sections (right) showing the cross-section and longitudinal view of the cords [2].
the load-carrying capacity of the air spring.Due to the dominant influence of the rubber matrix, the composite remains highly flexible against bending.Over the past years, many finite element models have been developed for cord-rubber composites.Frequently, these models use simplifications to deal with the complex structure of the cord, its anisotropic behavior, significant large deformations, or the challenging discretization of the composite [1].For instance, shell elements or rebar models are often employed to analyze the overall behavior of the component.However, these models are not capable of accurately determining local results such as stresses or strains.For such an evaluation, microstructural modelling is required, considering a representative volume element (RVE).These models have the advantage of providing a larger number of elements to resolve the large gradients near the interface.The challenge with these models lies in formulating the macroscopic boundary conditions in such a way that an efficient analysis of real loads on shell-like structures is enabled, while minimizing disruptions caused by boundary effects.Typical loadings on those structures are curvatures, membrane deformations, transversal shear, and pressure loads.The shell-like RVE can capture important features such as microscopic geometric structures and material properties that influence the overall behavior of the system.By analyzing such an RVE, valuable insights into the behavior of the larger system can be gained without the need to model the entire structure or material.Thus, the RVE can be replaced by a shell element of the same physical properties.This procedure is called homogenization.With the aid of a submodel, which is surrounded by an elastomer layer, the results in the interface region can be further enhanced.The boundary conditions of the submodel are derived from the macroscopic simulation of the RVE.

HOMOGENIZATION
Homogenization is efficiently applied in multiscale modelling when the complexity of the geometric and mechanical microstructure is too challenging to be handled by other conventional modelling techniques.In order to describe the behavior of composite materials, there are three methods of homogenization.The first method is analytical homogenization, followed by numerical homogenization techniques, such as fast Fourier transformation (FFT)-based homogenization for obtaining effective properties in linear elasticity and conductivity problems.Thus, this method is not suitable for arbitrary microscale material behavior, including physically nonlinear and time dependent behavior [3].The last method is computational homogenization, which represents an effective method for incorporating multiscale hierarchies into computational models.The basic idea behind computational homogenization is to couple two finite element models (FE 2 ) by solving a microscale boundary value problem at the integration points of the macroscale model.However, the main challenge of this approach is the high computational cost, which requires additional efforts to explore strategies for reducing the computational costs.The theory of shell computation homogenization is proposed by [4] and [5].The aim of this paper is to develop a shell-like RVE for highly loaded thin cord-elastomer composites that is capable of simulating not only membrane deformations, large curvatures, and internal pressure but also transverse shear.Two homogenization techniques are implemented and compared in order to identify the optimal approach for defining the periodic boundary conditions (PBC).The first technique, proposed by [4], represents a second-order computational homogenization framework of thin shells, as an extension of the classical first-order computational homogenization framework introduced by [6].This technique is presented in detail in Section 2.2 of the paper.The second technique, proposed by [1], extends the classical first-order PBC to allow for rigid body rotation of the material line elements at the periodic surfaces of the RVE.
In the context of this work, the second technique is further extended to incorporate transverse shear deformation.This technique is explained in Section 2.3.

Macroscale mindlin-reissner kinematics
The Mindlin-Reissner theory, also known as the first-order shear deformation theory, employs a linear approximation for the displacements of the cross-section, as illustrated in Figure 2. In contrast to the Kirchhoff-Love plate theory, in which a cross-section perpendicular to the plate's midplane remains straight and perpendicular after deformation, the Mindlin-Reissner concept allows for a rotation of the cross-section while maintaining its straightness.The number of bars under the tensor symbol represents the tensor's order.For instance, vectors are represented with a single bar, while secondorder tensors are denoted with two bars.Figure 2 shows the trajectory of an arbitrary material point.The position vectors r and  represent the point's location in the reference and current configurations, respectively.Although both vectors are expressed in the global coordinate system   , they can also be described using local curvilinear coordinates ( 1 ,  2 ,  3 ).
The coordinates   correspond to the mid-surface of the shell continuum with  = 1, 2, with the associated tensor being depicted in bold.Meanwhile, the coordinate  3 represents the direction transverse to the mid-surface.In the reference configuration, the position vector of any material point can be linearly approximated as follows [7].
The operators ( * ) , and ( * ) ,3 are defined as the partial derivatives of ( * ) with respect to   and  3 , respectively.In the current configuration, the director  is rotated with respect to the original director d due to two factors: transverse shear   and rotation of the reference plane by the angle   .The transverse shear deformation causes a rotation of the director along the thickness of the material, while the rotation of the reference plane is associated with the curvature of the structure.The transport of a material line element from the undeformed to the deformed configuration is described by the following operator: The macroscale deformation gradient tensor  of the three-dimensional continuum is approximated using the first-order Mindlin-Reissner theory with  =   ⊗ g as the in-plane deformation gradient tensor of the reference surface, which corresponds to the membrane deformation tensor in shell theories.Moreover,  = d, ⊗ g is the curvature tensor, and  represents the in-plane transverse shear.Additionally, in accordance with classical shell theories, an infinitesimal rotation description is employed for the in-plane deformation of the director, effectively restricting the main of curvature to a small value of  = 1 2  ⋅  −1 ≪ 1.The shell strains are obtained based on the assumed kinematics stated in assumption (2) and can be defined as follows [7] where  is a Hencky-tensor of the membrane deformation,  is denotes a symmetric relative curvature tensor,  represents the transverse shear strains, and ln((⋯)) denotes the tensorial natural logarithmic function.The membrane forces are quantified as   in [N mm] units, the bending moments are symbolized as   in [N mm 2 ] units, and the shear forces are indicated as   in [N mm] units, with their expressions as follows Here, ℎ represents the thickness of the shell.To define the nonlinear mechanical behavior of a shell section in terms of generalized section quantities using the ABAQUS subroutine UGENS, the force array [  ] and strain array [  ] are used.In these arrays, the direct membrane terms are stored first, followed by the shear membrane term, the direct and shear bending terms, and then the transverse shear components, as defined in Equation ( 5).The number of entries in these arrays depends on the specific element type, as only the active components are stored.
[ (5) The linearized constitutive relations that align with the macroscopic framework can be expressed as follows where the macroscopic constitutive tangents of a shell element are captured by an 8 × 8 stiffness matrix [  ].This matrix can be decomposed into four distinct parts.The first part [] is the so-called membrane stiffness, the second part [] represents the plate stiffness and the mixed derivatives [] correspond to coupling stiffnesses, accounting for the coupling between in-plane and out-of-plane deformations.The last part [] represent the transverse shear stiffness.

Microscale projection: Second-order homogenization
The transition from macro-to microscales involves averaging relationships that establish the boundary conditions for a microstructural RVE.These boundary conditions impose the macroscopic quantities, specifically the generalized strains , , and .This transition departs from the formulation (expressed in Equation ( 2)) for the material line element within a shell at the in-plane level.In order to consider the influence of a local fine-scale microstructure and its associated displacement field, a microfluctuation field   is introduced and superimposed onto the macroscopic displacement field.However, the transition is only applied to the surfaces of the microstructural RVE that have normal vectors ñ tangential to the midplane.Figure 3 provides a representation of these surfaces and their orientation relative to the midplane.On the other hand, the surfaces of the RVE that have a normal vector in the thickness direction remain unconstrained.Here, the displacement vector  and the average three-dimensional macroscale displacement gradient  are relevant.Consequently, the PBC couple the displacements of specific pairs of points.The Taylor expansion of the macroscopic field is The displacements of the cross-section are linearly approximated, which leads to deviations when the crosssection undergoes large rotations.

Microscale projection: Rotation-based homogenization
In the case of an air spring bellows, the curvatures are dominant and highly significant.Therefore, it is of great importance to develop an RVE-model with rotation-based PBC that is capable of simulating the significant rotations of the cross-section.Thereby, the modification of the PBC should be done without violating the assumptions of the Mindlin-Reissner shell theory, in order to maintain information exchange between the macro and micro FE-models.Henceforth, side surfaces of the rotation based RVE are divided into so-called "control surfaces" and "dependent surfaces", where each dependent surface lies opposite to a corresponding control surface.To aid in comprehension, the corresponding material line elements at the control and dependent surfaces are represented as d   and d   , respectively, as shown in the Figure 4.The PBCs ensure that two opposing surfaces are identical, with the exception of a rigid body rotation of the material line elements, which depends on membrane deformation  and curvature tensor .This leads to the constraint (8).By enforcing this constraint, the behavior and deformation of the dependent surfaces are directly linked to those of the corresponding control surfaces, preserving the periodicity and allowing for the analysis or modeling of periodic structures.However, to extend the PBC to include transverse shear , the director (only the edge between 3 and B) is rotated It is important to emphasize that only surfaces with a Gaussian curvature det[  ] = 0 are considered as so-called developable surfaces.This means that the curved surface can be flattened without distortion.Please note that the sphere is the only closed surface with a constant Gaussian curvature.In general, undesired surface effects increase with increasing Gaussian curvature.Therefore, either small RVEs of highly curved surfaces or large RVEs of slightly curved surfaces can be analyzed with negligible surface effects.In the case of an air spring, the curvature of the rolling lobe is large compared to the curvature in the circumferential direction.Hence, the highly loaded areas can be modeled as a parabolic RVE.

Implementation
ABAQUS provides a user-defined subroutine called the multi-point constraint (MPC), which allows the imposition of nonlinear constraints.In order to specify the macroscale displacement gradient in a local orthogonal coordinate system   , three additional pilot nodes   ( = 1, 2) in-plane and  3 out-of-plane are introduced (see Figure 3), which are not physically connected to the RVE.The displacement degrees of freedom    and  3  (with   3 =  3 3 = 0) are utilized to define the coefficients of the membrane tensor  and the transverse shear  =     , respectively.Conversely, the rotational degrees of freedom    (with   3 =  3  = 0) are employed to determine the coefficients of the curvature tensor .

Verification of surface effects and validation
To verify the surface effects, multiple homogeneous second-order RVEs are connected in series.These RVEs are defined by dimensions specified as  [ ×  × ℎ], with the number of these RVEs denoted as , and each having defined length (), width (), and thickness (ℎ).These RVEs are loaded with a high curvature of  11 = 0.4 mm −1 .The surface effects refer to disturbances in the peripheral regions of the RVE.On the left side of the illustration 5 (the first image from the left), it is clearly evident that changes in the von Mises stress distribution occur not only in the thickness direction but also in other coordinate directions, indicating disturbances in the peripheral areas.When the RVE is reduced in size (as seen in the third image from the left), the disturbances in the peripheral region decrease despite constant curvature.This becomes particularly evident when the stress distribution changes exclusively in the thickness direction.From this, the following conclusion can be drawn: When analyzing a small RVE with significant curvature, the curvature is indeed significant, but the RVEs are small enough to assume that surface effects are relatively negligible.The smaller the RVE, the more the surface effect is reduced until it converges to a free surface effect.In contrast, when examining large RVEs with minimal curvature, the curvature is less pronounced, but the RVEs are large enough to justify the assumption that surface effects are negligible.In the case of a rotation-based RVE, the stress distribution is only free from thickness any surface effects in the parabolic scenario, regardless of the size of the RVE, as depicted in Figure 6 on the left.This is due to the strict constraints in Equation ( 8).This emphasizes the statements of the theorems by Minding and Liebmann, which explain the surface effects in elliptical and hyperbolic structures.Only surfaces with a Gaussian curvature of zero are considered "developable", meaning that the curved surface can be flattened without distortion.
In order to prove the correct implementation of the constraints and the validity of the homogenization, a homogeneous FE-plate made of elastomer is considered.Initially, a FE-plate with dimensions of 10 × 40 × 1 mm 3 for length, width, and thickness is modeled, using volume elements (5 elements in the thickness direction).The plate is assigned a soft isotropic material with Hooke's parameters  = 42 MPa,  = 0.4, and a curvature of  11 = 0.4 mm −1 is applied.The results of the macroscopic simulation show at a specific location the maximum normal stress in the tensile region of  11max = 9.32 MPa (on the upper surface) and the minimum normal stress of  11min = −10.12MPa in the compressive region (on the lower surface).At the same location, an RVE with dimensions of 1 × 1 × 1 mm 3 is extracted, and the same material and boundary condition parameters are applied.The results of the microscopic simulation of the RVE using both methods are presented in Figure 6.It is shown that the rotation-based RVE provides more accurate results, and the small differences observed are due to numerical differentiation and can be disregarded when using a fine mesh.Instead, the second-order RVE exhibits some deviation, which diminishes as the size is reduced.It is shown that the rotation-based RVE provides more accurate results, and the small deviation in results between the FE plate and the RVE is attributed to numerical differentiation and can be disregarded when using a fine mesh.Instead, the second-order RVE exhibits some deviation, which diminishes as the size of RVE is reduced.

Cross-layer reinforced composites of an air spring
In this section, shell-like rotation-based RVEs for cord-elastomer composites with two intersecting reinforcement layers are considered.The objective is to analyze the locale stresses (localization) in the composite under realistic loading conditions that occur in a component such as an air spring.This aims to enable efficient optimization of cord-elastomer composites, including their durability, under realistic conditions.Figure 7 (left) shows the geometry of an undeformed air spring bellows with its dimensions in a top view.The RVEs represent sections of the bellows with different periodicity and loading directions.Although the initial geometry of many air spring bellows is cylindrical, the RVEs are assumed to be initially straight due to the large radius compared to the bellows thickness and the high number of cords .The parameters of the cross-layer reinforced composite are cord angle   , the thickness of the composite ℎ, the spacing of the cords   , and the relative position of the cords in the thickness direction.The RVE features a structured mesh with C3D8H hybrid elements, and the material model assigned to each integration point is assumed to be influenced by its corresponding material region.A hyperelastic model is chosen for the rubber matrix.In particular, the strain energy potential of Neo-Hookean type is used, which is split into deviatoric and volumetric contributions.The second material region that needs to be considered is the cord, which has a filament structure.The cord bundle is considered to possess transversely isotropic behavior with the preferred direction represented by   .The total strain energy potential is composed of an anisotropic part and an isotropic part that accounts for shear and volumetric deformations [1].The next step involves developing a global FE model in ABAQUS using the rebar formulation.In this context, the global model refers to an FE model of the air spring that does not rely on any assumed symmetries.The aim of the global model is to accurately capture all typical macroscopic deformation modes present in the air spring bellows, enabling their transfer to the RVE.These deformation modes specifically include membrane deformation, curvatures, and transverse shear effects occurring within the matrix between two cord layers (see, ref. [8]), as illustrated in Figure 8.A transfer of the macroscopic deformations to the RVE as PBC allows for a more precise analysis of the local stress distribution.However, the use of multiphase elements inevitably entails a loss of information concerning the interface.Since the interface between the cords and the elastomer often represents preferred locations for crack initiation and, consequently, for later component failure, the analysis of stress in this region is of particular interest.One way to significantly improve the quality of results near the interface is through the application of submodeling techniques.Figure 9 illustrates the distribution of maximum logarithmic strain in both the RVE and the submodel.Based on the simulations, it has been observed that the primary influence in changing the cord angle is a result of the combined effects of curvature and transverse shear.This interaction leads to a significant shearing of the elastomer layer located between the cord layers.The subsequent step F I G U R E 9 The distribution of maximum logarithmic strain in RVE (left), and in a submodel of a cord, which is coated with an elastomer layer (middle), and effective shell stiffness (right).RVE, representative volume elements.
involves calculating an effective shell stiffness using a second-order RVE, as depicted in Figure 9.This calculation is necessary for homogenization the effective mechanical behavior of the corresponding shell element and, thus, developing an appropriate constitutive model.
Additionally, all contributions highlighted in red appear as discretization errors.Besides the tension-torsion coupling, the coupling between bending and shearing is particularly interesting.Due to the curvature of the RVE, tension and compression sides are formed.Instead of simply stretching or compressing, the cords can also change their angle in response.

SUMMARY
The rotation-based RVE is expanded to incorporate transverse shear and is highly suitable for localization and accurately analyzing the local stress in cord-elastomer composites.Its effectiveness lies in its ability to simulate significant curvature effects at different sizes without violating the periodicity conditions.For its discretization, a simple structured meshing approach was used, combined with an integration point-wise assignment of material properties, resulting in a suitable FE model for cord-elastomer composites.The loss of information due to the implicitly modeled interface between the materials was compensated for by conducting a subsequent analysis of a submodel of a cord, which is coated with an elastomer layer.However, the rotation-based RVE is unable to simulate arbitrary curvatures in different directions.The curvature tensor must satisfy certain conditions, such as having zero Gaussian curvature, in order to avoid surface effects.
In case of non-zero Gaussian curvature, like in air springs, the second-order RVE is perfectly suited for homogenizing the effective properties and coupling with shell elements in FE 2 analyses.However, it is important to choose an appropriate size to minimize surface effects.

A C K N O W L E D G M E N T S
Open access funding enabled and organized by Projekt DEAL.

F I G U R E 2
Mindlin-Reissner shell: reference state (left) and deformed state (right).

F I G U R E 3 F I G U R E 4
Second-order RVE.RVE, representative volume elements.Rotation-based RVE.RVE, representative volume elements.

F I G U R E 5
Von Mises stress in RVEs connected in series of different sizes.RVE, representative volume elements.using a rotation tensor , which determines a rotation angle and a rotation axis using the Euler-Rodrigues formula.

F I G U R E 6
In-plane stress  11 distribution resulting from the simple bending of a rotation-based RVE (left), and of a second-order RVE (right).RVE, representative volume elements.

F I G U R E 7
Cross-ply RVEs as sections of air spring bellows.RVE, representative volume elements.F I G U R E 8 the local macroscopic deformation modes in a rebar model.