An efficient model for large deformations of nonlinear viscoelastic elastomeric membranes using predictor-corrector algorithm

Abstract

Many engineering systems incorporate viscoelastic membranes of different geometries and boundary conditions experiencing large deformations. This paper presents a formulation based on the theory of nonlinear viscoelasticity. First, the kinematics of membrane deformation is expressed in three-dimensional space, and then the viscoelastic formulation for membranes is obtained based on the multiplicative decomposition of the tensor of deformation gradient. Also, the right Cauchy−Green viscoelastic tensor is considered as an internal variable. To solve the integration of evolution equation, a predictor-corrector method is used. Finally, due to the nonlinearity of the equations governing the problem, a nonlinear finite element formulation is derived. To check the effectiveness of the obtained formulation, several problems are studied. The comparisons show that the results of this formulation are in good agreement with the analytical and experimental results in the literature. It is shown that the current simplified viscoelastic model can successfully predict the results in the literature with more complicated viscoelastic models. Moreover, it is proven that the present model can predict the experimental data with just four material parameters, while the previous models should employ 12 material parameters. Therefore, the model presented in this paper is capable of predicting the experimental results more accurately with fewer material parameters.

Graphical abstract

Keywords

elastomers / rubber-like materials / non-linear viscoelasticity theory / membranes / finite element analysis

1 Introduction

Structures have crucial role in several industries [1,2]. Membranes have many applications in medical [3,4], mechanical [5,6], aerospace [7,8] and civil sciences [9,10]. In general, membranes have very high tensile strength, while their endurance to bending loads is negligible [11–13]. In several applications, the membranes are made of rubber-like materials which are subject to the large deformations [14,15]. The analysis of membrane deformations is essential in engineering, particularly in studying rubber membranes that can undergo significant deformations over time [16]. While there has been extensive research on the behavior of elastic membranes, there is limited information available on the large deformations in viscoelastic membranes [17]. Some studies have focused on the finite deformation of liquid capsules enclosed by elastic membranes under simple shear flow [18,19]. These studies have emphasized the importance of nonlinear effects resulting from finite deformation and the continuous elongation of capsules when exposed to shear stress. Furthermore, the application of membranes in bioengineering is evident in the research on the tympanic membrane, where a viscoelastic model has been proposed to model large deformations [20]. Finite elastic deformations of membranes have been subject of many investigations. The finite element (FE) method is also a powerful tool in investigating large deformations of membranes and has been used in some studies such as Ref. [21]. To study viscoelastic membranes, Feng and Huang [22] investigated the behavior of circular membrane based on Christensen’s viscoelastic model, which is the simplest viscoelastic model for large deformations. The investigation of large deformations in membranes serves as the foundation for studying the viscoelastic model applied to the tympanic membrane [23]. Additionally, the utilization of viscoelastic models has been extended to simulate two-phase flows involving a deformable body with a viscoelastic membrane [24]. There has been also a development of a robust method for fluid-structure interaction involving membranes experiencing large displacement and added-mass effects [25]. The viscoelastic model on the fractional order viscoelastic model for stress relaxation of polyvinyl chloride geomembranes has been discussed in Mierunalan et al. [26], and Zhang et al. [27], highlighting the importance of exploring advanced viscoelastic models to understand membrane behavior under deformation. Furthermore, there is other study on Boltzmann simulation of two-phase flows with a deformable body and viscoelastic membrane [28], emphasizing the significance of numerical methods in studying membrane behavior. van Chung [29] used exact geometry to form the exact description of the circular defining curve and FE shape function to approximate the defining curve and verified the model with elasticity problems.

Bhandarkar et al. [30] developed a theory for viscoelastic behavior of polyimide membranes at both room and deep cryogenic temperatures. This enabled a quantitative specification of the behavior of these films. In the field of viscoelastic materials, Ghosh et al. [31] conducted a study that introduced a mathematical formulation utilizing tensor algebra and fractional calculus to determine the three dimensional (3D) complex modulus of various viscoelastic models, including Kelvin–Voigt and Maxwell arrangements. Their research presented advanced formulations that can enhance the modeling of viscoelastic materials in FE analysis. Additionally, Geith et al. [32] performed experimental and mathematical characterization of polyamide-12 balloon catheter membranes, revealing their semi-compliant, and viscoelastic properties and providing insights into their behavior in practical applications. Propagation and attenuation of the Rayleigh and pseudo surface waves in two types of viscoelastic seismic metamaterials was studied by Cai et al. [33]. Firouzi [34] studied behavior of nonlinear visco-hyper-hysteretic membranes. Nonlinear vibration of viscoelastic plates using isogeometric analysis was investigated by Shafei et al. [35]. Zhou et al. [36] conducted research on modeling the viscoelastic-viscoplastic behavior of glassy polymer membranes, which may contribute to the understanding of non-uniform sub-chains in entangled networks and its potential impact on overall viscoelastic behavior in membranes. Brighenti et al. [37] studied failure of viscoelastic elastomers using phase field method.

Recently, Khaniki and colleagues [38,39] provided two comprehensive review articles on large dynamics of hyperelastic structures as well as hyperelasticity with application on mechanics and biomechanics. Amabili [40] presented a comprehensive study on large nonlinear deformations of shells, composites and soft tissues in his book. In another study, Amabili [41] investigated nonlinear damping from viscoelasticity on vibrations using third-order harmonic balance method and compared with experimental data on plates and shells. Khorasani et al. [42] studied thermoelastic buckling of honeycomb sandwich microplates using quasi three-dimensional hyperbolic shear deformation theory. An isogeometric analysis for of anisotropic variable angle composite plates was done by Shafei et al. [43]. Study on coupled hemispherical-conical-conical shell structures using first-order shear deformation theory was performed by Sobhani et al. [44]. Mercan et al. [45] inspected free vibration of functionally graded annular composite plates. Franchini et al. [46] proposed a quasi-static and dynamic mechanical viscoelastic model experiment on 15 human thoracic aortas.

