Quasi-Static Deformations of Fibre-Reinforced Materials Based on Hyperelasticity

1. Introduction

Materials exhibiting inelastic behaviour are widely employed in engineering solutions within the field of civil engineering. Amongst the material types of interest are those sensitive to the rate of deformation and which, under standard service conditions, undergo finite elastic strains that exceed the range accurately captured by small-strain theory. These conditions apply in particular to elastomeric materials, utilised primarily as constituent components of bearings designed to attenuate vibrations. The work [1], for example, presents representative results of typical experimental tests performed on specimens of hard and soft elastomeric products.. The loading path on the nominal stress–strain diagram does not coincide with the unloading path, indicating energy dissipation. Furthermore, loading–unloading cycles performed at different displacement-controlled deformation rates yield quantitatively distinct results. This material characteristic is referred to as viscosity [2], which is a property intrinsic to fluids but may equally pertain to solid bodies. An example of a structure exhibiting comparatively large deformations is a roof comprising a textile membrane. A structure of this type was analysed in [3,4] employing fibre-reinforced hyperelastic material models, though without accounting for dissipative effects.

The rheological characterisation of materials, that is, their sensitivity to the strain rate (viscosity) [5] and the Mullins effect [6], has been presented in various forms within the literature on the mechanics of solids undergoing large deformations [7,8]. The Mullins effect is typically represented by introducing a parameter that defines the maximum elastic energy attained during loading [6,9], or alternatively by a parameter that specifies the maximum deformation intensity [1,10]. This constitutes a particular case of a model formulated within the framework of the internal variable theory [11,12], which provides a formal structure for incorporating inelastic effects in compliance with the second law of thermodynamics [13,14]. Another approach invokes the theory of materials with memory [15,16], in which approximations to the so-called response functional are introduced. One such approximation takes the integral form [17,18], leading to the class of quasi-linear viscoelastic (QLV) models [19,20]. Models of this type have also been developed for transversely isotropic and fibre-reinforced materials [21,22]. Since conventional models describing the Mullins effect do not account for strain-rate sensitivity, they can be generalised, for example, within the QLV framework [23,24].

In the present work, attention is directed towards finite element implementation and solutions to selected boundary value problems employing the fibre-reinforced material model described in the authors’ previous works [25,26]. It is emphasised that the model satisfies the requirements

1.

the requirement of objectivity associated with Galilean transformations and the given symmetry of the material in the initial configuration [27],

2.

thermodynamic admissibility through the restriction resulting from the Clausius-Duhem inequality [28],

3.

polyconvexity and growth conditions of the stored energy function [29].

The objectivity requirements in the context of constitutive relations pertain primarily to the function determining the stress tensor and the Helmholtz free energy function $\mathrm{\Psi }$. Adopting the material description, the following conditions must hold

$$\begin{array}{cc}\hfill \mathrm{\Psi }& =\widehat{\mathrm{\Psi }}(\mathbf{C},\alpha ),\hfill \\ \hfill \mathbf{T}& =\overline{\mathbf{T}}(\mathbf{X},t)=\widehat{\mathbf{T}}(\mathbf{C},\alpha ),\hfill \end{array}$$

(1)

where $\mathbf{F}$ is the deformation gradient tensor, $\mathbf{C}={\mathbf{F}}^T\mathbf{F}$ is the right Cauchy-Green tensor, $\alpha$ is the tensor internal variable, and $\mathbf{T}$ is the second Piola-Kirchhoff stress tensor.

In addition to the objectivity requirement (Euler’s) related to motion and the observer of the current configuration, the objectivity requirement (Lagrange’s) of the material symmetry given in the initial configuration is also imposed, therefore the rotation tensor $\mathbf{Q}$ here refers to the rotation of the material coordinate system. We then say that the group $\mathcal{G}$ is a group of material symmetry if

