Nonlinear analysis of plane frames considering hyperelastic models through the finite element positional method

Abstract

Computational mechanics has become an essential tool in engineering, just as the use of hyperelastic materials has seen remarkable growth in everyday applications. Therefore, it is fundamental to study hyperelastic models that represent the behavior of these materials, such as elastomers and polymers. With that in mind, the Mooney-Rivlin, Neo-Hookean, Ogden, and Yeoh models were implemented in a computational code in FORTRAN using the Positional Finite Element Method with Reissner kinematics and the Newton-Raphson method for nonlinear analysis of plane frames with samples of elastomers added with different percentages of carbon black. Ultimately, it was concluded that the Yeoh and Ogden models presented coherent values and that the use of the formulation for nonlinear analysis of plane frame performs well after the modifications proposed by this work. These modifications consisted of adding the first and second strain invariants of the simple shear formulation to include the consideration of distortion in the specific strain energy of hyperelastic models.

Keywords:
Computational mechanics; Hyperelastic models; Positional finite element method; Plane structures; Plane frames

Graphical Abstract

1 INTRODUCTION

The properties of materials and structural behavior are critical factors for a successful structural design, and understanding them is essential for designing safe and efficient structures. Computational mechanics, through numerical simulations to model and analyze complex systems, has become an essential tool for engineers in this regard (Wang et al., 2019; Navas et al., 2020).

By using computational mechanics, engineers can simulate how different materials will behave under various conditions, enabling them to predict the behavior of a structure before its construction. This can help identify potential problems early in the design process and lead to more efficient and cost-effective designs. However, while simulations can provide valuable information, there is always some degree of uncertainty involved, and it is essential to verify the accuracy of simulations through experimental testing (Samaniego et al., 2020).

The Finite Element Method (FEM) is a widely used numerical technique in structural analysis, mechanical engineering, and other fields. FEM is based on the principle of minimum potential energy, which states that the deformation of a structure under external loads will tend to minimize the potential energy of the system (David Müzel et al., 2020).

In positional formulations of the FEM, the nodal positions of a structure are taken as fundamental unknowns, which are used to derive the equations of motion and solve for stresses, strains, and other physical quantities of interest. This approach allows for flexible and efficient modeling of complex structures, as it divides the problem into smaller, simpler elements that can be analyzed individually and then assembled into a larger system. (Greco and Peixoto, 2021; Reis and Coda, 2014).

The use of highly deformable and elastic materials, such as elastomers and polymers, has increased in engineering applications, and hyperelastic models are commonly used to represent their behavior. They exhibit a highly nonlinear response to deformation, and therefore, traditional linear elasticity models are not suitable for describing their behavior. (Basheer, 2020; Jing and Jian, 2019).

Hyperelastic models are more appropriate because they can capture the nonlinear stress-strain relationships of these materials. These models are characterized by a strain energy density function, which describes the relationship between stored energy and material deformation. However, it should be noted that hyperelastic models are not universally applicable, and their accuracy may vary depending on the specific material properties and deformation regime under consideration. Therefore, it is important to select an appropriate model that fits the specific material and deformation conditions of interest (Motevalli et al., 2019; Taghizadeh and Darijani, 2018; Talebi and Darijani, 2021; Li et al., 2019).

There are some fundamental distinctions between hyperelastic and linear elastic materials. Specifically, it is noted that the stress-strain relationship of hyperelastic materials is derived from a strain energy density function, which is not the case for linear elastic materials. Additionally, hyperelastic materials exhibit physical nonlinearity, while the geometric nonlinearity of structures is due to displacement affecting the loading condition (Dastjerdi et al., 2022; Seguini and Nedjar, 2016; Darijani and Naghdabadi, 2010).

However, depending on the context, there may be additional nuances to consider, such as the specific behaviors of different types of hyperelastic materials or the impact of loading rate or temperature on material behavior.

2 POSITIONAL FEM FOR PLANE FRAMES

Ramos and Carrazedo (2020) describes that the main characteristic of positional FEM is that it adopts nodal positions as its primary variables instead of displacements. This approach allows for a more accurate nonlinear description of solid geometry. By using nodal positions as primary variables, the method can more precisely capture the behavior of the solid under deformation, which is important for many engineering applications (Paulino and Leonel, 2021).

In this topic, the formulation of the positional finite element method for plane frames developed in Maciel's (2008) studies will be presented.

2.1 Total potential energy functional

The potential energy functional of the system (Equation 1), expressed in the Lagrangian configuration of the object, is provided in the positional formulation by:

$\Pi =U_e+K_c+K_a-P$ (1)

Where $U_e$ = potential energy of elastic deformation; $P$ = potential energy of externally applied forces; $K_c$ = kinetic energy of the body (considered null in static or quasi-static situations); $K_a$ = energy loss due to damping.

The potential energy of elastic deformation ($U_e$) can be described in Equation 2 by the integral of the specific elastic strain energy of the initial volume ($V_0$) as follows:

$U_e={\int }_{V_0}{u}_{e}d{V}_0$ (2)

Where $u_e$ = specific elastic strain energy of the body.

The kinetic energy of the body ($K_c$​) can be described by the integral developed over the initial volume ($V_0$​), considering the density of the body (${\rho }_0$) in the Lagrangian reference frame, as indicated in Equation 3:

$K_c={\int }_{V_0}{\rho }_0\frac{\dot{x_i}\dot{x_i}}{2}d{V}_0$ (3)

Where $\dot{x_i}$ = vectorial velocity of the material point.

The vectorial velocity of the material point ($\dot{x_i}$​) is given in Equation 4 by:

$\dot{x_i}=\frac{d\left(x\right)}{dt}$ (4)

Where $x$ = position vector of the point.

The potential energy of the externally applied forces ($P$) for a system of conservative forces is given in Equation 5 by:

$P=F^{T}X$ (5)

Where $F$ = vector of external forces; $X$ = vector of positions of the applied external forces.

The energy loss due to damping ($K_a\mathit{}$) can be described by the integral resulting from the derivative with respect to positions and is given in Equation 6 by:

$\frac{\partial K_a}{\partial x_i}={\int }_{V_0}\frac{\partial k_a}{\partial x_i}d{V}_0={\int }_{V_0}{c}_m{\rho }_0\dot{x_i}d{V}_0$ (6)

Where $K_a$ = specific dissipated energy functional; $C_m$ = damping constant.

Thus, by substituting all the terms presented separately in Equations 2, 3, and 5 into Equation 1, the potential energy functional of the system can be described as:

$\Pi ={\int }_{V_0}{u}_{e}d{V}_0+{\int }_{V_0}{\rho }_0\frac{\dot{x_i}\dot{x_i}}{2}d{V}_0+K_a-F_j{X}_j$ (7)

Therefore, Maciel (2008) demonstrates up to Equation 7 how potential energy can be expressed as a function of the position of a body, including external forces and kinetic energy. Additionally, Maciel (2008) emphasizes the need to solve the elastic problem by minimizing the potential energy functional using the theorem of minimum total potential energy. This theorem is expressed in Equation 7 and must be applied to the unknowns of the problem at each instant in time, which in this case is the position vector at equilibrium. Thus, we have:

${\left.\frac{\partial \Pi }{\partial x_i}\right|}_t={\int }_{V_0}\frac{\partial u_e}{\partial x_i}d{V}_0+{\int }_{V_0}\frac{\partial }{\partial x_i}\left({\rho }_0\frac{\dot{x_i}\dot{x_i}}{2}\right)d{V}_0+\frac{\partial K_a}{\partial x_i}-\frac{\partial F_j{X}_j}{\partial x_i}$ (8)

Where $x$ = vector of unknown positional measure at body equilibrium.

By substituting Equation 6 into Equation 8, we have:

${\left.\frac{\partial \Pi }{\partial x_i}\right|}_{s+1}={\int }_{V_0}\frac{\partial u_e}{\partial x_i}d{V}_0+{\int }_{V_0}\frac{\partial }{\partial x_i}\left({\rho }_0\frac{\dot{x_i}\dot{x_i}}{2}\right)d{V}_0+{\int }_{V_0}{c}_m{\rho }_0\dot{x_i}d{V}_0-\frac{\partial F_j{X}_j}{\partial x_i}$ (9)

Where $s+1$ = current time; $s$ = immediately preceding time.

The variables of the discretized body, such as positions, accelerations, stresses, and velocities of a finite element, must be estimated using shape functions:

$x_i={Փ}_j{X}_i^j$ (10)

${\dot{x}}_i={Փ}_j{\dot{X}}_i^j$ (11)

${\dot{\ddot{x}}}_i={Փ}_j{\dot{\ddot{X}}}_i^j$ (12)

Where $X_i^j$ = nodal position value; ${\dot{X}}_i^j$ = nodal velocity value; ${\dot{\ddot{X}}}_i^j$ = nodal acceleration value; ${Փ}_j$ = shape functions corresponding to the number of nodes in the discretization.

Thus, from Equations 10, 11, and 12, Equation 9 can be rewritten to represent the equation of motion for the nonlinear problem as follows:

${\left.\frac{\partial \Pi }{\partial X}\right|}_{s+1}={\left.\frac{\partial U_e}{\partial X}\right|}_{s+1}+C{\dot{X}}_{s+1}+M{\ddot{X}}_{s+1}-F_{s+1}$ (13)

Where $M$ = mass matrix; $C$ = damping matrix proportional to mass; $F_{s+1}$ = nodal loads.

This mass matrix can be represented by:

$M={\int }_{V_0}{\rho }_0{Փ}_k{Փ}_{j}d{V}_0$ (14)

Similarly, the damping can be represented by:

$C=2c_{m}M$ (15)

For the nodal loads, we have:

$F_{s+1}=F_0\left[C_1+C_{2}t+C_3{t}^2+C_4{t}^3+C_{5}sin\left(C_{6}t\right)+C_{7}cos\left(C_{8}t\right)+C_9{e}^{c_{{10}^t}}\right]$ (16)

Where $C_1,C_2,C_3,C_4,C_5,C_6,C_7,C_8,C_9,C_{10}$ = define the variation of $F_0$​ over time.

In the scenario where there is no movement, the effects of inertia forces and damping forces balance each other out. As a result, Equation 13 transforms into Equation 17:

$\frac{\partial \Pi }{\partial X_i}=\frac{\partial U_e}{\partial X_i}-F=0$ (17)

The approach to solve Equation 17, which is inherently nonlinear, involves employing the iterative Newton-Raphson method for each time interval. This procedure is described in more detail in the following section.

2.2 Solution to the static problem

According to Maciel (2008), to solve Equation 17 corresponding to the static problem, an iterative technique is essential to ensure that the solution converges accurately. In the context of this study, the Newton-Raphson method is applied. Thus, Equation 17 can be expressed as Equation 18:

$\frac{\partial \Pi }{\partial X_i}=g_i\left(X\right)=f_i\left(X\right)-F_i=0$ (18)

Where $X$ = vector of nodal positions.

Equation 18 can be written in vector form in Equations 19 and 20 as follows:

$g\left(X,F\right)=0$ (19)

$f\left(X\right)-F=0$ (20)

It is worth noting that for the scope of this research, only conservative forces are used. However, Equation 20 encompasses the direct use of non-conservative forces if necessary. It is also emphasized that $g\left(X,F\right)$ demonstrates nonlinear characteristics regarding the parameters $X$ and $F$. According to Maciel (2008), to solve Equation 20, it is necessary to expand Taylor's series up to linear components for a value $\Delta X$, as shown in Equations 21 and 22:

$g\left(X\right)=0\simeq g\left(X_0\right)+\nabla g\left(X_0\right)\Delta X$ (21)

$\Delta X=-{\left[\nabla g\left(X_0\right)\right]}^{-1}g\left(X_0\right)$ (22)

Where $X$ = unknown position vector; $X_0$ = trial position vector.

It is worth noting that the trial position vector $X_0$ is usually determined in the previous iteration. Additionally, from the derivative relative to the nodal position $X_k$, one obtains the Hessian matrix $\nabla g\left(X_0\right)$, which can be expressed in Equation 23 as:

$\nabla g\left(X_0\right)=\frac{\partial f_i\left(X\right)}{\partial X_k}=\frac{{\partial }^2{U}_e}{\partial X_i\partial X_k}$ (23)

From Equation 23, one can observe the symmetry of the Hessian matrix. Consequently, the iterative Newton-Raphson process is described as follows:

1. Select the initial position $X_0$;

2. Using this provisional solution $X_0$, calculate the unbalanced vector $g\left(X_0\right)$ based on Equation 19;

3. Evaluate the Hessian matrix at $X_0$, using Equation 23;

4. Employ Equation 22 to calculate the change in the position vector $\Delta X$ to refine $X_0$;

5. Determine the convergence norm by dividing $\Delta X$ by $X_0$;

6. Check if the calculated norm is less than the specified tolerance;

7. If yes, the method has reached convergence;

