Analysis and Numerical Simulation of the Behavior of Composite Materials with Natural Fibers Under Quasi-Static Frictional Contact

Abstract

This paper analyzed the behavior of polymer composite materials reinforced with randomly oriented short natural fibers (hemp, flax, etc.) subjected to external stresses under quasistatic contact conditions with dry Coulomb friction. We presumed the composite body, a 2D flat rectangular plate, being in frictional contact with a rigid foundation for the quasistatic case. The manuscript proposes the finite element method approximation in space and the finite difference approximation in time. The problem of quasistatic frictional contact is described with a special finite element, which can analyze the state of the nodes in the contact area, and their modification, between open, sliding, and fixed contact states, in the analyzed time interval. This finite element also models the Coulomb friction law and controls the penetrability according to a power law. Moreover, the quasi-static case analyzed allows for the description of the load history using an incremental and iterative algorithm. The discrete problem will be a static and nonlinear one for each time increment, and in the case of sliding contact, the stiffness matrix becomes non-symmetric. The regularization of the non-differentiable term comes from the modulus of the normal contact stress, with a convex function and with the gradient in the sub-unit modulus. The non-penetration condition was achieved with the penalty method, and the linearization was conducted with the Newton–Raphson method.

1. Introduction

The contact between two solid bodies made of composite materials or between a composite body and a rigid foundation is a complex physical process that describes the opposing resistance to penetration, in the normal direction, and a relative tangential motion. The interaction between the two bodies made of composite materials in contact can be described by the laws of the contact phenomenon and the laws of friction, which describe the relationships between the displacements and the forces acting on the contact area.

In the analyzed time interval, the frictional states on the contact boundary can change among open, sliding, or fixed. This paper modeled the contact problem using the quasi-static case and used incrementally-iterative algorithms. In the specialized literature, there are many incremental approaches to elastic contact problems.

The manuscript analyzed a composite made up of 50% polymer and 50% short natural fibers, ±10%, randomly oriented in a homogeneous mixture. The properties of this composite material make it highly appreciated in industry, especially in the field of furniture and household appliances.

The roughness of the contact area and implicitly the behavior of the contact phenomenon depend on the nature of the composite material. There are studies in the specialized literature for this type of composite material that have addressed the behavior of contact problems with different friction laws and analyzed fiber-reinforced composites taking into account their tribological characteristics. The tribological characteristics of the composites are related to the fiber to polymer volumetric ratio and the roughness of the contact area. These relations can be expressed when analyzing the contact with friction behavior of the composite. In this regard, we mention [1,2,3,4].

The friction coefficient on the contact area depends on the proportions of the components in the composite material and its rigidity. The contact area in composite materials is much rougher than in the case of metals. In addition, the contact area will contain randomly distributed soft components (natural hemp fibers) and hard components (the polypropylene matrix).

In the contact problem, the anisotropy of the contact area is more important than the anisotropy of the composite material [5,6,7]. There are still many open problems in the study of composite materials with randomly oriented natural fibers.

We considered the composite as homogeneous (i.e., the polymer and the randomly oriented short natural fibers) well mixed. For the mathematical model, we assumed the composite body to be viscoelastic with linear elasticity. The mathematical model consisted of a series of variational inequalities that were approximated with the finite element method [8,9,10,11], in which the non-differentiable term given by the law of dry friction contact according to Coulomb’s law is replaced by a friction functional regularized with a convex function. The regularization method consists of the approximation of the non-differentiable term in the contact condition with a differentiable term.

2. Materials and Methods

2.1. Quasistatic Formulation of Contact Problems with Friction for Composite Materials

Contact problems for composite materials consist of determining the displacements, deformations, and stresses in two viscoelastic, isotropic linear bodies in contact, or of a composite body in contact with a rigid support, subjected to various friction laws.

In this manuscript, we analyzed the second variant when the traction forces act at a speed sufficiently low that the inertial effects can be neglected.

Let there be a viscoelastic composite body which at the moment t = 0 is set within the domain $\Omega \subset R^d$, where d = 2 or d = 3, with the Lipschitz boundary, $\partial \Omega ={\overline{\Gamma }}_u\cup {\overline{\Gamma }}_N\cup {\overline{\Gamma }}_C$ shared in three mutually disjunctive parts ${\Gamma }_u,{\Gamma }_N,{\Gamma }_C$ having the measure $\left({\Gamma }_u\right)>0.$

This is denoted by $Q=\Omega \times \left[0,T\right]$, and $T>0$, $S_u={\Gamma }_u\times \left[0,T\right]$, and $S_C={\Gamma }_C\times \left[0,T\right]$. The body is fixed on the boundary ${\Gamma }_u$, where the displacements are set to zero, that is, the body is embedded in that portion of the boundary, having tractions of intensity h acting on the boundary ${\Gamma }_N$, and the ${\Gamma }_C$ boundary being the area in contact with the rigid foundation. Here, we applied the Signorini condition and Coulomb friction law. Volume forces with the intensity f act in Ω. We have noted the following: $\mathit{u}=\left(u_1,\dots ,u_d\right)$ displacement vector, $\mathit{\sigma }=\left\{{\sigma }_{ij}\right\}$ tension tensor, $\mathit{\epsilon }\left(\mathit{u}\right)=\left\{{\epsilon }_{ij}\left(\mathit{u}\right)\right\}$, deformation tensor $i,j=1,\dots ,\mathrm{a}\mathrm{n}\mathrm{d}{\mathit{\tau }}_{\mathit{N}},{\mathit{\tau }}_{\mathit{T}},{\mathit{u}}_{\mathit{N}},{\mathit{u}}_{\mathit{T}}$ are the normal and tangential components of the tension tensor and displacement vector, respectively.

We considered the constitutive relation Kelvin–Voigt in the viscoelastic case.

${\sigma }_{ij}\equiv {{\sigma }_{ij}\left(\mathit{u},\dot{\mathit{u}}\right)=C}_{ijkl}^{\left(1\right)}{\epsilon }_{kl}\left(\mathit{u}\right)+C_{ijkl}^{\left(2\right)}{\epsilon }_{kl}\left(\dot{\mathit{u}}\right)$

where $C_{ijkl}^{\left(1\right)}$ și $C_{ijkl}^{\left(2\right)}$ with $i,j,k,l=1,\dots ,d$, are the components of the elastic and viscoelastic tensors, respectively, which satisfy the ellipticity condition.

The classical quasistatic equilibrium equations are:

(1)$-\mathit{div}\mathit{\sigma }=\mathbf{f}, \mathrm{in} \mathrm{Q}$

(2)$\mathit{\sigma }\left(\mathit{u},\dot{\mathit{u}}\right)=C^{\left(1\right)}\mathit{\epsilon }\left(\mathit{u}\right)+C^{\left(2\right)}\mathit{\epsilon }\left(\dot{\mathit{u}}\right), \mathrm{in} \mathrm{Q}$

(3)$\mathit{u}=0 \mathrm{on} S_u$

(4)$\mathit{\sigma }·\mathit{n}=\mathit{h}\mathrm{on}{S}_{\begin{array}{c}N\\ \end{array}}$

(5)${\mathit{u}}_N\le \mathit{g},{\mathit{\sigma }}_N\le 0,{\sigma }_N\left({\mathit{u}}_N-\mathit{g}\right)=0\mathrm{on}{S}_C$