$$\mathrm{\Psi }=\widehat{\mathrm{\Psi }}(\mathbf{C},\alpha )=\widehat{\mathrm{\Psi }}({\mathbf{Q}}^T\mathbf{C}\mathbf{Q},{\mathbf{Q}}^T\alpha \mathbf{Q}),\forall \mathbf{Q}\in \mathcal{G}\subset SO\left(3\right).$$

(2)

When $\mathcal{G}=SO\left(3\right)$ then the material is isotropic. Otherwise, we call the material anisotropic. We can describe the symmetry of the material through parametric tensors and, consequently, $\mathrm{\Psi }$ as a function of, e.g., anisotropic invariants of a symmetric second-order tensor [30,31]. If the parametric tensor is a symmetric second-order tensor $\overline{\mathbf{M}}$ [32], then the symmetry group

$$S_o=\left\{\mathbf{Q}\in O\left(3\right)\mid \overline{\mathbf{M}}=\mathbf{Q}\overline{\mathbf{M}}{\mathbf{Q}}^T\right\}$$

(3)

defines an orthotropic material in the initial configuration of the body in its natural state [31,33].

In the present work, attention is restricted to the special case of a transversely isotropic material within the framework of hyperelasticity. The model of a transversely isotropic material contains a parametric tensor in the form $\mathbf{M}=\mathbf{m}\otimes \mathbf{m}$, where the unit vector $\mathbf{m}$ is consistent with the direction of the distinguished fibres of the material in the initial configuration [34,35]. Therefore, the stored energy function can be represented in the form of five invariants

$$\widehat{W}\left(I_i\right)=\widehat{W}(tr\mathbf{C},trCof\mathbf{C},det\mathbf{C},tr\left(\mathbf{MC}\right),\mathrm{tr}\left(\mathbf{M}Cof\mathbf{C}\right)).$$

(4)

From now on, we restrict ourselves to a purely mechanical theory, i.e., neglecting all thermal effects. Fields such as temperature and entropy are omitted in the balance equations. Nevertheless, such an approach is consistent with the second law of thermodynamics by introducing the so-called mechanical dissipation [27]. Therefore, the requirements resulting from the principles of thermodynamics are reduced to satisfying the Clausius-Duhem inequality [28] in the form

$$D={\displaystyle \frac{1}{2}}\mathbf{T}·\dot{\mathbf{C}}-\mathrm{\Psi }\ge 0$$

(5)

Therefore, we determine the material time derivative of free energy with the introduction of one internal variable

$$\dot{\mathrm{\Psi }}={\displaystyle \frac{\partial \widehat{\mathrm{\Psi }}}{\partial \mathbf{C}}}·\dot{\mathbf{C}}+{\displaystyle \frac{\partial \widehat{\mathrm{\Psi }}}{\partial \alpha }}·\dot{\alpha }.$$

(6)

Substituting (6) into (5) we obtain

$$D=\left({\displaystyle \frac{1}{2}}\mathbf{T}-{\displaystyle \frac{\partial \widehat{\mathrm{\Psi }}}{\partial \mathbf{C}}}\right)·\dot{\mathbf{C}}-{\displaystyle \frac{\partial \widehat{\mathrm{\Psi }}}{\partial \alpha }}·\dot{\alpha }\ge 0.$$

(7)

The inequality should hold for all admissible $\dot{\mathbf{C}},\dot{\alpha }$. Therefore

$$\mathbf{T}-2{\displaystyle \frac{\partial \widehat{\mathrm{\Psi }}}{\partial \mathbf{C}}}{|}_{\mathbf{C}={\mathbf{C}}^T}=\mathbf{0},D=-{\displaystyle \frac{\partial \widehat{\mathrm{\Psi }}}{\partial \alpha }}·\dot{\alpha }\ge 0.$$

(8)