8. If not, update $X_0$ to $X$ and repeat the entire procedure until the norm is less than the tolerance.

According to Maciel (2008), depending on the load value and structural shape changes, it is necessary to partition the structure into load levels. This ensures that structural equilibrium is achieved at each level, progressing until the final load balance is achieved.

3 REISSNER KINEMATICS FOR PLANE FRAMES

In this section, Reissner kinematics for plane frames is discussed according to the work of Maciel (2008), where a refined approach is introduced for analyzing the deformation of structures.

3.1 Geometry Mapping

For the study of Reissner kinematics for plane frames according to Maciel (2008), it is essential to establish a mapping of the geometry of the analyzed body and determine the measures of deformation. From Figure 1, one can observe the initial configuration of the body (B₀) as well as the configuration of the body in its current state (B₁). Additionally, it is worth noting the presence of an auxiliary dimensionless space that serves as a facilitator of the connection between the initial configuration B₀ and the current state configuration B₁.

Figure 1
Body Deformation Representation – Adapted from Maciel (2008)

Furthermore, in Figure 1, one can also observe the gradient tensors A^I and A^II, which transition from the temporary configuration of the dimensionless state to the initial configuration of the body B₀, as well as from the temporary configuration of the dimensionless state to the configuration of the body in its current state B₁. It is worth noting that according to Maciel (2008), the gradient tensor A^I depends on the initial shape and the kinematics used, while the gradient tensor A^II depends on the current or final position of the body.

From this, assuming a generic point $p=(x,y)$, it can be represented in terms of dimensionless coordinates, as well as initial or current nodal positions in Equations 24 and 25 as follows:

$p^0\left(\xi ,\eta \right)=p_m^0\left(\xi \right)+\frac{h}{2}\eta {\overrightarrow{N}}^0\left(\xi \right)$ (24)

$p^1\left(\xi ,\eta \right)=p_m^1\left(\xi \right)+\frac{h}{2}\eta {\overrightarrow{N}}^1\left(\xi \right)$ (25)

Where $h$ = thickness; $p_m^i=\left(x_m^i,y_m^i\right)$ = point belonging to the midline; $\overrightarrow{N^i}$ = unit vector determining the position of the section.

Therefore, one can observe the generic cross-section of the element described in Figure 2:

Figure 2
Generic cross-sectional view of the element – Adapted from Maciel (2008)

Thus, according to Maciel (2008), by using the unit vector $\overrightarrow{N}$ in Equations 24 and 25, which describe a generic point in dimensionless coordinates, and converting to Cartesian coordinates, Equation 26 can be obtained:

$\left(x^i\left(\xi ,\eta \right),y^i\left(\xi ,\eta \right)\right)=\left(x_m^i\left(\xi \right),y_m^i\left(\xi \right)\right)+\frac{h}{2}\eta \left(-sen{\theta }^i\left(\xi \right),cos{\theta }^i\left(\xi \right)\right)$ (26)

With that, Equation 26 can be expressed in terms of $x$ and $y$ in Equations 27 and 28:

$x^i\left(\xi ,\eta \right)=x_m^i\left(\xi \right)-\frac{h}{2}\eta sen{\theta }^i\left(\xi \right)$ (27)

$y^i\left(\xi ,\eta \right)=y_m^i\left(\xi \right)+\frac{h}{2}\eta cos{\theta }^i\left(\xi \right)$ (28)

Where $i=0$ = initial configuration; $i=1$ = current configuration.

Furthermore, from Equations 27 and 28, one can derive in dimensionless coordinates to obtain the gradients A^I and A^II in Equation 29 as follows:

$A_i=\left[\begin{array}{cc}{A}_{11}^i& A_{12}^i\\ A_{21}^i& A_{22}^i\end{array}\right]=\left[\begin{array}{cc}\frac{d{x}^i}{d\xi }& \frac{d{x}^i}{d\eta }\\ \frac{d{y}^i}{d\xi }& \frac{d{y}^i}{d\eta }\end{array}\right]$ (29)

Thus, the elements of $A_i$ can be obtained in Equations 30, 31, 32, and 33:

$A_{11}^i=\frac{d{x}_m}{d\xi }-\frac{h}{2}\eta \frac{d\theta }{d\xi }cos\theta$ (30)

$A_{12}^i=\frac{-h}{2}sen\theta$ (31)

$A_{21}^i=\frac{d{x}_m}{d\xi }-\frac{h}{2}\eta \frac{d\theta }{d\xi }sen\theta$ (32)

$A_{22}^i=\frac{h}{2}cos\theta$ (33)

Where $i=I$ = dimensionless space transformation to configuration B₀; $i=II$= dimensionless space transformation to configuration B₁.

3.2 Deformation Measurement

After the proper mapping of the geometry and obtaining the gradients A^I and A^II, attention should be paid to the calculation of elongations in the dimensionless directions $\xi$ and $\eta$. For the transverse elongation ${\lambda }_t$, one has the unit vector $\overrightarrow{M_t}={\left[\mathrm{1,0}\right]}^T$, while for the longitudinal elongation ${\lambda }_n$, one has the unit vector $\overrightarrow{M_n}={\left[\mathrm{0,1}\right]}^T$ according to Maciel (2008). It is possible to observe the unit vectors $\overrightarrow{M_t}$ and $\overrightarrow{M_n}$ in Figure 3:

Figure 3
Auxiliary dimensionless space and geometry mapping - Adapted from Maciel (2008)

With that, elongations can be obtained through Equations 34 and 35:

${\lambda }_t^i=\lambda \left(\overrightarrow{M_t}\right)=\left|A^i\left[\begin{array}{c}1\\ 0\end{array}\right]\right|=\sqrt{{\left(A_{11}^i\right)}^2+{\left(A_{21}^i\right)}^2}$(34)

${\lambda }_n^i=\lambda \left(\overrightarrow{M_n}\right)=\left|A^i\left[\begin{array}{c}0\\ 1\end{array}\right]\right|=1$ (35)

From Equation 36, the calculation of the angle $\alpha$ represented in Figure 3 is obtained, as the angle between the vectors in configuration B₁.

$\alpha =arccos\left\{\frac{\left[10\right]\left[A_{II}^T{A}_{II}\right]\left[\begin{array}{c}0\\ 1\end{array}\right]}{A_t^{II}{A}_n^{II}}\right\}=arccos\left\{\frac{A_{11}^{II}{A}_{12}^{II}+A_{21}^{II}{A}_{22}^{II}}{\sqrt{{\left(A_{11}^{II}\right)}^2+{\left(A_{21}^{II}\right)}^2}}\right\}$ (36)

Thus, the distortion of configuration B₀ can be calculated from Equation 37:

${\gamma }_m=\alpha -\frac{\pi }{2}=arccos\left\{\frac{A_{11}^{II}{A}_{12}^{II}+A_{21}^{II}{A}_{22}^{II}}{\sqrt{{\left(A_{11}^{II}\right)}^2+{\left(A_{21}^{II}\right)}^2}}\right\}-\frac{\pi }{2}$ (37)

Similarly, the elongations of configuration B₀ can be calculated from Equations 38 and 39:

${\lambda }_t=\frac{{\lambda }_t^{II}}{{\lambda }_t^I}=\frac{\sqrt{{\left(A_{11}^{II}\right)}^2+{\left(A_{21}^{II}\right)}^2}}{\sqrt{{\left(A_{11}^I\right)}^2+{\left(A_{21}^I\right)}^2}}$ (38)

${\lambda }_n=\frac{{\lambda }_n^{II}}{{\lambda }_n^I}=1$ (39)

Additionally, the strains can be calculated from Equations 40 and 41:

${\epsilon }_t={\lambda }_t-1=\frac{\sqrt{{\left(A_{11}^{II}\right)}^2+{\left(A_{21}^{II}\right)}^2}}{\sqrt{{\left(A_{11}^I\right)}^2+{\left(A_{21}^I\right)}^2}}-1$ (40)

${\epsilon }_n={\lambda }_n-1=0$ (41)

Where ${\epsilon }_t$ = nonlinear engineering strain measurement in direction $\overrightarrow{T}$ in configuration B₀; ${\epsilon }_n$ = nonlinear engineering strain measurement in direction $\overrightarrow{N}$ in configuration B₀.

3.3 Static Equations

In the static equation, the terms related to the kinetic energy of the body and the damping energy are removed from the total potential energy functional of the system. Therefore, one has Equation 42:

$\Pi ={\int }_{V_0}\frac{1}{2}E\left({\epsilon }_t^2+\frac{{\gamma }_m^2}{2}\right)d{V}_0-F_j{X}_j$ (42)

Therefore, the discretization of the plane frame must be considered according to Figure 4:

Figure 4
Finite element with 'n' nodes and its nodal variables - Adapted from Maciel (2008)

From the consideration of the finite element shown in Figure 4, the unknowns $x_m$, $y_m$ e $\theta$ can be obtained as follows according to Maciel (2008) in Equations 43, 44, and 45:

$x_m\left(\xi ,\eta \right)={\Phi }_i{X}_i$ (43)

$y_m\left(\xi ,\eta \right)={\Phi }_i{Y}_i$ (44)

$\theta \left(\xi ,\eta \right)={\Phi }_i{\Theta }_i$ (45)

Where $X$, $Y$, $\Theta$ = nodal positions of the finite element; $\Phi$ = shape approximation functions.

With that, Equations 27 and 28 are resumed as follows in Equations 46 and 47:

$x^i={\left({\Phi }_k{X}_k\right)}^i-\frac{h}{2}\eta sen{\left({\Phi }_k{\Theta }_k\right)}^i$ (46)

$y^i={\left({\Phi }_k{Y}_k\right)}^i+\frac{h}{2}\eta cos{\left({\Phi }_k{\Theta }_k\right)}^i$ (47)

From this, Equations 30, 31, 32, and 33 are resumed as follows in Equations 48, 49, 50, and 51:

$A_{11}^i={\beta }_j{X}_j^i-\frac{h}{2}\eta \left({\beta }_j{\Theta }_j^i\right)cos\left({\Phi }_k{\Theta }_k^i\right)$ (48)

$A_{12}^i=\frac{-h}{2}\eta sen\left({\Phi }_k{\Theta }_k^i\right)$ (49)

$A_{21}^i={\beta }_j{Y}_j^i-\frac{h}{2}\eta \left({\beta }_j{\Theta }_j^i\right)sen\left({\Phi }_k{\Theta }_k^i\right)$ (50)

$A_{22}^i=\frac{h}{2}\eta cos\left({\Phi }_k{\Theta }_k^i\right)$ (51)

Where in Equation 52:

${\beta }_j=\frac{d{\Phi }_j}{d\xi }$ (52)

Therefore, one has in Equation 53 and 54:

$\Pi \to f\left(p\right)$ (53)

$p=\left(X_1,Y_1,{\Theta }_1,X_2,Y_2,{\Theta }_2,X_3,Y_3,{\Theta }_3,...,X_{\left(3n-1\right)},Y_{\left(2n-1\right)},{\Theta }_{\left(3n\right)}\right)$ (54)

With this, the method consists of minimizing the total potential energy functional for the plane frames (Equation 42) in terms of nodal positions. According to Maciel (2008), the energy functional is derived with respect to the degree of freedom p_i in Equations 55 and 56:

$\frac{\partial \Pi }{\partial p_i}={\int }_{V_0}\frac{1}{2}E\frac{\partial }{\partial p_i}\left({\epsilon }_t^2+\frac{{\gamma }_m^2}{2}\right)d{V}_0-F_i$ (55)

$f_i^m\left(p\right)={\int }_{V_0}\frac{1}{2}E\frac{\partial }{\partial p_i}\left({\epsilon }_t^2+\frac{{\gamma }_m^2}{2}\right)d{V}_0$ (56)

Where $f_i^m\left(p\right)$ = vector of internal forces referring to position p_i.

4 HYPERELASTIC MODELS

Hyperelastic materials are a particular type of material that exhibit large deformations and are commonly used in applications such as rubber and biological tissues. The notable difference between hyperelastic materials and linear elastic materials is that to accurately model the behavior of hyperelastic materials using the positional finite element method, it is necessary to express their constitutive equation in terms of the derivative of the strain energy density, as can be observed in Figure 5. This is because hyperelastic materials store and release energy during deformation, and this energy is related to the material's stress state. However, there are challenges associated with deriving the strain energy density for a hyperelastic material, which can be a complex and time-consuming process. (Perônica et al., 2022; Anssari-Benam and Horgan, 2022; Khaniki et al., 2022)

Figure 5
Stress-strain curve for linear elastic (a) and hyperelastic (b) models respectively - Adapted from Perônica et al. (2022)

For a more in-depth study of hyperelastic models, according to Perônica et al. (2022), one should address the work of Ogden (1984) on the description of the motion of bodies with the theory of continuum mechanics. According to Ogden (1984), one has Equation 57:

$X=\chi \left(x,t\right)$ (57)

Where $X$ = position vector of the body with respect to the deformed configuration at time t; $x$ = position vector of the body with respect to the reference configuration at time t.

With this, one has the deformation gradient A in Equations 58 and 59:

$A_{ij}=\frac{\partial X_i}{\partial x_j}$ (58)

$A_{ij}=\left[\begin{array}{ccc}{\lambda }_1& 0& 0\\ 0& {\lambda }_2& 0\\ 0& 0& {\lambda }_3\end{array}\right]$ (59)

Where ${\lambda }_i$ = principal stretches (eigenvalues of the elastic tensors).

From this, one has the right and left Cauchy-Green deformation tensors respectively in Equations 60 and 61 according to Perônica et al. (2022):

$C=A^{T}A$ (60)

$B=A{A}^T$ (61)

So, one has Equation 62:

$C_{ij}=B_{ij}=\left[\begin{array}{ccc}{\lambda }_1^2& 0& 0\\ 0& {\lambda }_2^2& 0\\ 0& 0& {\lambda }_3^2\end{array}\right]$ (62)

Where the first principal invariant I₁ can be obtained through the trace of the Cauchy-Green deformation tensor in Equation 63:

$I_1={tr\left(BouC\right)=\lambda }_1^2+{\lambda }_2^2+{\lambda }_3^2$ (63)

Additionally, the second principal invariant I₂ can be obtained through Equation 64:

$I_2=\frac{I_1^2-tr\left(B^2\mathrm{}ou{\mathrm{}C}^2\right)}{2}={\lambda }_1^2{\lambda }_2^2+{\lambda }_2^2{\lambda }_3^2+{\lambda }_1^2{\lambda }_3^2$ (64)

Moreover, the third principal invariant I₃ can also be obtained through the determinant of the Cauchy-Green deformation tensor in Equation 65:

$I_3=det\left(B\mathrm{}ou\mathrm{}C\right)={\lambda }_1^2{\lambda }_2^2{\lambda }_3^2=J^2$ (65)

Where $J$ = determinant of A (volume change).

It is worth noting that according to Perônica et al. (2022), considering the assumption of incompressibility, the second and third directions are equal, therefore ${\lambda }_2=\frac{1}{\sqrt{{\lambda }_1}}$ and ${\lambda }_3=\frac{1}{\sqrt{{\lambda }_1}}$. Thus, one has Equations 66, 67, and 68:

$I_1={\lambda }_1^2+{\left(\frac{1}{\sqrt{{\lambda }_1}}\right)}^2+{\left(\frac{1}{\sqrt{{\lambda }_1}}\right)}^2={\lambda }_1^2+\frac{2}{{\lambda }_1}$ (66)

$I_2={\lambda }_1^2{\left(\frac{1}{\sqrt{{\lambda }_1}}\right)}^2+{\left(\frac{1}{\sqrt{{\lambda }_1}}\right)}^2{\left(\frac{1}{\sqrt{{\lambda }_1}}\right)}^2+{\lambda }_1^2{\left(\frac{1}{\sqrt{{\lambda }_1}}\right)}^2=\frac{1}{{\lambda }_1^2}+2{\lambda }_1$ (67)

$I_3={\lambda }_1^2{\left(\frac{1}{\sqrt{{\lambda }_1}}\right)}^2{\left(\frac{1}{\sqrt{{\lambda }_1}}\right)}^2=1$ (68)

Therefore, as an essential aspect of a hyperelastic model, there is the energy density function (often denoted as $W$ or $\Psi$), which is a mathematical function that describes how the stored energy within a material changes as it deforms. (Jerábek and Écsi, 2019)

This function is expressed in terms of deformation measures (Equation 69), which quantify how much a material has been deformed. The most common deformation measures include the use of strain invariants, as well as the principal stretches.

$\psi =\psi \left(I_1,I_2,I_3\right)=\psi \left({\lambda }_1,{\lambda }_2,{\lambda }_3\right)$ (69)

With that, according to Perônica et al. (2022), the Cauchy tensor component can be represented in terms of the principal stretches in Equations 70 and 71 as follows:

${\sigma }_i={\lambda }_i\frac{\partial \psi }{\partial {\lambda }_i}$ (70)

$\sigma =2\left(\lambda -\frac{1}{{\lambda }^2}\right)\left(\frac{\partial \psi }{\partial I_1}+\frac{1}{\lambda }\frac{\partial \psi }{\partial I_2}\right)$ (71)

4.1 Mooney-Rivlin Model

The Mooney-Rivlin hyperelastic model is a constitutive model used to describe the deformation behavior of rubber-like materials such as elastomers. In this model, the strain energy density function is represented as a polynomial function of the components of the strain tensor. (Mooney, 1940)

The mathematical expression for the strain energy density function in the Mooney-Rivlin hyperelastic model is given by Equation 72 as:

${\psi }_{MR}\left(I_1,I_2\right)=C_{10}\left(I_1-3\right)+C_{01}\left(I_2-3\right)$ (72)

Where $C_{10}$ and $C_{01}$ are material constants that depend on the specific material being modeled.

Substituting the first and second strain invariants obtained previously, one has Equation 73:

${\psi }_{MR}=C_{10}\left({\lambda }_1^2+\frac{2}{{\lambda }_1}-3\right)+C_{01}\left(\frac{1}{{\lambda }_1^2}+2{\lambda }_1-3\right)$ (73)

From this, we can obtain the Cauchy stress tensor component using Equation 74:

${\sigma }_{MR}=2\left({\lambda }_1-\frac{1}{{\lambda }_1^2}\right)\left(C_{10}+\frac{C_{01}}{{\lambda }_1}\right)$ (74)

4.2 Neo-Hookean Model

The Neo-Hookean model is a type of hyperelastic model obtained through a particularization of the Mooney-Rivlin model and is used to describe the mechanical behavior of materials undergoing large deformations, such as rubber or soft biological tissues. This model assumes that the material is isotropic, meaning that its properties are the same in all directions.

The strain energy function for the Neo-Hookean model can be expressed as Equation 75:

${\psi }_{NH}\left(I_1\right)=C_{10}\left(I_1-3\right)$ (75)

Substituting the first strain invariant obtained previously, one has in Equation 76:

${\psi }_{NH}=C_{10}\left({\lambda }_1^2+\frac{2}{{\lambda }_1}-3\right)$ (76)

From this, we can obtain the Cauchy stress tensor component in Equation 77:

${\sigma }_{NH}=2\left({\lambda }_1-\frac{1}{{\lambda }_1^2}\right)\left(C_{10}\right)$ (77)

4.3 Ogden Model

The Ogden hyperelastic model has a mathematical structure, more complex than previous methods, that is primarily used to describe the mechanical behavior of elastomers and predict the stress-strain relationship of elastomers under different loading conditions and for very high deformations. (Ogden, 1972; Ogden, 1984)

The strain energy density function for the Ogden model is expressed in Equation 78 as follows:

${\psi }_{OG}\left({\lambda }_1,{\lambda }_2,{\lambda }_3\right)={\sum }_{p=1}^N\frac{2{\mu }_p}{{\alpha }_p}\left({\lambda }_1^{\left({\alpha }_p\right)}+{\lambda }_2^{\left({\alpha }_p\right)}+{\lambda }_3^{\left({\alpha }_p\right)}-3\right)$ (78)

Where ${\mu }_p,{\alpha }_p$= material constants depending on the specific material being modeled.

It is worth noting that according to Perônica et al. (2022), considering the assumption of incompressibility made in the previous chapter, the second and third directions are equal. Therefore, ${\lambda }_2=\frac{1}{\sqrt{{\lambda }_1}}$ and ${\lambda }_3=\frac{1}{\sqrt{{\lambda }_1}}$. Thus, Equation 79 is obtained.

${\psi }_{OG}\left({\lambda }_1,{\lambda }_2,{\lambda }_3\right)={\sum }_{p=1}^N\frac{2{\mu }_p}{{\alpha }_p}\left({\lambda }_1^{\left({\alpha }_p\right)}+{\left(\frac{2}{\sqrt{{\lambda }_1}}\right)}^{\left({\alpha }_p\right)}-3\right)$ (79)

To demonstrate the mathematical formulation of the Ogden hyperelastic model, we can start with the definition of the model for $N=3$ in Equation 80.

${\psi }_{OG}=\frac{2{\mu }_1}{{\alpha }_1}\left({\lambda }_1^{\left({\alpha }_1\right)}+{\left(\frac{2}{\sqrt{{\lambda }_1}}\right)}^{\left({\alpha }_1\right)}-3\right)+\frac{2{\mu }_2}{{\alpha }_2}\left({\lambda }_1^{\left({\alpha }_2\right)}+{\left(\frac{2}{\sqrt{{\lambda }_1}}\right)}^{\left({\alpha }_2\right)}-3\right)+\frac{2{\mu }_3}{{\alpha }_3}\left({\lambda }_1^{\left({\alpha }_3\right)}+{\left(\frac{2}{\sqrt{{\lambda }_1}}\right)}^{\left({\alpha }_3\right)}-3\right)$ (80)

From this, the Cauchy stress tensor component can be obtained, for $N=3$, in Equation 81.

${\sigma }_{OG}=2{\mu }_1\left({\lambda }_1^{\left({\alpha }_1-1\right)}-{{\lambda }_1}^{\left(-\frac{{\alpha }_1}{2}-1\right)}\right)+2{\mu }_2\left({\lambda }_1^{\left({\alpha }_2-1\right)}-{{\lambda }_1}^{\left(-\frac{{\alpha }_2}{2}-1\right)}\right)+2{\mu }_3\left({\lambda }_1^{\left({\alpha }_3-1\right)}-{{\lambda }_1}^{\left(-\frac{{\alpha }_3}{2}-1\right)}\right)$ (81)

4.4 Yeoh Model

The Yeoh (1990) hyperelastic model is a mathematical framework derived from the polynomial model and an extension of the Neo-Hookean model. The strain energy density function $\psi$ in the Yeoh hyperelastic model is typically represented as follows in Equation 82:

${\psi }_{YE}\left(I_1\right)={\sum }_{p=1}^N{C}_{p0}{\left(I_1-3\right)}^p$ (82)

To demonstrate the mathematical formulation of the Yeoh model, one can start with the definition of the model for $N=3$ in Equation 83:

${\psi }_{YE}=C_{10}\left(I_1-3\right)+C_{20}{\left(I_1-3\right)}^2+C_{30}{\left(I_1-3\right)}^3$ (83)

Substituting the first invariant of deformation obtained previously for $N=3$, one has in Equation 84:

${\psi }_{YE}=C_{10}\left({\lambda }_1^2+\frac{2}{{\lambda }_1}-3\right)+C_{20}{\left({\lambda }_1^2+\frac{2}{{\lambda }_1}-3\right)}^2+C_{30}{\left({\lambda }_1^2+\frac{2}{{\lambda }_1}-3\right)}^3$ (84)

From this, the Cauchy tensor component can be obtained, for $N=3$, in Equation 85:

${\sigma }_{YE}=2\left({\lambda }_1-\frac{1}{{\lambda }_1^2}\right)\left(C_{10}+{2C}_{20}\left(I_1-3\right)+3C_{30}\left(I_1-3\right)²\right)$ (85)

4.5 Internal force vector and Hessian matrix

For the context of this research, it is necessary to calculate the parameters internal force vector $q_i$ and Hessian matrix $K_{ij}$ for plane frames according to Reissner's kinematics and the formulation of hyperelastic models based on the stretch $\lambda$ addressed in previous chapters.

4.5.1 Linear Elastic Model

According to Maciel (2008), for the linear elastic model, the specific strain energy is given by Equation 86:

$u_e=\frac{1}{2}E\left({\epsilon }_t^2+\frac{{\gamma }_{tn}^2}{2}\right)$ (86)

With that, Equation 86 is substituted into Equation 87 to obtain the elastic potential energy in Equation 88:

$U_e={\int }_{V_0}{u}_{e}d{V}_0$ (87)

$U_e={\int }_{V_0}\frac{1}{2}E\left({\epsilon }_t^2+\frac{{\gamma }_{tn}^2}{2}\right)d{V}_0$ (88)

Thus, the internal force vector $q_i$ is obtained from the gradient of the elastic potential energy represented by Equation 89:

$q_i=\nabla U_e=\frac{\partial U_e}{\partial X_i}=\frac{\partial }{\partial X_i}\left[{\int }_{V_0}\frac{1}{2}E\left({\epsilon }_t^2+\frac{{\gamma }_{tn}^2}{2}\right)d{V}_0\right]$ (89)

$q_i=\frac{1}{2}E\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}\frac{\partial }{\partial X_i}\left({\epsilon }_t^2+\frac{{\gamma }_{tn}^2}{2}\right)b{J}_{0}d\xi d\eta$ (90)