Study through the literature shows that despite the extensive study of the behavior of elastic membranes, the analysis of large deformations of viscoelastic membranes is limited to very few articles and in special cases such as the axially symmetric state. In this paper, the hysteretic effect is ignored, and a simplified efficient formulation is developed to gain a comprehensive understanding of viscoelastic membrane behavior under large deformations. The research contributes to the existing knowledge gap by introducing an advanced nonlinear viscoelastic model and utilizing the FE method. By integrating these approaches, a thorough comprehension of membrane behavior is achieved, overcoming limitations imposed by geometry and loading conditions. It is also noteworthy that in this work, a predictor-corrector method is employed to deal with evolution equation resulting from internal variable, in contrast to the Euler-Backward method employed in Ref. [47].

The rest of the paper is organized as follows. The kinematic quantities of thin elastomers are described in Section 2. Afterwards, in Section 3, the nonlinear viscoelastic formulation for elastomeric membrane is derived. Section 4 is allocated to obtaining FE formulation for addressing highly nonlinear formulations. The predictor-corrector method for dealing with evolution equation is discussed in Section 5. In Section 6, some demonstrative problems are solved to investigate the performance of derived formulation. Conclusion is discussed in Section 7.

2 Kinematics of thin elastomeric materials

The kinematical quantities for thin elastomeric materials are briefly described in this section. As depicted in Fig.1, a Cartesian system \{x_1,x_2,x_3\} is defined with orthonormal basis vectors \{{\mathit{e}}_1,{\mathit{e}}_2,{\mathit{e}}_3\}. The location X of each point in pre-deformation membrane may be defined as

(1)\mathit{X}({\mathit{\eta }}^1,{\mathit{\eta }}^2,{\mathit{\eta }}^3)=\mathit{R}({\mathit{\eta }}^1,{\mathit{\eta }}^2)+{\mathit{\eta }}^3\mathit{N},

where \mathit{R} is the location vector of a point in the middle plane, \{{\mathit{\eta }}^1,{\mathit{\eta }}^2,{\mathit{\eta }}^3\} are the coordinate system of the mid-surface of the membrane, and \mathit{N} is the normal vector.

When the membrane is deformed, its new position is obtained as follows

(2)\mathit{x}({\mathit{\eta }}^1,{\mathit{\eta }}^2,{\mathit{\eta }}^3)=\mathit{r}({\mathit{\eta }}^1,{\mathit{\eta }}^2)+{\mathit{\eta }}^3{\lambda }_3({\mathit{\eta }}^1,{\mathit{\eta }}^2)\mathit{n}({\mathit{\eta }}^1,{\mathit{\eta }}^2),

where r is the location vector in the middle plane of after deformation and n is the vector perpendicular to it. Also, {\lambda }_3 is the extension in thickness direction. Moreover, the displacement vector is defined as \mathit{q}=\mathit{r}-\mathit{R} and the tensor of deformation gradient \mathit{F} in the membrane will be calculated as follows

(3)\mathit{F}={\mathit{g}}_{\alpha }\otimes {\mathit{G}}^{\alpha }+{\lambda }_3\mathit{n}\otimes \mathit{N},

where {\mathit{g}}_{\alpha } and {\mathit{G}}^{\alpha }, \alpha =1,2, are the covariant and contravariant basis vectors in spatial and material configurations, respectively. The right Cauchy–Green tensor \mathit{C} in middle plane is expressed as

(4)\mathit{C}={\mathit{F}}^{\mathrm{T}}\mathit{F}={\mathit{C}}_{\alpha \beta }{\mathit{G}}^{\alpha }\otimes {\mathit{G}}^{\beta }+{\mathit{C}}_{33}\mathit{N}\otimes \mathit{N},

where {\mathit{C}}_{\alpha \beta } are the components of right Cauhy–Green tensor in the middle plane of the thin elastomer.

3 Derivation of viscoelastic formula for elastomeric membrane

In this part, an efficient formula for finite deformations of viscoelastic thin elastomers is obtained. The main basis of this theory is the multiplicative decomposition of the deformation gradient tensor into elastic and viscous tensors [48]

(5)\mathit{F}={\mathit{F}}^{\mathrm{e}}{\mathit{F}}^{\mathrm{v}}.

Based on this decomposition, an internal variable named H is defined as follows

(6)\mathit{H}={\mathit{F}}_{\mathrm{v}}^{\mathrm{T}}{\mathit{F}}_{\mathrm{v}},{\mathit{H}}^{\mathrm{T}}=\mathit{H}.

Moreover, incompressibility is assumed for both elastic and viscous states, that is

(7)\mathrm{d}\mathrm{e}\mathrm{t}(\mathit{F})=\mathrm{d}\mathrm{e}\mathrm{t}({\mathit{F}}^{\mathrm{e}})=\mathrm{d}\mathrm{e}\mathrm{t}({\mathit{F}}^{\mathrm{v}})=1.

The strain energy density function is defined as follows

(8)U=U^{\text{e}}(\mathit{C})+U^{\text{v}}({\mathit{C}}^{\text{e}}),{\mathit{C}}^{\text{e}}={\mathit{F}}_{}^{{\text{v}}^{-\text{T}}}{\mathit{C}\mathit{F}}_{}^{{\text{v}}^{-1}}.

