Introduction

Engineering materials usually include a certain level of porosity, which is ignored most of the time in performing mechanical simulations due to their negligible size and density. However, it is unavoidable to consider such pores in certain circumstances so that these voids/pores may have prominent effects on the mechanical performance of the structures in macro scale [1, 2]. What is more, the porous micro-structures are desirable in special applications, e.g., wastewater treatment systems utilize porous filters to remove pollutants and extract usable water.

It is of great importance to characterize the mechanical behaviour of porous structures in wastewater treatment systems as they may be subjected to high working pressures up to 10 bars in an ultra-filtration process [3]. Wang et al. [4] presented an extended literature report on the evaluation of mechanical properties for desalination and wastewater treatment membranes. de Wit et al. [5] discussed the candidate methods, namely, three-point bending, four-point bending and diametrical compression tests to assess the mechanical strength of ceramic hollow fibre membranes. The same authors also examined the impact of the production methods on the mechanical performance of alumina porous hollow fibre membranes [6]. Lee et al. [7] investigated potential use of alumina hollow fibre membranes for wastewater treatment in terms of different micro-structures, which may affect the permeability and mechanical performance of membranes.

Aside to the specific works addressing the mechanical response of porous membrane materials, the simplified formulations can be adopted for estimating the porosity dependence of mechanical properties of solids in general. Several decades ago, elastic constants of solids containing spherical voids for low porosity ratios were estimated by Mackenzie [8], while the effective material properties of solids with high porosity were predicted by Gibson and Ashby formulae [9]. Coble and Kingery [10] carried out extensive experimental works to establish the relationship between porosity, temperature and physical properties, i.e., strength, elastic modulus, rigidity and coefficient of thermal expansion in porous ceramics. Brown et al. [11] estimated the mechanical strength of brittle porous materials by a theoretical approach by taking both pore density and geometry into account. Walsh et al. [12] reported the impacts of porosity on the dynamic response of the glass under compressive loads. Phani and Niyogi [13] proposed a semi-empirical formulation to predict the elastic modulus of porous brittle solids, which was tested up to porosity level of 0.5. Boccaccini et al. [14] estimated effective Young’s modulus of porous media by volume fraction ratio of closed porosity and their micro-structural characteristics, such as shape and orientation. The predicted results were compared with the experimental measurements for metals, ceramics and glasses. Ji and Gu [15] estimated mechanical properties of porous media by taking them as a special case of two-phase composites, in which pores are dispersed within a solid body. Accordingly, they employed generalized mixture rule and compared the performance of their estimations with other approximate methods. Recently, Manoylov et al. [16] improved the porosity–elastic modulus relations for isolated spherical holes proposed by Vavakin and Salganik [17, 18]. Furthermore, these expressions were modified for merged and open pores.

The above mentioned works have mostly employed theoretical or empirical approximations to construct a relationship between the mechanical and physical properties of porous solids instead of direct simulation. However, peridynamics (PD) [19] offer a great potential to assess the mechanical response of porous bodies. In this study, we adopt Ordinary State-Based (OSB)—PD [20] reformulated by Madenci and Oterkus [21]. Unlike the bond-based PD [22], the OSB-PD formulation eliminates the restrictions on Poisson’s ratio of the material by defining the force state of a material point (particle) in terms of deformation states of other material points within the neighbourhood, so called horizon, in addition to the extension of the PD bonds between two particles. In the recent decade, there has been a remarkable progress in the PD and relevant nonlocal methods. Ren et al. [23] proposed dual horizon PD, which eliminates the ghost force effect due to the variable horizon sizes. Dual horizon PD naturally allows to handle complex problems with relatively high computational efficiency by means of reducing the horizon sizes as well as the particle distance in desired locations, e.g., crack tip area, while keeping the particle distance large in other parts of a model. Madenci et al. [24] proposed the PD differential operator to obtain partial derivatives in the form of nonlocal integrations. Likewise, a nonlocal operator method for solving partial differential equations was proposed by Rabczuk and his colleagues [25, 26]. Ren et al. [27] later presented a higher order nonlocal operator method, which brings additional advantages to the original form of nonlocal operator method and it can be adopted to obtain partial derivatives of higher orders.

The wave propagation in a solid medium can be assumed as an indication of the material stiffness [12]. The apparent material properties of a porous medium is estimated by simplified approaches most of the time [11, 13,14,15]. These approaches rely on certain assumptions and simplifications, and basically do not take into account the variability of the porosity ratios and randomness of the micro-structure. PD is capable of representing the complicated micro-structures by means of the pre-damage [28,29,30]; we thus aim to capture dynamic wave propagation in porous media, which allows to characterize material properties in a more reliable manner by considering non-uniform and random micro-structure. In this respect, we utilize the assets of OSB-PD formulation to model porosity as well as the micro-cracks in a porous medium to assess their influence on the wave propagation speed for the first time in the literature. The ultimate motivation of the present work is to establish a numerical framework to evaluate the mechanical response of porous solids that can be used as water treatment membranes.

