Investigation of Hydrodynamically Dominated Membrane Rupture, Using Smoothed Particle Hydrodynamics–Finite Element Method

Abstract

The rupturing process of a membrane, located between two fluids at the center of a three-dimensional channel, is numerically investigated. The smoothed particle hydrodynamics (SPH) and the finite element method (FEM) are used, respectively, for modeling the fluid and solid phases. A range of pressure differences and membrane thicknesses are studied and two different rupturing processes are identified. These processes differ in the time scale of the rupture, the location of the rupture initiation, the level of destruction and the driving mechanism.

1. Introduction

Predicting the exact values of impact loads exerted from fluids on structures is of great importance due to its many applications, for example, in cardiovascular systems [1,2,3,4,5,6,7] and industrial flows [8,9]. In particular, simulation of fluid-solid interaction (FSI) between a membrane and two fluids have many applications in petrochemical [10] and aerospace industries. Better understanding of membrane-fluid interactions helps mitigate failures of heat exchangers [9], fuel cells, ship tankers and aircraft fuel containers, which hinge on fluid pressure difference on the sides of membranes. Due to complex nature of FSI, progress in developing effective numerical methods have been slow. However, with the advances in meso-scale methods, for example the Lattice Boltzmann method (LBM) and smoothed particle hydrodynamics (SPH), satisfactory computational schemes have been emerged for handling the complex physics associated with FSI.

The SPH method was first introduced by Monaghan and Gingold [11] for solving astrophysical problems. In this method, the fluid is discretized to many particles, which are tracked in both time and space. The most important advantage of the SPH is the absence of typical errors of Eulerian methods, such as numerical diffusion occurring in high gradient regions of the flow.

The SPH method was first employed in 1990 for simulation of solid materials [12]. It was also used to study the failure of structures. Subsequently, several important developments in SPH were achieved. Swegle et al. [13] analyzed various numerical features of the SPH and explored its stability margins, zero energy modes, interface modeling and artificial viscosity to overcome the numerical challenges of the SPH method. Further developments were obtained by optimizing kernel functions [14] and symmetrical boundaries with extra particles, that is, ghost particles [15].

In early twenty first century, the SPH method was widely used in computational fluid dynamics (CFD), including simulation of compressible flows [16], incompressible flows [17], two-phase flows [18], flows in porous media [19], explosion modeling [20], and non-Newtonian flows [21]. The SPH method has been used as a suitable context to simulate viscoelastic fluid flow by implementing Maxwell [22] and cross models [23]. In addition, Hashemi et al. [24] studied the effect of fluid elasticity on the migration mechanism of solid bodies in Oldroyd-B shear flows. Xu and Ouyang [25] used SPH to simulate the challenging jet buckling and rod-climbing effects of viscoelastic fluids. SPH has shown good results for molding flows of Non-Newtonian fluids for both isothermal [26] and non-isothermal [27] conditions. Recently, Venkatesan and Ganesan [28] employed SPH to investigate viscoelasticity on droplet spreading dynamics on hydrophilic surfaces.

With advancements in simulation of fluid flow and growing interest in FSI, the smoothed particle hydrodynamics was used to simultaneously simulate the fluids and structures. In 2008, Maurel and Combescure [29] provided a formulation for SPH, which enabled it to simulate behaviors of shell structures. This formulation was capable of simulating elastic and plastic deformations of solids until failure. Concurrently, Potapov et al. [30] made an effort to simulate FSI in a problem in which random collision of fast-moving fluid to a structure caused major deformation of the solid.

In 2013, Caleyron et al. [31] modeled the rupture of a shell structure using the SPH method. He considered a cylinder and a piston system, in which a shell was placed at one end of the cylinder. The shell underwent large deformations after a weight was collided with the piston, generating extreme pressure waves. Recently, Faucher et al. [32] implemented a high resolution adaptive framework for fast transient FSI with interfaces and structural failure for application to failing tanks under impact. However, this work was focused on accurate evaluation of crack propagation under fluid pressure but did not analyze the structure failure modes. The present work aims to systematically identify the hydrodynamically dominated membrane rupture processes.