Apart from the requirements resulting from the second law of thermodynamics, objectivity and material symmetry, additional requirements on the constitutive relations of hyperelasticity should be introduced so that a well-defined initial-boundary value problem of equations of motion has a solution [36]. Here, we primarily require the stored-energy function to be polyconvex and satisfy appropriate growth conditions [29,37]. Since the invariant $I_4=\mathrm{tr}\left(\mathbf{C}\mathbf{M}\right)$ is a convex function with respect to $\mathbf{F}\in Lin$, while the invariant $I_5=\mathrm{tr}\left(\mathbf{M}Cof\mathbf{C}\right)$ is a convex function with respect to $Cof\mathbf{F}\in Lin$ then the function in the form

$$\widehat{W}\left(I_i\right)=W_1\left(I_1\right)+W_2\left(I_2\right)+W_3\left(I_3\right)+W_4\left(I_4\right)+W_5\left(I_5\right)$$

(9)

is polyconvex if the functions $W_1\left(I_1\right)$ and $W_4\left(I_4\right)$ are convex functions with respect to the tensor $\mathbf{F}\in Lin$, the functions $W_2\left(I_2\right)$ and $W_5\left(I_5\right)$ are convex functions with respect to $Cof\mathbf{F}\in Lin$, while the function $W_3\left(I_3\right)$ is convex with respect to $det\mathbf{F}>0$. It is known, however, that the condition is too strong for several reasons and should be replaced by a weaker condition than polyconvexity, namely quasiconvexiy [38]. However, there is no known useful description of quasi-convexity under the local growth condition, i.e. $W\left(\mathbf{F}\right)\to \infty$ as $det\mathbf{F}\to 0^{+}$.

The article is structured as follows: Section 2 develops the constitutive framework for fibre-reinforced hyperelastic materials with a scalar internal variable, including the material model and its quasi-static characterisation. Section 3 describes the finite element implementation in Abaqus/Standard via a UMAT formulation and the associated rate-form constitutive equations. Section 4 presents numerical results for homogeneous tests and a non-uniform inflation-like problem, illustrating dissipation, rate effects and instability-like responses. Section 5 summarises the main conclusions and outlines directions for future work, particularly regarding permanent deformations and dynamic applications.

2. Fibre-Reinforced Material Models

2.1. General Framework

Fibre-reinforced material models are formulated by an additive decomposition of the Helmholtz energy [39,40] in the form

$$\mathrm{\Psi }=(1-\sum _{n=1}^N{v}_n){\mathrm{\Psi }}_M+\sum _{n=1}^N{v}_n{\mathrm{\Psi }}_{Zn},N\le 6,$$

(10)

The inelastic behaviour is introduced through an internal variable represented by a scalar quantity, analogous to the constitutive relations of material models formulated within the framework of continuum damage mechanics (CDM)[41]. The principal distinction between standard CDM models and the model presented here lies in the ability to describe the reversibility of the material damage process[42]. We propose a Helmholtz free-energy function, ${\mathrm{\Psi }}_M$ , analogous to that adopted in [43], where wave-propagation problems were studied using the isotropic Murnaghan constitutive model. If we have

$${\mathrm{\Psi }}_M=\widehat{\mathrm{\Psi }}(\mathbf{C},g)={\varphi }_1\left(g\right)\widehat{W}\left(\mathbf{C}\right)+{\varphi }_2\left(g\right),\mathbf{C}=\widehat{\mathbf{C}}(\mathbf{X},t),g=\widehat{g}(\mathbf{C},t),$$

(11)

then according to (7) we obtain

$$D=\left({\displaystyle \frac{1}{2}}\mathbf{T}-{\displaystyle \frac{\partial \widehat{\mathrm{\Psi }}}{\partial \mathbf{C}}}\right)·\dot{\mathbf{C}}+P\dot{g}\ge 0,$$

(12)

and

$$\mathbf{T}=2{\varphi }_1\left(g\right){\displaystyle \frac{\partial \widehat{W}\left(\mathbf{C}\right)}{\partial \mathbf{C}}}{|}_{\mathbf{C}={\mathbf{C}}^T},P=-{\displaystyle \frac{\partial \widehat{\mathrm{\Psi }}(\mathbf{C},g)}{\partial g}}=-{\displaystyle \frac{\partial {\varphi }_1\left(g\right)}{\partial g}}\widehat{W}\left(\mathbf{C}\right)-{\displaystyle \frac{\partial {\varphi }_2\left(g\right)}{\partial g}}.$$