Instead of modelling every single detail of pores, these pores can be imposed by means of the pre-defined damage into the body within the PD perspective [28, 29]. Chen and Bobaru [28] proposed concentration dependent damage technique to model pitting corrosions in the PD framework. Then, this technique was implemented for porosity modelling by intermediate homogenization algorithm [29]. Chen et al. [29] employed the bond-based PD with conical micro-modulus function to simulate elastic wave propagation in porous media under an impact load. Apparent modulus of the material was then predicted by means of the wave propagation speed. Furthermore, crack propagation simulation in porous media was also performed by Ref. [29]. De Meo et al. [31, 32] developed a PD model, based on the concentration dependent damage method [28], to simulate the crack onset and propagation from pitting corrosions, and the model was implemented in a commercial software framework. Oterkus et al. [33] developed a coupled model for fluid driven fracture in porous media. Xia et al. [34, 35] proposed a homogenization scheme using a representative volume element approach for periodically micro-structured materials, e.g., composites, porous media so as to predict effective material properties from the PD displacement gradient tensor. Li et al. [36] investigated the influence of porosity on the fracture behaviour of brittle granular materials. It was reported that as inter-granular strength increases, the influence of porosity declines considerably. Karpenko et al. [37] recently addressed the porosity effects on the fatigue crack nucleation of additively manufactured alloys using the bond-based PD.

It has been reported so far that porosity [12, 29] and small-sized defects [38,39,40,41] have an suppressing impact on the dynamic response in the form of waves and crack propagation. To the best knowledge of the authors, the combination of porosity and micro-cracks has not been studied comprehensively yet. Accordingly, the present work is aimed to fill this gap in the literature.

1 A brief recap of Peridynamic Theory

1.1 Fundamental equations

The basic formulation utilized in this paper is the reformulated OSB-PD by Madenci and Oterkus [21]. The equation of motion for a particle, at position vector \({\mathbf {x}}\), is expressed as

$$\begin{aligned} \rho ({\mathbf {x}}){\ddot{{\mathbf {u}}}}({\mathbf {x}},t)&=\int _{{{\mathcal {H}}}}\left[ {\mathbf {t}}\left( {\mathbf {u}}^\prime -{\mathbf {u}}, {\mathbf {x}}^\prime -{\mathbf {x}}, t\right) -{\mathbf {t}}^\prime \left( {\mathbf {u}}-{\mathbf {u}}^\prime , {\mathbf {x}}-{\mathbf {x}}^\prime , t\right) \right] d{{\mathcal {H}}}+{\mathbf {b}}({\mathbf {x}},t), \end{aligned}$$
(1)

where \(\rho \) is the mass density at \({\mathbf {x}}\) and the displacement vector is represented by \({\mathbf {u}}\). The force interactions between the particle at \({\mathbf {x}}\) and it’s neighbour at \({\mathbf {x}}^\prime \) are denoted by \({\mathbf {t}}\) and \({\mathbf {t}}^\prime \), respectively. \({{\mathcal {H}}}\) defines the neighbourhood, i.e., horizon. In the OSB-PD theory, the interaction forces are not only related with the relative deformation of particles each other, but also the deformation states of the other particles within the horizon. Accordingly, the force vectors defining the interactions between the particles are aligned with the relative position vector in the deformed configuration, but may have unequal magnitudes. These vectors can be written as follows [21]:

$$\begin{aligned} {\mathbf {t}}\left( {\mathbf {u}}^\prime -{\mathbf {u}}, {\mathbf {x}}^\prime -{\mathbf {x}}, t\right) =\frac{1}{2}C_1\frac{{\mathbf {y}}^\prime -{\mathbf {y}}}{|{\mathbf {y}}^\prime -{\mathbf {y}}|}, \end{aligned}$$
(2)

and

$$\begin{aligned} {\mathbf {t}}\left( {\mathbf {u}}^\prime -{\mathbf {u}}, {\mathbf {x}}^\prime -{\mathbf {x}}, t\right) =-\frac{1}{2}C_2\frac{{\mathbf {y}}^\prime -{\mathbf {y}}}{|{\mathbf {y}}^\prime -{\mathbf {y}}|}. \end{aligned}$$
(3)

In Eqs. (2) and (3), \(C_1\) and \(C_2\) represent parameters depending on the material properties as well as the deformation states. These parameters can be derived from the strain energy functional by neglecting the temperature effects as below:

$$\begin{aligned}&C_1=4w\left[ d\frac{{\mathbf {y}}^\prime -{\mathbf {y}}}{|{\mathbf {y}}^\prime -{\mathbf {y}}|}\cdot \frac{{\mathbf {x}}^\prime -{\mathbf {x}}}{|{\mathbf {x}}^\prime -{\mathbf {x}}|}\left( a\theta _{{\mathbf {x}}}\right) +b\left( |{\mathbf {y}}^\prime -{\mathbf {y}}|-|{\mathbf {x}}^\prime -{\mathbf {x}}| \right) \right] , \end{aligned}$$
(4)
$$\begin{aligned}&C_2=4w\left[ d\frac{{\mathbf {y}}-{\mathbf {y}}^\prime }{|{\mathbf {y}}-{\mathbf {y}}^\prime |}\cdot \frac{{\mathbf {x}}-{\mathbf {x}}^\prime }{|{\mathbf {x}}-{\mathbf {x}}^\prime |}\left( a\theta _{{\mathbf {x}}^\prime }\right) +b\left( |{\mathbf {y}}-{\mathbf {y}}^\prime |-|{\mathbf {x}}-{\mathbf {x}}^\prime | \right) \right] . \end{aligned}$$
(5)