$q_i=Eb{J}_0\left({\mathrm{\epsilon }}_{\mathrm{t}}\frac{\partial {\mathrm{\epsilon }}_t}{\partial X_i}+\frac{{\mathrm{\gamma }}_{tn}}{2}\frac{\partial {\mathrm{\gamma }}_{tn}}{\partial X_i}\right){\omega }_1{\omega }_2$ (91)

Where $E$ = Young's modulus; $b$ = width of the element; $J_0$ = Jacobian of the numerical integration transformation by Gauss-Legendre method; ${\omega }_1{,\omega }_2$ = weights associated with numerical integration by Gauss-Legendre method.

Furthermore, by applying the gradient to the internal force vector in Equation 89, the Hessian matrix $K_{ij}$ can be found from Equation 92.

$K_{ij}=\nabla q_i=\frac{\partial }{\partial X_j}\left[\frac{\partial }{\partial X_i}\left[{\int }_{V_0}\frac{1}{2}E\left({\epsilon }_t^2+\frac{{\gamma }_{tn}^2}{2}\right)d{V}_0\right]\right]$ (92)

$K_{ij}=\frac{1}{2}E\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}\frac{\partial }{\partial X_j}\frac{\partial }{\partial X_i}\left({\epsilon }_t^2+\frac{{\gamma }_{tn}^2}{2}\right)b{J}_{0}d\xi d\eta$ (93)

$K_{ij}=Eb{J}_0\left(\left(\frac{\partial {\mathrm{\epsilon }}_t^2}{\partial X_{ij}}+{\mathrm{\epsilon }}_{\mathrm{t}}\frac{{\partial }^2{\mathrm{\epsilon }}_t}{\partial X_{ij}^2}\right)+\frac{1}{2}\frac{\partial {\mathrm{\gamma }}_{tn}^2}{\partial X_{ij}}+\frac{1}{2}{\mathrm{\gamma }}_{tn}\frac{{\partial }^2{\mathrm{\gamma }}_{tn}}{\partial X_{ij}^2}\right){\omega }_1{\omega }_2$ (94)

4.5.2 Neo-Hookean Model

In a similar manner to what was done for the linear elastic model, in hyperelastic models, the specific free deformation energy is used as the material constitutive law. By substituting the specific free energy of the Neo-Hookean model given in Equation 76 into Equation 87, one has from Equation 95:

$U_e={\int }_{V_0}{C}_{10}\left({\lambda }_1^2+\frac{2}{{\lambda }_1}-3\right)d{V}_0$ (95)

However, when conducting preliminary tests on implementing this formulation, it was noted that the equations for the strain invariants from uniaxial tensile tests are not sufficient for the approach of planar frames with Reissner kinematics. Therefore, the first strain invariant from simple shear tests was added to consider distortion in the deformation energy of the model.

Therefore, one has the following equation from Equation 96:

$U_e={\int }_{V_0}{C}_{10}\left(I_{1_{TU}}+I_{1_{CS}}-3\right)d{V}_0$ (96)

$U_e={\int }_{V_0}{C}_{10}\left({(\lambda }_1^2+\frac{2}{{\lambda }_1})+({\mathrm{\gamma }}_{tn}^2+3)-3\right)d{V}_0$ (97)

$U_e={\int }_{V_0}{C}_{10}\left({\lambda }_1^2+\frac{2}{{\lambda }_1}+{\mathrm{\gamma }}_{tn}^2\right)d{V}_0$ (98)

Where $I_{1_{TU}}$ is the first strain invariant for uniaxial tensile test, and $I_{1_{CS}}$ is the first strain invariant for simple shear test.

Thus, the internal force vector $q_i$ is obtained by the gradient of the elastic potential energy represented by Equation 99:

$q_i=\nabla U_e=\frac{\partial U_e}{\partial X_i}=\frac{\partial }{\partial X_i}\left[{\int }_{V_0}{C}_{10}\left({\lambda }_1^2+\frac{2}{{\lambda }_1}+{\mathrm{\gamma }}_{tn}^2\right)d{V}_0\right]$ (99)

$q_i=C_{10}\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}\frac{\partial }{\partial X_i}\left({\lambda }_1^2+\frac{2}{{\lambda }_1}+{\mathrm{\gamma }}_{tn}^2\right)b{J}_{0}d\xi d\eta$ (100)

$q_i=C_{10}b{J}_0\left(2{\lambda }_1\frac{\partial {\lambda }_1}{\partial X_i}-\frac{2\frac{\partial {\lambda }_1}{\partial X_i}}{{\lambda }_1^2}+2{\mathrm{\gamma }}_{tn}\frac{\partial {\mathrm{\gamma }}_{tn}}{\partial X_i}\right){\omega }_1{\omega }_2$ (101)

Furthermore, by applying the gradient to the internal force vector in Equation 99, one can find the Hessian matrix $K_{ij}$ from Equation 102:

$K_{ij}=\nabla q_i=\frac{\partial }{\partial X_j}\left[\frac{\partial }{\partial X_i}\left[{\int }_{V_0}{C}_{10}\left({\lambda }_1^2+\frac{2}{{\lambda }_1}+{\mathrm{\gamma }}_{tn}^2\right)d{V}_0\right]\right]$ (102)

$K_{ij}=C_{10}b{J}_0\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}\frac{\partial }{\partial X_j}\frac{\partial }{\partial X_i}\left({\lambda }_1^2+\frac{2}{{\lambda }_1}+{\mathrm{\gamma }}_{tn}^2\right)d\xi d\eta$ (103)

$K_{ij}=C_{10}b{J}_0\left(2\left(\frac{\partial {\lambda }_1^2}{\partial X_{ij}}+{\lambda }_1\frac{{\partial }^2{\lambda }_1}{\partial X_{ij}^2}\right)-\left(\frac{2\left({\lambda }_1\frac{{\partial }^2{\lambda }_1}{\partial X_{ij}^2}-2\frac{\partial {\lambda }_1^2}{\partial X_{ij}}\right)}{{\lambda }_1^3}\right)+2\left(\frac{\partial {\mathrm{\gamma }}_{tn}^2}{\partial X_{ij}}+{\mathrm{\gamma }}_{tn}\frac{{\partial }^2{\mathrm{\gamma }}_{tn}}{\partial X_{ij}^2}\right)\right){\omega }_1{\omega }_2$ (104)

4.5.3 Mooney-Rivlin Model

By substituting the specific free energy of the Mooney-Rivlin model as shown in Equation 73 into Equation 87, one has from Equation 105:

$U_e={\int }_{V_0}{C}_{10}\left({\lambda }_1^2+\frac{2}{{\lambda }_1}-3\right)+C_{01}\left(\frac{1}{{\lambda }_1^2}+2{\lambda }_1-3\right)d{V}_0$ (105)

However, upon conducting preliminary implementation tests of this formulation, it was noticed that the equation for the strain invariants of uniaxial traction test is not sufficient for the approach of planar frames with Reissner kinematics. Therefore, the first and second strain invariants of simple shear test were added so that distortion is considered in the deformation energy of the model.

Therefore, the following equation is obtained from Equation 106:

$U_e={\int }_{V_0}{C}_{10}\left(I_{1_{TU}}+I_{1_{CS}}-3\right)+C_{01}\left(I_{2_{TU}}+I_{2_{CS}}-3\right)d{V}_0$ (106)

$U_e={\int }_{V_0}{C}_{10}\left({(\lambda }_1^2+\frac{2}{{\lambda }_1})+({\mathrm{\gamma }}_{tn}^2+3)-3\right)+C_{01}\left(\left(\frac{1}{{\lambda }_1^2}+2{\lambda }_1\right)+({\mathrm{\gamma }}_{tn}^2+3)-3\right)d{V}_0$ (107)

$U_e={\int }_{V_0}{C}_{10}\left({\lambda }_1^2+\frac{2}{{\lambda }_1}+{\mathrm{\gamma }}_{tn}^2\right)+C_{01}\left(\frac{1}{{\lambda }_1^2}+2{\lambda }_1+{\mathrm{\gamma }}_{tn}^2\right)d{V}_0$ (108)

Thus, the internal force vector $q_i$ is obtained from the gradient of the elastic potential energy represented by Equation 109:

$q_i=\nabla U_e=\frac{\partial U_e}{\partial X_i}=\frac{\partial }{\partial X_i}\left[{\int }_{V_0}{C}_{10}\left({\lambda }_1^2+\frac{2}{{\lambda }_1}+{\mathrm{\gamma }}_{tn}^2\right)+C_{01}\left(\frac{1}{{\lambda }_1^2}+2{\lambda }_1+{\mathrm{\gamma }}_{tn}^2\right)d{V}_0\right]$ (109)

$q_i=\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}\left(C_{10}\frac{\partial }{\partial X_i}\left({\lambda }_1^2+\frac{2}{{\lambda }_1}+{\mathrm{\gamma }}_{tn}^2\right)+C_{01}\frac{\partial }{\partial X_i}\left(\frac{1}{{\lambda }_1^2}+2{\lambda }_1+{\mathrm{\gamma }}_{tn}^2\right)\right)b{J}_{0}d\xi d\eta$ (110)

$q_i=b{J}_0\left(C_{10}\left(2{\lambda }_1\frac{\partial {\lambda }_1}{\partial X_i}-\frac{2\frac{\partial {\lambda }_1}{\partial X_i}}{{\lambda }_1^2}+2{\mathrm{\gamma }}_{tn}\frac{\partial {\mathrm{\gamma }}_{tn}}{\partial X_i}\right)+C_{01}\left(-\frac{2\frac{\partial {\lambda }_1}{\partial X_i}}{{\lambda }_1^3}+2\frac{\partial {\lambda }_1}{\partial X_i}+2{\mathrm{\gamma }}_{tn}\frac{\partial {\mathrm{\gamma }}_{tn}}{\partial X_i}\right)\right){\omega }_1{\omega }_2$ (111)

Furthermore, by applying the gradient to the internal force vector in Equation 109, we can find the Hessian matrix $K_{ij}$ from Equation 112:

$K_{ij}=\nabla q_i=\frac{\partial }{\partial X_j}\left[\frac{\partial }{\partial X_i}\left[{\int }_{V_0}{C}_{10}\left({\lambda }_1^2+\frac{2}{{\lambda }_1}+{\mathrm{\gamma }}_{tn}^2\right)+C_{01}\left(\frac{1}{{\lambda }_1^2}+2{\lambda }_1+{\mathrm{\gamma }}_{tn}^2\right)d{V}_0\right]\right]$ (112)

$K_{ij}=\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}\left(C_{10}\frac{\partial }{\partial X_j}\frac{\partial }{\partial X_i}\left({\lambda }_1^2+\frac{2}{{\lambda }_1}+{\mathrm{\gamma }}_{tn}^2\right)+C_{01}\frac{\partial }{\partial X_j}\frac{\partial }{\partial X_i}\left(\frac{1}{{\lambda }_1^2}+2{\lambda }_1+{\mathrm{\gamma }}_{tn}^2\right)\right)b{J}_{0}d\xi d\eta$ (113)

$K_{ij}=b{J}_0\left(C_{10}\left(2\left(\frac{\partial {\lambda }_1^2}{\partial X_{ij}}+{\lambda }_1\frac{{\partial }^2{\lambda }_1}{\partial X_{ij}^2}\right)-\left(\frac{2\left({\lambda }_1\frac{{\partial }^2{\lambda }_1}{\partial X_{ij}^2}-2\frac{\partial {\lambda }_1^2}{\partial X_{ij}}\right)}{{\lambda }_1^3}\right)+2\left(\frac{\partial {\mathrm{\gamma }}_{tn}^2}{\partial X_{ij}}+{\mathrm{\gamma }}_{tn}\frac{{\partial }^2{\mathrm{\gamma }}_{tn}}{\partial X_{ij}^2}\right)\right)+C_{01}\left(-\left(\frac{2\left({\lambda }_1\frac{{\partial }^2{\lambda }_1}{\partial X_{ij}^2}-3\frac{\partial {\lambda }_1^2}{\partial X_{ij}}\right)}{{\lambda }_1^4}\right)+2\frac{{\partial }^2{\lambda }_1}{\partial X_{ij}^2}+2\left(\frac{\partial {\mathrm{\gamma }}_{tn}^2}{\partial X_{ij}}+{\mathrm{\gamma }}_{tn}\frac{{\partial }^2{\mathrm{\gamma }}_{tn}}{\partial X_{ij}^2}\right)\right)\right){\omega }_1{\omega }_2$ (114)

4.5.4 Ogden Model

By substituting the specific strain energy of the Ogden model as exposed in Equation 79 into Equation 87, one has, for $N=1$, Equation 115:

$U_e={\int }_{V_0}\frac{2{\mu }_1}{{\alpha }_1}\left({\lambda }_1^{\left({\alpha }_1\right)}+{\left(\frac{2}{\sqrt{{\lambda }_1}}\right)}^{\left({\alpha }_1\right)}-3\right)d{V}_0$ (115)