In this study, the transient FSI problem of membrane ruputre is investigated, in which the flow and the solid deformations are computed using, respectively, the SPH and finite element method (FEM). Initially, a membrane is located at the middle of a vertical three-dimensional rectangular channel between two weakly compressible liquids with different properties. By increasing the pressure difference, the membrane experiences large deformation and finally ruptures when the pressure difference becomes sufficiently large. The main contribution of the present study is in identifying different membrane failure modes, which are governed by different mechanisms, namely, the pressure waves and steady pressure force field. The rupture initiation time, the affected area, pressure difference sensitivity, thickness sensitivity and the level of destruction are the factors that influence different modes of membrane rupture.

2. Governing Equations

In this section, the governing equations of fluid flow and the solid phase are described. In addition, the damage model of the solid phase is briefly discussed. Note that the isothermal condition is assumed throughout this study and therefore, there was no need to sue the energy equation in fluid and solid phases. In the SPH method, the fluid dynamics are derived in a Lagrangian framework, similar to the lattice Boltzmann and MD [33], which obviates the need for a mesh to solve the governing equations, in contrast to the traditional CFD methods [34,35,36,37].

2.1. Fluid Phase

In this work, the fluid phase is modeled using the SPH method that is capable of handling free-surface flows involving ruptures and multiphase flows with small density differences [38]. The corresponding governing equations are the Lagrangian form of the Navier-Stokes equations, which in SPH context are given as,

$$\begin{array}{c}\frac{D{\rho }_i}{Dt}=\sum _{j=1}^N{m}_j{v}_{ij}^{\alpha }\frac{\partial W_{ij}}{\partial x_i^{\alpha }},\hfill \\ \frac{D{v}_i^{\alpha }}{Dt}=-\sum _{j=1}^N{m}_j(\frac{{\sigma }_i^{\alpha \beta }}{{\rho }_i^2}+\frac{{\sigma }_j^{\alpha \beta }}{{\rho }_j^2})\frac{\partial W_{ij}}{\partial x_i^{\beta }}+F_i^{\alpha },\hfill \end{array}$$

(1)

where, the superscripts $\alpha$ and $\beta$ are coordinate direction indicators and the subscripts i and j denote the central particle and its neighboring particles, respectively. In addition, $\rho$ is density, m is mass, W is kernel function, v is velocity, $\sigma$ is total stress tensor, F is external force and $v_{ij}=v_i-v_j$ [39,40]. The continuity equation for a weakly compressible fluid is given as,

$$\frac{\partial \rho }{\partial t}+\rho (\nabla .u)=0.$$

(2)

Note, in Equation (2), $u.\nabla \rho$ that appears in the compressible continuity equation has been ignored. This assumption is valid when density variation is orthogonal to the velocity vector. This occurs when the maximum velocity of the fluid is negligible in comparison with the speed of sound [17]. Using Equation (2), for calculation of pressure p, one can use the following equation,

$$dp=c^{2}d\rho$$

(3)

Since the flow is considered to be inviscid, the viscous effects are also ignored throughout this study. Thus, the momentum equation becomes:

$$\frac{D{v}^{\alpha }}{Dt}=-\frac{1}{\rho }\frac{\partial p}{\partial x^{\alpha }}+F^{\alpha }$$

(4)

Substituting Equations (2) and (4) into Equation (1) leads to:

$$\begin{array}{c}{(\frac{D\rho }{Dt})}_i={\rho }_i\sum _{j=1}^N\frac{m_j}{{\rho }_j}(v_{ij}^{\alpha })\frac{\partial W_{ij}}{\partial x_i^{\alpha }},\hfill \\ \frac{D{v}_i^{\alpha }}{Dt}=-\sum _{j=1}^N{m}_j(\frac{p_i}{{\rho }_i^2}+\frac{p_j}{{\rho }_j^2}+\prod _{ij}{\delta }^{\alpha \beta })\frac{\partial W_{ij}}{\partial x_i^{\beta }}+F_i^{\alpha },\hfill \end{array}$$