Based on the Clausius–Planck inequality, it can be written [49]

(9){\mathcal{D}}_{\mathrm{i}\mathrm{n}\mathrm{t}}=\frac{1}{2}\mathit{S}:\dot{\mathit{C}}-\dot{U}⩾0.

The superscript ‘·’ indicates the time derivative and {\mathcal{D}}_{\mathrm{i}\mathrm{n}\mathrm{t}}⩾0 is the internal loss and is defined as follows

(10){\mathcal{D}}_{\mathrm{i}\mathrm{n}\mathrm{t}}=\frac{\mathrm{d}\psi }{\mathrm{d}\dot{\mathit{H}}}:\dot{\mathit{H}}⩾0,

where \psi is the dissipation potential. Finally, the Clausius–Planck inequality can be re-arranged as

(11)\left(\frac{1}{2}\mathit{S}-\frac{\mathrm{\partial }U}{\mathrm{\partial }\mathit{C}}\right):\dot{\mathit{C}}-\left(\frac{\mathrm{\partial }U}{\mathrm{\partial }\mathit{H}}+\frac{\text{d}\psi }{\text{d}\mathit{H}}\right):\dot{\mathit{H}}⩾0.

As a result, the second Piola–Kirchhoff stress tensor and the inequality related to the growth of deformation are extracted as follows

(12)\mathit{S}=2\frac{\mathrm{\partial }U}{\mathrm{\partial }\mathit{C}}=P{\mathit{C}}^{-1}+2\frac{\mathrm{\partial }{\hat{U}}^{\text{e}}(\mathit{C})}{\mathrm{\partial }\mathit{C}}+2\frac{\mathrm{\partial }{\hat{U}}^{\text{v}}({\mathit{C}}^{\text{e}})}{\mathrm{\partial }\mathit{C}},

(13)\left(\frac{\mathrm{\partial }U}{\mathrm{\partial }\dot{\mathit{H}}}+\frac{\text{d}\psi }{\text{d}\dot{\mathit{H}}}\right):\dot{\mathit{H}}⩾0.

By performing calculations on Eq. (13), finally the main form of the evolution equation is obtained as follows [47]

(14)\text{dev}\left\{\left(\frac{\mathrm{\partial }{\hat{U}}^v({\mathit{C}}^{\text{e}})}{\mathrm{\partial }\mathit{H}}+\frac{\text{d}\psi }{\text{d}\dot{\mathit{H}}}\right)\mathit{H}\right\}=0,

where ‘dev’ represents the deviatoric part and is defined as follows [50]

(15)\text{dev}(\cdot )=(\cdot )-\frac{1}{3}\text{tr}(\cdot )\mathit{I},

where I is the identity tensor. Assuming incompressibility constraint, the general state of the evolution equation is defined as follows

(16)\begin{array}{rl}& \text{dev}\{-\left(\frac{\mathrm{\partial }{\hat{U}}^{\text{v}}}{\mathrm{\partial }{I}_1^{\text{e}}}+I_{\text{1e}}\frac{\mathrm{\partial }{\hat{U}}^{\text{v}}}{\mathrm{\partial }{I}_2^{\text{e}}}\right){\mathit{H}}^{-1}\mathit{C}\\ +& \frac{\mathrm{\partial }{\hat{U}}^{\text{v}}}{\mathrm{\partial }{I}_2^{\text{e}}}{\mathit{H}}^{-1}\mathit{C}{\mathit{H}}^{-1}\mathit{C}+\frac{\text{d}\psi }{\text{d}\dot{\mathit{H}}}\mathit{H}\}=0.\end{array}

In most cases, the Newtonian fluid is assumed for the viscous part, and in this case, the dissipation potential is defined as follows

(17)\psi =\kappa \text{tr(}{\mathit{D}}^{{\text{v}}^{\text{2}}}),

where \kappa is the viscosity and {\mathit{D}}^{\mathrm{v}} is the viscous deformation rate tensor. Here, assuming the neo-Hookean energy function for the viscous part, that is

(18){\hat{U}}^{\text{v}}({\mathit{C}}^{\text{e}})=\frac{{\mu }^{\text{v}}}{2}(I_1^{\text{e}}-3).

The evolution function can be obtained as follows

(19)\dot{\mathit{H}}+\frac{1}{\tau }\left\{\frac{1}{3}\text{tr}({\mathit{H}}^{-1}\mathit{C})\mathit{H}-\mathit{C}\right\}=0,\mathit{H}(0)=\mathit{I},

where \tau =\kappa /{\mu }^{\text{v}} is called the relaxation time. It should be noted that this is the simplest evolution equation that can be calculated, it will be a nonlinear equation in terms of internal variable H.

4 Finite element analysis

In this section, a Total Lagrangian FE formula is obtained. The basis of the formulation is the use of the virtual work principle, which is expressed as follows [51]

(20)\delta W^{\mathrm{i}\mathrm{n}\mathrm{t}}=\delta W^{\text{ext}},

where \delta W^{\mathrm{i}\mathrm{n}\mathrm{t}} and \delta W^{\text{ext}} are virtual internal energy and virtual external work, respectively, which can be calculated as follows

(21)\begin{array}{rl}& \delta W^{\mathrm{i}\mathrm{n}\mathrm{t}}={\int }_{V_0}\mathit{S}\mathbf{:}\delta \mathit{L}\text{d}{V}_0={\int }_{V_0}\delta {\overrightarrow{\mathit{L}}}^{\text{T}}\overrightarrow{\mathit{S}}\text{d}{V}_0,\\ & \delta W^{\text{ext}}={\int }_{V_0}\stackrel{⌢}{\mathit{b}}\cdot \delta \mathit{q}\text{d}{V}_0+{\oint }_{A_0}\stackrel{⌢}{\mathit{T}}\cdot \delta \mathit{q}\text{d}{A}_0,\end{array}