In the above equations, w is the scalar influence function, which is related to the initial distance between the particles. \(\theta _{{\mathbf {x}}}\) and \(\theta _{{\mathbf {x}}^\prime }\) denote the dilations for the particles at \({\mathbf {x}}\) and \({\mathbf {x}}^\prime \), respectively. By the definition of \({\varvec{\xi }}={\mathbf {x}}^\prime -{\mathbf {x}}\) and \({\varvec{\eta }}={\mathbf {y}}^\prime -{\mathbf {y}}\), the unit vectors along the relative position vectors in the deformed and initial configurations are \({\mathbf {m}}={\varvec{\eta }}/|{\varvec{\eta }}|\) and \({\mathbf {n}}={\varvec{\xi }}/|{\varvec{\xi }}|\), respectively. In case of small deformation analysis, the scalar product of these vectors would become \({\mathbf {m}}\cdot {\mathbf {n}}\approx 1.0\). Since the OSB-PD is capable of dealing with the large deflections as well as the plastic deformation problems [42], the original form of the force state expressions are kept throughout the paper. Accordingly, the PD dilatation \(\theta _{\mathbf {x}}\) for a particle at \({\mathbf {x}}\) is expressed as

$$\begin{aligned} \theta _{\mathbf {x}}=d\int _{{\mathcal {H}}}\left( wS{\varvec{\xi }} \cdot {\mathbf {m}}\right) d{{\mathcal {H}}}. \end{aligned}$$
(6)

In a similar manner, the strain energy density for the same particle can be given as

$$\begin{aligned} W_{\mathbf {x}}=a\theta _{\mathbf {x}}^2+b\int _{{\mathcal {H}}}w \left( |{\varvec{\eta }}|-|{\varvec{\xi }}|\right) ^2 d{{\mathcal {H}}}. \end{aligned}$$
(7)

For 2D OSB-PD, the influence function, based on the dimensional analysis, can be obtained as \(w=\delta /|{\varvec{\xi }}|\) [21]. \(\delta \) denotes the horizon size. S is the stretch between two particles within the neighbourhood, which is given as

$$\begin{aligned} S=\dfrac{|{\varvec{\eta }}|-|{\varvec{\xi }}|}{|{\varvec{\xi }}|}. \end{aligned}$$
(8)

The other PD material parameters, a, b and d in Eqs. (4)–(7) are related to the material properties and horizon size, \(\delta \) [21]:

$$\begin{aligned} a=\dfrac{1}{2}\left( \kappa -2\mu \right) , \quad b=\dfrac{6\mu }{h\pi \delta ^4}, \quad d=\dfrac{2}{h\pi \delta ^3}, \end{aligned}$$
(9)

where \(\kappa \) and \(\mu \), respectively, stand for the bulk and shear moduli of the material. h is the thickness of a 2D body. The plane stress condition is employed throughout the work.

1.2 Discrete form of the fundamental equations

Equations (1)–(8) have been defined in the continuous form so far. As PD is a particle based method in essence, these equations can be discretized accordingly. The equation of motion in discrete form is written for a particle (i) as

$$\begin{aligned} \rho _{(i)}{\ddot{{\mathbf {u}}}}_{(i)}=\sum _{j=1}^{{{\mathcal {N}}}_{(i)}}\left[ {\mathbf {t}}_{(i)(j)}-{\mathbf {t}}_{(j)(i)}\right] V_{(j)}+{\varvec{b}}_{(i)}, \end{aligned}$$
(10)

where \({{\mathcal {N}}}_{(i)}\) stands for the number of particles within the neighbourhood of particle (i). \(V_{(j)}\) is the finite volume assigned to particle (j) within the horizon of (i) in the discretized body. Similarly, the force interaction vectors in discrete form can be written as

$$\begin{aligned} {\mathbf {t}}_{(i)(j)}=\frac{1}{2}C_1\frac{{\varvec{\eta }}_{(i)(j)}}{|{\varvec{\eta }}_{(i)(j)}|}, \end{aligned}$$
(11)

and

$$\begin{aligned} {\mathbf {t}}_{(j)(i)}=-\frac{1}{2}C_2\frac{{\varvec{\eta }}_{(i)(j)}}{|{\varvec{\eta }}_{(i)(j)}|}, \end{aligned}$$
(12)

where \({\varvec{\eta }}_{(i)(j)}={\mathbf {y}}_{(j)}-{\mathbf {y}}_{(i)}\). \(C_1\) and \(C_2\) parameters can be expressed in discretized form as below:

$$\begin{aligned} C_1= & {} 4w_{(i)(j)}\Big [d\frac{{\varvec{\eta }}_{(i)(j)}}{|{\varvec{\eta }}_{(i)(j)}|}\cdot \frac{{\varvec{\xi }}_{(i)(j)}}{|{\varvec{\xi }}_{(i)(j)}|}\left( a\theta _{(i)}\right) \nonumber \\&\quad +b\left( |{\varvec{\eta }}_{(i)(j)}|-|{\varvec{\xi }}_{(i)(j)}|\right) \Big ], \end{aligned}$$
(13)

and