(13)

We assume the evolution equation of the variable (g) to satisfy the inequality (12)

$${\displaystyle \frac{\partial {\varphi }_1\left(g\right)}{\partial g}}\widehat{W}\left(\mathbf{C}\right)+{\displaystyle \frac{\partial {\varphi }_2\left(g\right)}{\partial g}}+\alpha \dot{g}=0.$$

(14)

Indeed, if $\alpha >0$ then

$$D=P\dot{g}=\alpha {\dot{g}}^2\ge 0.$$

(15)

In the case of $\alpha =0$ or $\alpha \to +\infty$ the model does not predict dissipation, i.e., it describes a thermodynamically reversible process. From the thermodynamic requirements it additionally follows that

$$\left({\displaystyle \frac{\partial {\varphi }_1\left(g\right)}{\partial g}}\widehat{W}\left(\mathbf{C}\right)+{\displaystyle \frac{\partial {\varphi }_2\left(g\right)}{\partial g}}\right){|}_{\dot{g}=0}=0.$$

(16)

We assume that ${\varphi }_1>0$ and for (g=0) we obtain a hyperelastic relation, therefore ${\varphi }_1\left(0\right)=1,{\varphi }_2\left(0\right)=0$. Given the assumption of a natural state with zero dissipation as an initial condition $t=0\to \mathbf{C}=\mathbf{I}$ and (g=0), then (14) leads to

$$\widehat{W}\left(\mathbf{I}\right)=0\Rightarrow {\displaystyle \frac{\partial {\varphi }_2\left(g\right)}{\partial g}}{|}_{g=\widehat{g}(\mathbf{I},0)=0}=0.$$

(17)

The function associated with the anisotropic part ${\mathrm{\Psi }}_{Zn}$ of the energy is assumed in the form

$${\mathrm{\Psi }}_{Zn}={\tilde{W}}_Z\left(I_{4n}\right),$$

(18)

where $I_{4n}=\mathrm{tr}\left({\mathbf{M}}_n\mathbf{C}\right)$. As explained in the previous section of the article, the parametric tensors ${\mathbf{M}}_n\equiv {\mathbf{m}}_n\otimes {\mathbf{m}}_n$ are defined by the unit vectors ${\mathbf{m}}_n$ associated with the distinguished directions of the fibre families. For $n=1$, function (10) leads to a simplified model of a transversely isotropic material [44,45,46].

The constitutive relations of the presented model classes can be interpreted as a composite material composed of a viscoelastic matrix and elastic fibres. It should be noted that in an analogous way we could formulate equations describing a model with the interpretation of an elastic matrix and viscoelastic fibres.

2.2. Material Model

The function $W=\widehat{W}\left(\mathbf{C}\right)$ appearing in the previous equations pertains to hyperelasticity. Here, we adopt an isotropic function $W=\tilde{W}\left(I_1,I_2,J\right)$ in the form of a generalised Ciarlet model [47,48]

$$\widehat{W}\left(\mathbf{C}\right)=\tilde{W}\left(I_1,I_2,J\right)=a\left(I_1-3\right)+b\left(I_2-3\right)+c{\left(I_1-3\right)}^2+\beta \left(J\right)$$

(19)

where $\beta$ is a convex function. From the assumption of a natural state it follows that

$${\displaystyle \frac{\partial \beta \left(J\right)}{\partial J}}{|}_{J=1}=-2\left(a+2b\right).$$

(20)

In relation to the basic Ciarlet model, i.e., when $c=0$[49], function (19) predicts, above all, a different qualitative response to a given deformation associated with shear. We assume the function $\beta$ in the form

$$\beta \left(J\right)=d\left(J^2-1\right)-2(a+2b+d)lnJ.$$

(21)