(6)$\parallel {\mathit{\sigma }}_{\mathit{T}}\parallel \le \mu \left|{\mathit{\sigma }}_N\right|\mathrm{and}\parallel {\mathit{\sigma }}_{\mathit{T}}\parallel <\mu \left|{\mathit{\sigma }}_N\right|\Rightarrow {\dot{\mathit{u}}}_T=0\mathrm{on}{S}_C$

(7)$\parallel {\mathit{\sigma }}_{\mathit{T}}\parallel =\mu \left|{\mathit{\sigma }}_N\right|\Rightarrow \exists \lambda \ge 0,\mathrm{s}.\mathrm{t}.{\dot{\mathit{u}}}_{\mathit{T}}=-\lambda {\mathbf{u}}_T\mathrm{on}{S}_C$

(8)$\mathit{u}\left(\mathit{x},0\right)={\mathit{u}}_{\mathbf{0}},\dot{\mathit{u}}\left(\mathit{x},0\right)={\mathit{u}}_1\mathrm{on}{S}_C$

where n is the external unit normal on $\partial \Omega ,\mu$ is the friction coefficient that depends on the tangential velocity, and $\dot{u}={\displaystyle \frac{\partial u}{\partial t}}$ is the deformation velocity.

Relation (1) is the equilibrium equation, (3) and (4) are the Dirichlet and Neumann boundary conditions, respectively, (5) is the Signorini contact condition, (6) and (7) are the Coulomb type friction laws, and (8) are the initial conditions.

Because of the Coulomb friction law, the problem becomes time dependent, and all the values can be time dependent. This makes the studied process irreversible and dependent on the loading way. Thus, the final displacements and tensions depend on the way the volume forces f and tractions h change in time, and not only their final values. The numerical algorithm proposed in this work can describe this by performing incremental loadings, in small steps, followed by iterative solutions and testing the convergence of the chosen method.

The quasistatic problem with dry friction following the Coulomb’s law for homogeneous and viscoelastic composite materials can be formulated as follows:

Problem P1: Find out the displacements and tensions that satisfy Equations (1)–(8) for any given moment t ∈ [0,T] and the initial condition u(x,0) = u₀ on $S_C$.

Physical arguments and mathematical difficulties regarding the local friction laws (6) and (7) made Duvant [12], Oden and Pires [13], and Cocu [14] take into consideration the nonlocal, even nonlinear, friction laws derived from (6) and (7) by replacing |${\sigma }_n$| with a convenient regularized R(${\sigma }_n$). These models fit even better, given that in the case of composite materials reinforced with fibers, the contact area has a micro-structure consisting of pronounced roughness, no longer being a perfect plain surface, and as the normal force N on the contact area gradually increases, the asperities gradually flatten to an equilibrium with the normal forces. However, the normal force is not single pointed but distributed over the real contact area, and thus the tangential tension ${\sigma }_T$ is also distributed over the real contact area, whose shape is influenced by surfaces with deformed asperities in contact.

In the accepted hypothesis, these asperities can be modeled as smooth periodic perturbations of the contact area with a period ε, where ε > 0; in this case, Coulomb’s law becomes:

$\parallel {\mathit{\sigma }}_T\left(\mathit{x}\right)\parallel \le \mu \mathit{N}\left(\mathit{y}\right)·\delta \left(\left|\mathbf{x}-\mathbf{y}\right|\right)$

where δ is Dirac’s distribution for a point x.

However, because the contact area is finite, we can assume that:

(9)$\parallel {\mathit{\sigma }}_{\mathit{T}}\left(\mathit{x}\right)\parallel \le \mu \mathit{N}\left(\mathit{y}\right){\omega }_{\delta 0}\left(\left|\mathbf{x}-\mathbf{y}\right|\right)$

where ${\left({\omega }_{\rho }\right)}_{\rho \ge 0}$ is a series of test functions that usually approximate $\delta$, and ${\rho }_o$ is the radius of the deforming contact area.

However, (9) only applies in the case of a normal force acting only in point y, while in our case, where N is distributed, we obtain:

(10)$\parallel {\mathit{\sigma }}_{\mathit{T}}\left(\mathit{x},t\right)\parallel \le \mu \left|\mathit{N}\left(\mathit{y}\right)\right|\ast {\omega }_{{\rho }_0}\left(\left|\mathbf{x}-\mathbf{y}\right|\right)$

where ∗ stands for the convolution operation.

Generalizing, a convenient regularization might be defined as:

(11)$\mathcal{R}\left({\mathit{\sigma }}_{\mathit{N}}\right)\left(\mathit{x}\right)={\int }_{{\Gamma }_C}{\omega }_{{\rho }_0}(\left|\mathbf{x}-\mathbf{y}\right|){\mathit{\sigma }}_N\left(\mathit{y}\right)\mathit{d}\mathit{y}$

where