where \stackrel{⌢}{\mathit{b}} is the noncontact body force and \stackrel{⌢}{\mathit{T}} is the surface traction acting on area A_0. Besides, V_0 is the initial volume, and L is the Green−Lagrange strain. The middle surface of S_0^{} is discretized as S_0^{}\approx {\cup }_{e=1}^{\tilde{m}}{S}_0^{\mathrm{e}}, where \tilde{m} represents the total number of elements in the discretized geometry. Then, the virtual strain approximation is obtained as follows

(22)\delta \overrightarrow{\mathit{L}}=\sum _{I=1}^{n_{\text{e}}}{\mathit{B}}_I\delta {\mathit{q}}_I=\sum _{I=1}^{n_{\text{e}}}\left({\mathit{B}}_{0I}+{\mathit{B}}_{1I}\right)\delta {\mathit{q}}_I,

where {\mathit{q}}_I is the displacement vector and the strain–displacement matrices are obtained as follows

(23)\begin{array}{rl}& {\mathit{B}}_{0I}=\left[\begin{array}{cc}{N}_{I,1}& 0\\ 0& N_{I,2}\\ N_{I,2}& N_{I,1}\end{array}\right],\\ & {\mathit{B}}_{1I}=\left[\begin{array}{cc}{N}_{I,1}{q}_{1,1}& N_{I,1}{q}_{2,1}\\ N_{I,2}{q}_{1,2}& N_{I,2}{q}_{2,2}\\ N_{I,1}{q}_{1,2}+N_{I,2}{q}_{1,1}& N_{I,1}{q}_{2,2}+N_{I,2}{q}_{2,1}\end{array}\right],\end{array}

where N_I is the shape function. As a result, the discrete form of virtual internal energy will be as follows

(24)\delta W^{\mathrm{i}\mathrm{n}\mathrm{t}}={\mathit{A}}_{e=1}^{\text{m}}\sum _{I=1}^{n_e}\delta {\mathit{q}}_I^{\text{T}}{\mathit{f}}_I^{\mathrm{i}\mathrm{n}\mathrm{t}},{\mathit{f}}_I^{\mathrm{i}\mathrm{n}\mathrm{t}}={\int }_{V_0^e}{\mathit{B}}_I^{\text{T}}\overrightarrow{\mathit{S}}\text{d}{V}_0^e,

where {\mathit{f}}_I^{\mathrm{i}\mathrm{n}\mathrm{t}} is the internal force vector and {\mathit{A}}_{e=1}^{\text{m}} is an operator that considers the effect of all elements. Also, the discrete form of virtual work of external forces is obtained as follows

(25)\begin{array}{rl}& \delta W^{\text{ext}}={\mathit{A}}_{e=1}^{\text{m}}\sum _{I=1}^{n_e}\delta {\mathit{q}}_I^{\text{T}}{\mathit{f}}_I^{\text{ext}}{=\mathit{A}}_{e=1}^{\text{m}}\delta {\mathit{q}}_e^{\text{T}}{\mathit{f}}_e^{\text{ext}},\\ & {\mathit{f}}_I^{\text{ext}}={\int }_{V_0^e}{N}_I\stackrel{⌢}{\mathit{b}}\text{d}{V}_0^e+{\int }_{A_0^e}{N}_I\stackrel{⌢}{\mathit{T}}\text{d}{A}_0^e,\end{array}

where {\mathit{f}}_e^{\text{ext}} is the external load vector per element. Finally, by placing the obtained values for internal energy and external work in the principle of virtual work, the following relationship is obtained

(26){\hat{\mathit{R}}=\mathit{A}}_{e=1}^{\text{m}}{\hat{\mathit{R}}}_e^{}={\mathit{A}}_{e=1}^m({\mathit{f}}_e^{\mathrm{i}\mathrm{n}\mathrm{t}}-{\mathit{f}}_e^{\text{ext}})=\mathbf{0}.

This relationship is a set of nonlinear equations based on displacements, which is treated via the Newton–Raphson method. The linearization format of this relationship can be obtained by calculating the growth of virtual energy