The existence theorem of a solution in hyperelasticity requires that the parameter values satisfy $a,b,d>0$ and $c\ge 0$. The stored energy function is then polyconvex and satisfies the growth requirements.

If we consider moderately large deformations, then the function (19) can be considered in a quadratic approximation with respect to the strain tensor $\mathbf{E}$[49]. Then, the initial bulk modulus $K_0$ and Lamé’s constant ${\lambda }_0$ can be expressed through parameters such that

$$K_0={\displaystyle \frac{4}{3}}(a+4b+6c+3d),{\lambda }_0=4(b+2c+d).$$

(22)

If we additionally assume that $c=0$ and

$$a={\displaystyle \frac{1}{2}}f{\mu }_0,b={\displaystyle \frac{1}{2}}(1-f){\mu }_0,f\in (0,1),$$

(23)

then we obtain

$$d={\displaystyle \frac{1}{4}}{K}_0-{\displaystyle \frac{1}{3}}\left(a+4b\right)>0.$$

(24)

We can also assume the parameter (c) in the form $c=h{\mu }_0>0$. Therefore, the model parameters can be represented by elastic constants. If we additionally assume that $K_0={\overline{k}}_0{\mu }_0$ then the stored energy function is scalable by ${\mu }_0$ and we can use a set of three dimensionless parameters $f,h,{\overline{k}}_0$, because

$${\displaystyle \frac{a}{{\mu }_0}}={\displaystyle \frac{1}{2}}f,{\displaystyle \frac{b}{{\mu }_0}}={\displaystyle \frac{1}{2}}(1-f),{\displaystyle \frac{c}{{\mu }_0}}=h,{\displaystyle \frac{d}{{\mu }_0}}={\displaystyle \frac{1}{12}}(6f+3{\overline{k}}_0-24h-8).$$

(25)

The parameter values must satisfy the constraints

$$f\in (0,1),{\overline{k}}_0>0,h\ge 0,6f+3{\overline{k}}_0-24h-8>0.$$

(26)

From the point of view of polyconvexity [32], we consider the function ${\tilde{W}}_Z$ as follows:

$${\tilde{W}}_Z\left(I_4\right)={\displaystyle \frac{1}{4}}{E}_Z{(I_4-1)}^2,$$

(27)

where $E_Z>0$ can be interpreted as the initial Young’s modulus of the fibre family.

In view of the requirements established above, the functions ${\varphi }_1,{\varphi }_2$ pertaining to the dissipative part of the material model may be taken [42] in the form

$${\varphi }_1\left(g\right)=(1-g)\vartheta ,{\varphi }_2\left(g\right)=-{\displaystyle \frac{1}{2}}\gamma ln(1-g^2),g\in (-1,1),$$

(28)

where the parameters $\gamma ,\vartheta$ are positive scalars and $\gamma$ has the dimension of energy per unit volume in the initial configuration.

In the remainder of this work, the term ,,FRD model” shall denote the constitutive relation defined by the Helmholtz function (10) together with the fibre description (27). The dissipative part of the Helmholtz function within this model is understood as (11) in conjunction with the generalised Ciarlet model (19) and (28).

2.3. Quasi-Static Problem with Constant Deformation Rate

To illustrate some fundamental characteristics of the model under consideration in the isotropic part, in this subsection we examine a quasi-static deformation problem with a constant deformation rate over time-scale intervals. The deformation type considered is triaxial compression/tension, such that $\mathbf{F}\left(t\right)={\lambda }_1\left(t\right)\mathbf{I}$. In this case, the Cauchy stress tensor is also isotropic and is denoted as $\sigma \left(t\right)={\sigma }_{11}\left(t\right)\mathbf{I}$. The stretch ratio ${\lambda }_1$ is defined according to the plot shown in Figure 1.

The problem is solved assuming the parameter values $f=0.35$, $h=0.1$, and ${\overline{k}}_0=2$ or ${\overline{k}}_0=20$. These values were selected to obtain less conventional solutions than those reported in the literature on damage mechanics [50,51]. Since the evolution equation for the variable g is non-linear, the solution is determined using the procedure NDSolve in the Wolfram Mathematica environment [52].