However, upon conducting preliminary implementation tests of this formulation, it was noted that the equation is not sufficient for the approach of planar frames with Reissner kinematics. Therefore, the first invariant of simple shear deformation test was added so that distortion is considered in the deformation energy of the model.

Therefore, one has the following equation starting from Equation 116:

$U_e={\int }_{V_0}\frac{2{\mu }_1}{{\alpha }_1}\left({\lambda }_1^{\left({\alpha }_1\right)}+{\left(\frac{2}{\sqrt{{\lambda }_1}}\right)}^{\left({\alpha }_1\right)}+I_{1_{CS}}-3\right)d{V}_0$ (116)

$U_e={\int }_{V_0}\frac{2{\mu }_1}{{\alpha }_1}\left({\lambda }_1^{\left({\alpha }_1\right)}+{\left(\frac{2}{\sqrt{{\lambda }_1}}\right)}^{\left({\alpha }_1\right)}+({\mathrm{\gamma }}_{tn}^2+3)-3\right)d{V}_0$ (117)

$U_e={\int }_{V_0}\frac{2{\mu }_1}{{\alpha }_1}\left({\lambda }_1^{\left({\alpha }_1\right)}+{\left(\frac{2}{\sqrt{{\lambda }_1}}\right)}^{\left({\alpha }_1\right)}+{\mathrm{\gamma }}_{tn}^2\right)d{V}_0$ (118)

Thus, the internal force vector $q_i$ is obtained from the gradient of the elastic potential energy represented by Equation 119:

$q_i=\nabla U_e=\frac{\partial U_e}{\partial X_i}=\frac{\partial }{\partial X_i}\left[{\int }_{V_0}\frac{2{\mu }_1}{{\alpha }_1}\left({\lambda }_1^{\left({\alpha }_1\right)}+{\left(\frac{2}{\sqrt{{\lambda }_1}}\right)}^{\left({\alpha }_1\right)}+{\mathrm{\gamma }}_{tn}^2\right)d{V}_0\right]$ (119)

$q_i=\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}\frac{2{\mu }_1}{{\alpha }_1}\frac{\partial }{\partial X_i}\left({\lambda }_1^{\left({\alpha }_1\right)}+{\left(\frac{2}{\sqrt{{\lambda }_1}}\right)}^{\left({\alpha }_1\right)}+{\mathrm{\gamma }}_{tn}^2\right)b{J}_{0}d\xi d\eta$ (120)

$q_i=b{J}_0\left(\frac{2{\mu }_1}{{\alpha }_1}\left({\alpha }_1{\lambda }_1^{\left({\alpha }_1-1\right)}\frac{\partial {\lambda }_1}{\partial X_i}-{\alpha }_1\frac{\partial {\lambda }_1}{\partial X_i}{\lambda }_1^{\left(\frac{-{\alpha }_1-2}{2}\right)}+2{\mathrm{\gamma }}_{tn}\frac{\partial {\mathrm{\gamma }}_{tn}}{\partial X_i}\right)\right){\omega }_1{\omega }_2$ (121)

Furthermore, by applying the gradient to the internal force vector in Equation 119, one can find the Hessian matrix $K_{ij}$ from Equation 122:

$K_{ij}=\nabla q_i=\frac{\partial }{\partial X_j}\left[\frac{\partial }{\partial X_i}\left[{\int }_{V_0}\frac{2{\mu }_1}{{\alpha }_1}\left({\lambda }_1^{\left({\alpha }_1\right)}+{\left(\frac{2}{\sqrt{{\lambda }_1}}\right)}^{\left({\alpha }_1\right)}+{\mathrm{\gamma }}_{tn}^2\right)d{V}_0\right]\right]$ (122)

$K_{ij}=\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}\frac{2{\mu }_1}{{\alpha }_1}\frac{\partial }{\partial X_j}\frac{\partial }{\partial X_i}\left({\lambda }_1^{\left({\alpha }_1\right)}+{\left(\frac{2}{\sqrt{{\lambda }_1}}\right)}^{\left({\alpha }_1\right)}+{\mathrm{\gamma }}_{tn}^2\right)b{J}_{0}d\xi d\eta$ (123)

$K_{ij}=b{J}_0\left(\frac{2{\mu }_1}{{\alpha }_1}\left({\alpha }_1\left(\left({\alpha }_1-1\right)\frac{\partial {\lambda }_1^2}{\partial X_{ij}}{\lambda }_1^{\left({\alpha }_1-2\right)}+\frac{{\partial }^2{\lambda }_1}{\partial X_{ij}^2}{\lambda }_1^{\left({\alpha }_1-1\right)}\right)-{\alpha }_1\left(\frac{\left({-\alpha }_1-2\right)\frac{\partial {\lambda }_1^2}{\partial X_{ij}}{\lambda }_1^{\frac{\left(-{\alpha }_1-4\right)}{2}}}{2}+\frac{{\partial }^2{\lambda }_1}{\partial X_{ij}^2}{\lambda }_1^{\frac{\left(-{\alpha }_1-2\right)}{2}}\right)+2\left(\frac{\partial {\mathrm{\gamma }}_{tn}^2}{\partial X_{ij}}+{\mathrm{\gamma }}_{tn}\frac{{\partial }^2{\mathrm{\gamma }}_{tn}}{\partial X_{ij}^2}\right)\right)\right){\omega }_1{\omega }_2$ (124)

It is worth noting that for $N=2$, $N=3$, and $N=4$, the values ​​of the internal force vector $q_i$ and Hessian matrix $K_{ij}$ can be obtained by adding cumulative terms of the model for $N=1$, differentiating only the constitutive coefficients $\mu$ and $\alpha$.

4.5.5 Yeoh Model

Upon substituting the specific free energy of the Yeoh model as shown in Equation 84 into Equation 87, one has from Equation 125:

$U_e={\int }_{V_0}{C}_{10}\left({\lambda }_1^2+\frac{2}{{\lambda }_1}-3\right)+C_{20}{\left({\lambda }_1^2+\frac{2}{{\lambda }_1}-3\right)}^2+C_{30}{\left({\lambda }_1^2+\frac{2}{{\lambda }_1}-3\right)}^{3}d{V}_0$ (125)

However, upon conducting preliminary tests of implementing this formulation, it was noticed that the equation for the strain invariants from the uniaxial tensile test is not sufficient for the approach of plane frames with Reissner kinematics. Therefore, the first strain invariant from the simple shear test was added so that distortion is considered in the deformation energy of the model.

Therefore, one has the following equations starting from Equation 126:

$U_e={\int }_{V_0}{C}_{10}\left(I_{1_{TU}}+I_{1_{CS}}-3\right)+C_{20}{\left(I_{1_{TU}}+I_{1_{CS}}-3\right)}^2+C_{30}{\left(I_{1_{TU}}+I_{1_{CS}}-3\right)}^{3}d{V}_0$ (126)

$U_e={\int }_{V_0}{C}_{10}\left(({\lambda }_1^2+\frac{2}{{\lambda }_1})+{(\mathrm{\gamma }}_{tn}^2+3)-3\right)+C_{20}{\left(({\lambda }_1^2+\frac{2}{{\lambda }_1})+{(\mathrm{\gamma }}_{tn}^2+3)-3\right)}^2+C_{30}{\left(({\lambda }_1^2+\frac{2}{{\lambda }_1})+{(\mathrm{\gamma }}_{tn}^2+3)-3\right)}^{3}d{V}_0$ (127)

$U_e={\int }_{V_0}{C}_{10}\left({\lambda }_1^2+\frac{2}{{\lambda }_1}+{\mathrm{\gamma }}_{tn}^2\right)+C_{20}{\left({\lambda }_1^2+\frac{2}{{\lambda }_1}+{\mathrm{\gamma }}_{tn}^2\right)}^2+C_{30}{\left({\lambda }_1^2+\frac{2}{{\lambda }_1}+{\mathrm{\gamma }}_{tn}^2\right)}^{3}d{V}_0$ (128)

Thus, the internal force vector $q_i$ is obtained by the gradient of the elastic potential energy represented from Equation 129:

$q_i=\nabla U_e=\frac{\partial U_e}{\partial X_i}=\frac{\partial }{\partial X_i}\left[{\int }_{V_0}{C}_{10}\left({\lambda }_1^2+\frac{2}{{\lambda }_1}+{\mathrm{\gamma }}_{tn}^2\right)+C_{20}{\left({\lambda }_1^2+\frac{2}{{\lambda }_1}+{\mathrm{\gamma }}_{tn}^2\right)}^2+C_{30}{\left({\lambda }_1^2+\frac{2}{{\lambda }_1}+{\mathrm{\gamma }}_{tn}^2\right)}^{3}d{V}_0\right]$ (129)

$q_i=\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}\left(C_{10}\frac{\partial }{\partial X_i}\left({\lambda }_1^2+\frac{2}{{\lambda }_1}+{\mathrm{\gamma }}_{tn}^2\right)+C_{20}\frac{\partial }{\partial X_i}{\left({\lambda }_1^2+\frac{2}{{\lambda }_1}+{\mathrm{\gamma }}_{tn}^2\right)}^2+C_{30}\frac{\partial }{\partial X_i}{\left({\lambda }_1^2+\frac{2}{{\lambda }_1}+{\mathrm{\gamma }}_{tn}^2\right)}^3\right)b{J}_{0}d\xi d\eta$ (130)

$q_i=b{J}_0\left(C_{10}\left(2{\lambda }_1\frac{\partial {\lambda }_1}{\partial X_i}-\frac{2\frac{\partial {\lambda }_1}{\partial X_i}}{{\lambda }_1^2}+2{\mathrm{\gamma }}_{tn}\frac{\partial {\mathrm{\gamma }}_{tn}}{\partial X_i}\right)+C_{20}\left(4{\lambda }_1^3\frac{\partial {\lambda }_1}{\partial X_i}+4\frac{\partial {\lambda }_1}{\partial X_i}+4{\lambda }_1^2{\mathrm{\gamma }}_{tn}\frac{\partial {\mathrm{\gamma }}_{tn}}{\partial X_i}-\frac{8\frac{\partial {\lambda }_1}{\partial X_i}}{{\lambda }_1^3}+\frac{8{\mathrm{\gamma }}_{tn}\frac{\partial {\mathrm{\gamma }}_{tn}}{\partial X_i}}{{\lambda }_1}+4{\mathrm{\gamma }}_{tn}^2{\lambda }_1\frac{\partial {\lambda }_1}{\partial X_i}-\frac{4{\mathrm{\gamma }}_{tn}^2\frac{\partial {\lambda }_1}{\partial X_i}}{{\lambda }_1^2}+4{\mathrm{\gamma }}_{tn}^3\frac{\partial {\mathrm{\gamma }}_{tn}}{\partial X_i}\right)+C_{30}\left(\frac{24{\mathrm{\gamma }}_{tn}\left({\lambda }_1\frac{\partial {\mathrm{\gamma }}_{tn}}{\partial X_i}-{\mathrm{\gamma }}_{tn}\frac{\partial {\lambda }_1}{\partial X_i}\right)}{{\lambda }_1^3}+\frac{6\left(-{\mathrm{\gamma }}_{tn}^4\frac{\partial {\lambda }_1}{\partial X_i}+4{\mathrm{\gamma }}_{tn}^3{\lambda }_1\frac{\partial {\mathrm{\gamma }}_{tn}}{\partial X_i}\right)}{{\lambda }_1^2}+6{\lambda }_1^5\frac{\partial {\lambda }_1}{\partial X_i}+6{\mathrm{\gamma }}_{tn}^5\frac{\partial {\mathrm{\gamma }}_{tn}}{\partial X_i}+18{\lambda }_1^2\frac{\partial {\lambda }_1}{\partial X_i}+12\left({\mathrm{\gamma }}_{tn}^2\frac{\partial {\lambda }_1}{\partial X_i}+2{\mathrm{\gamma }}_{tn}{\lambda }_1\frac{\partial {\mathrm{\gamma }}_{tn}}{\partial X_i}\right)+3\left(2{\lambda }_1^4{\mathrm{\gamma }}_{tn}\frac{\partial {\mathrm{\gamma }}_{tn}}{\partial X_i}+4{\lambda }_1^3{\mathrm{\gamma }}_{tn}^2\frac{\partial {\lambda }_1}{\partial X_i}\right)+3\left(4{\mathrm{\gamma }}_{tn}^3{\lambda }_1^2\frac{\partial {\mathrm{\gamma }}_{tn}}{\partial X_i}+2{\mathrm{\gamma }}_{tn}^4{\lambda }_1\frac{\partial {\lambda }_1}{\partial X_i}\right)-\frac{24\frac{\partial {\lambda }_1}{\partial X_i}}{{\lambda }_1^4}\right)\right){\omega }_1{\omega }_2$ (131)