(27)\begin{array}{rl}\mathrm{\Delta }\delta W^{\mathrm{i}\mathrm{n}\mathrm{t}}=& {\mathit{A}}_{e=1}^{\text{m}}{\int }_{V_0^e}[(\delta {\overrightarrow{\mathit{L}}}^{\text{T}}\mathrm{\Delta }\overrightarrow{\mathit{S}}+\mathrm{\Delta }\delta {\overrightarrow{\mathit{L}}}^{\text{T}}\overrightarrow{\mathit{S}})\text{d}{V}_0^e\\ =& {\mathit{A}}_{e=1}^{\text{m}}\sum _{I=1}^{n_e}\sum _{J=1}^{n_e}\delta {\mathit{q}}_I^{\text{T}}{\mathit{K}}_{IJ}\mathrm{\Delta }{\mathit{q}}_J={\mathit{A}}_{e=1}^m\delta {\mathit{q}}_e^{\text{T}}{\mathit{K}}_e\mathrm{\Delta }{\mathit{q}}_e\end{array},

where {\mathit{K}}_e is the stiffness matrix of the element and {\mathit{K}}_{IJ} is the stiffness matrix in the IJ block, which can be calculated from the following equation

(28){\mathit{K}}_{IJ}={\int }_{V_0^e}({{\mathit{B}}_I}^{\mathrm{T}}\tilde{\mathit{C}}{\mathit{B}}_J+{\mathit{G}}_{IJ})\text{d}{V}_0^e,

where \tilde{\mathit{C}} is the material matrix and {\mathit{G}}_{IJ} is the geometry part of the stiffness matrix. Finally, the linearized form of the algebraic equations system will be \mathit{K}\mathrm{\Delta }\mathit{Q}={\mathit{A}}_{e=1}^m{\mathit{K}}_e\mathrm{\Delta }{\mathit{q}}_e=-\hat{\mathit{R}}.

5 Time integration method

In this section, to solve the evolution equation for internal variable H, a predictor-corrector method is used. To do so, it is assumed that the parameters are known in time step t_n and the method is employed to obtain the variable in time t_{n+1}. At predictor stage, the internal variable H is predicted by Euler-backward method as follows

(29){\mathit{H}}_{n+1}^{(0)}={\mathit{H}}_n+\frac{\mathrm{\Delta }t}{\tau }\left\{{\mathit{C}}_n-\frac{1}{3}\text{tr}({\mathit{H}}_n^{-1}{\mathit{C}}_n){\mathit{H}}_n\right\}.

Next, based on this prediction of H, the FE formulation is solved and the approximation for total displacement vector Q is obtained. In the corrector stage, the midpoint rule is used as follows

(30)\begin{array}{rl}{\mathit{H}}_{n+1}^{(1)}=& {\mathit{H}}_n+\frac{\mathrm{\Delta }t}{2\tau }\{{\mathit{C}}_n+{\mathit{C}}_{n+1}^{(1)}\\ & -\frac{1}{3}\text{tr}({\mathit{H}}_n^{-1}{\mathit{C}}_n){\mathit{H}}_n-\frac{1}{3}\text{tr}({\mathit{H}}_{n+1}^{{(1)}^{-1}}{\mathit{C}}_{n+1}^{(1)})\}{\mathit{H}}_{n+1}^{(1)}.\end{array}

In this case, the FE formulation is solved by new corrected form of H, and new values for Q are derived.

6 Demonstrative examples

In this section, in order to validate the performance of the formulas obtained in the previous sections, some numerical examples are presented. All programming is done in MATLAB software and four-noded elements with 12 degrees of freedom are used to discretize the geometry. Also, 2 × 2 Gauss–Logendre quadrature is used to calculate all integrals.

6.1 Equibiaxial tension of a viscoelastic rubber membrane

In the first problem, the deformation of a rubber subjected to the equibiaxial tension is studied. Neo–Hookean energy function for the elastic part and the viscous branch are considered

(31){\hat{U}}^e=c_1(I_1-3),{\hat{U}}^{\text{v}}=({\mu }^{\text{v}}/2)(I_1^{\text{e}}-3).

Also, relaxation time \tau =1s and material parameters are referred to as c_1=0.5\text{MPa} and {\mu }^{\text{v}}=1\text{MPa} and due to symmetry, only one-fourth of the geometry will be examined. The rubber is suddenly stretched to the tension {\lambda }_1={\lambda }_2=1.5 and then kept in this position until it relaxes. Convergence analysis shows that a 4 × 4 meshing is enough to reach a converged response. The stress parameters P_{11}/\mu and P_{22}/\mu are shown in Fig.2(a), and the inverse of the components for internal variable levels of {\mathit{H}}_{}^{-1} are shown in Fig.2(b). As can be seen, the results of this modeling match well with the results of Refs. [50, 52]. It is necessary to explain that the model considered in Ref. [52] is totally different from the model proposed in this article. Moreover, in Ref. [50], the authors used general nonlinear viscoelastic model, and used trapezoidal implicit model for time integration. Besides, in Ref. [52], the internal variable was the inverse of H which led to completely different evolution equation for the viscoelastic model. However, in this work, H is considered as the internal variable and there is no need to calculate the inverse of it. Hence, it decreases the computational cost.

6.2 Viscoelastic rubber membrane subjected to simple shear

In this example, shear of a viscoelastic rubber will be analyzed. Following Ref. [53], a viscoelastic rubber made of bromobutyl (BIIR) is considered. The rubber is extended by simple shear method as far as {\lambda }_1=2 via two strain rates, {\dot{\lambda }}_1= 20%/min and {\dot{\lambda }}_1= 5%/min which means that for the first rate, the deformation occurs between and for the latter one is . Following Ref. [53], a hyperelastic Mooney–Rivlin model of which material parameters are c_1=79.11\text{kPa} and c_2=8.73\text{kPa} for elastic branch is considered. For viscous branch, a neo–Hookean model is utilized as follows

(32){\hat{U}}^{\text{e}}=c_1(I_1-3)+c_2(I_2-3),{\hat{U}}^{\text{v}}=({\mu }^{\text{v}}/2)(I_1^{\text{e}}-3).

Because the strain rate is different, consequently, the relaxation time will be different, and \tau =110\mathrm{s} for the first strain rate and \tau =450\text{s} for the latter one are adjusted. The results of the present approach are compared against four integral-based viscoelastic models including as Bernstein−Kearsley−Zapas (BKZ), Christensen, Simo, Fosdick–Yu [54–57]). The curves of stress versus strain values for integral-based viscoelastic models and current differential-based nonlinear viscoelastic theory are displayed in Fig.3(a) for \dot{\epsilon }= 20%/min and Fig.3(b) for \dot{\epsilon }= 5%/min. As can be seen, the current model with only one viscoelastic strain energy is able to predict the laboratory results better than the previous models that considered four viscoelastic branches. This example clearly shows the superiority of nonlinear viscoelastic theory over linear viscoelastic theory. To be more specific, the nominal stress values at strains 0.2, 0.4, 0.6, 0.1, 1.0 are shown in Tab.1 for \dot{\epsilon }= 20%/min and Tab.2 for \dot{\epsilon }= 5%/min. As can be deduced from the results in the Tab.1 and Tab.2, the values of nominal stress from the current study are closer to the experimental results in comparison to other models, although just four material parameters are used in the proposed model in comparison to 12 material parameters for other ones.