The plots in Figure 2 show the obtained values of the rescaled stress ${\sigma }_{11}/{\mu }_0=\widehat{\sigma }\left({\lambda }_1\right)$ for ${\overline{k}}_0=2$ and ${\overline{k}}_0=20$, respectively, in the case of an stretch defined on two intervals according to Figure 1(a), that is, a single loading/unloading cycle. The pronounced difference in the results arises not only from the fact that, with increasing ${\overline{k}}_0$, a stiffer material response is modelled, but also because the larger values of elastic energy significantly affect the solution for the internal variable g.

The obtained values of the internal variable g and the rescaled stress ${\sigma }_{11}/{\mu }_0$ as functions of ${\lambda }_1$ for the prescribed stretch shown in Figure 1(b) are presented in Figure 3 and Figure 4. From the form of the evolution equation for the variable g, it is known that the solution depends directly on the value of the stored-energy function, which can be observed in the behaviour of the plot of g. For ${\overline{k}}_0=20$, the values of elastic energy are significantly larger than for ${\overline{k}}_0=2$, which results in a much greater reduction of the obtained stress between the values in the time intervals $L_1$ and $L_2$. This may be interpreted as a more pronounced damage of the modelled material.

3. Material Model Implementation

3.1. Abaqus/Standard User Subroutines

Abaqus provides multiple means of implementing constitutive relations for user-defined material models [53]. The primary tools are the so-called user subroutines, written in Fortran or C/C++, which interact with the core Abaqus module through well-defined interfaces. These subroutines may either modify standard procedures or introduce, amongst other things, custom material models. From the standpoint of hyperelasticity, the subroutines UHYPER and UANISOHYPER_INV are of particular relevance, as they enable the definition of constitutive relations for isotropic and anisotropic materials by specifying the stored energy density function together with its derivatives with respect to the appropriate deformation invariants [54]. On the basis of these quantities, the programme determines the components of the Cauchy stress tensor $\sigma$ and the components of the tangent stiffness tensor ${\mathfrak{C}}_B$, accounting for the objectivity requirement.

$$\dot{\tau }=J\left(\mathbf{W}\sigma +\sigma {\mathbf{W}}^T+{\mathfrak{C}}_B·\mathbf{D}\right).$$

(29)

In Abaqus/Standard, the Zaremba–Jaumann rate is employed [55]. For the Kirchhoff stress tensor, one has

$$\dot{\tau }=\stackrel{\circ }{\tau }+\mathbf{W}\tau -\tau \mathbf{W}=\stackrel{\circ }{\tau }+\mathbf{L}\tau +\tau {\mathbf{L}}^T-(\mathbf{D}\tau +\tau \mathbf{D}).$$

(30)

For constitutive relations incorporating dissipative effects, the UMAT subroutine interface is employed. It requires the definition of the components of the Cauchy stress tensor and the components of the tangent stiffness tensor $\partial \Delta \sigma /\partial \Delta ϵ$ arising from the constitutive relation. If the tangent stiffness tensor is understood as ${\mathfrak{C}}_B$ in (29), the implementation is consistent with the hyperelastic equations, and such an approach is referred to as the total formulation [55]. In this case, finite elements with a mixed formulation may be employed, though the definition of appropriate third-order derivatives is required in a manner analogous to the UHYPER procedure.

3.2. Rate Form of Constitutive Equations

The incremental form of the constitutive relation of FRD model, incorporating the Zaremba–Jaumann objective rate, may be expressed as

$$\stackrel{\circ }{\tau }=\mathbf{C}·\mathbf{D}={\varphi }_1\left(g\right){\stackrel{\circ }{\tau }}_M+{\dot{\varphi }}_1\left(g\right){\tau }_M+{\stackrel{\circ }{\tau }}_Z,$$

(31)

where