Furthermore, by applying the gradient to the internal force vector in Equation 129, the Hessian matrix $K_{ij}$ can be obtained from Equation 132:

$K_{ij}=\nabla q_i=\frac{\partial }{\partial X_j}\left[\frac{\partial }{\partial X_i}\left[{\int }_{V_0}{C}_{10}\left({\lambda }_1^2+\frac{2}{{\lambda }_1}+{\mathrm{\gamma }}_{tn}^2\right)+C_{20}{\left({\lambda }_1^2+\frac{2}{{\lambda }_1}+{\mathrm{\gamma }}_{tn}^2\right)}^2+C_{30}{\left({\lambda }_1^2+\frac{2}{{\lambda }_1}+{\mathrm{\gamma }}_{tn}^2\right)}^{3}d{V}_0\right]\right]$ (132)

$K_{ij}=\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}\left(C_{10}\frac{\partial }{\partial X_j}\frac{\partial }{\partial X_i}\left({\lambda }_1^2+\frac{2}{{\lambda }_1}+{\mathrm{\gamma }}_{tn}^2\right)+C_{20}\frac{\partial }{\partial X_j}\frac{\partial }{\partial X_i}{\left({\lambda }_1^2+\frac{2}{{\lambda }_1}+{\mathrm{\gamma }}_{tn}^2\right)}^2+C_{30}\frac{\partial }{\partial X_j}\frac{\partial }{\partial X_i}{\left({\lambda }_1^2+\frac{2}{{\lambda }_1}+{\mathrm{\gamma }}_{tn}^2\right)}^3\right)b{J}_{0}d\xi d\eta$ (133)

$K_{ij}=b{J}_0\left(C_{10}\left(2\left(\frac{\partial {\lambda }_1^2}{\partial X_{ij}}+{\lambda }_1\frac{{\partial }^2{\lambda }_1}{\partial X_{ij}^2}\right)-\frac{2\left({\lambda }_1\frac{{\partial }^2{\lambda }_1}{\partial X_{ij}^2}-2\frac{\partial {\lambda }_1^2}{\partial X_{ij}}\right)}{{\lambda }_1^3}+2\left(\frac{\partial {\mathrm{\gamma }}_{tn}^2}{\partial X_{ij}}+{\mathrm{\gamma }}_{tn}\frac{{\partial }^2{\mathrm{\gamma }}_{tn}}{\partial X_{ij}^2}\right)\right)+C_{20}\left(4\left(3{\lambda }_1^2\frac{\partial {\lambda }_1^2}{\partial X_{ij}}+{\lambda }_1^3\frac{{\partial }^2{\lambda }_1}{\partial X_{ij}^2}\right)+4\frac{{\partial }^2{\lambda }_1}{\partial X_{ij}^2}+4\left(2{\lambda }_1\frac{\partial {\lambda }_1}{\partial X_{ij}}{\mathrm{\gamma }}_{tn}\frac{\partial {\mathrm{\gamma }}_{tn}}{\partial X_{ij}}+\left(\frac{\partial {\mathrm{\gamma }}_{tn}^2}{\partial X_{ij}}+{\mathrm{\gamma }}_{tn}\frac{{\partial }^2{\mathrm{\gamma }}_{tn}}{\partial X_{ij}^2}\right){\lambda }_1^2\right)-\frac{8\left({\lambda }_1\frac{{\partial }^2{\lambda }_1}{\partial X_{ij}^2}-3\frac{\partial {\lambda }_1^2}{\partial X_{ij}}\right)}{{\lambda }_1^4}+\frac{8\left(\left(\frac{\partial {\mathrm{\gamma }}_{tn}^2}{\partial X_{ij}}+{\mathrm{\gamma }}_{tn}\frac{{\partial }^2{\mathrm{\gamma }}_{tn}}{\partial X_{ij}^2}\right){\lambda }_1-\frac{\partial {\lambda }_1}{\partial X_{ij}}{\mathrm{\gamma }}_{tn}\frac{\partial {\mathrm{\gamma }}_{tn}}{\partial X_{ij}}\right)}{{\lambda }_1^2}+4\left(2{\mathrm{\gamma }}_{tn}\frac{\partial {\mathrm{\gamma }}_{tn}}{\partial X_{ij}}{\lambda }_1\frac{\partial {\lambda }_1}{\partial X_{ij}}+\left(\frac{\partial {\lambda }_1^2}{\partial X_{ij}}+{\lambda }_1\frac{{\partial }^2{\lambda }_1}{\partial X_{ij}^2}\right){\mathrm{\gamma }}_{tn}^2\right)-\frac{4{\mathrm{\gamma }}_{tn}\left(\left(\frac{{\partial }^2{\lambda }_1}{\partial X_{ij}^2}{\mathrm{\gamma }}_{tn}+2\frac{\partial {\mathrm{\gamma }}_{tn}}{\partial X_{ij}}\frac{\partial {\lambda }_1}{\partial X_{ij}}\right){\lambda }_1-2{\mathrm{\gamma }}_{tn}\frac{\partial {\lambda }_1^2}{\partial X_{ij}}\right)}{{\lambda }_1^3}+4\left(3{\mathrm{\gamma }}_{tn}^2\frac{\partial {\mathrm{\gamma }}_{tn}^2}{\partial X_{ij}}+{\mathrm{\gamma }}_{tn}^3\frac{{\partial }^2{\mathrm{\gamma }}_{tn}}{\partial X_{ij}^2}\right)\right)+C_{30}\left(\frac{24\left({\mathrm{\gamma }}_{tn}\frac{{\partial }^2{\mathrm{\gamma }}_{tn}}{\partial X_{ij}^2}{\lambda }_1^2-{\lambda }_1\frac{{\partial }^2{\lambda }_1}{\partial X_{ij}^2}{\mathrm{\gamma }}_{tn}^2+\frac{\partial {\mathrm{\gamma }}_{tn}^2}{\partial X_{ij}}{\lambda }_1^2-4\frac{\partial {\lambda }_1}{\partial X_{ij}}{\mathrm{\gamma }}_{tn}\frac{\partial {\mathrm{\gamma }}_{tn}}{\partial X_{ij}}{\lambda }_1+3\frac{\partial {\lambda }_1^2}{\partial X_{ij}}{\mathrm{\gamma }}_{tn}^2\right)}{{\lambda }_1^4}+\frac{6{\mathrm{\gamma }}_{tn}^2\left(\left(4\frac{{\partial }^2{\mathrm{\gamma }}_{tn}}{\partial X_{ij}^2}{\mathrm{\gamma }}_{tn}{\lambda }_1-{\mathrm{\gamma }}_{tn}^2\frac{{\partial }^2{\lambda }_1}{\partial X_{ij}^2}+12\frac{\partial {\mathrm{\gamma }}_{tn}^2}{\partial X_{ij}}{\lambda }_1\right){\lambda }_1-2\frac{\partial {\lambda }_1}{\partial X_{ij}}{\mathrm{\gamma }}_{tn}\left(4\frac{\partial {\mathrm{\gamma }}_{tn}}{\partial X_{ij}}{\lambda }_1-\frac{\partial {\lambda }_1}{\partial X_{ij}}{\mathrm{\gamma }}_{tn}\right)\right)}{{\lambda }_1^3}+6\left(5{\lambda }_1^4\frac{\partial {\lambda }_1^2}{\partial X_{ij}}+\frac{{\partial }^2{\lambda }_1}{\partial X_{ij}^2}{\lambda }_1^5\right)+6\left(5{\mathrm{\gamma }}_{tn}^4\frac{\partial {\mathrm{\gamma }}_{tn}^2}{\partial X_{ij}}+\frac{{\partial }^2{\mathrm{\gamma }}_{tn}}{\partial X_{ij}^2}{\mathrm{\gamma }}_{tn}^5\right)+18\left(2{\lambda }_1\frac{\partial {\lambda }_1^2}{\partial X_{ij}}+\frac{{\partial }^2{\lambda }_1}{\partial X_{ij}^2}{\lambda }_1^2\right)+12\left({\mathrm{\gamma }}_{tn}^2\frac{{\partial }^2{\lambda }_1}{\partial X_{ij}^2}+2{\mathrm{\gamma }}_{tn}{\lambda }_1\frac{{\partial }^2{\mathrm{\gamma }}_{tn}}{\partial X_{ij}^2}+4{\mathrm{\gamma }}_{tn}\frac{\partial {\lambda }_1}{\partial X_{ij}}\frac{\partial {\mathrm{\gamma }}_{tn}}{\partial X_{ij}}+2{\lambda }_1\frac{\partial {\mathrm{\gamma }}_{tn}^2}{\partial X_{ij}}\right)+3\left(2{\lambda }_1^4\frac{\partial {\mathrm{\gamma }}_{tn}^2}{\partial X_{ij}}+2{\lambda }_1^4{\mathrm{\gamma }}_{tn}\frac{{\partial }^2{\mathrm{\gamma }}_{tn}}{\partial X_{ij}^2}+4{\lambda }_1^3{\mathrm{\gamma }}_{tn}^2\frac{{\partial }^2{\lambda }_1}{\partial X_{ij}^2}+16{\lambda }_1^3{\mathrm{\gamma }}_{tn}\frac{\partial {\mathrm{\gamma }}_{tn}}{\partial X_{ij}}\frac{\partial {\lambda }_1}{\partial X_{ij}}+12{\lambda }_1^2{\mathrm{\gamma }}_{tn}^2\frac{\partial {\lambda }_1^2}{\partial X_{ij}}\right)+3\left(2{\mathrm{\gamma }}_{tn}^4\frac{\partial {\lambda }_1^2}{\partial X_{ij}}+2{\mathrm{\gamma }}_{tn}^4{\lambda }_1\frac{{\partial }^2{\lambda }_1}{\partial X_{ij}^2}+4{\mathrm{\gamma }}_{tn}^3{\lambda }_1^2\frac{{\partial }^2{\mathrm{\gamma }}_{tn}}{\partial X_{ij}^2}+16{\mathrm{\gamma }}_{tn}^3{\lambda }_1\frac{\partial {\mathrm{\gamma }}_{tn}}{\partial X_{ij}}\frac{\partial {\lambda }_1}{\partial X_{ij}}+12{\mathrm{\gamma }}_{tn}^2{\lambda }_1^2\frac{\partial {\mathrm{\gamma }}_{tn}^2}{\partial X_{ij}}\right)-\frac{24\left({\lambda }_1\frac{{\partial }^2{\lambda }_1}{\partial X_{ij}^2}-4\frac{\partial {\lambda }_1^2}{\partial X_{ij}}\right)}{{\lambda }_1^5}\right)\right){\omega }_1{\omega }_2$ (134)

5 RESULTS AND DISCUSSIONS

In this chapter, numerical examples were developed based on the formulation discussed in the previous chapters, employing the positional finite element method for plane frames using the Reissner kinematics. Subsequently, comparisons were made with theoretical results and experimental data available in the structural engineering literature, aiming to validate the developed software.

5.1 Uniaxial traction in hyperelastic bar

Based on the work of Pascon (2008) and Perônica et al. (2022), Figure 6 depicts a hyperelastic bar with a second-kind support at node 1 and a first-kind support at node 3, under uniaxial traction with horizontal force applied from left to right at node 3.

Figure 6
Hyperelastic bar under uniaxial traction - Adapted from Perônica et al. (2022)

As observed in Figure 6, the length is L = 10 cm, the area is A = 1.0 cm², and the inertia is I_z = 8.33 x 10^-2 cm⁴. Additionally, 50 load steps were used with 1 finite element using quadratic approximation and a convergence tolerance in positions of 10^-6.

The purpose of this numerical example is to compare and validate the results obtained from the developed software with the formulations described in the previous chapters for hyperelastic models, using experimental data from the study by Paula et al. (2019). Paula et al. (2019) conducted uniaxial tensile tests on elastomer samples with different percentages of carbon black, namely 0%, 10%, 30%, and 60%, and concluded that the higher the percentage of carbon black added to the elastomers, the greater the stiffness of the samples.

Furthermore, by comparing the developed formulation with these tests, it is possible to adjust the calibration of the constitutive material constants using the least squares method. This improves the accuracy of the models and allows their use in subsequent numerical examples. As mentioned earlier, this research analyzed four hyperelastic models: Neo-Hookean (Mod. NH), Mooney-Rivlin (Mod. MR), Yeoh (Mod. YE), and Ogden (Mod. OG), in addition to the linear elastic model (Mod. EL).

Therefore, the least squares method can be applied to determine the material constitutive coefficients by defining an objective function that quantifies the difference between the strain energy density function and the stress vs. strain data observed experimentally by Paula et al. (2019). With this, based on the equations obtained in the previous sections, it is possible to determine the constitutive coefficients for points (${\epsilon }_i$,${\sigma }_i$). In this way, the sum of the squares of the differences between the experimental and estimated points is minimized, resulting in a linear system where the unknowns are the constitutive coefficients of the hyperelastic models.