$$\begin{aligned} C_2= & {} 4w_{(j)(i)}\Big [d\frac{{\varvec{\eta }}_{(j)(i)}}{|{\varvec{\eta }}_{(j)(i)}|}\cdot \frac{{\varvec{\xi }}_{(j)(i)}}{|{\varvec{\xi }}_{(j)(i)}|}\left( a\theta _{(j)}\right) \nonumber \\&\quad +b\left( |{\varvec{\eta }}_{(j)(i)}|-|{\varvec{\xi }}_{(j)(i)}|\right) \Big ], \end{aligned}$$
(14)

where \({\varvec{\xi }}_{(i)(j)}={\mathbf {x}}_{(j)}-{\mathbf {x}}_{(i)}\). The scalar influence function is discretized as \(w_{(i)(j)}=\delta /|{\varvec{\xi }}_{(i)(j)}|\). The PD dilatation term \(\theta \) for the particle (i) is expressed in discrete form as

$$\begin{aligned} \theta _{(i)}=d\sum _{j=1}^{{{\mathcal {N}}}_{(i)}}\left( w_{(i)(j)}S_{(i)(j)}{\varvec{\xi }}_{(i)(j)} \cdot \frac{{\varvec{\eta }}_{(i)(j)}}{|{\varvec{\eta }}_{(i)(j)}|}\right) V_{(j)}. \end{aligned}$$
(15)

In a similar way, the strain energy density is written in discrete form as

$$\begin{aligned} W_{(i)}=a\theta ^2_{(i)}+b\delta \sum _{j=1}^{{{\mathcal {N}}}_{(i)}}\left( S_{(i)(j)}^2 |{\varvec{\xi }}_{(i)(j)}| \right) V_{(j)}. \end{aligned}$$
(16)

Finally, the stretch between the particles (i) and (j) in the discrete form is given below:

$$\begin{aligned} S_{(i)(j)}=\dfrac{|{\varvec{\eta }}_{(i)(j)}|-|{\varvec{\xi }}_{(i)(j)}|}{|{\varvec{\xi }}_{(i)(j)}|}. \end{aligned}$$
(17)

In the numerical solution of the problems, both surface correction and volume correction procedures are considered. Since the present formulation is OSB-PD, the dilatation and deviatoric deformation terms are decoupled. In this respect, the dilatation term is corrected for isotropic expansion case, while the strain energy terms are corrected for uniform shear deformation. The details regarding both the surface and volume correction procedures have been provided by Ref. [21].

2 Porosity implementation

Porous micro-structure can be implemented in the PD through several ways. The first one is the direct modelling of pores [37]. This approach can be effective for low porosity ratios. However, if the variable porosity is concerned with high ratios, direct modelling might be challenging. Second, the homogenization approach, in which the material is assumed to be homogeneous with reduced material modulus [34, 35], can be employed for the regular micro-structure cases. This technique, however, omits the random or complex micro-structure effects. Finally, porosity can be implemented as pre-defined damage by the intermediate homogenization method as proposed by Chen et al. [29]. This approach can be easily tailored for a micro-structure including complex and variable porosities. In the present work, the main focus is uniform porosity ratio, yet it involves a certain level of randomness in terms of the number and orientation of broken bonds for each particle.

Fig. 1
figure 1

Porosity representation by means of PD damage for the low porosity ratios: a \(\delta =4dx\), b \(\delta =8dx\)

Full size image

The intermediate homogenization approach [29] is adopted in this work. In the PD framework, the damage parameter for the particle (i) can be estimated by the ratio of the broken bonds to the total number of bonds associated with the particle, i.e., \(d_{(i)}={{\mathcal {N}}}_b/{{\mathcal {N}}}_{(i)}\). The number of broken bonds associated with the particle (i) is represented by \({{\mathcal {N}}}_b\), while \({{\mathcal {N}}}_{(i)}\) stands for the total number of bonds connected to the particle (i). Based on the definition of damage parameter, porosity can be represented by the pre-damage index \(d_{\varphi (i)}=\varphi _{(i)}/\varphi _c\). The critical porosity level, beyond which the materials can exist only in the form of suspension, is \(\varphi _c\). This value is taken as 1.0 (\(\varphi _c=1.0\)) in the present work, which is consistent with the work of Chen et al. [29].

Fig. 2
figure 2

Porosity representation by means of PD damage for the high porosity ratios: a \(\delta =4dx\), b \(\delta =8dx\)

Full size image

The numerical implementation of uniform porosity by the intermediate homogenization technique is quite straightforward. The algorithm for generating porosity for the particle (i) is simply defined as

figure a

Detailed description of the intermediate homogenization method can be found in Ref. [29]. Being dx discretization size and m horizon factor, the horizon size can be expressed as \(\delta =mdx\). For efficient modelling of the porosity, horizon factor m must be adjusted to cover enough number of bonds. For a non-porous body and the low porosity ratios, \(m=4\) can be taken. As for the higher porosity ratios, e.g., \(d_{\varphi }\ge 0.5\), m must be set as a large number, e.g., \(m=8\). These values were adopted by Chen et al. [29] as well. The above algorithm is applied for generating porous media considering both \(m=4\) and \(m=8\) cases for the same discretization size, dx.