$${\stackrel{\circ }{\tau }}_M=\left(1-\sum _{n=1}^N{v}_n\right){\overline{\mathfrak{C}}}_M·\mathbf{D},{\stackrel{\circ }{\tau }}_Z=\sum _{n=1}^N{v}_n{\mathfrak{C}}_{Z_n}·\mathbf{D},N\le 6.$$

(32)

The tangent stiffness tensors ${\overline{\mathfrak{C}}}_M,{\mathfrak{C}}_{Zn}$, associated with the generalised Ciarlet model and the polynomial of the invariant $I_{4n}$, were determined, for example, in [40]. Consequently, the evaluation of tensor $\mathbf{C}$ in (31) additionally requires the determination of the expression ${\dot{\varphi }}_1\left(g\right){\tau }_M$. From the adopted form of the evolution equation for the internal variable, it follows directly that $\dot{g}$ may be treated as a function of the variable $W_M$. It therefore follows that

$$\begin{array}{cc}\hfill & {\dot{\varphi }}_1\left(g\right){\tau }_M={\displaystyle \frac{\partial {\varphi }_1}{\partial g}}{\displaystyle \frac{\partial g}{\partial W_M}}{\dot{W}}_M{\tau }_M=\left[g_W{\varphi }_1^{\prime }\left(g\right){\tau }_M\otimes {\tau }_M\right]·\mathbf{D},\hfill \\ \hfill & g_W\equiv {\displaystyle \frac{\partial g}{\partial W_M}},{\varphi }_1^{\prime }\left(g\right)\equiv {\displaystyle \frac{\partial {\varphi }_1\left(g\right)}{\partial g}},\hfill \end{array}$$

(33)

together with

$${\varphi }_1^{\prime }\left(g\right)+\alpha {\dot{g}}_W=0,W_M=(1-v)W_0,v=\sum _{n=1}^N{v}_n.$$

(34)

The implementation of equation (34) is carried out by applying the forward Euler method over integration steps $s_{n+1}=s_n+\Delta s$, employing

$$g_{W,k+1}=g_{W,k}-{\displaystyle \frac{\Delta s}{\alpha }}{\varphi }_{1,k}^{\prime }.$$

(35)

If ${\varphi }_1\left(g\right)=1-g$, the derivative ${\varphi }_1^{\prime }$ is constant, namely ${\varphi }_1^{\prime }\left(g\right)=-1$. In this case, a straightforward relation is obtained

$$g_{W,k+1}=g_{W,k}+{\displaystyle \frac{\Delta s}{\alpha }}.$$

(36)

4. Results

4.1. Numerical Validation

4.1.1. Hyperelasticity

The correctness of the implementation of the constitutive relation of FRD model via the UMAT subroutine was initially verified through uniaxial strain and simple shear tests on a single C3D8R finite element, without accounting for inelastic effects. All boundary conditions were specified as displacement-controlled. The parameters of the hyperelastic model are adopted as ${\mu }_0=0.5\mathrm{MPa},k_0=2{\mu }_0,f=0.7,h=0.1,v=0.05,E_z/{\mu }_0=10$.

The computed stress values are in agreement with those obtained from the analytical expressions, within the limits of the prescribed numerical accuracy. Figure 5 and Figure 6 present plots of the rescaled components of the Cauchy stress tensor.

4.1.2. Fibre-Reinforced Model with Dissipative Effects

The next step in the verification of the implementation of FRD model is the inclusion of viscous effects. To this end, relaxation and creep tests were carried out, adopting in the former the parameter values $\alpha =1,\gamma =1$, and in the latter $\alpha =0.5,\gamma =0.2$.