6.3 Inflation and creep of clamped-clamped cylindrical elastomeric membrane

In this example, inflation of clamped-clamped cylinder is investigated. It is assumed that the elastic branch follows the Mooney–Rivlin model, while the viscoelastic branch follows the Neo–Hookean model as follows

(33){\hat{U}}^{\text{e}}=c_1(I_1-3)+c_2(I_2-3),{\hat{U}}^{\text{v}}=({\mu }^{\text{v}}/2)(I_1^{\text{e}}-3).

Following Ref. [58] for elastic problem, it is assumed that the relation between material parameters and the relation between radius and length are c_2=0.1c_1 and L=6r_0, respectively. In this example, the radius is r_0=1 and the relaxation time is \tau =20\text{s}. Because of the symmetry property, one-quarter of membrane is analyzed. At first stage, it is assumed that the inflation procedure is performed between as far as the maximum deformation which occurs in the middle of cylinder at point A reaches 4. Our numerical study shows that 6\times 6 elements suffice to get converged results. However, a mesh of 10\times 10 elements is employed in this problem for more accuracy. Fig.4(a) shows the curves of pressure parameter p_{{}_{\ast }}=p{r}_0/c_1{h}_0 versus maximum deflection of cylindrical membrane for different values of {\mu }^{\text{v}}. As shown, the result is coincident with the elastic result reported by Ref. [58] for {\mu }^{\text{v}}=0. In the next stage, the creep phenomenon relaxation time \tau =40\text{s} is studied. It is assumed that the membrane nondimensional pressure increases linearly during until p_{{}_{\ast }}=2.1, and then is held constant to permit the membrane to creep during . The curves of creep of clamped-clamped cylindrical membrane are displayed in Fig.4(b) for various values of {\mu }^{\text{v}}. Furthermore, the deformation profiles of membrane for miscellaneous time steps of creep are portrayed in Fig.5. Finally, the convergence analyses are done for inflation and creep procedures. For inflation, the maximum deflection is considered to be 1, and the convergence analysis shows that considering 10\times 10 is enough for achieving converged results for p_{{}_{\ast }}=p{r}_0/c_1{h}_0 for all values of {\mu }^{\text{v}} (Fig.6(a)). Moreover, for creep phenomenon, the maximum deflection at point A is considered, and the convergence is analyzed at time t=300\text{s}. In this analysis, as displayed in Fig.6(b), it is concluded that a mesh of 100 elements is enough for reaching converged results.

7 Conclusions

In this paper, an efficient nonlinear viscoelastic formula for large deformation of thin elastomeric membranes is derived. The basis of the derivation is multiplicative decomposition of the tensor of deformation gradient into elastic and viscous partitions. The right Cauchy–Green deformation tensor in the viscous branch is considered as internal variable, and the evolution equation are obtained. To solve the evolution equation, a predictor-corrector method is employed. Moreover, to deal with nonlinear viscoelastic formulation, a FE method in Total Lagrangian framework is derived. Finally, three examples are solved to show the power of the derived model. The findings of the present model can be mentioned as follows.

1) In the first example, it is shown that the present model is in excellent agreement with the results in literature. Indeed, using predictor-corrector algorithm gives faster convergence in comparison to common Euler-backward and trapezoidal method.

2) In the second example, it is proven that the present model can overcome the drawbacks of results obtained from previous linear viscoelasticity theory and anticipate the experimental results more accurately with fewer material parameters. In fact, the present model which is a simplified nonlinear model can predict the experimental data with four material parameters by higher accuracy in comparison to finite linear models with 12 material parameters.

3) Moreover, in the third example, inflation, and inflation-creep of a cylindrical elastomeric membrane is studied, and it is shown that incorporating viscoelasticity leads to requirement more pressure the reach the same amount of maximum deflection.

In general, the proposed model in this paper which is a simplified nonlinear viscoelastic model considering neo–Hookean model for viscoelastic branch takes the advantage of simplicity as well as predicting the experimenal results with higher accuracy while uing fewer material parameter. Moreover, using predictor-corrector algorithm results in faster convergence in comparison to common Euler-backward implicit method.

References

[1] Demir Ç, Civalek Ö. On the analysis of microbeams. International Journal of Engineering Science, 2017, 121: 14–33

[2] Abouelregal A E, Ersoy H, Civalek Ö. Solution of Moore–Gibson–Thompson equation of an unbounded medium with a cylindrical hole. Mathematics, 2021, 9(13): 1536

[3] Adiga S P, Jin C, Curtiss L A, Monteiro-Riviere N A, Narayan R J. Nanoporous membranes for medical and biological applications. Wiley Interdisciplinary Reviews. Nanomedicine and Nanobiotechnology, 2009, 1(5): 568–581

[4] Radu E R, Voicu S I, Thakur V K. Polymeric membranes for biomedical applications. Polymers, 2023, 15(3): 619

[5] Firouzi N, Misra A. New insight into large deformation analysis of stretch-based and invariant-based rubber-like hyperelastic elastomers. Thin-walled Structures, 2023, 192: 111162