As can be seen from Figs. 1 and 2, it is possible to generate porous media by employing both small (\(m=4\)) and large horizon parameters (\(m=8\)). However, as the porosity ratio increases, the numerical model with small horizon parameter poorly represents the exact porosity. For instance, let’s compare \(d_{\varphi }=0.1\) cases in Fig. 1; the obtained porosity representations can satisfy the average value of 0.100 with the standard deviation of 0.047 and 0.023 for small and large horizon parameters, respectively. As for the \(d_{\varphi }=0.7\) cases in Fig. 2, average porosity value obtained by the small horizon parameter is 0.699 with the standard deviation of 0.071, while these values are 0.701 and 0.035, respectively, for the larger horizon parameter case. In summary, even though small horizon parameters can represent the porous micro-structure in some extent, the larger horizon parameter implementation is much meaningful to represent the porosity in terms of the PD damage.

3 Dynamic simulation of bodies with constant porosity

Our numerical implementation is first verified by the reference results presented by Chen et al. [29]. Accordingly, case studies are adopted from the reference work. Problem set up is schematically illustrated in Fig. 3. Both the material properties, loading condition and PD discretization are provided in Fig. 3. The impulse load is applied during the first 5 \({\upmu }\)s, and the velocity waves propagating through the porous body are captured and average wave speeds are calculated. Total simulation time is considered as \(t_{\mathrm{total}}=200\) \({\upmu }\)s. Non-porous body as well as the constant porosity cases, \(d_{\varphi }=0.1\), 0.3, 0.5 and 0.7, are considered as the numerical examples. Time step size is adopted as \(t_{\mathrm{inc}}=0.25\) \({\upmu }\)s, which is far below the allowable limit, 1.55 \(\upmu \text{s}\), based on Refs. [21, 22].

Fig. 3
figure 3

Problem setup

Full size image

3.1 Wave propagation results

The vertical velocity (\(v_y\)) wave patterns for the non-porous and low porosity bodies are obtained and presented in Fig. 4. The wave tip locations at the instance of \(t=40\) \(\upmu \text{s}\) seem to be very close to each other, the differences cannot be distinguished at the first glance; however, the influence of porosity on the location of wave tip can be distinguished easily at the instance of \(t=200\) \(\upmu \text{s}\). In addition to the wave tip locations, magnitude of the velocity bands are also suppressed by the porosity.

Fig. 4
figure 4

Vertical velocity contours for the non-porous body and the low porosity ratios: a \(t=40\) \({\upmu }\)s, b \(t=200\) \({\upmu }\)s

Full size image

Velocity wave contours for the high porosity ratios are demonstrated in Fig. 5. As it is obvious from the figure, the high porosity influences the wave propagation behaviour remarkably. Not only the wave tip location is altered by the porous micro-structure, but also wave patterns, i.e., velocity magnitude and number of waves are reduced by the high porosity ratios.

Fig. 5
figure 5

Vertical velocity contours for the high porosity ratios: a \(t=40\) \({\upmu }\)s, b \(t=200\) \({\upmu }\)s

Full size image

This behaviour is due to the fact that the impact load is compensated mainly by the deformation of the particles in the vicinity of the loading area. The disconnected bonds as a result of the porous micro-structure prevent the transmission of physical deformations, which in turn excessive wave oscillations can be observed in the localized zone, see Fig. 5b.

Fig. 6
figure 6

Wave tip locations along the mid-span (\(x=L/2\))

Full size image

The wave tip locations are captured along the mid-span (\(x=L/2\)) for the instances of \(t=40\), 100, 160 and 200 \({\upmu }\)s, and are given in Fig. 6. This figure sheds light on the almost linear change of the wave tip locations along the mid-span, which allows to obtain average wave speeds efficiently. Moreover, it is obvious that for the non-porous and low porosity (\(d_{\varphi }=0.1\)) cases, impact of the horizon parameter on the average wave speed is insignificant.

Fig. 7
figure 7

Wave patterns of the non-porous body along the mid-span (\(x=L/2\)) for the horizon parameter, \(m=8\) at different instances: a \(t=40\) \({\upmu }\)s, (b) \(t=200\) \({\upmu }\)s

Full size image

Based on the insights from Fig. 6, average wave speeds are obtained for the interval between \(t=40\) \(\mu\text{s}\) and \(t=200\) \(\upmu \text{s}\). Therefore, as to explain the procedure more clearly, Fig. 7 gives the wave patterns along the mid-span of the non-porous body with \(m=8\) for \(t=40\) \(\upmu \text{s}\) and \(t=200\) \(\upmu \text{s}\) instances. In Fig. 7, being \(\Delta y\) the distance between the wave tips, the average speed is obtained as \(V_{\mathrm{avg}}=\Delta y/\Delta t\) with \(\Delta t=160\) \(\upmu \text{s}\).

Employing the previously defined approach, the average wave speeds are evaluated and compared with the reference results presented by Chen et al. [29] for the different porosity ratios in Fig. 8. In this figure, both present and reference results are provided for horizon factor, \(m=8\). It must be noted that the reference results are digitized from the graphs in [29].

Overall agreement between the results presented in Fig. 8 is satisfactory. As the porosity ratio increases, the difference between the wave speeds obtained by the present method and the reference values tends to increase slightly. It is worth noting that the reference results have been obtained by the bond-based PD with the conical micro-modulus function. Furthermore, the reference paper [29] does not report whether surface and volume correction techniques were implemented or not. In summary, these small differences between the present and reference results may be attributed to above mentioned factors.