(5)

where W denotes a cubic B-spline kernel function and ${\delta }^{\alpha \beta }$ is Kronecker’s delta function [41,42]. In Equation (5), artificial viscosity terms (${\prod }_{ij}$) are added to the momentum equation to enhance the stability of the method for the computation of discontinuities [43]. Note that the reason for including the fluid phase in this study, rather than a simple pressure loading on the membrane, is that we are interested in evaluation of flow after the rupture of the membrane which is considered for the future work.

2.2. Solid Phase

The behavior of the materials used in the solid phase are assumed to be independent of strain rate and thus von Mises’ elastic-plastic laws are utilized. The Lemaitre-Chaboche damage laws are also coupled with Von Mises’ elastic-plastic laws to model the rupturing mechanism of the membrane. The Lemaitre-Chaboche failure criterion ($f_D$) is presented in the following form [44],

$$f_D=[\frac{2}{3}(1+\nu )+3(1-2\nu ){(\frac{{\sigma }_h}{{\sigma }_{eq}})}^2]P-{ϵ}_p^s$$

(6)

where $\nu$ is Poisson’s ratio, ${\sigma }_h$ is the hydrostatic stress, ${\sigma }_{eq}$ is Von Mises’ equivalent stress and p is the isotropic pressure. Also, ${ϵ}_p^s$ denotes the threshold equivalent plastic strain, beyond which microcrackes start to grow. This criterion gives the damage evolution (${\dot{D}}_{nc}$) in the following form,

$${\dot{D}}_{nc}=\frac{D_c}{{ϵ}_p^c-{ϵ}_p^s}[\frac{2}{3}(1+\nu )+3(1-2\nu ){(\frac{{\sigma }_h}{{\sigma }_{eq}})}^2]\dot{P}$$

(7)

where ${ϵ}_p^c$ is the critical equivalent plastic strain, beyond which fracture of material occurs. Also, $D_c$ is the volume fraction of defects at the fracture point [31].

Bazant and Belytschko [45] and Needleman [46] showed that the damage localization problem for models similar to Equation (7) reduces to transformation of the original hyperbolic problem to an ill-posed elliptic one, when material softens. Therefore, the majority of proposed solutions emphasizes on preserving the hyperbolic nature of equations. Mostly suggested solutions utilize a characteristic time or length to moderate the interactions between the elements. similar to recent work in Reference [31], we compute damage rate ($\dot{D}$) as [47],

$$\dot{D}=\frac{1}{{\tau }_c}(1-exp(-a|{\dot{D}}_{nc}-D|)),$$

(8)

where ${\dot{D}}_{nc}$ is calculated form Equation (7) and ${\tau }_c$ represents the inverse of maximum damage rate. Also, a denotes a coefficient related to the material properties of the membrane.

3. Numerical Modeling

In this section, the material characteristics and element topologies of the fluid phase are first described and then, the features of the employed materials and elements of solid phase are presented. Finally, the initial and boundary conditions that includes the utilized FSI interface are discussed. The entire numerical modeling is performed using the open version of EUROPLEXUS, which is a transient FSI computer code [48].

3.1. Fluid Phase

For fluid discretization, spherical fluid particles with a radius of 650 $\mathsf{\mu }$m are employed. To decide on the size of the particles, dependency studies are performed, which are presented in the next section. Each of these particles have 6 degrees of freedom and a total of 211,146 particles are utilized for both high and low pressure fluids. In Figure 1, the blue spheres show high pressure liquid (water) and gold spheres represent low pressure liquid (oil). Also, the pink plane placed between the two liquids is the membrane under study. As it is shown in this figure, the lower bottom of the channel is closed by a rigid wall, while the upper end is open and thus, a free surface of low pressure fluid is created. The characteristics of the fluid phase are presented in