[6] Cooley C G, Lowe R L. Nonlinear vibration of dielectric elastomer membranes with axial inertia effects. International Journal of Mechanical Sciences, 2023, 248: 108205

[7] Zhao B, Hu J, Chen W, Chen J, Jing Z. A nonlinear uniaxial stress-strain constitutive model for viscoelastic membrane materials. Polymer Testing, 2020, 90: 106633

[8] Li E H, Li Y Z, Xie J Y, Sun Y H, Yang L Z, Ning X W. A fuzzy coordination control of a water membrane evaporator cooling system for aerospace electronics. Applied Thermal Engineering, 2021, 191: 116872

[9] Huang H, Li X, Xue S, Luo Y, Shi D, Hou X, Liu Y, Li N. Performance and measurement devices for membrane buildings in civil Engineering: A review. Applied Sciences, 2022, 12(17): 8648

[10] Tian L, Chen F, Zhu L, Yu T. Large deformation of square plates under pulse loading by combined saturation analysis and membrane factor methods. International Journal of Impact Engineering, 2020, 140: 103546

[11] Hou X, Mao Y, Zhang R, Fang D. Super-flexible polyimide nanofiber cross-linked polyimide aerogel membranes for high efficient flexible thermal protection. Chemical Engineering Journal, 2021, 417: 129341

[12] Khmelinskii I, Makarov V I. On the effects of mechanical stress of biological membranes in modeling of swelling dynamics of biological systems. Scientific Reports, 2020, 10(1): 8395

[13] Yin Y, Pu D, Xiong J. Analysis of the comprehensive tensile relationship in electrospun silk fibroin/polycaprolactone nanofiber membranes. Membranes, 2017, 7(4): 67

[14] Guerra N B, Pegorin G S A, Boratto M H, de Barros N R, de Oliveira Graeff C F, Herculano R D. Biomedical applications of natural rubber latex from the rubber tree Hevea brasiliensis. Materials Science and Engineering C, 2021, 126: 112126

[15] Lima L R, Caetano G F, Soares V O, Santos R J, Malmonge J A, Silva M J, Yarin A L. Evaluation of tensile, thermal, and biological properties of natural rubber-based biocomposite with biosilicate and 45SS-K bioglass. Journal of Applied Polymer Science, 2023, 140(22): e53894

[16] Tamadapu G, Dasgupta A. Finite inflation of a hyperelastic toroidal membrane over a cylindrical rim. International Journal of Solids and Structures, 2014, 51(2): 430–439

[17] Farsad K, de Camilli P. Mechanisms of membrane deformation. Current Opinion in Cell Biology, 2003, 15(4): 372–381

[18] Foster B, Knobloch E. Elastic fingering in a rotating hele-shaw cell. Physical Review E, 2023, 107(6): 065104

[19] Véron A R, Martins A F. Tensorial form of leslie-ericksen equations and applications. Molecular Crystals and Liquid Crystals, 2009, 508(1): 309–336

[20] MarkmanMSimanovskyL. Method of making a parabolic-cylindrical concentrator and device for its implementation. 1982

[21] Gruttmann F, Taylor R. Theory and finite element formulation of rubberlike membrane shells using principal stretches. International Journal for Numerical Methods in Engineering, 1992, 35(5): 1111–1126

[22] Feng X, Huang R Y. Liquid separation by membrane pervaporation: A Review. Industrial & Engineering Chemistry Research, 1997, 36(4): 1048–1066

[23] Cheng T, Dai C, Gan R Z. Viscoelastic properties of human tympanic membrane. Annals of Biomedical Engineering, 2007, 35(2): 305–314

[24] Murayama T, Yoshino M, Hirata T. Three-dimensional lattice boltzmann simulation of two-phase flow containing a deformable body with a viscoelastic membrane. Communications in Computational Physics, 2011, 9(5): 1397–1413

[25] LauberMWeymouthG DLimbertG. Immersed-boundary fluid-structure interaction of membranes and shells. Journal of Physics Conference Series, 2023, 2647, 052002

[26] Mierunalan S, Dassanayake S P, Mallikarachchi H M Y C, Upadhyay S H. Simulation of ultra-thin membranes with creases. International Journal of Mechanics and Materials in Design, 2022, 19(1): 73–94

[27] Zhang X, Ma Z Q, Wu Y, Liu J. Response of mechanical properties of polyvinyl chloride geomembrane to ambient temperature in axial tension. Applied Sciences, 2021, 11(22): 10864

[28] Sarkezi-SelskyPRastedtMLatzAJahnkeT. Lattice-boltzmann simulation of two-phase flows in the GDL and MPL of polymer electrolyte fuel cells. In: EFCF 2019 Proceedings (PEY). Lucerne: EFCF, 2019

[29] van Nguyen C. Numerical investigation of circle defining curve for two-dimensional problem with general boundaries using the scaled boundary finite element method. Frontiers of Structural and Civil Engineering, 2019, 13(1): 92–102

[30] Bhandarkar S D, Betcher J, Smith R, Lairson B M, Ayers T. Constitutive models for the viscoelastic behavior of polyimide membranes at room and deep cryogenic temperatures. Fusion Science and Technology, 2016, 70(2): 332–340

[31] Ghosh A, Both A K, Cheung C L. Characterization of three-dimensional fractional viscoelastic models through complex modulus analysis and polar decomposition. Physics of Fluids, 2022, 34(7): 077115