It is a well-known fact that wave propagation in a solid body can be affiliated with the characterization of the micro-structure and approximate prediction of apparent material properties. In this regard, the reliable modelling and simulation of wave propagation in porous solids is rather essential. Ravi-Chandar [43] briefly presented the relationship between the bulk waves and the material properties within the linear elastodynamic perspective.

Fig. 8
figure 8

Average wave speeds obtained by the present OSB-PD implementation and the reference work [29] for the horizon parameter, \(m=8.0\)

Full size image

3.2 Parametric sensitivity analysis

After verifying our implementation with a reference work and investigating the influence of the porosity ratios on the wave propagation speed in a porous medium, a series of parametric work is carried out to assess the sensitivity of the wave patterns and prospectives of porous media with respect to the PD parameters, i.e., the horizon size, \(\delta \), and the horizon factor m. As being the lowest and highest porosity ratios, \(d_{\varphi }=0.1\) and 0.7 are considered for parametric studies, respectively. For the purpose of conducting a comprehensive parametric study, the particle distance is varied as \(dx=L/150, L/200\) and L/250. In addition, the horizon parameter, m, is varied as 4, 8 and 12 for each discretization size, dx. With these assumptions, the horizon sizes for \(m=4\) become \(\delta =26.67\), 20.0 and 16.0 mm, respectively. Likewise, the horizon sizes for \(m=8\) are \(\delta =53.33\), 40.0 and 32.0 mm, and for \(m=12\) they are \(\delta =80.0\), 60.0 and 48.0 mm, respectively. In total, 18 cases are considered. It must be noted that some of the above combinations have been studied in the previous section; these cases are: \(dx=L/200\) with \(m=4\) and 8 for \(d_{\varphi }=0.1\); and \(dx=L/200\) with \(m=8\) for \(d_{\varphi }=0.7\) porosity ratio. For the sake of completeness, the results for these cases will be included in the following as well.

First, the vertical velocity contours at \(t=200\) \(\upmu \text{s}\) are captured. The contours for the lowest porosity ratio are given in Fig. 9. This figure suggests that the location of the wave tip is not influenced considerably by the horizon sizes as long as the number of particles within the horizon is the same, i.e., constant m values. However, it is clear that the wave numbers (inversely proportional to the wave length for a unit length) decline as the horizon size becomes smaller for the same m values. What is more, the wave patterns become more coherent as the discretization becomes finer.

Fig. 9
figure 9

Vertical velocity contours for \(d_{\varphi }=0.1\) at 200 \(\mu \)s: a \(dx=L/150\), b \(dx=L/200\), c \(dx=L/250\)

Full size image

Let us consider the cases with the same particle discretization but varying m values. As obviously seen in Fig. 9, the wave numbers decline as increase of the m values for the same dx assumptions. It is also worth noting that the velocity oscillations in the vicinity of the loading region can be reduced significantly by increasing the number of particles within the horizon. The wave patterns become more coherent for the smaller values of m, which is mainly because of the limited (short range) interactions between the particles. Overall, Fig. 9 indicates that the wave length in solids can be increased substantially by the increase of the horizon size in the PD perspective, also see Eq. (9) for the relationship between the PD material constants and the horizon size.

The vertical velocity contours for the highest porosity ratio are given in Fig.  10. The similar interpretations can be obtained from the wave patterns in this figure. The most obvious difference between the wave patterns of the lowest and highest porosity ratios is observed for \(m=4\) case. The excessive oscillations of the velocity component are apparent in Fig. 10 for the mentioned case. The main reason for the such unstable behaviour is the limited interactions between the particles within the same horizon, as the number of particles within the horizon is small for \(m=4\) case. When many of the PD bonds are broken to represent the porosity by means of pre-damage, the remaining PD bonds can not recover the necessary PD interaction forces under the suddenly applied impact loading.

Fig. 10
figure 10

Vertical velocity contours for \(d_{\varphi }=0.7\) at 200 \(\mu \)s: a \(dx=L/150\), b \(dx=L/200\), c \(dx=L/250\)

Full size image

The wave patterns for \(m=4\) in Fig. 10 also suggest that the wave tip may not propagate along the mid-span. Since the bonds are randomly broken to generate porosity, the directions of the intact bonds may play an influential role in the direction of wave propagation for lower m values. This effect can be avoided by increasing the number of particles within the horizon, as shown in Fig. 10 for \(m=8\) and 12 cases.

Fig. 11
figure 11

Wave-crest locations in the vertical direction for the lowest (\(d_{\varphi }=0.1\)) and the highest (\(d_{\varphi }=0.7\)) porosity ratios: a \(dx=L/150\), b \(dx=L/200\), c \(dx=L/250\)

Full size image

To examine the influence of assumed PD parameters on the wave propagation speed quantitatively, the wave tip locations are captured for \(t=40\), 100, 160 and 200 \(\mu \)s, then presented in Fig. 11 for the lowest and highest porosity ratios. For the lowest porosity case, the impacts of the discretization size as well as the horizon parameter are found to be very limited. The average wave speeds for these cases are around 4100–4200 m/s. As the numerical discretization becomes finer, their impact is lesser. On the other hand, the impacts of the particle discretization and the horizon parameters become visible in terms of the wave tip locations for the highest porosity case. This is mainly because of the limited number of interactions between the particles for higher porosity ratios. For \(d_{\varphi }=0.7\) case, one can predict the average speed excluding Fig. 11a. The average speeds for the remaining discretization sizes are evaluated between 2200 and 2500 m/s. As can be inferred from Fig. 11 and the average wave speeds, the higher porosity ratios are more sensitive to the PD parameters. For reliable modelling of higher porosity ratios in terms of pre-damage in PD framework, the number of particles within horizon must be sufficiently large, e.g., \(m=8\) or higher. Otherwise, unstable velocity fluctuations may occur.