From obtaining the constitutive constants for the different carbon black samples, it was possible to use the formulations of the hyperelastic models, develop the computational code in FORTRAN, and obtain the results of the displacements of the bar considering 0%, 10%, 30%, and 60% of carbon black, as per the experiment conducted by Paula et al. (2019), illustrated in Figure 7:

Figure 7
Force vs. displacement plot for elastomer with carbon black additions: a) 0% b) 10% c) 30% d) 60%

It is worth noting that for the sample with 0% carbon black, only the OG and YE models achieved a correlation coefficient R² greater than 0.99, while the others obtained R² slightly above 0.90. Additionally, the NH and MR models reached such R² after reducing points from the stress-strain curve. The OG and YE models were able to capture all 25 points from the experimental curves of Paula et al. (2019), and the NH and MR models captured 18 and 20 points respectively, making the models coherent only in the initial portion of the curve where there is predominance of a linear behavior. Such behavior was expected given the formulation with fewer constitutive constants of the NH and MR models.

For the sample with 10% carbon black, only the OG and YE models achieved an R² greater than 0.99, while the others obtained an R² slightly above 0.90. Additionally, the NH and MR models reached such R² after reducing points from the stress-strain curve. The OG and YE models were able to capture all 25 points from the experimental curves of Paula et al. (2019), and the NH and MR models captured 14 and 19 points respectively, making the models coherent only in the initial portion of the curve, similar to the case with 0% carbon black.

For the sample with 30% carbon black, only the OG and YE models achieved an R² greater than 0.99, while the others obtained an R² slightly above 0.90. Additionally, the NH model reached such R² after reducing points from the stress-strain curve. The OG, YE, and MR models were able to capture all 21 points from the experimental curves of Paula et al. (2019), and the NH model captured only 9 points, making the model coherent only in the initial portion of the curve where there is predominance of linear behavior. Although the MR model captured all points of the curve, the NH and MR models continued to exhibit the expected linear behavior due to the formulation with fewer constitutive constants. Furthermore, it can be observed that with the increase in the percentage to 30% carbon black, and consequently, an increase in the stiffness of the elastomer, the OG and YE models maintain high accuracy, but their curves show a subtle deviation from the curves of the NH, MR, and EL models.

For the sample with 60% carbon black, only the OG, YE, and MR models achieved an R² greater than 0.99, while the NH model obtained an R² slightly above 0.90. All models were able to capture the 17 points from the experimental curves of Paula et al. (2019), and following the subtle trend observed in the previous case, it can be noted that with the increase in the percentage to 60% carbon black, and consequently, an increase in the stiffness of the elastomer, the OG and YE models remain the models with the highest observed accuracy. However, the curves of these models show a tendency to approximate the curves of the NH, MR, and EL models, demonstrating a trend towards linearization with the increase in material stiffness.

Such results for elastomers with carbon black additions, obtained through the developed computational code, were consistent with the experimental study by Paula et al. (2019) and with the computational code development study for truss analysis by Perônica et al. (2022). However, it can be observed that the most accurate models in all elastomer samples regardless of stiffness were the Ogden and Yeoh hyperelastic models, which directly corroborates with the aforementioned literature.

5.2 Arch plane frame instability

Based on the work of Maciel (2008), Figure 8 depicts a plane frame in an arch configuration, fixed at one end and supported at the other, with a transverse load applied from top to bottom along the axis of symmetry of the arch.

Figure 8
Arch plane frame - Adapted from Maciel (2008)

As observed in Figure 8, one has a radius R = 100, EA = 100EI, EI = 10⁶, and angle α = 145°. Additionally, we used 50 load steps and 20 finite elements with cubic approximation and convergence tolerance of 10^-6. For the validation of the computational code in plane frames, Figure 9 presents the comparison between longitudinal, transverse displacements, and rotation from the study by Maciel (2008), which, in turn, corroborated with the analytical solutions of critical load from the studies by DaDeppo and Schmidt (1975) and Ibrahimbegovic (1995).

Figure 9
Force vs. displacement and rotation graphs for linear elastic model

It is worth highlighting the high accuracy of the developed computational code with the study by Maciel (2008) and the aforementioned literature, considering that the maximum critical load convergence value of Maciel (2008), Ibrahimbegovic (1995), and the present study was 897.3 when using the Reissner kinematics, as well as that of the study by DaDeppo and Schmidt (1975) was 897 with Euler-Bernoulli kinematics. Such validation is necessary to evaluate the same numerical example for the hyperelastic models. Therefore, for the analysis of the hyperelastic models with the elastomer samples from the study by Paula et al. (2019), the values ​​of the constitutive constants obtained for 0%, 10%, 30%, and 60% carbon black were used respectively.

With that, one has in Figure 10 the graphs of F vs. δx, F vs. δy, and F vs. ω for the circular arch elastomer plane frame with 0% carbon black:

Figure 10
Force vs. displacement and rotation graphs for elastomer with 0% carbon black.

For the sample with 0% carbon black, the differences between the hyperelastic models are accentuated, especially with the OG model compared to the others, given that the MR, NH, and YE models showed similar results in longitudinal, transverse, and rotational displacements. This difference partially corroborates with the results of the previous numerical examples, considering that only the OG model showed discrepancies, whereas in the previous numerical examples, the NH and MR models were the most discrepant due to the lower number of constants and consequently, the lower complexity of the formulations, making the models less accurate.

With that, one has in Figure 11 the graphs of F vs. δx, F vs. δy, and F vs. ω for the circular arch plane frame made of elastomer with 10% carbon black.

Figure 11
Force vs. displacement and rotation graphs for elastomer with 10% carbon black.

For the sample with 10% carbon black, the differences between the hyperelastic models are pronounced, especially between the models with greater differences in the number of constants and formulation complexity. The OG and YE models showed similar results to each other, while the NH and MR models also showed similar results for longitudinal displacement, transverse displacement, and rotation.

With that, one has in Figure 12 the graphs of F vs. δx, F vs. δy, and F vs. ω for the circular arch plane frame of elastomer with 30% carbon black.

Figure 12
Force vs. displacement and rotation graphs for elastomer with 30% carbon black.

For the sample with 30% carbon black, the differences between the hyperelastic models are accentuated, mainly between the MR model and the others. The OG, NH, and YE models showed similar results in terms of longitudinal, transverse displacement, and rotation. This difference partially corroborates with the results of previous numerical examples, considering that only the MR model showed a significant difference, whereas in previous numerical examples, the NH and MR models were the most discrepant due to the lower number of constants and consequently, the lower complexity of the formulations, making the models less accurate.

With that, one has in Figure 13 the graphs of F vs. δx, F vs. δy, and F vs. ω for the circular arch elastomer plane frame with 60% carbon black.

Figure 13
Force vs. displacement and rotation graphs for elastomer with 60% carbon black.

For the 60% carbon black sample, the trend of results directly corroborated with the analysis of the results for the 10% carbon black sample.

5.3 Articulated diamond-shaped plane frame

From the work of Maciel (2008), one has in Figure 14 the representation of an articulated diamond-shaped plane frame with applied symmetry.

Figure 14
Articulated diamond-shaped plane frame - Adapted from Maciel (2008)

As observed in Figure 14, the length is L=1, E=1, I=1, and area A=10000. Additionally, 50 load steps were used along with 20 finite elements with quadratic approximation and a convergence tolerance of 10⁻⁶. It is important to highlight that whenever the nomenclatures "(T)" and "(C)" appear in this section, they refer to tension and compression forces in the structure, respectively.

For the validation of the computational code in plane frames, a comparison is presented in Figure 15 between the longitudinal and transverse displacements, as well as the rotation for traction and compression in the plane frame from the study by Maciel (2008), which in turn corroborated with the analytical solutions from the studies by Jenkins et al. (1966) and Mattiasson (1981).

Figure 15
Force vs. displacement and rotation graphs for linear elastic model

It is worth noting the high accuracy of the developed computational code with the study by Maciel (2008) and the aforementioned literature. Such validation is necessary to evaluate the same numerical example for hyperelastic models. Therefore, for the analysis of hyperelastic models with elastomer samples from the study by Paula et al. (2019), a load value 10⁶ times higher than the load in the study by Maciel (2008) and a cross-sectional area 10 times smaller were used to allow the applicability of the obtained constitutive constant values for 0%, 10%, 30%, and 60% carbon black content, respectively.

With that, one has in Figure 16 the graphs of F vs. δx, F vs. δy, and F vs. ω for the articulated diamond-shaped plane frame made of elastomer with 0% carbon black:

Figure 16
Force vs. displacement and rotation graphs for elastomer with 0% carbon black.

For the sample with 0% carbon black, the differences between the hyperelastic models are very subtle both in traction and compression. In traction, the largest difference between the models was 4.35%, while in compression, it was 6.19%. Additionally, the smallest difference in traction was 0.43%, while in compression, it was 0.87%. This difference in accuracy between traction and compression results corroborates with the studies of Marczak and Iturrioz (2006), as well as with the studies of Perônica et al. (2022), which state that although compression follows the same theoretical equations in hyperelastic models, the difference arises because the experimental adjustment is made based on an experimental curve for uniaxial traction, thus, the compression behavior is predicted and not adjusted.

With that, one has in Figure 17 the graphs of F vs. δx, F vs. δy, and F vs. ω for the articulated diamond-shaped elastomer plane frame with 10% carbon black.

Figure 17
Force vs. displacement and rotation graphs for elastomer with 10% carbon black.

For the sample with 10% carbon black, the greatest difference between the hyperelastic models in traction and compression occurs between the OG and MR models, while the smallest difference in traction and compression occurs between the YE and OG models. The differences between the hyperelastic models are more pronounced in compression.

With this, one has in Figure 18 the graphs of F vs. δx, F vs. δy, and F vs. ω for the articulated diamond-shaped elastomer plane frame with 30% carbon black.

Figure 18
Force vs. displacement and rotation graphs for elastomer with 30% carbon black.

For the 30% carbon black sample, the largest difference between the hyperelastic models in traction and compression occurs between the YE and MR models, while the smallest difference in traction and compression occurs between the NH and OG models. The differences between the hyperelastic models are pronounced, both in traction and compression, mainly from the MR model compared to the others, considering that the OG, NH, and YE models showed similar results in longitudinal, transverse displacement, and rotation.

With that, one has in Figure 19 the graphs of F vs. δx, F vs. δy, and F vs. ω for the articulated diamond-shaped elastomer plane frame with 60% carbon black.

Figure 19
Force vs. displacements and rotation graphs for elastomer with 60% carbon black.

For the sample with 60% carbon black, the greatest difference between the hyperelastic models in traction and compression occurs between the OG and NH models, while the smallest difference in traction and compression occurs between the YE and OG models. The differences between the hyperelastic models are pronounced in both compression and traction, especially among the models with the greatest difference in the number of constants and formulation complexity. The OG and YE models showed similar results, while the NH and MR models also exhibited similar longitudinal, transverse displacement, and rotation.

6 CONCLUSION

Finally, it was successfully developed a FORTRAN computational code for the nonlinear analysis of plane frames considering hyperelastic models using the positional finite element method combined with the use of Reissner kinematics. It is worth noting the difficulty of directly comparing the developed formulation for plane frames using hyperelastic models, given that there aren’t numerical or experimental studies in the literature for plane frames with hyperelastic models.

From this, it is worth highlighting the importance of experimental tests and appropriate stress vs. strain curves in the study and accuracy of hyperelastic models, given that the formulation of these models depends directly on the adjustment of the constitutive constants from experimental test data. Such statement was evidenced in the numerical example where there was an analysis of traction and compression behavior in the plane frame, where for traction, the accuracy of the computational code was much more evident than in compression, starting from the fact that the models were adjusted with experimental data from the uniaxial traction study by Paula et al. (2019), and thus, the behavior in compression was predicted by theoretical equations, and the constitutive constants were not obtained for both traction and compression behaviors. This difference in accuracy between traction and compression results corroborates with the studies of Marczak and Iturrioz (2006), as well as with the studies of Perônica et al. (2022).

Regarding the accuracy of the hyperelastic models, it is noted that the Yeoh and Ogden models presented the most adjusted and coherent values, while the models that deviated the most were the Neo-Hookean and Mooney-Rivlin models. Such observation can be explained by the complexity and number of constants present in the formulation of these models, with the endorsement of studies by Paula et al. (2019) and Perônica et al. (2022), which also found higher precision in models with a greater number of constants. The differences between the hyperelastic models and the linear elastic model were significant but expected, considering that hyperelastic models can better capture the nonlinearity of the stress-strain curve of hyperelastic materials.

It can be stated that the use of the positional finite element method formulation for the nonlinear analysis of plane frames with the use of hyperelastic models performed well, especially after the modifications proposed by this work in the formulations of the hyperelastic models. These modifications consisted of adding the first and second invariants of deformation from the simple shear formulation to consider the distortion in the specific strain energy of the hyperelastic models using Reissner kinematics.