[32] Geith M A, Eckmann J D, Haspinger D C, Agrafiotis E, Maier D, Szabo P, Sommer G, Schratzenstaller T G, Holzapfel G A. Experimental and mathematical characterization of coronary polyamide-12 balloon catheter membranes. PLoS One, 2020, 15(6): e0234340

[33] Cai R, Jin Y, Rabczuk T, Zhuang X, Djafari-Rouhani B. Propagation and attenuation of rayleigh and pseudo surface waves in viscoelastic metamaterials. Journal of Applied Physics, 2021, 129(12): 124903

[34] Firouzi N. Mechanics of nonlinear visco-hyperelastic-hysteresis membranes. International Journal of Non-linear Mechanics, 2022, 147: 104231

[35] Shafei E, Faroughi S, Rabczuk T. Nonlinear transient vibration of viscoelastic plates: A NURBS-based isogeometric HSDT approach. Computers & Mathematics with Applications, 2021, 84: 1–15

[36] Zhou X, Qiu D, Peng L, Lai X. Modeling the viscoelastic-viscoplastic behavior of glassy polymer membrane with consideration of non-uniform sub-chains in entangled networks. International Journal of Solids and Structures, 2023, 282: 112469

[37] Brighenti R, Rabczuk T, Zhuang X. Phase field approach for simulating failure of viscoelastic elastomers. European Journal of Mechanics. A, Solids, 2021, 85: 104092

[38] Khaniki H B, Ghayesh M H, Chin R, Amabili M. A review on the nonlinear dynamics of hyperelastic structures. Nonlinear Dynamics, 2022, 110(2): 963–994

[39] Khaniki H B, Ghayesh M H, Chin R, Amabili M. Hyperelastic structures: A review on the mechanics and biomechanics. International Journal of Non-linear Mechanics, 2023, 148: 104275

[40] AmabiliM. Nonlinear mechanics of shells and plates: Composite, Soft and Biological Materials. New York, NY: Cambridge University Press, 2018

[41] Amabili M. Derivation of nonlinear damping from viscoelasticity in case of nonlinear vibrations. Nonlinear Dynamics, 2019, 97(3): 1785–1797

[42] Khorasani M, Soleimani-Javid Z, Arshid E, Lampani L, Civalek Ö. Thermo-elastic buckling of honeycomb micro plates integrated with FG-GNPs reinforced Epoxy skins with stretching effect. Composite Structures, 2021, 258: 113430

[43] Shafei E, Faroughi S, Rabczuk T. Multi-patch NURBS formulation for anisotropic variable angle tow composite plates. Composite Structures, 2020, 241: 111964

[44] Sobhani E, Masoodi A R, Civalek Ö, Ahmadi-Pari A R. Agglomerated impact of CNT vs. GNP nanofillers on hybridization of polymer matrix for vibration of coupled hemispherical-conical-conical shells. Aerospace Science and Technology, 2022, 120: 107257

[45] Mercan K, Baltacıoglu A K, Civalek Ö. Free vibration of laminated and FGM/CNT composites annular thick plates with shear deformation by discrete singular convolution method. Composite Structures, 2018, 186: 139–153

[46] Franchini G, Breslavsky I D, Holzapfel G A, Amabili M. Viscoelastic characterization of human descending thoracic aortas uncer cyclic load. Acta Biomaterialia, 2021, 130: 291–307

[47] Le Tallec P, Rahier C, Kaiss A. C, Rahier, A Kaiss, Three-dimensional incompressible viscoelasticity in large strains: Formulation and numerical approximation. Computer Methods in Applied Mechanics and Engineering, 1993, 109(3–4): 233–258

[48] Reese S, Govindjee S. A theory of finite viscoelasticity and numerical aspects. International Journal of Solids and Structures, 1998, 35(26–27): 3455–3482

[49] HolzapfelG A. Nonlinear Solid Mechanics, A Continuum Approach for Engineering. New York: Wiley, 2000

[50] Firouzi N, Żur K K, Amabili M, Rabczuk T. On the time-dependent mechanics of membranes via the nonlinear finite element method. Computer Methods in Applied Mechanics and Engineering, 2023, 407: 115903

[51] ReddyJ N. An Introduction to Nonlinear Finite Element Analysis. Oxford: Oxford University Press, 2015

[52] Kroon M. A constitutive framework for modelling thin incompressible viscoelastic materials under plane stress in the finite strain regime. Mechanics of Time-Dependent Materials, 2011, 15(4): 389–406

[53] Jridi N, Arfaoui M, Hamdi A, Salvia A, Bareille O, Ichchou M, Ben Abdallah J. Separable finite viscoelasticity: Integral-based models vs. experiments. Mechanics of Time-Dependent Materials, 2019, 23(3): 295–325

[54] Bernstein B, Kearsley E, Zapas L. A study of stress relaxation with finite strain. Transactions of the Society of Rheology, 1963, 7(1): 391–410

[55] Christensen R M. A nonlinear theory of viscoelasticity for application to elastomers. Journal of Applied Mechanics, 1980, 47(4): 762–768

[56] Simo J C. On a fully three-dimensional finite-strain viscoelastic damage model: formulation and computational aspects. Computer Methods in Applied Mechanics and Engineering, 1987, 60(2): 153–173

[57] Fosdick R, Yu J H. Thermodynamics stability and non-linear oscillations of viscoelastic solids-II history type solids. International Journal of Non-linear Mechanics, 1998, 33(1): 165–188

[58] Patil A, Nordmark A, Eriksson A. Free and constrained inflation of a pre-stretched cylindrical membrane. Proceedings of the Royal Society A-Mathematical Physical and Engineering Sciences, 2015, 470: 20140282

RIGHTS & PERMISSIONS

Higher Education Press