Table 1. Fluid phase characterestics.

Fluid	Density (Kg/m${}^3$)	Speed of Sound (m/s)
Water	1000	1450
Diesel Oil	900	1250

.

3.2. Solid Phase

The solid phase is composed of two main parts, namely: channel walls and the membrane placed in the middle of the channel, both of which are modeled using FEM. The three-dimensional channel has a length of 20 cm, which is placed vertically. It is closed from the bottom but its top end is open. It has a cross section of 5 cm × 5 cm square. In simulating the behavior of the structural components, 16,200 two-dimensional plane elements are used with the optimal dimensions of 1.6 cm × 1.6 mm. Each of these elements has 20 Gaussian points for stress calculations and every 4 Gaussian points are placed in 5 different layers of its thickness. Three different thickness length are considered for the membrane, ranging from 1 to 3 mm. Given the fact that the components of the structure are expected to show different behaviors, we utilize two different materials. The channel walls are assumed to be made of iron with perfectly-elastic behavior. The features of this component are given in

Table 2. Solid phase characteristics.

Material	Application	Density [Kg/m${}^3$]	Young’s Modulus [GPa]	Poisson’s Ratio	Elastic Limit [MPa]
Iron	Channel wall	7860	200	0.33	Not applicable
XC 38 steel	Membrane	7860	200	0.33	166.5

. For the membrane, von Mises’ elastic-plastic laws are employed, coupled by the Lemaitre-Chaboche damage law with XC 38 steel material characteristics. Features of this component are illustrated in Table 2 and

Table 3. Required parameters for damage modeling of XC38 steel [44,49].

Parameter	${\mathit{ϵ}}_{\mathit{p}}^{\mathit{c}}$	${\mathit{ϵ}}_{\mathit{p}}^{\mathit{s}}$	${\mathit{D}}_{\mathit{c}}$	${\mathit{\tau }}_{\mathit{c}}[\mathit{ns}]$	a
Value	0.0	0.56	0.22	10	1

and also in Figure 2.

3.3. Initial and Boundary Conditions

As noted before, Figure 1 displays the initial conditions of the problem, wherein the pressure on both sides of the membrane is fixed at 1 KPa. For the region placed beneath the membrane, a high pressure is suddenly applied to the fluid. The magnitude of applied pressure varies depending on the test cases as will be discussed later. The free surface of the low-pressure fluid is modeled by periodic boundary condition, in order to let the high-pressure fluid to flow towards the low-pressure side after rupturing the membrane. All the membrane elements, except the ones attached to the wall, are treated using the master-slave method [50] for detection of interpenetration between the fluid and solid parts. Whenever interpenetration of two bodies is detected, proper contact forces are applied to them by the Lagrangian multiplier method [51].

4. Grid Resolution Study

The size of both SPH fluid particles and FEM solid elements are altered to study their effect on the displacement of the membrane central node (see

Table 4. Displacement of membrane central node for different particle and element sizes (all sizes are in mm).

Element Side/Particle Radius	5	2.5	1.25	0.83	0.65	Relative Error (%) (Effect of Element Size)
10	4.72	6.00	5.74	5.62	5.37	4.4
5	4.61	5.99	6.38	6.73	6.53	3
2.5	4.33	5.42	6.66	7.21	7.10	1.5
1.66	4.36	5.09	6.4	7.16	7.19	0.4
Relative error (%) (effect of particle radius)	0.7	6.1	4	0.7	1.3	-

). In all of these cases, a pressure difference of 250 MPa is exerted to a membrane and the final displacement of the central node of the membrane is measured.

As shown in Table 4, by reducing the size of the elements, relative error, which is defined as the ratio of the changes in central node displacement to its initial displacement, is decreased from 4.4% to 0.4%. In this table, only the last relative errors of the last two columns (for the study of element size) and the last two rows (for the study of particle radius) are presented for the sake of brevity. To be conservative, a 1.66 mm × 1.66 mm element was used.