The relaxation problem concerns a uniaxial strain state with displacement boundary conditions. The deformation gradient tensor is taken in the form $\mathbf{F}=2{\mathbf{b}}_1\otimes {\mathbf{b}}_{1}H\left(t\right)$, and the fibre direction as $\mathbf{m}={\displaystyle \frac{1}{\sqrt{2}}}({\mathbf{b}}_1+{\mathbf{b}}_2$. In the creep problem, a plane stress is considered, i.e. the stress state is defined by the tensor $\sigma ={\sigma }_{11}{\mathbf{b}}_1\otimes {\mathbf{b}}_1$, and the deformation state by $\mathbf{C}={\left({\lambda }_{11}\right)}^2{\mathbf{b}}_1\otimes {\mathbf{b}}_1+{\left({\lambda }_{22}\right)}^2({\mathbf{b}}_2\otimes {\mathbf{b}}_2+{\mathbf{b}}_3\otimes {\mathbf{b}}_3)$. It is assumed that the “11” stress component is controlled by prescribing the traction vector on the undeformed surface of the finite element, ${\mathbf{t}}_s=2{\mu }_{0}H\left(t\right){\mathbf{b}}_1$.

Figure 7 presents plots of the internal variable and the rescaled, non-zero components of the Cauchy stress tensor in the relaxation problem. It is worth noting that the value of the ${\sigma }_{12}$ component does not change over the considered time scale, since the constitutive relation of model 1af in the anisotropic part does not provide for relaxation.

In the creep problem, the fibre direction is taken to be aligned with the vector $\mathbf{m}={\mathbf{b}}_1$. The resulting values of the internal variable and the stretches as functions of time are shown in Figure 8 and Figure 9. The results obtained with Abaqus are consistent with those computed by means of the NDSolve procedure in Mathematica.

4.2. Non-Uniform Deformations

An inflation-like problem is considered for a geometry in the form of a circle with radius R and thickness $R/20$, subject to the boundary conditions $u_{2\left(y\right)}=0,u_{3\left(z\right)}=0$. The analysis is carried out for the FRD model with the parameter values $\alpha =0.2,\gamma /{\mu }_0=0.01,{\mu }_0=0.5\mathrm{MPa},k_0=2{\mu }_0,f=0.7,h=0.1,v=0.05,E_z/{\mu }_0=10$. The fibre family direction is assumed to be aligned with the x-direction. Since the implemented UMAT procedure is dedicated to a three-dimensional stress state, a rectangular, eight-node, reduced-integration brick element of type C3D8R is employed. This element type is not associated with a shell-theory formulation [56]. The loading and unloading are applied by prescribing a pressure on the upper surface of the circle, as illustrated in Figure 10. The discretised domain consists of 44 finite elements.

Figure 11 shows plots of the normalised pressure as a function of the displacement of the region midpoint. It is worth noting that, in the final unloading phase, a “jump” of the displacement value to positive values occurs, where the next cycle of loading is also presented. It may be interpreted as a form of instability; see [57]. The rescaled stress values ${\sigma }_{11}/{\mu }_0$ as a function of stretch in the loading–unloading process are presented in Figure 11. These values were extracted from one of the central elements and clearly indicate energy dissipation. The internal variable, evaluated for several selected elements, attains maximum values between 0.15 and 0.30; see Figure 12.

5. Conclusions

A proposed generalisation of the constitutive relations of hyperelastic models, incorporating fibre reinforcement and viscous effects, involves the introduction of an internal variable in the form of a scalar quantity, analogous to the constitutive relations of material models formulated within the framework of continuum damage mechanics (CDM). These models belong to the class of fibre-reinforced material models [35], formulated additionally within the framework of a simple mixture theory [39]. The form of the isotropic component of the Helmholtz free-energy function was adopted from [42] and subsequently extended to encompass the generalised Ciarlet model and fibre reinforcement, yielding a transversely isotropic material model incorporating dissipative effects [MMS]. The principal distinction between CDM models and the class of models presented in this work lies in the ability to describe the reversibility of the material damage process. As demonstrated in the boundary-value problems discussed herein, besides modelling the phenomenon of material damage, the constitutive relations within this class also capture strain-rate sensitivity as well as creep and stress-relaxation phenomena. The principal limitation of the models presented lies in their inability to capture permanent deformations. This aspect, together with the investigation of dynamic effects in initial–boundary value problems based on the proposed material model, will be addressed in future work.