4 Wave propagation in solids with constant porosity including micro-cracks

Once our porosity implementation has been validated through several numerical examples, and has been discussed in detail; we can proceed for further applications. In this section, random and regular micro-crack configurations are going to be introduced into the porous media, which were considered in the previous section.

As for the modelling of micro-cracks, discretization size, dx, and horizon parameter, m, have to be adjusted. Hence, the discretization size is set as \(dx=L/250\). Numerical trials suggest that if the horizon size becomes larger than the micro-crack size, the micro-cracks cannot be modelled as a straight line crack, since they become like a damage zone. Therefore, the horizon factor is taken as \(m=4\) in micro-crack models by keeping the porosity ratios small, i.e., \(d_{\varphi }=0.1\), 0.2 and 0.3.

4.1 Random micro-crack configuration

In random micro-crack configurations, both crack location and orientation are set arbitrarily, while the micro-crack size, l, is set to be varied as \(l=6\mathtt {\sim }10\cdot dx\). The number of micro-cracks is considered as \(nc=20\) and 50. Figure 12 demonstrates the random micro-crack configurations for the various porosity ratios.

Fig. 12
figure 12

Random micro-crack configurations: a \(d_{\varphi }=0.1\), b \(d_{\varphi }=0.2\), c \(d_{\varphi }=0.3\)

Full size image

As it can be seen from Fig. 12, the micro-cracks are allocated throughout the 2D body except for the boundary regions. There are no other limitations on either location or orientation of the cracks. The micro-cracks also interact, if they intersect to each other. Then, the same loading condition with the previous section is implemented for the specimens involving micro-cracks. The material properties as well as the solution procedure are identical with the conditions assigned in the previous section. Since we adopt smaller horizon parameter, i.e, \(m=4\), lower porosity ratios, \(d_{\varphi }=0.1\), 0.2 and 0.3, are considered here. The vertical velocity patterns are obtained and presented in Figs. 13 and 14 for \(nc=20\) and 50 cases, respectively.

Fig. 13
figure 13

Vertical velocity contours for random micro-crack configurations (\(nc=20\)): a non-porous, b \(d_{\varphi }=0.1\), c \(d_{\varphi }=0.2\), d \(d_{\varphi }=0.3\)

Full size image
Fig. 14
figure 14

Vertical velocity contours for random micro-crack configurations (\(nc=50\)): a non-porous, b \(d_{\varphi }=0.1\), c \(d_{\varphi }=0.2\), d \(d_{\varphi }=0.3\)

Full size image

It must be noted that only the number of micro-cracks are taken as constant parameters; however, their size, location and orientation are unique to each run of the simulation. Even though this implementation causes some level of uncertainty, it is still possible to make judgement on the wave patterns. The micro-cracks affect the coherence of the wave patterns, while the wave-tip locations do not change significantly by the existence of micro-cracks. The most obvious behaviour, which can be deduced from Figs. 13 and 14, is that the energy introduced by the impulse loading is mainly absorbed by the excessive dynamic response in the vicinity of the loading area, see dark blue and red areas in Figs. 13 and 14. This behaviour may be explained by two phenomena: (1) the reduced horizon parameter is the main reason, which confines the interactions between the particles to a localized area, which in turn impulse cannot be easily transmitted to remote particles; and the PD forces cannot be recovered by the remaining intact bonds to compensate the impact load, (2) the waves are inhibited to the vicinity of the loading area partly because of the micro-crack barriers.

Fig. 15
figure 15

Wave speeds obtained along the mid-span (\(x=L/2\)) for random micro-crack configurations

Full size image

The average wave speeds for vertical velocity contours are calculated with the same procedure described in the previous section. Then, these values are plotted with respect to the porosity ratios for porous bodies with and without micro-cracks in Fig. 15. As discussed earlier, the micro-cracks do not influence the wave propagation speed significantly but their magnitudes.

4.2 Regular micro-crack configuration

In the previous part, it has been shown that the random micro-cracks can reduce the wave propagation speed just marginally. However, their influence on the velocity magnitude and the wave coherence is significant. Their impact is mainly governed by the orientation of micro-cracks. In this part, we, therefore, decided to investigate the influence of micro-cracks in regular patterns, which are allocated in equally spaced manner and certain orientations. The micro-cracks are allocated in a regular pattern on the central region of the 2D body. Number of crack arrays in x- and y-directions is varied, and certain crack angles, \(0^{\circ }\), \(45^{\circ }\) and \(90^{\circ }\) with respect to x-axis are considered. The crack configurations to be considered are: \(5 \times 1\), \(5 \times 5\), \(10 \times 1\) and \(10 \times 10\). In these configurations, the first terms are the number of micro-cracks in x-direction, while the second terms stand for the number of micro-cracks in y-direction. Figure 16 illustrates a porous body with regular micro-crack configurations.