Reducing the particle radius does not show similar effects due to the immobility of the particles before rupturing. For this reason, we resorted to employ equivalent plastic strain to study the effects of particle radius. According to Figure 3, it is clear that for particle sizes of 5 and 2.5 mm equivalent plastic strain patterns are highly asymmetric and thus incorrect, as far as the particles’ impacting area to the membrane is seen. These figures also show that by altering the particle size from 0.83 to 0.65 mm equivalent plastic strains do not change. Consequently, a relatively conservative choice is SPH particles of radius 0.65 mm.

5. Code Validation

Three related test cases were previously studied to validate the EUROPLEXUS code. First, acoustic oscillations of a fluid contained in a cylinder with a piston was modeled using the SPH method and the results are compared with the FE-FV method [52]. The second case investigated the bending of a plate, which was modeled using FEM, under the impact of a column of water modeled with the SPH method [30]. In this case, the capabilities of the SPH-FEM, included in the code for FSI simulations, were studied. Third, the SPH results were compared to experimental data for simulation of the leakage rate of a metal shell at the bottom of water column, impacted by a weight at its top piston [52]. All these test cases gave satisfactory results for SPH-FEM implemented in the EUROPLEXUS code for simulating different FSI problems.

Here, it is necessary to satisfy two critical conditions. First, momentum transport from fluid particles to the membrane elements must be modeled properly. Also, the structural modeling of the membrane must be of high accuracy to simulate material behavior under external forces. Due to the lack of experimental data for von Mises’ elastic-plastic material interacting with fluids, we resorted to empirical data of materials with hyperelastic behavior for validation purposes. For code validation, the experimental data of a FEM modeling of a rubber membrane under multi-axial stress is used [53]. For modeling natural rubber, we use the Mooney-Rivlin model. The fluid is modeled using the SPH method. Also, the Master-Slave method is employed for the FSI interfaces. The comparison of SPH-FEM results with the experimental data and a full FEM method is presented in Figure 4, showing satisfactory results (maximum error of 4%).

6. Results and Discussion

In this section, we investigated the occurrences of the start of membrane deformation until complete failure in two stages:

• membrane deformation (at this stage, instantaneous rise in water pressure causes membrane deformation) and,
• membrane rupture (the mechanism which controls the initiation and propagation of cracks).

Note that throughout this work, we used three different membrane thicknesses, namely: thin membrane (1 mm thickness), membrane with intermediate thickness (2 mm thickness) and thick membrane (3 mm thickness). Therefore, the deformation and rupture results are presented for these membrane thicknesses.

6.1. Membrane Deformation

The effects of of pressure difference on the membrane deflection is first examined. Then, the time variation of deformation of the membrane subjected to the imposed pressure difference is studied. Clearly, by applying higher pressure differences on the membrane, its deflection increases, until micro cracks start to emerge at a critical pressure. In order to provide quantitative results, the behavior of a particular point on the membrane is examined. Here the deflection of the central point is studied, since it is more sensitive to pressure differences. Figure 5 shows the variation of membrane central node deflection versus imposed pressure difference. It is from this figure that the central point deflection increases linearly with pressure difference.

The equivalent plastic strain contours of a thin membrane subjected to a pressure difference was also studied and the results are shown in Figure 6. It is important to note that the maximum plastic strain occur at the center and at the edges of the membrane. Strain concentration on the edges of the membrane is because in those areas one side of the membrane is fixed to the wall. Therefore, the membrane deformation in the vicinity of the edges leads to large amount of strain. At the center of the membrane, since maximum displacement occurs, strain is large. Consequently, these areas are the potential points for onset of membrane rupture.