${\omega }_{{\rho }_0}\left(r\right)=\left\{\begin{array}{c}{C}_{0}exp\left({\displaystyle \frac{{\rho }_0^2}{r^2-{\rho }_0^2}}\right)\\ 0ifr>{\rho }_0\end{array}\right.ifr\le {\rho }_0$

C₀ is constant and ${\rho }_0$ is a measure of the characteristic radius of the deformed asperities on the contact area.

Thus, for ${\delta }_0$ > 0, the slipping tangential tension is proportional to R(${\sigma }_N$), which is the weighted average of ${\sigma }_N$ in a vicinity of radius 2${\rho }_0$ on the contact area.

Oden and Pires [15] considered the problem of contact between an elastic body and a rigid support with a nonlocal and nonlinear friction law for the static case on ${\Gamma }_C$ and d = 3, ($\Omega \subset R^3$).

(12)${\mathit{\sigma }}_{T_i}\left(\mathit{u}\right)=-\mu \left|\mathcal{R}\left({\mathit{\sigma }}_N\right)\left(\mathit{u}\right)\right|{\phi }_{\epsilon }^{\prime }\left(\parallel {\mathit{u}}_{\mathit{T}}\parallel \right){\displaystyle \frac{{\mathit{u}}_{T_i}}{\parallel {\mathit{u}}_{T_i}\parallel }}$

For the quasistatic case, it will be

(13)${\mathit{\sigma }}_{T_i}\left(\mathit{u}\right)=-\mu \left|\mathcal{R}\left({\mathit{\sigma }}_{\mathit{N}}\right)\left(\mathit{u}\right)\right|{\phi }_{\epsilon }^{\prime }\left(\parallel {\dot{\mathit{u}}}_{\mathit{T}}\parallel \right){\displaystyle \frac{{\dot{\mathit{u}}}_{T_i}}{\parallel {\dot{\mathit{u}}}_{T_i}\parallel }}$

where $i=\mathrm{1,2}$, ${\phi }_{\epsilon }:\left(0,\infty \right)\to R$ is an increasing function, piecewise differentiable so that $\underset{\epsilon \to 0}{\mathit{lim}}{\phi }_{\epsilon }^{\prime }\left(r\right)=1$ and $\underset{r\to \infty }{\mathit{lim}}{\phi }_{\epsilon }^{\prime }\left(r\right)=1$, and ${\displaystyle \frac{1}{ℇ}}$ is a measurement of the stiffness of the asperities’ junctions on the contact area of the composite material. If there is no reversal of the loading direction, the following form can be chosen for the function ${\phi }_{\epsilon }$:

${\phi }_{\epsilon }\left(r\right)=\left\{\begin{array}{c}r-{\displaystyle \frac{\epsilon }{2}}ifr>\epsilon \\ {\displaystyle \frac{r^2}{2\epsilon }}ifr\le \epsilon .\end{array}\right.$

The friction laws (12) and (13) are general laws containing particular cases (i.e., for the static and quasi-static cases, respectively):

(a). if $\delta \to 0$ and $\epsilon >0$, we obtain local but nonlinear laws;

(b). if $\delta >0$ and $\epsilon \to 0$, we obtain nonlocal laws;

(c). if $\delta \to 0$ and $\epsilon \to 0$, we obtain classical Coulomb’s laws;

(d). if $\delta \to 0$ and $\epsilon \to 0$, if ${\mathit{\sigma }}_{\mathit{N}}\left(\mathit{u}\right)$ is known, the contact area is known, a priori;

(e). if $\mu =0$ we obtain the classical problem of frictionless unilateral contact.

The unilateral Signorini-type contact conditions (5) represent impenetrability conditions, valid for smooth contact surfaces, and remain smooth throughout the process. These can be represented by functions in C¹(Ω) as continuous with continuous derivatives.

Both the theoretical and experimental data led to the conclusion that a power law should be adopted to model the normal contact tension variation ${\sigma }_N$ in accordance with the penetration degree for sufficiently small loadings.

These aspects led Oden and Martins [16] to consider the following power law relationship:

(14)${\mathit{\sigma }}_{\mathit{N}}=c_N({\mathit{u}}_{\mathit{N}}-\mathit{g}{)}_{+}^{m_N}on{S}_C$

In the case of a deformable body in contact with a rigid support, $c_N$ > 0 and $m_N$ > 0 are experimentally determined constants representing the contact physical characteristics, the notation ${(·)}_{+}$ represents the positive part (i.e., ${\left(x\right)}_{+}$ = max (x, 0)), and g is the distance between the rigid surface and the highest asperities of the body.

Note that the Signorini law (5) can be formally obtained from (14), making $c_N\to \infty$.

The contact with friction problem for a viscoelastic body in contact with a rigid support can be therefore formulated as:

Problem P2: Find out the displacements and tensions that satisfy Equations (1)–(8), (13) and (14) and the initial condition $\mathit{u}\left(\mathit{x},0\right)={\mathit{u}}_0\in S_C.$

The proposed friction law is a generalization of the friction law considered so far to be isotropic.

For the case d = 3 (i.e., $\Omega \subset R^3)$, we denote by ${\nu }_1$, ${\nu }_2$ the main orthotropic axes defined on the surface portion ${\Gamma }_C$, with $\parallel {\nu }_1\parallel =\parallel {\nu }_2\parallel$ and ${\nu }_1·{\nu }_2=0$, then instead of Relations (7) and (8), we may consider, without essential changes, the following orthotropic friction law, suitable for composite materials reinforced with natural fibers oriented in two perpendicular directions:

$\parallel {\mathit{\sigma }}_{\mathit{T}}·{\mathit{\nu }}_{\mathit{i}}\parallel \le {\mu }_i\left|{\mathit{\sigma }}_{\mathit{N}}\right| \mathrm{and}$

$\parallel {\mathit{\sigma }}_{\mathit{T}}·{\mathit{\nu }}_{\mathit{i}}\parallel <{\mu }_i\left|{\mathit{\sigma }}_{\mathit{N}}\right|\Rightarrow {\dot{\mathit{u}}}_{\mathit{T}}·{\mathit{\nu }}_{\mathit{i}}=0on{\Gamma }_C,$

$\parallel {\mathit{\sigma }}_{\mathit{T}}·{\mathit{\nu }}_{\mathit{i}}\parallel ={\mu }_i\left|{\mathit{\sigma }}_{\mathit{N}}\right|\Rightarrow \exists {\lambda }_i\ge 0s.t.{\dot{\mathit{u}}}_{\mathit{T}}·{\mathit{\nu }}_{\mathit{i}}=-\lambda {\mathit{\sigma }}_{\mathit{T}}{\mathit{\nu }}_i,i=\mathrm{1,2}.$

2.2. Variational Formulation of Quasistatic Contact Problems for Composite Materials

In order to numerically solve the quasistatic contact problem using the finite element method, we need to express it variationally. We start by introducing the following function spaces:

(15)$V=\left\{\mathit{v}={(v}_1,\dots ,v_d\right)\in \left[H^1{\left(Q\right)]}^d∣\mathit{v}=0pe{S}_u\right\}$

(16)$K=\left\{\mathit{v}\in V∣{\mathit{v}}_N\le \mathit{g} a.e. on S_u\right\}$

with $\mathit{g}\in L^2\left(S_C\right) \mathrm{a}\mathrm{n}\mathrm{d} \mathit{g}\ge 0 a.e. on S_C$, we define

(17)$a\left(\mathit{u},\dot{\mathit{u}},\mathit{v}\right)={\int }_{\Omega }{\mathit{C}}_{ijkl}^{\left(1\right)}{\mathit{\epsilon }}_{kl}\left(\mathit{u}\right){\mathit{\epsilon }}_{kl}\left(\mathit{v}\right)d\mathit{x}+{\int }_{\Omega }{\mathit{C}}_{ijkl}^{\left(2\right)}{\mathit{\epsilon }}_{kl}\left(\dot{\mathit{u}}\right){\mathit{\epsilon }}_{kl}\left(\mathit{v}\right)d\mathit{x}$

(18)$j\left(\mathit{u},\mathit{v}\right)={\int }_{{\Gamma }_C}\mathsf{\mu }\mid \mathcal{R}\left({\mathit{\sigma }}_N\left(\mathit{u}\right)\right)\mid \parallel {\mathit{v}}_{\mathit{T}}\parallel ds$

(19)$j_T\left(\mathit{u},\mathit{v}\right)={\int }_{{\Gamma }_C}{C}_T{\left({\mathit{u}}_{\mathit{T}}-\mathit{g}\right)}_{+}^{m_T}{\parallel \mathit{v}}_T\parallel ds$

(20)$j_N\left(\mathit{u},\mathit{v}\right)={\int }_{{\Gamma }_C}{C}_N{\left({\mathit{u}}_{\mathit{N}}-\mathit{g}\right)}_{+}^{m_N}{\mathit{v}}_{N}ds$

(21)$F\left(\mathit{u},\mathit{v}\right)={\int }_{\Omega }\mathit{f}\mathit{v}d\mathit{x}+{\int }_{{\Gamma }_N}\mathit{h}\mathit{v}ds$

In order to obtain variational formulations for P₂, we used Green type relations, where for u, v ∈ V we obtained:

(22)${\int }_{\Omega }{\mathit{C}}_{ijkl}^{\left(1\right)}{\mathit{\epsilon }}_{kl}\left(\mathit{u}\right){\mathit{\epsilon }}_{kl}\left(\mathit{v}\right)d\mathit{x}+{\int }_{\Omega }{\mathit{C}}_{ijkl}^{\left(2\right)}{\mathit{\epsilon }}_{kl}\left(\dot{\mathit{u}}\right){\mathit{\epsilon }}_{kl}\left(\mathit{v}\right)dx={\int }_{{\Gamma }_N}(\mathit{u},\dot{\mathit{u})}·\mathit{n}\mathit{v}d\mathit{s}$

The Signorini condition (5) and the local friction laws (6) and (7) will accordingly become:

(23)${\int }_{{\Gamma }_C}\left[{\mathit{\sigma }}_T\left(\mathit{u}\right)\left({\mathit{v}}_{\mathit{T}}-{\dot{\mathit{u}}}_{\mathit{T}}\right)+\mu ∣{\mathit{\sigma }}_N\left(\mathit{u}\right)∣\left(\parallel {\mathit{v}}_T\parallel -\parallel {\dot{\mathit{u}}}_T\parallel \right)\right]\cdot ds\ge 0,\forall \mathit{v}\in K$

(24)${\int }_{{\Gamma }_C}[{\mathit{\sigma }}_T\left(\mathit{u}\right)({\mathit{v}}_{\mathit{T}}-{\dot{\mathit{u}}}_{\mathit{T}})+\mu \mid \mathcal{R}{\mathit{\sigma }}_N\left(\mathbf{u}\right)\mid (\parallel {\mathit{v}}_T\parallel -\parallel {\dot{\mathit{u}}}_T\parallel )]\cdot ds\ge 0$

By using Green’s formula on Ω, we obtain the following variational expression of problem P₁:

Problem P^’₁: Find out u: Q → K so that $\mathit{u}\left(x,0\right)={\mathit{u}}_0$ and

(25)$a\left(\mathit{u},\dot{\mathit{u}},\mathit{v}-\dot{\mathit{u}}\right)+j\left(\mathit{u},\mathit{v}\right)-j\left(\mathit{u},\dot{\mathit{u}}\right)\ge F\left(\mathit{v}-\dot{\mathit{u}}\right),\forall \mathit{v}\in K.$

Analogous, a variational expression for problem P₂ is:

Problem P^‘₂: Find out u: Q →V so that $\mathit{u}\left(\mathit{x},0\right)={\mathit{u}}_0$ and

(26)$$\begin{gathered} a\left(\mathit{u},\dot{\mathit{u}},\mathit{v}-\dot{\mathit{u}}\right)+j_N\left(\mathit{u},\mathit{v}-\dot{\mathit{u}}\right)+j_T\left(\mathit{u},\mathit{v}\right)-j_T\left(\mathit{u},\dot{\mathit{u}}\right)\ge F\left(\mathit{v}-\dot{\mathit{u}}\right), \\ \forall \mathit{v}\in V \end{gathered}$$

Please note that the inequalities (23) and (24) can be written in an elegant synthetic mathematical form using the notion of sub-differentials.

In order to solve these problems numerically, we need to take into account approximations of these problems, obtained by discretization using the finite element method in space and finite differences in time. This procedure leads to incremental formulations, to whom we may apply, in principle, the same methods as for static problem cases.

2.3. Incremental Formulations Obtained by Temporal Discretization

We obtained the incremental formulations from the variational formulations by replacing the temporal derivatives of displacements with finite differences. The quasistatic problems were solved step by step: at each step, we calculated the small deformations and displacements, obtained as a result of the small incremental changes of the applied forces, and we added them to the values calculated in the previous step.

Thus, the dependence on the way of loading is neglected in every increment; this sequence of steps also takes into account how the applied forces vary in time, which is very important for describing the contact state of the points on the contact boundary.

In order to obtain incremental formulations of the quasistatic contact problems, we introduced a partition $\left(t^0,t^1,\dots ,t^p\right)$ of the time interval $\left[0,T\right]$ and the following notations.

$\mathit{u}\left(t^m\right)={\mathit{u}}^m,∆{\mathit{u}}^m={\mathit{u}}^{m+1}-{\mathit{u}}^m,∆t^m=t^{m+1}-t^m,{\mathit{F}}^m=\mathit{F}{(t}^m),$

${∆\mathit{F}}^m={\mathit{F}}^{m+1}-{\mathit{F}}^m,m=\mathrm{0,1},\dots ,p-1$

Derivative, by the time t, $\dot{\mathit{u}}\left(\mathit{x},t\right)={\displaystyle \frac{\partial \mathit{u}\left(\mathit{x},t\right)}{\partial t}}$ can be approximated using the finite differences formula:

$\dot{\mathit{u}}\left(t^m\right)\cong {\displaystyle \frac{\mathit{u}\left(t^{m+1}\right)-\mathit{u}\left(t^m\right)}{{t^{m+1}-t}^m}}={\displaystyle \frac{{∆\mathit{u}}^m}{{∆t}^m}},m=\mathrm{0,1},\dots ,p-1$

Using the above notations, the time discretized form of the problem P^‘₁ becomes:

(27)$a\left({\mathit{u}}^m+{∆\mathit{u}}^m,{\displaystyle \frac{{∆\mathit{u}}^m}{{∆t}^m}}\mathit{v}-{\displaystyle \frac{{∆\mathit{u}}^m}{{∆t}^m}}\right)+j\left({\mathit{u}}^m+{∆\mathit{u}}^m,\mathit{v}\right)-j\left({\mathit{u}}^m+{∆\mathit{u}}^m,{\displaystyle \frac{{∆u}^m}{{∆t}^m}}\right)\ge \left({\mathit{F}}^m+{∆\mathit{F}}^m\right)\left(\mathit{v}-{\displaystyle \frac{{∆\mathit{u}}^m}{{∆t}^m}}\right)$

Denoting by $K^m$ the set:

$K^m=\left\{{\mathit{w}}^{\mathit{m}}\in V∣{\mathit{w}}_{\mathit{N}}^{\mathit{m}}\le \mathit{g}-{\mathit{u}}_{\mathit{N}}^{\mathit{m}},a.e.on{S}_C\right\}$, we see that:

${∆\mathit{u}}^m\in K^m⟺{\mathit{u}}^{m+1}={\mathit{u}}^m+{∆\mathit{u}}^m$, and from (27), we obtain:

Problem P^’’₁: We looked for ${∆\mathit{u}}^m\in K$ so that

(28)$$\begin{gathered} a\left({∆\mathit{u}}^m,{∆\mathit{u}}^m,{\mathit{w}-∆\mathit{u}}^m\right)+j\left({\mathit{u}}^m+{∆\mathit{u}}^m,\mathit{w}\right)-j\left({\mathit{u}}^m+{∆\mathit{u}}^m,{∆\mathit{u}}^m\right)\ge \\ ∆F^m\left(\mathit{w}-{∆\mathit{u}}^m\right)-G\left({\mathit{u}}^m,\mathit{w}-{∆\mathit{u}}^m\right),\forall \mathit{w}\in K^m, \end{gathered}$$

where we denoted:

$G\left({\mathit{u}}^m,\mathit{w}-{∆\mathit{u}}^m\right)=a\left({\mathit{u}}^m,{∆\mathit{u}}^m,{\mathit{w}-∆\mathit{u}}^m\right)-F^m\left(\mathit{w}-{∆\mathit{u}}^m\right),$

which, using Green’s formula and considering the equilibrium in step m due to the displacement ${\mathit{u}}^m$, can also be written as:

$G\left({\mathit{u}}^m,\mathit{w}-{∆\mathit{u}}^m\right)={\int }_{{\mathsf{\Gamma }}_C}\left[{\mathit{\sigma }}_N\left({\mathit{u}}^m\right)\left({\mathit{w}}_N-∆{\mathit{u}}_N^m\right)+{\mathit{\sigma }}_T\left({\mathit{u}}^m\right)\left({\mathit{w}}_T-∆{\mathit{u}}_T^m\right)\right]ds.$

One can note that the incremental inequality (28) has a similar form with a quasi-variational inequality describing a static problem, and the condition ${∆\mathit{u}}^m\in K^m$ ensuring the impenetrability condition ${\mathit{u}}^{m+1}\in K^m$, m = 0, 1, …, p − 1.

Similarly, problem P^′₂ leads to the incremental one.

Problem P^’’₂: Find out the sequence ${\mathit{u}}^1,\dots ,{\mathit{u}}^p$ with ${\mathit{u}}^{m+1}={\mathit{u}}^m+{∆\mathit{u}}^m\in V,\left({\mathit{u}}^0={\mathit{u}}_0\right)$, where ${∆\mathit{u}}^m\in V$ is the inequality solution.

(29)$$\begin{gathered} a\left({∆\mathit{u}}^m,{∆\mathit{u}}^m,{\mathit{w}-∆\mathit{u}}^m\right)+j_N\left({\mathit{u}}^m+{∆\mathit{u}}^m,\mathit{w}-{∆\mathit{u}}^m\right)+j_T\left({\mathit{u}}^m+{∆\mathit{u}}^m,\mathit{w}\right)-j_T\left({\mathit{u}}^m+{∆\mathit{u}}^m,{∆\mathit{u}}^m\right) \\ \ge F^m\left(\mathit{w}-{∆\mathit{u}}^m\right)-G\left({\mathit{u}}^m,\mathit{w}-{∆\mathit{u}}^m\right)\forall \mathit{w}\in V. \end{gathered}$$

It can be seen that the quasi-variational inequality (29) is equivalent with the following incremental, time independent (static) problem.

It is required that ${∆\mathit{u}}^m\in V$, so that

${\sigma }_{ij,j}\left({∆\mathit{u}}^m\right)+f_i=0 \mathrm{in} \mathrm{Q}$

${\sigma }_{ij}\left({∆\mathit{u}}^m\right)n_j=∆h_i\mathrm{on}{S}_N$

$∆{\mathit{u}}^m=0\mathrm{on}{S}_u$

${\mathit{\sigma }}_N\left({\mathit{u}}^m+{∆\mathit{u}}^m\right)=-C_N({\mathit{u}}_N^m+∆{\mathit{u}}_N^m-g{)}_{+}^{m_N}on{S}_C$

$\parallel {\mathit{\sigma }}_T\left({\mathit{u}}^m+{∆\mathit{u}}^m\right)\parallel \le C_T({\mathit{u}}_N^m+∆{\mathit{u}}_N^m-\mathit{g}{)}_{+}^{m_N}and$

$\parallel {\sigma }_T\left({\mathit{u}}^m+{∆\mathit{u}}^m\right)\parallel <C_T({\mathit{u}}_N^m+∆{\mathit{u}}_N^m-\mathit{g}{)}_{+}^{m_T}\Rightarrow ∆{\mathit{u}}_T^m=0$

${\mathit{\sigma }}_T\left({\mathit{u}}^m+{∆\mathit{u}}^m\right)\parallel =C_T({\mathit{u}}_N^m+∆{\mathit{u}}_N^m-\mathit{g}{)}_{+}^{m_T}\Rightarrow \exists \lambda \ge 0s.t.$

$∆{\mathit{u}}_T^m=-\lambda {\mathit{\sigma }}_{\mathit{T}}\left({\mathit{u}}^m+{∆\mathit{u}}^m\right)$

The existence and uniqueness of the solutions of the quasi-variational inequations have been analyzed by Anderson [16], Klarbing [17], and Cocu [14]. The necessary conditions for this are as follows: the friction coefficient should be sufficiently small, slipping directions on the contact area should be predictable, and Hooke’s law should be dependent on the velocity of the deformation.

All of these restrictions are met in the case of quasistatic contact in viscoelastic materials such as natural fiber-reinforced composites.

To correctly model these characteristics, the incremental and iterative Newton-Raphson method was used.

2.4. Approximation Using the Finite Elements Method and the Description of the Incremental and Iterative Algorithms

2.4.1. The Existence and the Uniqueness of the Solution of the Incremental Problem: Contact Finite Element with Friction

Using the standard procedure for the finite element method, the approximation of the variational Equation (26) can be made in the finite dimensional space $V_h\subset V.$ For a certain (h > 0) that defines the fineness of the partition and the approximations of the displacements for the time $t\in \left[0,T\right]$ are elements, ${\mathit{u}}^h\left(x,t\right)\in V_h$. In each finite element ${\Omega }_h^e\left(e=1,\dots ,N_h\right)$ in $\Omega ,$ the displacement fields from within a finite element are interpolated between the nodal values of the element as:

(30)${\mathit{u}}_d\left(x,t\right)={\sum }_{i=1}^{N_e}{\mathit{u}}_d^i\left(t\right)N_d\left(\mathit{x}\right),$

where $d=2,3,$ $N_e$ is the node number of the finite element ${\Omega }_h^e$, and $N_d$ is the form function associated with the nodal point i of the finite element ${\Omega }_h^e$.

For every incremental step of Problem P^″₂, we took an internal approximation $\left(V_h,K_h,J_{Nh},J_{Th}\right)$ of the finite element type, with the discretization parameter $h>0$ [18] as:

Problem P^″_2h: Find out ${\mathit{u}}_h^m$ so that

(31)$a\left({\mathit{u}}_h^m,{\mathit{u}}_h^m,{{\mathit{w}}_h-\mathit{u}}_h^m\right)+J_{Nh}\left({\mathit{w}}_h^m,{\mathit{v}}_h-{\mathit{u}}_h^m\right)+J_{Th}\left({\mathit{u}}_h^m,{\mathit{v}}_{\mathit{h}}\right)-J_{Nh}\left({\mathit{u}}_h^m,{\mathit{u}}_h^m\right)\ge F_h^m\left({{\mathit{w}}_h-\mathit{u}}_h^m\right)-G\left({{{\mathit{u}}_h^m,\mathit{w}}_h-\mathit{u}}_h^m\right)$

Given that for every time increment $∆t$ the problem P₂ is a static one, the existence of the finite element approximate solution develops into proving the existence of the solution ${\mathit{u}}_h^m$ of inequation (31), which is a fixed point for the application $f_h:K_h\to K_h$, which associates $\forall {\mathit{z}}_h^m\in K_h$, where the element $f_h{\mathit{z}}_h^m\in K_h$ is defined as

(32)$a\left({\mathit{u}}_h^m,{\mathit{u}}_h^m,{{\mathit{w}}_h-\mathit{u}}_h^m\right)+J_{Nh}\left({\mathit{w}}_h^m,{\mathit{v}}_h-{\mathit{u}}_h^m\right)+J_{Th}\left({\mathit{u}}_h^m,{\mathit{v}}_{\mathit{h}}\right)-J_{Nh}\left({\mathit{u}}_h^m,{\mathit{u}}_h^m\right)\ge F_h^m\left({{\mathit{w}}_h-\mathit{u}}_h^m\right)-G\left({{{\mathit{u}}_h^m,\mathit{w}}_h-\mathit{u}}_h^m\right)$

so that the series $\left\{{\mathit{w}}_{hn}^m\right\}$ approaches ${\mathit{w}}_h^m$, where ${\mathit{w}}_{hn}^m={f_h\mathit{w}}_{h\left(n-1\right)}^m,n\ge 1,$ because the map $f_h$ is a contraction.

Proof. The proof can be found in [14] for a sufficiently small frictional coefficient.

The matrix form of the equilibrium equation approximated with the finite element method for Problem P^″_2h is:

Problem P^″_2hm: Find a function $\mathit{U}:\left[0,t\right]\to R^{dx{N}_t^{\Omega }}$ so that

(33)$KU\left(t\right)-P\left(U\left(t\right)\right)+J\left(U\left(t\right),\dot{U}\left(t\right)\right)=F\left(t\right),\mathrm{t}\in \left[0,T\right]$

using the approximation with the finite element method $\left(V_h,K_h,J_{Nh},J_{Th}\right)$, according to the notations (16)–(22). The use of the finite element method in the study of composites can be seen in [19,20,21,22,23,24]. Where K is the stiffness matrix, and $F\left(t\right)$ is the loadings vector, obtained from (18), and (22), respectively, by approximating and assembling with finite elements. The P and J matrices, obtained via (20) and (21), respectively, represent the nodal forces vector and the friction forces on the contact boundary ${\Gamma }_C$, respectively, obtained by assembling the contact finite elements.

System (33) is an algebraic system of equations, resulting from the approximation of the spatial variables with the finite element method and the approximation of the temporal variables with finite differences, resulting from the boundary conditions that also depend on the time t.

The contact finite element in 2D (see Figure 1) consists of three nodes on the contact boundary, two “master” nodes belonging to a body (or the rigid support) and a “slave” node belonging to the second body. This contact finite element is a node-segment type and has linear, first degree, interpolation functions and models the friction law and the non-penetration or penetration condition by a power type law [25,26].

2.4.2. The Newton–Raphson Algorithm in Iterative and Incremental Solving of Nonlinear Algebraic Systems

Let $K_N=K-P+J,$ the stiffness matrix containing also the contribution of the frictional contact; for this reason, the matrix $K_N$ will contain both nonlinearities and non-symmetries, caused by the slipping contact and the anisotropicity of the composite material. We considered a partition of the time interval [0,T] = $\bigcup _{k=1}^{Nt}[t_{k-1},t_k]$ with 0 = $t_0$ < ... < $t_{Nt}=T$ and denoted $∆t_{k+1}=t_{k+1}-t_K$, a time sub-interval where the system can be taken as static. We assumed that we found the solution ${\mathit{U}}_k$ on the sub-interval $∆t_k$, corresponding to this time step, and searched for a solution on the sub-interval $∆t_{k+1}$, applying iteratively the Newton–Raphson algorithm to obtain the solution ${\mathit{U}}_{k+1}^{n+1}$, which satisfies the condition ${|\mathit{U}}_{k+1}^{n+1}-{\mathit{U}}_k^n|$ < TOL, where index n defines the iteration on this sub-interval and TOL is a positive constant, sufficiently small, stop condition for the iteration. Considering (33), we denoted the matrix function as follows:

${T\left(\mathit{U}\left(t\right)\right)≔K}_N\mathit{U}\left(t\right)-F\left(t\right)=0$

or in the incremental and iterative form ${T(U}_k^n)≔K_N{U}_k^n-F_k^n$ = 0, n $\ge 1,$ and knowing it is differentiable in $U_k^n$ and $U_k^n$ is a solution on the sub-interval $∆t_k$ (i.e., $T\left(U_k^n\right)=0),$ we developed the Taylor series for function T and retained only the first term, thus resulting in the Newton method:

(34)${\mathit{U}}_{k+1}^{n+1}={\mathit{U}}_k^n+{\left[{\displaystyle \frac{\partial T\left({\mathit{U}}_k^n\right)}{\partial {\mathit{U}}_k^n}}\right]}^{-1}{T(\mathit{U}}_{k+1}^{n+1})$

Finally, on interval [0,T], the solution for the quasistatic problem will contain the solutions on every time sub-interval $∆t_k$, meaning that {U₀, U₁,..., U_T}, which correspond to the incremental loadings on each sub-interval.

A synthetic presentation of the algorithm of the quasistatic contact problem is given below:

Initialize data t₀ = 0, k = 0, n = 0, ${\mathit{U}}_k^n=0;$

Calculate the stiffness matrix K of the composite;

Calculate the contact nodal forces P and the frictional contact forces J, assemble the global stiffness matrix K_N and the global incremental loading vector $∆F_k$ on the sub-interval $∆t_k$, followed by resolving with the Newton–Raphson algorithm {${\mathit{U}}_{k+1}^n$};

Check the convergence: if $|{\mathit{U}}_{k+1}^{n+1}-{\mathit{U}}_k^n|$ < TOL, go to (vii);

Else go to (vi);

Update the displacement fields, the condition ${\mathit{U}}_{k+1}^{n+1}={\mathit{U}}_K^n$ + ${∆\mathit{U}}_k^n$, and the increments of the loading forces $F_{k+1}=F_k+∆F_k$;

n = n + 1 and go to (ii);

Exit.

2.4.3. Numerical Simulation

For the numerical simulation with finite elements, we chose a classical test used by many authors [27,28]. The test consists of a plane composite plate made of 50% polypropylene and 50% hemp fibers, ±10%, whose elasticity module E_c was determined through measurements in [29] using the resonance method. The plate is in dry friction contact with a rigid support, the contact area being along the length of the plate. This test has the advantage of being very simple geometrically and symmetric by the vertical axis in the middle of the rectangular plate. The stress loadings manifest in the plate plane (Figure 2), so that on the contact area with the rigid support, all of the possible contact states can appear: open, sliding, or fixed. For symmetry reasons, working with only half of the plate will suffice. Table 1 presents the zone distributions on the contact boundary, with the contact states open, sliding, or fixed, for different fiber to composite ratios, different stress loadings, and different friction coefficients, in the hypothesis of plane state deformation. Half of the plane plate was discretized into 128 quadrilateral finite elements (plane state deformation), and the contact area ABCD was discretized into 32 node-segment contact finite elements. The contact area ABCD had 17 nodes.

The elasticity module (E_c) formula for a composite material reinforced with random oriented fibers derives from composite theory as it is an extension of the law of the mixture, adjusted for the random orientation of the fibers. The most frequent form of this is:

E_c = η V_f E_f + V_m E_m

where:

-. V_f, V_m are the volume ratios of the fibers and matrix (V_f + V_m = 1), respectively;

-. E_f, E_m are the elasticity modules of the fibers and matrix, respectively,

-. η is an orientation factor that depends on the fiber length, orientation, adhesion, etc. (for random orientation in 3D, η ≈ 1/5; for random 2D, η ≈ 3/8).

E_c = (3/8)V_f E_f + V_m E_m, in the case of 2D.

This formula is accepted in the literature [28,30,31] for composites reinforced with random oriented short fibers, mainly in studies on thermoplastic materials reinforced with natural or synthetic fibers. An extrapolation of the mechanical properties of the composite with its proportions was searched, but for this, the elasticity module E of every constituent is needed, (random oriented fibers + melted matrix), and not only the E of the composite at a 50–50% mixture.

However, it is known that the law of mixture is only approximately valid and does not take into account the shape or appearance of the fibers, their distribution or orientation, and most importantly, the fiber–matrix adhesion. To verify these aspects, it was necessary to find out the Young’s modulus E for hemp fiber and for molten polypropylene, in our case, from the specialized literature [32].

From the numerical simulation with the finite element method, we obtained the following results, which characterized with a good approximation, the mechanical properties of the composite material.

In Table 1, the distribution of contact states for E_C = 130 GPa (see [32]) on the ABCD boundary is represented.

The ordinates of the deformation of the contact nodes, for three different values of E_c, is presented in Figure 3.

The main goal of this numerical simulation was to emphasize the importance of the ratios of the two constituents of the composite by analyzing the behavior of the contact states on the contact boundary $\overline{ABCD}$ in the case of the plane state deformation in Figure 2. The variation in the volume ratio of the fibers had an important influence on the distribution of the contact areas and states and on the normal and tangential contact loadings (see the results in Table 1 and Table 2).

The contact states for a composite with different values of the friction coefficient are presented in Table 3.

It was found that the composite properties were sufficiently well-characterized by the boundary area in contact with the rigid solid through harshness, friction coefficient, and the distribution of the contact areas states, which can be open, sliding, or fixed, and through the normal and tangential stress loadings in the contact area.

In practical problems, nonlinear problems often arise, as in our case, in the study of the elastic contact problem with friction for composites reinforced with randomly oriented short natural fibers. The implemented algorithms are restrictive and users have to continuously modify them to solve certain problems. For this reason, it is desirable to use modular software, whose modules can be selectively used to conveniently solve the desired problem with only the help of input data, without modifying the software. This was achieved and presented in this manuscript with the help of an interpretive language using macro-instructions that allowed us to build a variable algorithm, depending on the type of problem to be solved (see [33]).

For example, if in the numerical simulation the composite material is considered homogeneous and isotropic, the global stiffness matrix will be symmetric, and the TANG macro-instruction will be used; in the case of an anisotropic composite or if sliding contact also occurs on the contact area, the stiffness matrix will be non-symmetric, in which case the UTAN macro-instruction will be used.

If we want to form the right member of the equations, modified in the quasistatic case with time dependent boundary conditions, we will use the FORM macro-instruction; the macro-instruction SOLV will give the system solution U, the printing of the displacements is made with DISP, and of the stresses with STRE.

LOOP n and NEXT can be used to repeat n times the group of instructions between LOOP n and NEXT.

The macro-instructions set in the first column of Table 4 implements the Newton–Raphson linearization method, which at each step recalculates both the stiffness matrix and the free terms vector, thus requiring a high computational effort. This algorithm consists of 10 iterations and the displacements are printed every two iterations.

If the nonlinearities are not very pronounced, the modified Newton–Raphson method can be used, which forms the stiffness matrix only once and modifies the vector of free terms at each iteration; an example of this case is in the second column of Table 4.

This method is faster at computing time, but the convergence speed to the solution is lower than in the case of the classical method. For this reason, there are times when a combination of the two methods is a better approach: for example, case C in Table 4, when the stiffness matrix is formed only once every five iterations and the free term is updated in each iteration.

In the case where we want to solve an incremental and iterative nonlinear problem in which the load force is given incrementally, case D in Table 4 should be used, where:

DT = initializes iterations;

PROP = performs the proportional incremental loading;

CV = performs the convergence test of the displacements.

3. Conclusions

One of the main results of this manuscript is that the properties of the composite materials were sufficiently well-characterized by the boundary area in contact with the rigid solid through the harshness, friction coefficient, and the distribution of the contact areas states, which could be open, sliding, or fixed, and through the normal and tangential stress loadings in the contact area. All of these are considered in the contact with friction phenomenon.

In addition, it was observed that composite materials exhibit a certain shape memory, pronounced damping, and a delayed return to the initial state.

The main contribution of the algorithm presented in this manuscript and proven by examples from the numerical simulations is that the mechanical behavior of the new composite obtained with different fiber volume proportions can be predicted with a good approximation.

Author Contributions

Conceptualization, M.R.A. and N.P.; Methodology, A.M.M. and T.S.; Software, N.P. and V.M.M.; Validation, M.R.A., A.M. and N.P.; Formal analysis T.S.; Data curation, T.S. and A.M.M.; Writing—original draft preparation, N.P. and M.R.A.; Writing—review and editing, A.M.M. and A.M.; Supervision, N.P. and T.S. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The original contributions presented in this study are included in the article. Further inquiries can be directed to the corresponding authors.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1 Contact finite element in 2D. 1 and 2 are the master nodes, S is the slave node and g is the gap.

Figure 2 Geometry (h = 40 mm) and stress loading. AB segment (boundary) is open contact, BC segment in a sliding contact, CD is stick contact. F stands for the horizontal force applied on the AG segment, and f for the vertical force applied on the GE segment of the plate.

Figure 3 The three deformations of the contact boundary for η = 0.2, μ = 0.01, V_f = 0.3 with E_c = 27 MPa, V_f = 0.5 with E_c = 3650 MPa, and V_f = 0.7 and E_c = 4590 MPa [31].

Table 1

The contact states on the ABCD boundary.

$\mathit{\mu }$	F[N]	f[N]	Gap Status AB[mm]	Slip Status BC[mm]	Stick Status CD[mm]
0.2	100,000	−50,000	0.	2.5	37.5
				(2 nodes)	(15 nodes)
0.2	100,000	−150,000	0.	10	30
				(5 nodes)	(12 nodes)
0.2	100,000	−250,000	0.	35	5
				(14)	(3 nodes)

Table 2

Contact state for the composite material for Ec =3.65 GPa with variable η $\in$ {0.25, 0.2, 0.15} and V_f = 0.5, ν = 0.39 with 17 contact nodes, where E_f = 30 GPa and E_m =1.3 GPa [31].

$\mathit{\mu }$	F[N]	f[N]	Gap Status AB[mm]	Slip Status BC[mm]	Stick Status CD[mm]
0.2	100,000	−50,000	0.	15	25
0.2	100,000	−150,000	0.	32.5	7.5
0.2	100,000	−250,000	0.	37.5	2.5

Table 3

Contact state for the composite material with E_c = 130 GPa, ν = 0.3, V_f = 0.5, $\eta =0.2$ with 17 contact nodes, for different values of $\mu$ [31,32].

$\mathit{\mu }$	F[N]	F[N]	Gap Status AB[mm]	Slip Status BC[mm]	Stick Status CD[mm]
0.0	200,000	−100,000	0	40	0
0.2	200,000	−100,000	0	35	5

Table 4

Software usage examples.

CaseA	CaseB	CaseC	CaseD	Meaning
LOOP 10	TANG (UTAN)	LOOP 2	DT 1	$∆t=1$
TANG (UTAN)	LOOP 10	TANG (UTAN)	PROP 1	$F=t{F}_0$
FORM	FORM	LOOP 5	LOOP 10
SOLV	SOLV	FORM	TIME	$t=t+∆t$
DISP 2	DISP 2	SOLV	TANG (UTAN)	Calculate the K matrix
NEXT	NEXT	DISP	FORM	Calculate the F updated free term
STRE	STRE	NEXT	SOLV	$\mathrm{Calculate}U$

References

1. Rodriguez-Tembleque, L.; Buroni, F.; Abascal, R.; Sáez, A. Analysis of FRP composites under frictional contact condition. Int. J. Solids Struct.; 2013; 50, pp. 3947-3959. [DOI: https://dx.doi.org/10.1016/j.ijsolstr.2013.08.007]

2. Alzweighi, M.; Tryding, J.; Mansour, R.; Borgqvist, E.; Kulachenko, A. Anisotropic damage behavior in fiber-based materials: Modeling and experimental validation. J. Mech. Phys. Solids; 2023; 181, 105430. [DOI: https://dx.doi.org/10.1016/j.jmps.2023.105430]

3. Pickering, K.L.; Efendy, M.G.A.; Le, T.M. A review of recent developments in natural fibre composites and their mechanical performance. Compos. Part A Appl. Sci.; 2016; 83, pp. 98-112. [DOI: https://dx.doi.org/10.1016/j.compositesa.2015.08.038]

4. Abdollahiparsa, H.; Shahmirzaloo, A.; Teuffel, P.; Blok, R. A review of recent developments in structural applications of natural fiber-Reinforced composites (NFRCs). Compos. Adv. Mater.; 2023; 32, pp. 1-18. [DOI: https://dx.doi.org/10.1177/26349833221147540]

5. Vlassak, J.J.; Ciavarella, M.; Barber, J.R.; Wang, X. The indentation modulus of elastically anisotropic materials for indenters of arbitrary shape. J. Mech. Phys. Solids; 2003; 51, pp. 1701-1721. [DOI: https://dx.doi.org/10.1016/S0022-5096(03)00066-8]

6. Ning, X.; Lovell, M.R.; Slaughter, W.S. Asymptotic solutions for axisymmetric contact of a thin, transversely isotropic elastic layer. Wear; 2006; 260, pp. 693-698. [DOI: https://dx.doi.org/10.1016/j.wear.2005.03.024]

7. Bagault, C.; Nélias, D.; Baietto, M.C.; Ovaert, T.C. Contact analyses for anisotropic half-space coated with an anisotropic layer: Effect of the anisotropy on the pressure distribution and contact area. Int. J. Solids Struct.; 2013; 50, pp. 743-754. [DOI: https://dx.doi.org/10.1016/j.ijsolstr.2012.11.002]

8. Birleanu, C.; Udroiu, R.; Cioaza, M.; Pustan, M.; Paul, B.; Vilau, C. The Effect of fiber Weight Fraction on Tribological Behavior for Glass Fiber Reinforced Polymer. Polymers; 2025; 17, 720. [DOI: https://dx.doi.org/10.3390/polym17060720]

9. Tserpes, K.; Sioutis, I. Advances in Composite Materials for Space Applications: A Comprehensive Literature Review. Aerospace; 2025; 12, 215. [DOI: https://dx.doi.org/10.3390/aerospace12030215]

10. Goda, T.; Vàradi, K.; Wetzel, B.; Friedrich, K. Finite element simulation of the fiber–matrix debonding in polymer composites produced by a sliding indentor: Part II: Parallel and anti-parallel fiber orientation. J. Compos. Mater.; 2004; 38, pp. 1607-1618. [DOI: https://dx.doi.org/10.1177/0021998304043760]

11. Vàradi, K.; Nèder, Z.; Friedrich, K.; Flöck, J. Finite-element analysis of a polymer composite subjected to a ball indentation. Compos. Sci. Technol.; 1999; 59, pp. 271-281.

12. Duvaut, G. Equilibre d’un solide elastique avec contact unilateral et frottement de Coulomb. C. R. Acad. Sc. Paris; 1980; 290, pp. 263-265.

13. Oden, J.T.; Pires, E. Contact Problems in Elastostatics with Non-Local Friction Laws; TICOM Report The University of Texas at Austin: Austin, TX, USA, 1981; pp. 81-112.

14. Cocu, M. Existence of solutions of Signorini problems with friction. Int. J. Eng. Sci.; 1984; 22, pp. 567-575. [DOI: https://dx.doi.org/10.1016/0020-7225(84)90058-2]

15. Oden, J.T.; Pires, E. Numerical analysis of certain contact problems in elasticity with nonlocal friction laws. Comp. Struct.; 1983; 16, pp. 481-485. [DOI: https://dx.doi.org/10.1016/0045-7949(83)90187-6]

16. Andersson, L.E. A global existence result for a quasistatic contact problem with friction. Adv. Math. Sci. Appl.; 1995; 5, pp. 249-286.

17. Klarbring, A.; Mikelic, A.; Schillor, A. A global existence result for the quasistatic frictional contact problem with normal compliance. Unilateral Problems in Structural Analysis IV (Capri 1989); Birkhauser: Basel, Switzerland, 1991; pp. 85-111.

18. Pop, N. Preconditioning Uzawa algorithm for contact problems. PAMM Proc. Appl. Math. Mech.; 2008; 8, pp. 10985-10986. [DOI: https://dx.doi.org/10.1002/pamm.200810985]

19. Bonari, J.; Paggi, M.; Dini, D. A new finite element paradigm to solve contact problems with roughness. Int. J. Solids Struct.; 2022; 253, 111643. [DOI: https://dx.doi.org/10.1016/j.ijsolstr.2022.111643]

20. Jiang, L.; Wu, M.; Yu, Q.; Shan, Y.; Zhang, Y. Investigations on the Adhesive Contact Behaviors between a Viscoelastic Stamp and a Transferred Element in Microtransfer Printing. Coatings; 2021; 11, 1201. [DOI: https://dx.doi.org/10.3390/coatings11101201]

21. Bensaada, A.; Essoufi, E.-H.; Zafrar, A. Primal-dual formulation for parameter estimation in elastic contact problem with friction. Appl. Math. Sci. Eng.; 2024; 32, 2367025. [DOI: https://dx.doi.org/10.1080/27690911.2024.2367025]

22. Katouzian, M.; Vlase, S.; Marin, M.; Pop, N. Temperature Influence on the Creep Behavior of a Carbon Fiber/Polyetheretherketone Composite. J. Appl. Comput. Mech.; 2025; pp. 1-10. [DOI: https://dx.doi.org/10.22055/jacm.2025.48137.5007]

23. Katouzian, M.; Vlase, S. A Model to Study the Creep Behavior of Carbon Fiber/Epoxy Resin Composites Under Temperature. Appl. Sci.; 2025; 15, 4206. [DOI: https://dx.doi.org/10.3390/app15084206]

24. Vasile, A.; Coropeţchi, I.C.; Constantinescu, D.M.; Sorohan, S.; Picu, C.R. Simulated annealing algorithms used for microstructural design of composites. Mater. Today Proc.; 2023; 93, pp. 680-684. [DOI: https://dx.doi.org/10.1016/j.matpr.2023.05.109]

25. Wriggers, P.; Simo, J.C. A note on tangent stiffness for fully nonlinear contact problems. Commun. Appl. Numer. Methods; 1985; 1, pp. 199-203. [DOI: https://dx.doi.org/10.1002/cnm.1630010503]

26. Simo, J.C.; Wriggers, P.; Taylor, R.L. A perturbed Lagrangian formulation for the finite element solution of contact problems. Comput. Methods Appl. Mech. Eng.; 1985; 50, pp. 163-180. [DOI: https://dx.doi.org/10.1016/0045-7825(85)90088-X]

27. Raous, M.; Chabrand, P.; Lebon, F. Numerical method for frictional contact problems and applications. J. Theor. Appl. Mech.; 1988; 7, pp. 111-128.

28. Hull, D.; Clyne, T.W. On Introduction to Composite Materials; 2nd ed. Cambridge University Press: Cambridge, UK, 1996.

29. Apsan, M.R.; Mitu, A.M.; Neagoe, C.N.; Pop, N.; Sireteanu, T. Non-Destructive Testing for Evaluation of Young’s Modulus by Using Free Vibration Response of Composite Materials. Appl. Sci.; 2025; 15, 189. [DOI: https://dx.doi.org/10.3390/app15010189]

30. Halpin, J.C.; Kardos, J.L. The Halpin-Tsai Equations: A Review. Polym. Eng. Sci.; 1976; 16, pp. 344-352.

31. Mina, J.H.; González, A.V.; Muñoz-Vélez, M.F. Micro- and Macromechanical Properties of a Composite with a Ternary PLA–PCL–TPS Matrix Reinforced with Short Fique Fibers. Polymers; 2020; 12, 58. [DOI: https://dx.doi.org/10.3390/polym12010058]

32. Alart, P.; Lebon, F. Numerical Study of a Stratified Composite Coupling Homogenization and Frictional Contact. Math. Comput. Model.; 1988; 28, pp. 273-286. [DOI: https://dx.doi.org/10.1016/S0895-7177(98)00122-8]

33. Zienkiewicz, O.C.; Taylor, R.L. The Finite Element Method; 4th ed. McGraw-Hill: London, UK, 1991; Volume 2.