Fig. 16
figure 16

Regular micro-crack configurations: a \(5 \times 1\), b \(5 \times 5\), c \(10 \times 1\), d \(10 \times 10\)

Full size image

In Fig. 16, only the case for crack orientation angle, \(\theta \)=\(0^{\circ }\) is illustrated for saving the space; however, we consider other orientations, \(45^{\circ }\) and \(90^{\circ }\), in the computations as well.

We employ the same modelling techniques implemented in the random micro-crack configurations. Figure 17 demonstrates the wave patterns for the non-porous body with various crack configurations and orientations. It is the most obvious deduction from the figure that the velocity waves are transmitted to entire body by means of the intact bonds representing the non-porous micro-structure. In addition, the waves are reflected back when the crack orientation is \(\theta \)=\(0^{\circ }\) for all configurations. Reflected waves are also seen in \(\theta \)=\(45^{\circ }\) case. The wave patterns are mostly influenced by the number of crack arrays in y-direction for \(\theta \)=\(0^{\circ }\) and \(45^{\circ }\) cases. If the cracks are aligned with the impulse direction, the cracks play an accelerating role for the propagating wave, especially in the \(10 \times 10\) configurations.

Fig. 17
figure 17

Vertical velocity contours at \(t=200\) \(\upmu \text{s}\) for regular micro-crack configurations in non-porous media: a \(5 \times 1\), b \(5 \times 5\), c \(10 \times 1\), d \(10 \times 10\)

Full size image

We present the wave patterns for the porosity ratio of \(d_{\varphi }=0.3\) case including regular micro-cracks in Fig. 18. This figure clearly illustrates effect of the smaller horizon parameter in the form of localized velocities with very high magnitudes in the vicinity of the loading region as similar to the case of random micro-cracks in Figs. 13 and 14. Aside from the localized velocities, wave propagation is suppressed along the mid-span by the existence of micro-cracks depending on the number of crack arrays in y-direction and their orientation. When the number of crack arrays in y-direction is one, their influence on the wave tip location can be neglected regardless of the crack orientation, see Fig. 18a, c. On the other hand, wave propagation is suppressed substantially along the mid-span for \(\theta \)=\(0^{\circ }\) case with \(5 \times 5\) and \(10 \times 10\) configurations. If the orientation of micro-cracks is set as \(\theta \)=\(45^{\circ }\), the wave shapes are disturbed in a way of asymmetrical pattern.

Fig. 18
figure 18

Vertical velocity contours at \(t=200\) \(\upmu \text{s}\) for regular micro-crack configurations in porous media (\(d_{\varphi }=0.3\)): a \(5 \times 1\), b \(5 \times 5\), c \(10 \times 1\), d \(10 \times 10\)

Full size image
Fig. 19
figure 19

Average wave speeds for regular micro-crack configurations: a non-porous, b \(d_{\varphi }\)=0.1, c \(d_{\varphi }\)=0.2, d \(d_{\varphi }\)=0.3

Full size image

By the previously defined approach for measuring the average wave speeds, influence of the regular micro-cracks on the wave propagation speed is investigated. It must be noted that the wave tips are traced along the mid-span and their locations are captured at the instances of \(t=40\) \(\upmu \text{s}\) and \(t=200\) \(\upmu \text{s}\). These wave speeds are obtained for various configurations of micro-cracks in certain orientations and are presented in Fig. 19. At the first glance, it can be noticed that the regular micro-cracks alter the wave propagation speeds for both non-porous and porous bodies in a similar way. As already been reported in the previous section, wave propagations become slower with the porous micro-structure. On the other hand, both non-cracked body, and the micro-cracked bodies with \(5 \times 1\) and \(10 \times 1\) configurations behave similarly in terms of the wave patterns and wave propagation speeds. In case of the micro-crack configurations of \(5 \times 5\) and \(10 \times 10\), the wave propagation speed drops considerably, if the cracks are aligned perpendicularly to the impulse direction (\(\theta \)=\(0^{\circ }\)). If the cracks are aligned with the impulse direction (\(\theta \)=\(90^{\circ }\)), the wave propagation speeds become much closer to the non-cracked case.

5 Concluding remarks

Porous micro-structure has been implemented through OSB-PD formulation using the intermediate homogenization approach. First, the wave propagation problem in porous media has been addressed and it has been observed that as the porosity increases, the wave propagation speed drops. The present results were verified against the available reference works and a good agreement was achieved between them. Parametric sensitivity analysis with respect to discretization size and horizon parameters have been carried out. It was found that the higher porosity ratios, the more sensitive to the PD parameters. Once the wave propagation problem has been addressed appropriately, the micro-cracks have been introduced into the porous body then. The micro-cracks have been considered in the random and regular configurations. In summary, the random micro-cracks have altered the wave patterns but the wave propagation speeds have remained more or less the same. However, the regular micro-cracks have had significant influences on the wave propagation speed and patterns depending on the crack orientation and the number of crack arrays in y-direction. In conclusion, the present OSB-PD implementation holds great potential in simulating the mechanical behaviour of porous bodies. Since the main aim of the present work is to examine the wave propagation in porous media, which is relevant to characterization of mechanical response of porous materials; evaluation of the crack propagation in porous materials remains as a future work.