Another important parameter is the time variation of membrane deformation after imposing a pressure difference. The time histories of the membrane central point displacement for different thicknesses under various imposed pressure differences are examined. Figure 7a shows the corresponding time evolutions of the central point deflection for different thickness subject to a pressure difference of 450 MPa. It is seen that the deflection increases with time and reaches a saturation level. Also as the membrane thickness increases, the slope of the time evolution curve and the peak deformation decrease. It is also seen that all curves reach to their maximum displacement almost at the same time of about 0.12 ms.

Similarly, Figure 7b shows similar trends for a pressure difference of 250 MPa. Note, for this pressure difference, the steady state reaches about 0.15 ms for all membrane thicknesses. It seems that the time required for reaching the membrane steady state deflection is a function of the applied pressure difference.

In order to investigate the time history of membrane deformation under a given pressure difference, the equivalent plastic strain contours for a 1 mm thickness membrane under a 450 MPa pressure difference are shown in Figure 8. These contour plots show that initially, the edges of the membrane experience the maximum strain, that is, the pressure pulse exerted on the membrane generates higher stress around fixed boundaries of the membrane. The strain in the central regions increases with time and eventually exceeds the levels at the edges. These trends are important for properly predicting the onset of membrane rupture. That is, at high pressure differences, the rupture is initiated at the membrane edges, while at low pressure difference, the rupture is initiated from the central point. These points will be further discussed in the next subsection.

6.2. Membrane Rupture

Different modes of membrane rupture are discussed in this section. As mentioned before, membrane failure is simulated using the Lemaitre-Chaboche method. First, the failure initiation time for different pressure differences is investigated. Then, two different rupturing modes are introduced. One objective of this study is to predict the minimum pressure difference that leads to the membrane rupture. In addition, the failure time as a function of pressure difference is also evaluated.

The failure time of several membranes are presented in Figure 9. This figure illustrates the first sign of rupture for several membrane thicknesses and pressure differences. The criterion utilized to interpret the first sign of membrane rupture is the failure of all 20 Gaussian points of a single membrane element. The simulation data in Figure 9 are fitted by power laws that are shown by solid lines in this figure. For a given membrane thickness, it is seen that there is a lower limit for pressure difference below which the rupture does not occur. The time required to reach the minimum pressure difference needed for rupture, appears to be roughly the same for all membranes with different thickness. We refer to the time for failure under minimum pressure difference as the “critical time”, which is roughly independent of membrane thickness but perhaps is a function of the size and shape of the membrane. For the present membrane, the critical time is about 0.21 ms (varies between 0.204 to 0.218) and is illustrated in Figure 9 by a dashed line. Thus, if the membrane does not fail before the critical time, rupturing will not occur, regardless of the thickness of the membrane or the exerted pressure difference.

Figure 9 shows that the failure time for low values of pressure difference is highly sensitive to the value of membrane thickness. On the other hand, for high pressure differences, the failure time is not as sensitive to changes in the membrane thickness. In this context, membrane rupture initiation can be categorized into slow rupture (thickness/pressure-sensitive) and fast rupture (thickness/pressure-less sensitive) and also maybe a mixed mode between the two, namely: medium speed rupture (mildly sensitive to thickness/pressure). Membrane rupture evolution and the amount of failure over the surface of the membrane are addressed next. We first discuss slow rupturing process and then consider the fast one.

6.2.1. Slow Rupture

Figure 10 shows top-view of membrane failure by slow rupture for a pressure difference of 700 MPa (at six time instants). It is seen that for slow rupture mode, the equivalent plastic strain profiles vary smoothly. The pattern of plastic strain distribution at time 30 $\mathsf{\mu }$s (Figure 10b) is similar to the last strain patterns of the fast rupturing mode before rupturing occurs (see Figure 11b). However, in the slow mode, strain magnitudes do not reach their critical values and thus, rupturing does not occur. Consequently, peak strain has sufficient time to develop towards central part of the membrane. In other word, the effective destructive mechanism in this mode is the pressure force magnitude (and not the pressure pulse). As shown in Figure 10c, at time 80 $\mathsf{\mu }$s, the central parts of the membrane experience the growth of plastic strain. As the intensity of strain in the bulk of the membrane increases with time, the magnitude of strain in the central areas becomes larger than the values of the stain at edges that are relaxing time at 120 $\mathsf{\mu }$s (Figure 10d). Then, at 190 $\mathsf{\mu }$s, the first element at the center of the membrane fails (Figure 10e), causing its neighboring elements to fail and reach the final rupturing state at 250 $\mathsf{\mu }$s (Figure 10f).

