A mixed variational‐based dynamic simulation method for fiber‐reinforced continua in non‐isothermal rotordynamical systems

Three‐dimensional woven carbon‐fiber reinforced polymers (3D‐CFRP) are being increasingly used for blades in rotordynamical systems as turbine fans, shafts as well as disc brake rotors. By using dynamic mixed finite‐element methods, these materials can be simulated as anisotropic continua with non‐isothermal behaviour in the matrix as well as the fiber material. Therefore, we present a novel dynamic mixed finite‐element method which is variational‐based and capable to simulate these materials. This dynamic mixed finite‐element method preserves each basic balance law emanating from the variational principle exact in a discrete setting. We demonstrate this physically‐consistent simulation behaviour by using dynamic simulations of a disc brake rotor and a rotating pipe subject to thermal and mechanical loads.


Introduction
The motivation of this paper is the increasing application of 3D-CFRP as meta-materials in rotating machinery. We find woven composites in turbine and pump rotors, disc brake rotors and shafts. A woven composite is a multiscale problem in space. Zooming in a preform of a composite part, we see a representative volume element, in which yarns are connected according to a weaving pattern. A yarn is made of thousands of carbon fibers, and therefore possesses a bending stiffness in contrast to one carbon fiber. We aim at dynamic finite element (FE) simulations of a body B 0 , which take into account this fiber bending stiffness, and also a micro inertia to obtain dispersive waves. Therefore, we design dynamic variational principles for continua with these multiscale effects, and derive from that new energy-momentum schemes. Our strategy is to consider generalized continuum mechanics with an unsymmetric Kirchhoff stress tensor. The reason is that anisotropic Cauchy continua, formulated with structural tensors A 0 , assume that fibers are strings with no bending stiffness [1]. But, a three-point bending test of a 3D-woven composite reveals, that the yarn bending stiffness has a large influence on the macroscopic scale. Beam ends are stiff and beam curvatures are large, in contrast to FE simulations with an anisotropic Cauchy continuum [2]. These curvature influences can be taken into account by tensor invariants based on the second-gradient of the deformation.

Variational setting
As basis of the variational formulation of a generalized Cauchy continuum, we consider the principle of virtual power for a Cauchy continuum in Ref. [3]. But, in the kinetic power, we add a rotational energy functional which introduces a micro-inertia based on length scale parameters l 0 and l F , given bẏ The vectors α, ω and π denote the independent fields of continuum rotation, angular velocity and rotational momentum, and I the unity tensor. In the external power, we prescibe Dirichlet rotationsᾱ and Neumann torquesW by adding the functional Here, we denote by τ t skw the transpose of the skew-symmetric Kirchhoff stress tensor, and by Z andẐ generalized reactions on the Dirichlet rotation boundary ∂ α B 0 . The set ∂ W B 0 indicates the Neumann torque boundary. The permutation tensor is denoted by ǫ. In the internal power, we introduce the deformation gradient and the first Piola-Kirchhoff stress tensor as independent fieldsF and P , respectively. This provides the application of a second-gradient free energy function. Hence, we obtain the internal poweṙ with the skew-symmetric projection tensor 2 I skw = ǫ · ǫ. We denote by S the second Piola-Kirchhoff stress tensor, and by S V and S F the independent volumetric and fiber stress, respectively. The last term in Eq. (3) defines the continuum rotation vector α, which is not present in Reference [3]. For more details about the remainder quantities, we may refer to Ref. [3].
After the variation with respect to temporally continuous time rate fields and temporally discontinuous Lagrange multipliers (cf. Ref. [3]), the considered generalized Cauchy continuum leads to additional weak forms in comparison to Ref. [3]: We obtain the weak forms for the continuum rotation vector in the first row of Eq. (4), and the weak balance of rotational momentum in the second row. Further, we obtain the weak forms for the rotational momentum and angular velocity in the third row, and the first Piola-Kirchhoff stress tensor and the deformation gradient in the last row. These additional weak forms leads to additional balance laws in comparison to Ref. [3]. First, we obtain the balance law of rotational momentum by setting the variation of the continuum rotation vector constant. This balance law is part of the angular momentum balance law. Second, we obtain the following balance law of rotational energy, which is finally part of the total energy balance law:

Numerical examples
In this paper, we show our first step: the application of the transversely isotropic first-gradient material law in Ref. [3] in the variational setting of a generalized Cauchy continuum defined above. We recognize that this formulation provides a Neumann boundary torque load. Therefore, we first demostrate this boundary condition by a dynamic torsion of a long composite square pipe with two cells, discretized by 20-node hexahedral finite elements (see Fig. 1).

Dynamic torsion of a composite square pipe
The fibers in the cell corners lie along the pipe axis, but, in the cell walls, the fibers are diagonal. On the pipe ends, we prescribe the time evolutions Θ A =Θ f (t) of nodal temperatures and nodal torque vectors W A x = ±Ŵ x f (t) along with the pipe axis. The time evolution of the used loads are given by the wave function f (t) in Ref. [3]. The initial conditions are homogeneous. In Fig. 1, red arrows indicate the nodal torque vectors, and colours the current temperature. Note the different temperature values on the top of the pipe. 3.2 Driven rotation of the composite square pipe subject to pressure and cooling As second example, we show a driven rotation of the same pipe. The pipe is free falling due to gravity. At the ends, we prescribe again the nodal temperatures and torques, but we prescribe at x = 0 a high drive torqueŴ drive x = 3.5 · 10 8 and at x = 100 a low friction torqueŴ friction x = 3.5 · 10 3 . The inner boundaries (blue patches in Fig. 2) are loaded by a  transient outward heat flux and a pressure follower load (cf. Ref. [3]). The initial conditions are again homogeneous. On the green patches, no boundary conditions are prescribed. By means of the angular velocity vector (see Fig. 2, first row), we recognize that the drive torque accelerate the pipe in steps until the torque is switched off. We also switched off the prescribed temperature evolution and heat flux (see Fig. 2, first row). The time evolutions of displacements and velocities are comparable with analytical solutions of accelerated rotors (see Fig. 2, last row).

Driven rotation of a composite disc brake rotor subject to traction and heating
Finally, we show a driven rotation of a composite disc brake rotor, in which the fibers lie in radial direction. We load the disc from below with a traction load, fix the inner ring in vertical direction, and apply a drive torque on the middle ring with the holes. On the cyan patches, we apply an inward heat flux load (heat flux during braking), and the inner ring represents a heat sink. The initial conditions are also homogeneous. In Fig. 3, we show the time evolutions of a node near a spoke of the disc. Inspite of the heat flux, the temperature is decreasing due to the heat conduction in the heat sink. The angular velocity is increasing due to the drive torque, and the time evolutions of displacements and velocities fits to an accelerated rotor.