Furthermore, it was observed that there is a significant increase in material stiffness with the increase in the percentage of carbon black addition, which consequently caused a decrease in displacements in the numerical examples. Moreover, with the increase in the percentage of carbon black and the consequent decrease in displacements, there was no reduction in computational accuracy. On the contrary, with the increase in the percentage of carbon black and the consequent increase in material stiffness, there was a tendency for approximation and linearization of the behavior of the hyperelastic models with the linear elastic model.

In summary, the research analysis objectives were successfully achieved, considering that the study could contribute to various gaps in the nonlinear analysis of plane frames using hyperelastic materials with the positional finite element method. Furthermore, such a study was conducted without the need for extensive computational resources, demonstrating the efficiency of the computational code and the FORTRAN language.

References

• W. Yi Wang, J. Li, W. Liu e Z.-K. Liu, “Integrated computational materials engineering for advanced materials: A brief review”,Comput. Mater. Sci., vol. 158, pp. 42–48, fevereiro, 2019, https://doi.org/10.1016/j.commatsci.2018.11.001
  » https://doi.org/10.1016/j.commatsci.2018.11.001
• J. Navas, S. Bonnet, J. Voirin e G. Journaux, “Models as enablers of agility in complex systems engineering”,INCOSE Int. Symp., vol. 30, n.º 1, pp. 339–355, julho, 2020, https://doi.org/10.1002/j.2334-5837.2020.00726.x
  » https://doi.org/10.1002/j.2334-5837.2020.00726.x
• E. Samaniego et al., “An energy approach to the solution of partial differential equations in computational mechanics via machine learning: Concepts, implementation and applications”,Comput. Methods Appl. Mechanics Eng., vol. 362, p. 112790, abril, 2020, https://doi.org/10.1016/j.cma.2019.112790
  » https://doi.org/10.1016/j.cma.2019.112790
• S. David Müzel, E. P. Bonhin, N. M. Guimarães e E. S. Guidi, “Application of the Finite Element Method in the Analysis of Composite Materials: A Review”,Polymers, vol. 12, n.º 4, p. 818, abril, 2020, https://doi.org/10.3390/polym12040818
  » https://doi.org/10.3390/polym12040818
• M. Greco e D. H. N. Peixoto, “Comparative assessments of strain measures for nonlinear analysis of truss structures at large deformations”,Eng. Computations, novembro, 2021, https://doi.org/10.1108/ec-01-2021-0056
  » https://doi.org/10.1108/ec-01-2021-0056
• M. C. J. Reis e H. B. Coda, “Physical and geometrical non-linear analysis of plane frames considering elastoplastic semi-rigid connections by the positional FEM”,Latin Amer. J. Solids Struct., vol. 11, n.º 7, pp. 1163–1189, dezembro, 2014, https://doi.org/10.1590/s1679-78252014000700006
  » https://doi.org/10.1590/s1679-78252014000700006
• A. A. Basheer, “Advances in the smart materials applications in the aerospace industries”,Aircr. Eng. Aerosp. Technol., vol. 92, n.º 7, pp. 1027–1035, junho, 2020, https://doi.org/10.1108/aeat-02-2020-0040
  » https://doi.org/10.1108/aeat-02-2020-0040
• L. Jing e B. Jian, “Study on Factors Affecting Sealing Performance of Proton Exchange Membrane Fuel Cell End Plate”, in2019 Int. Conf. Advances Construction Machinery Vehicle Eng. (ICACMVE), Changsha, China, 2019-05-14–16. IEEE, 2019, https://doi.org/10.1109/icacmve.2019.00023
  » https://doi.org/10.1109/icacmve.2019.00023
• M. Motevalli, J. Uhlemann, N. Stranghöner e D. Balzani, “Geometrically nonlinear simulation of textile membrane structures based on orthotropic hyperelastic energy functions”,Composite Struct., vol. 223, p. 110908, setembro, 2019, https://doi.org/10.1016/j.compstruct.2019.110908
  » https://doi.org/10.1016/j.compstruct.2019.110908
• D. M. Taghizadeh e H. Darijani, “Mechanical Behavior Modeling of Hyperelastic Transversely Isotropic Materials Based on a New Polyconvex Strain Energy Function”,Int. J. Appl. Mechanics, vol. 10, n.º 09, p. 1850104, novembro, 2018, https://doi.org/10.1142/s1758825118501041
  » https://doi.org/10.1142/s1758825118501041
• S. Talebi e H. Darijani, “A pseudo-strain energy density function for mechanical behavior modeling of visco-hyperelastic materials”,Int. J. Mech. Sci., vol. 208, p. 106652, outubro, 2021, https://doi.org/10.1016/j.ijmecsci.2021.106652
  » https://doi.org/10.1016/j.ijmecsci.2021.106652
• H. Li, P. Xue, T. Zhang e B. Wen, “Nonlinear vibration study of fiber-reinforced composite thin plate with strain-dependent property based on strain energy density function method”,Mechanics Adv. Mater. Struct., vol. 27, n.º 9, pp. 761–773, Janeiro, 2019, https://doi.org/10.1080/15376494.2018.1495792
  » https://doi.org/10.1080/15376494.2018.1495792
• S. Dastjerdi, A. Alibakhshi, B. Akgöz e Ö. Civalek, “A Novel Nonlinear Elasticity Approach for Analysis of Nonlinear and Hyperelastic Structures”,Eng. Anal. with Boundary Elements, vol. 143, pp. 219–236, outubro, 2022, https://doi.org/10.1016/j.enganabound.2022.06.015
  » https://doi.org/10.1016/j.enganabound.2022.06.015
• M. Seguini e D. Nedjar, “Modelling of soil–structure interaction behaviour: geometric nonlinearity of buried structures combined to spatial variability of soil”,Eur. J. Environmental Civil Eng., vol. 21, n.º 10, pp. 1217–1236, fevereiro, 2016, https://doi.org/10.1080/19648189.2016.1153525
  » https://doi.org/10.1080/19648189.2016.1153525
• H. Darijani e R. Naghdabadi, “Hyperelastic materials behavior modeling using consistent strain energy density functions”,Acta Mech., vol. 213, n.º 3-4, pp. 235–254, Janeiro, 2010, https://doi.org/10.1007/s00707-009-0239-3
  » https://doi.org/10.1007/s00707-009-0239-3
• É. S. Ramos e R. Carrazedo, “Cross-section modeling of the non-uniform corrosion due to chloride ingress using the positional finite element method”,J. Brazilian Soc. Mech. Sci. Eng., vol. 42, n.º 10, setembro, 2020, https://doi.org/10.1007/s40430-020-02627-5
  » https://doi.org/10.1007/s40430-020-02627-5
• D. M. S. Paulino e E. D. Leonel, “Topology optimization and geometric nonlinear modeling using positional finite elements”,Optim. Eng., julho, 2021, https://doi.org/10.1007/s11081-021-09661-9
  » https://doi.org/10.1007/s11081-021-09661-9
• D. N. Maciel, “Análise de problemas elásticos não lineares geométricos empregando o método dos elementos finitos posicional”, Univ. Sao Paulo, 2008, http://www.teses.usp.br/teses/disponiveis/18/18134/tde-08052008-090039/
  » http://www.teses.usp.br/teses/disponiveis/18/18134/tde-08052008-090039
• D. S. Perônica, D. N. Maciel, R. Barros, J. A. d. Nascimento Neto e J. N. d. Silva Filho, “Análise não linear de treliças considerando modelos hiperelásticos pelo método dos elementos finitos posicional”,Res., Soc. Develop., vol. 11, n.º 10, Agosto, 2022, art. n.º e449111032684, https://doi.org/10.33448/rsd-v11i10.32684
  » https://doi.org/10.33448/rsd-v11i10.32684
• A. Anssari-Benam e C. O. Horgan, “New results in the theory of plane strain flexure of incompressible isotropic hyperelastic materials”,Proc. Roy. Soc. A: Math., Physical Eng. Sci., vol. 478, n.º 2258, fevereiro, 2022, https://doi.org/10.1098/rspa.2021.0773
  » https://doi.org/10.1098/rspa.2021.0773
• H. B. Khaniki, M. H. Ghayesh, R. Chin e M. Amabili, “A review on the nonlinear dynamics of hyperelastic structures”,Nonlinear Dyn., Agosto, 2022, https://doi.org/10.1007/s11071-022-07700-3
  » https://doi.org/10.1007/s11071-022-07700-3
• R. W. Ogden, “Non-linear elastic deformations”,Eng. Anal., vol. 1, n.º 2, p. 119, junho, 1984, https://doi.org/10.1016/0264-682x(84)90061-3
  » https://doi.org/10.1016/0264-682x(84)90061-3
• R. Jerábek e L. Écsi, “Numerical Study of Material Degradation of a Silicone Cross-Shaped Specimen Using a Thermodynamically Consistent Mooney-Rivlin Material Model”,Mater. Sci. Forum, vol. 952, pp. 258–266, abril, 2019, https://doi.org/10.4028/www.scientific.net/msf.952.258
  » https://doi.org/10.4028/www.scientific.net/msf.952.258
• M. Mooney, “A Theory of Large Elastic Deformation”,J. Appl. Phys., vol. 11, n.º 9, pp. 582–592, setembro, 1940, https://doi.org/10.1063/1.1712836
  » https://doi.org/10.1063/1.1712836
• R. W. Ogden “Large deformation isotropic elasticity – on the correlation of theory and experiment for incompressible rubberlike solids”,Proc. Roy. Soc. London. A. Math. Physical Sci., vol. 326, n.º 1567, pp. 565–584, fevereiro, 1972, https://doi.org/10.1098/rspa.1972.0026
  » https://doi.org/10.1098/rspa.1972.0026
• O. H. Yeoh, “Characterization of Elastic Properties of Carbon-Black-Filled Rubber Vulcanizates”,Rubber Chemistry Technol., vol. 63, n.º 5, pp. 792–805, novembro, 1990, https://doi.org/10.5254/1.3538289
  » https://doi.org/10.5254/1.3538289
• J. P. Pascon, “Modelos constitutivos para materiais hiperelásticos: estudo e implementação computacional”, Univ. Sao Paulo, 2008. http://www.teses.usp.br/teses/disponiveis/18/18134/tde-17042008-084851/
  » http://www.teses.usp.br/teses/disponiveis/18/18134/tde-17042008-084851/
• J. P. A. Paula, D. F. Lalo, Á. M. Almeida e M. Greco, “Comparative curve fitting through hyperelastic constitutive models for several formulations of carbonbrack-filled rubber vulcanizates”, inIbero-latin-amer. congr. comput. methods eng., Natal, Brasil, 2019-11-11–14. Belo Horizonte: ABMEC, 2019, pp. 1–18.
• D. A. DaDeppo e R. Schmidt, “Instability of Clamped-Hinged Circular Arches Subjected to a Point Load”,J. Appl. Mechanics, vol. 42, n.º 4, pp. 894–896, dezembro, 1975, https://doi.org/10.1115/1.3423734
  » https://doi.org/10.1115/1.3423734
• A. Ibrahimbegović, “On finite element implementation of geometrically nonlinear Reissner's beam theory: three-dimensional curved beam elements”,Comput. Methods Appl. Mechanics Eng., vol. 122, n.º 1-2, pp. 11–26, abril, 1995, https://doi.org/10.1016/0045-7825(95)00724-f
  » https://doi.org/10.1016/0045-7825(95)00724-f
• J. A. Jenkins, T. B. Seitz e J. S. Przemieniecki, “Large deflections of diamond-shaped frames”,Int. J. Solids Struct., vol. 2, n.º 4, pp. 591–603, outubro, 1966, https://doi.org/10.1016/0020-7683(66)90041-2
  » https://doi.org/10.1016/0020-7683(66)90041-2
• K. Mattiasson, “Numerical results from large deflection beam and frame problems analysed by means of elliptic integrals”,Int. J. Numer. Methods Eng., vol. 17, n.º 1, pp. 145–153, Janeiro, 1981, https://doi.org/10.1002/nme.1620170113
  » https://doi.org/10.1002/nme.1620170113
• R. J. Marczak e I. Iturrioz, in Caracterização do Comportamento de Materiais Hiperelásticos. SENAI – Dep. Reg. Rio Gd. Sul, 2006, pp. 1–137. https://www.senairs.org.br/sites/default/files/documents/caracterizacao-de-comportamento-de-materiais-hiperelasticos.pdf
  » https://www.senairs.org.br/sites/default/files/documents/caracterizacao-de-comportamento-de-materiais-hiperelasticos.pdf

Edited by

Editor:

Marco L. Bittencourt

Publication Dates

• Publication in this collection
  13 Sept 2024
• Date of issue
  2024

History

• Received
  15 Apr 2024
• Reviewed
  24 May 2024
• Accepted
  15 July 2024
• Published
  17 July 2024