6.2.2. Fast Rupture

As mentioned before, the fast rupturing mode takes place for high pressure differences. In order to model fast rupturing, a pressure difference of 4 GPa is imposed suddenly on a 1 mm thick membrane. Figure 11 shows the initiation, progress and complete destruction of the membrane by the fast rupturing model at six different times. The non-smoothness of equivalent plastic strain contours particularly for the initial times is clearly seen in this figure, which is due to the highly unsteady nature of this kind of rupturing. The first signs of fast failure is observed in about 5 $\mathsf{\mu }$s after imposing the pressure difference (Figure 11b). Figure 11c illustrates that the edges of the membrane are the first areas that the rupturing is initiated and progresses. The displacement of the central parts of the membrane in this mode is not significant. That is, due to the fast nature of this rupturing mode, the membrane does not have sufficient time to propagate the strain to its central area. At 12 $\mathsf{\mu }$s, the membrane remains attached to the duct walls only through its corners (Figure 11d). As the rupturing process continuing, as shown in Figure 11e, at 18 $\mathsf{\mu }$s the membrane completely detached from the duct walls with only some small pieces of the membrane which have more constraints remaining on the duct walls.

To summarize the findings of this section, different features of the membrane rupturing modes are presented in

Table 5. Features of fast and slow rupturing modes.

Feature	Fast Rupturing	Slow Rupturing
Initiation Time	Order of $\mathsf{\mu }$s	Order of ms
Initiation Area	Edges	Center
Thickness/Pressure Sensitivity	low sensitivity	Sensitive
Destruction Level	Smaller Destruction per Area with Complete Separation from Channel Walls	High Destruction per Area without Separation from Channel Walls
Destructive Mechanism	Pressure pulse	Steady Pressure Field Forces

. It is clear that magnitude of the pressure difference is a major factor that controls the rupturing mode. As noted before, there is a mixed rupturing mode which occurs between the slow and fast modes. However, because of the transitional nature of the mixed mode, its detailed study is left for a future work.

7. Concluding Remarks

This study utilized the SPH and FEM for modeling of fluid and solid phases, respectively. Based on the presented results, the following conclusions are drawn:

• The maximum deflection of the membrane occurs at its center. The deflection is also a linear function of the applied pressure difference.
• The time required to reach the final deformation of the membrane is only a function of the applied pressure difference and the effect of membrane thickness is negligible.
• The maximum plastic strain on the membrane, which is first seen at the constrained edges, move to the center of the membrane with time.
• The minimum pressure difference required for initiation of membrane rupture varies with membrane thickness but the rupturing time (referred to as critical time) is roughly the same for all membrane thicknesses.
• Simulations of the rupturing time as a function of pressure difference suggest three different rupturing modes, namely, thickness/pressure-sensitive mode (slow rupturing), mild thickness/pressure-sensitive mode (medium speed rupturing) and low thickness/pressure-sensitive mode (fast rupturing).
• Slow rupturing causes cracks in the central parts of the membrane that spreads to other areas.
• In fast rupturing, failure first occurs at the edges of the membrane and quickly detaches the membrane from the duct walls.

Given the potential of the SPH method for simulating FSI and multiphase flows, as well as its applicability to handling different types of Non-Newtonian flows, for future work, we plan to study membrane rupturing in Non-Newtonian flows and investigate the effect of shear thinning and shear thickening as well as the viscoelastic effects on the process, particularly on the onset of rupture.