The Use of Foot of the Green-Lipped Mussel is <em>Perna viridis</em> as an Alternative Method to Reduce the Gender Effect on the Bioaccumulation of Cu and Zn in the Mussel

Asghar Farhadi; Asghar Farhadi

ISSN: 2641-3078

Annals of Limnology and Oceanography

Research Article Open aaacess Peer-Reviewed

The Use of Foot of the Green-Lipped Mussel is Perna viridis as an Alternative Method to Reduce the Gender Effect on the Bioaccumulation of Cu and Zn in the Mussel

Asghar Farhadi*

Author and article information

Young Researchers and Elite Club, Shiraz Branch, Islamic Azad University, Shiraz, Iran

*Corresponding author: Asghar Farhadi, Young Researchers and Elite Club, Shiraz Branch, Islamic Azad University, Shiraz, Iran, E-mail: farhadiasghar@shirazu.ac.ir

DOI: 10.17352/alo.000002

Received: 01 December, 2016 | Accepted: 13 December, 2016 | Published: 15 December, 2016

Keywords: ISPH; Kernel function; Outskirts decay coefficient; Solitary wave; Wave celerity

Cite this as

Farhadi A (2016) Simulating Solitary Wave Generation Using Incompressible SPH. Ann Limnol Oceanogr 1(1): 013-021. DOI: 10.17352/alo.000002

Abstract

The current study examines the generation and propagation of a Third order solitary water wave along the channel. Surface displacement and wave profile prediction challenges are interesting subjects in the field of marine engineering and many researchers have tried to investigate these parameters. To study the wave propagation problem, here, firstly the meshless Incompressible Smoothed Particle Hydrodynamics (ISPH) numerical method is described. Secondly, the boundary condition handling method, discretization, timestep selection and geometry provitions are presented. The numerical model is then used to simulate solitary wave propagation along the fixed depth channel. Here two still water depths of h = 0.2 m and h = 0.3 m are assumed and the dimensionless height of desired wave ranging from ε = 0.1 to ε = 0.6 are simulated. The numerical results show that studied Grimshaw Third order method can track the wave profile and it has acceptable relative variation for 5 seconds after the wave propagation, about 10%. In general, the numerical model gives satisfactory results for the wave kinematics.

Main article text

Introduction

Over the past few decades, offshore structures, such as oil platforms, offshore wind-power plants, have been in rapid growth in coastal and deep ocean regions, and wave-structure interaction has long been a strong interest in coastal and offshore engineering. A thorough understanding of the interaction of waves with offshore structures is vital in the safe and design of such structures. In addition, the flow field near the structures is helpful to understand the scour, sediment transport process in the coastal regions. In designing these structures, it is critical to be able to calculate wave forces acting on each individual structure.

Information on wave forces can be obtained by means of laboratory experiments or numerical simulations. Since laboratory experiments are usually constrained by the physical dimensions of laboratory facilities, it is not very often feasible to perform extensive parameter studies (e.g., variation of water depth, wave parameters, breaker type, etc.) even if the costs are of no concern. The alternative is to use numerical simulations as supplements to laboratory experiments, where accurate numerical simulations will also provide much more detailed insights into the physical processes that could not be achieved by experimental approach. In other words, a limited numbers of experiments can be designed so that the laboratory data can be effectively used to validate numerical models. The validated numerical models are then used to simulate scenarios with much wider range of physical parameters of interest.

So far, wide range of the numerical simulation models developed for wave propagation as the key concept in the maritime engineering which have been built upon the Navier-Stokes equations. One of these methods is the SPH method that has gradually matured over time into a suitable tool for computational fluid dynamics because of its flexibility to simulate complex problems such as flow through porous medium [1], multi-phase flows [2-4], heat conduction [5], free surface problems [6-8], fluid structure interactions [9,10], fuel cell [11], etc. However, compared with the Finite Difference (FD) or Finite Volume (FV) methods, SPH is still a relatively novel method in computational fluid dynamics and its shortcomings are still being improved. Shao and Lo [12], introduced ISPH algorithm based on the projection scheme. Numerical results have shown that ISPH produces reasonable accurate predictions of velocity and forces on solids.

Solitary wave generation is a traditional benchmark for numerical wave model tests. It is a permanent progressing wave form consisting of a single elevation above the undisturbed surface that propagates without the change of form on a constant still water depth over a flat bottom. A wide variety of analytical theories have been introduced for the solitary wave generation. It was first reported by Russell [13], who made remarkable experiments and gave an empirical relationship for the wave speed, which was later established theoretically, to the lowest order by Boussinesq [14] and Rayleigh [15], as part of an overall approximate solution. Since then, there have been several attempts to improve upon this solution, e.g. see [16-18]. Of the approximate solutions, several methods for the solitary wave have obtained series expansions in terms of wave amplitude, these being taken as far as the third order Grimshaw. He considered the one-dimensional modulations formed on the Boussinesq solitary wave and obtained third order equations analogous to those used by Boussinesq for the case of constant undisturbed depth in the higher order form [19]. Furthermore, an asymptotic solution was presented which described a slowly varying solitary wave [20]. An early work on WCSPH for solitary wave was done by Monaghan [21], while Lo and Shao used ISPH to generate solitary waves [22]. In this study, an ISPH method will be used to simulate the solitary wave generation and propagation in constant water depth.

Governing equations

The motion of a continuum in the Lagrangian description subjected to the action of body force in the isothermal condition, is represented by the continuity equation:

$\frac{D ρ}{D t} = - ρ \frac{\partial v_{i}}{\partial x_{i}} (1)$

and the momentum equation:

$\frac{D v_{i}}{D t} = f_{i} + \frac{1}{ρ} \frac{\partial σ_{i j}}{\partial x_{j}} (2)$

Where p is density, t is time, vi is the velocity vector, xi is the position vector, fi is the body force vector, σij is the stress tensor and the notation implies summation over repeated indices. The stress tensor can be decomposed into deviatoric viscous stress tensor τij and isotropic pressure p, according to the following equation:

$σ_{i j} = - p δ_{i j} + τ_{i j} (3)$

Where δij is the Kronecker delta. Pressure can be formally defined by the equation of state in the compressible flows, while for incompressible flows, it is derived from the divergence free condition of the velocity field. Assuming an incompressible Newtonian fluid, the continuity Equation (1) reduces to:

$\nabla . \vec{u} = 0 (4)$

and the momentum Equation (2) will be:

$\frac{D \vec{u}}{Dt} = - \frac{1}{ρ} \nabla p + \nabla . (v \nabla \vec{u}) + \vec{F} (5)$

Where ν is the kinematic viscosity. The conventional incompressible approach deals with pressure and velocity as primitive variables. The classical projection method [23] is used to calculate the pressure field and enforce incompressibility, simultaneously. The discretized form of the momentum equation is split into two parts. The first being the prediction step and is based on viscous and body forces. In this step, the intermediate velocity field ${\vec{u}}^{*}$ is obtained from velocity at (n)th time step:

$\frac{{\vec{u}}^{*} - {\vec{u}}^{(n)}}{Δ t} = {[\nabla . (ν \nabla \vec{u}) + \vec{F}]}^{(n)} (6)$

In each time step, the intermediate velocity field is calculated for fluid and boundary particles. In the second step, correction step, pressure force is included:

$\frac{{\vec{u}}^{(n + 1)} - {\vec{u}}^{*}}{Δ t} = {[- \frac{1}{ρ} \nabla p]}^{(n + 1)} (7)$

The intermediate velocity field is usually not divergence free but this is imposed upon ${\vec{u}}^{(n + 1)}$ . Hence, the intermediate velocity is projected on the divergence free space by taking the divergence of Equation (7) as:

$\frac{ρ}{Δ t} \nabla . ({\vec{u}}^{*}) = \nabla^{2} p^{(n + 1)} (8)$

where the ∇2 is the Laplacian operator. Once the pressure is obtained from pressure Poisson Equation (8), the velocity vector is updated by using the computed new pressure gradient:

${\vec{u}}^{(n + 1)} = {\vec{u}}^{*} - (\frac{1}{ρ} \nabla p^{(n + 1)}) Δ t (9)$

Finally, particles are moved according to this corrected velocity as:

$r^{(n + 1)} = r^{(n)} + {\vec{u}}^{(n + 1)} Δ t (10)$

SPH interpolation

The foundation of mesh free SPH method is based on integral interpolants which represents that any field variable X can be calculated over a set of SPH particles on domain of interest in terms of its values by taking a good interpolation kernel function. The exact integral representation of field variable X is:

$X (r) = \int_{Ω}^{} X (r^{'}) δ (r - r^{'}) d r^{'} (11)$

Where δ(r – r′) is Dirac delta function and Ω represents the computational domain. Equation (11) can be represented by defining a proper kernel function, W, with effective smoothing length h as:

$X (r) = \int_{Ω}^{} X (r^{'}) W (r - r^{'}) d r^{'} (12)$

In discrete notation, this approximation leads to the following approximation of the function at a interpolation particle a:

$X (\vec{r}) = \sum_{b} m_{b} \frac{X_{b}}{ρ_{b}} W_{a b} (13)$

where b is all the particles within the kernel function’s support domain. mb and pb are the mass and density of particle b, respectively, and weight function or kernel is denoted by. $W_{a b} = W ({\vec{r}}_{a} - {\vec{r}}_{b}, h)$

The parameter h is influence domain or smoothing domain, and controls the size of the area around particle a where contribution from the rest of the particles cannot be neglected. Considering the computational accuracy and efficiency [24], the following kernel function based on the cubic spline function and normalized in two dimensional is adopted:

$W (q, h) = α_{D} {\begin{matrix} 1 - \frac{3}{2} q^{2} + \frac{3}{4} q^{3} 0 \leq q \leq 1 \\ \frac{1}{4} {(2 - q)}^{3} 1 \leq q \leq 2 \\ 0 2 \leq q \end{matrix} (14)$

where for two dimensional cases $α_{D} = 10 / 7 π h^{2}$ . The gradient, divergence and Laplacian operators need to be formulated in ISPH algorithm. In the current work, the following commonly used forms are employed for gradient of a scalar A [25]:

$\nabla_{a} A \approx ρ_{a} \sum_{b} m_{b} (\frac{A_{a}}{ρ_{a}^{2}} + \frac{A_{b}}{ρ_{b}^{2}}) {\vec{\nabla}}_{a} W_{a b} (15)$

and divergence of a vector u* [26]:

$\nabla . {\vec{u^{*}}}_{a} \approx - \frac{1}{ρ_{a}} \sum_{b} m_{b} {\vec{u^{*}}}_{a b} . {\vec{\nabla}}_{a} W_{a b} (16)$

where ${\vec{\nabla}}_{a} W_{a b}$ is the gradient of the kernel function with respect to particle a and calculated as:

${\vec{\nabla}}_{a} W_{h} (r_{a b}) = \frac{{\vec{x}}_{a} - {\vec{x}}_{b}}{r_{a b}} \frac{\partial W_{h} (r_{a b})}{\partial r_{a b}} (17)$

Viscous term is discretized according to the relation given in [5]:

${(υ_{0} \nabla^{2} \vec{v})}_{a} \approx \sum_{b} m_{b} (\frac{16 υ_{0}}{(ρ_{a} + ρ_{b})} \frac{{\vec{v}}_{a b} . {\vec{r}}_{a b}}{{| {\vec{r}}_{a b} |}^{2} + η^{2}}) {\vec{\nabla}}_{a} W_{a b} (18)$

where ${\vec{A}}_{a b} = {\vec{A}}_{a} - {\vec{A}}_{b}$ and $η^{2} = 0.01 h^{2}$ is a parameter to avoid a zero denominator. Also Laplacian equation is discretized according to the relation given in [27]:

${(\nabla^{2} p)}_{a} = \sum_{b} m_{b} (\frac{ρ_{a} + ρ_{b}}{ρ_{a} ρ_{b}} \frac{{\vec{r}}_{a b} . {\vec{\nabla}}_{a} W_{a b}}{{| {\vec{r}}_{a b} |}^{2} + η^{2}}) p_{a b} (19)$

Resolution of the linear systems are widely studied by mathematicians as the demand for an efficient and smoothly-converging solver increases from numerical simulations. There are numerous iterative methods that are widely used in academic and commercial codes to solve the Pressure Poisson Equation (PPE). Here two solvers, Conjugate Gradient (CG) and Bi-Conjugate Gradient (Bi-CG) [28,29], can be applied to solve PPE.

The computational domain is divided into square cells of side 2h. Thus, for a particle located inside a cell, only the interactions with the particles of the same cell and its neighbors need to be considered (only 9 cells in 2-D). The searching algorithm is applied at the beginning of each time step updating the particle’s neighbors and the corresponding kernel derivatives. The time step limit for this method is the minimum of three conditions, the CFL, the mass and the viscous force conditions such that [26]:

$Δ t = \underset{a}{m i n} (\frac{h}{| u_{r e f} |}, \frac{h^{2}}{ν_{a}}, \sqrt{\frac{h}{f_{a}}}) (20)$

where fa is the force per unit mass, equivalent to the magnitude of particle acceleration and u_ref is the maximum fluid velocity in the domain [26].

Boundary condition handling

The Lagrangian nature of SPH method will cause the implementation of the boundary conditions less straightforward than in common mesh based methods. Different boundary conditions are used in the SPH method. In the present study, moving and stationary solid wall boundary conditions are used. There are different boundary types in SPH to simulate solid walls, namely the repulsive force [30], ghost or mirror particles [31] and dummy particles [12,32]. The repulsive force boundary condition, first proposed by Monaghan [30], uses forces similar to inter-molecular interactions. A force is exerted on a fluid particle having a distance r from a boundary particle, which has the form of Lennard-Jones potential. This force is increased as the distance r between a boundary and a fluid particle is decreased, preventing the fluid particles from penetrating the wall, Figure 1.

The mirror particle method is used to enforce the no-slip as well as the Neumann boundary conditions. In this method the particles whose support domain is truncated by a solid boundary are reflected on the other side of the wall. The mirror or ghost particles have the same pressure as their corresponding fluid particles but have velocities extrapolated from the fluid and wall velocities.

One of the sources of inaccuracy in the SPH method is the truncation of the boundary particles. This means that, not enough particles might be present in the support domain of a fluid particle. The other method to model solid walls is the use of dummy particles. In this method several layers of dummy particles are placed parallel to the boundary particles. So, the support domain of the particles located close to the solid wall will not be truncated any longer. These layers of dummy particles are linked to their corresponding boundary particles and have the same pressure and velocity as their linked particles. In the present work, dummy particles are used to model solid walls.

The number of dummy particle layers are decided from the radius of the compact support. In the following simulations, three layers of dummy particles are used. Governing equations are solved only for fluid and boundary particles and the pressure and velocity (also intermediate velocity) of the dummy particles are updated to their corresponding boundary wall particles.

The velocity and the intermediate velocities are held constant on the boundary particles, but their pressure is calculated from the PPE equation. Afterwards, the pressure of the dummy particles are updated to their corresponding boundary particles. In this way the Neumann boundary condition on the walls is approximated.

Wave characteristics

Three key parameters to identify waves are their lengths and heights, and the water depth over which they are propagating. All other parameters can be calculated from these quantities, e.g. wave induced water accelerations and velocities. A two dimensional schematic of a wave propagating in the x direction is shown in the Figure 2. The wave length, λ, is the horizontal distance between two successive wave troughs. This length is related to the water depth,h, and wave period, T, which is the time required for two successive troughs to pass a particular point. As the wave moves a distance λ, in time T, its speed called celerity is defined as C =λ/ T.

Solitary waves in constant water depth

Two dimensional dam break flow is chosen as the first suitable validation test case. Dam break flows over dry and wet beds have attracted wide research areas due to their theoretical, engineering and scientific considerations. If a dam break occurred over a dry bed, the generated wave is described by a tongue of water extended rapidly along the dry bed and if dam flows toward downstream over the wet bed, the attributed fluid flow features become remarkably different and some vorticity is developed at the front of the dam break. As a result, characteristics of the fluid flow are represented by the wave generation, wave crest development, wave breaking, and its impact with the downstream calm water and strew of water that generates some splash-up flow. Due to these features, dam break over the wet bed is an interesting benchmark to validate the numerical methods [33,34].

Solitary wave propagation in constant water depth is a classical benchmark problem for numerical wave model test. It propagates without the change of form in constant water depth over a flat bottom. Different analytical theories are already developed for the solitary wave. Therefore, it is well suited to evaluate the accuracy of the numerical model. For example, by checking the free surface location, the quality of volume tracking algorithm can be evaluated. In the numerical simulation, the free surface elevation and velocity distribution prescribed on the incident wave boundary are calculated by Third-order Grimshaw solitary wave solution [20,35]:

$η (x, t) = h [ε s^{2} - \frac{3}{4} ε^{2} s^{2} q^{2} + ε^{3} (\frac{5}{8} s^{2} q^{2} - \frac{101}{80} s^{4} q^{2})] (21)$

where ε = H/h; H is the wave height; h is the still water depth; $s = s e c h (α X / h)$ ; $q = t a n h (α X / h)$ ; $X = x - C t$ in which C is the wave speed; the coefficient α:

$α = \sqrt{\frac{3}{4} ε} (1 - \frac{5}{8} ε + \frac{71}{128} ε^{2}) (22)$

and the wave speed C is:

$C = \sqrt{g h} {(1 + ε - \frac{1}{20} ε^{2} - \frac{3}{70} ε^{3})}^{0.5} (23)$

and the velocity distribution is:

$\begin{array}{l} \frac{u}{\sqrt{g h}} = ε s^{2} + ε^{2} \\ {- \frac{3}{4} s^{2} + s^{2} q^{2} + {(\frac{z}{h})}^{2} (\frac{3}{4} s^{2} - \frac{9}{4} s^{2} q^{2})} + \\ ε^{3} {\begin{cases} \frac{21}{40} s^{2} - s^{2} q^{2} - \frac{6}{5} s^{4} q^{2} + {(\frac{z}{h})}^{2} (- \frac{9}{4} s^{2} + \frac{15}{4} s^{2} q^{2} + \frac{15}{2} s^{4} q^{2}) + \\ {(\frac{z}{h})}^{4} (\frac{3}{8} s^{2} - \frac{45}{16} s^{4} q^{2}) \end{cases}} \end{array}$

$\begin{array}{l} \frac{w}{\sqrt{g h}} = \sqrt{3 ε} (\frac{z}{h}) s^{2} q \\ [\begin{array}{l} ε + ε^{2} {- \frac{3}{8} - 2 s^{2} + {(\frac{z}{h})}^{2} (- \frac{1}{2} + \frac{3}{2} s^{2})} + \\ ε^{3} {\begin{cases} - \frac{49}{640} - \frac{17}{20} s^{2} - \frac{18}{5} s^{4} + {(\frac{z}{h})}^{2} (- \frac{13}{16} - \frac{25}{16} s^{2} + \frac{15}{2} s^{4}) + \\ {(\frac{z}{h})}^{4} (- \frac{3}{40} + \frac{9}{8} s^{2} - \frac{27}{16} s^{4}) \end{cases}} \end{array}] (24) \end{array}$

The coordinate system is defined in Figure 3. The Grimshaw solution is well suited for solitary wave of $ε < 0.5$ .

The Third order Grimshaw solitary waves is simulated to be compared with the analytical theory. Here two still water depths of h=0.2 m and h=0.3 m are assumed (Figure 4) and the dimensionless height of desired wave ranging from $ε = 0.1$ to $ε = 0.6$ h = 2.0 m are simulated, where particle resolution of 37000 and 77000 are used for the water depths of h=0.2 m and h=0.3 m, respectively. The computational domain is 50h in the stream wise direction x and 2h in the vertical direction y. In the stream wise direction, 1260 and 1860 particles with uniform grid size of $δ r = 0.008 m$ are used, while we use 50 and 75 particles in the vertical direction in two cases, respectively. The fluid has a density of $1000 k g / m^{3}$ and kinematic viscosity of $ν = 1 \times 10^{- 6} m^{2} / s$ .Three layers of dummy particles are used to handle the wall boundary condition.

Figure 5 shows the free surface elevation profiles at three different times (1, 3 and 5 seconds from beginning of the simulation) corresponding to solitary waves generated using the Third order Grimshaw method for a water column with $h = 0.3 m$ and $ϵ = 0.4$ . It is noticeable that compared with the analytical solution, the wave profile heights are decreased as it follows throughout the channel. The discrepancy between the predicted and analytical heights are 10% when the wave propagates along the channel for first 5 seconds.

On Figure 6(a), the paddle laws of motion corresponding to the Third order wave generation is plotted for $ϵ = 0.4$ for a water column with $h = 0.3 m$ . Greatest accelerations occur somewhere between the beginning of motion (zero velocity) and mid-stroke (maximum velocity). Qualitatively, the larger the maximum velocity and the shorter the duration of motion, the greater the acceleration. For the same wave we also plot on Figure 6 (b) the dimensionless paddle velocity corresponding to this law of motion. The maximum nondimensionalized velocity of 0.27 is observed.

Measurements of the free surface elevation at different distances from the paddle for a solitary wave of desired dimensionless amplitude $ε = 0.2$ in a water column of $h = 0.2 m$ where was generated using Third order numerical integration are plotted on Figure 7. The solitary wave amplitude is decreases slowly as it travels along the channel.

Figure 8 shows the nondimensionalized velocity component variations versus nondimensionalized vertical coordinate of the wave at $t = 1.5 s$ from the beginning of the simulation for two locations ( $x = 2 m$ and $x = 3 m$ regarding the initial paddle location). The velocity profiles are close to the analytic ones. Near the water surface, under both wave surfaces, the measured particle horizontal velocity is somewhat greater than the theoretical values, while the vertical velocity components is underestimated compared with the analytic ones, the reason of which may be due to the particle resolution.

The flow field contours at three different times are illustrated in Figures (9-11) for the $ϵ = 0.3$ case for a water column of $h = 0.2 m$ . h = 2.0m As seen, the pressures are captured precisely and the pressure distributions under the wave crests are accurately predicted, Figure 9. The wave speed at this case is $C = 1.59 m / s$ , $C = 1.59 m / s$ that proves the wave crest propagation at these selected times, where the wave crest transferred to x = 2.61 m, x = 5.83 m and x = 9.05 m, respectively, Figure 10. Furthermore, the symmetric behavior of the vertical velocity at the wave is noticeable at these selected times, Figure 11.

In order to investigate the solitary wave generated behavior on the pureness of the generated wave on other wave heights, this problem is solved for a wide range of the desired wave heights, $ϵ = 0.1 \to 0.6$ . As expected, greater the desired wave height, the greater amount of the water volume must be pushed. Figure 12(a) illustrates the depth averaged net mass displacement L of a solitary wave. This net mass displacement is the total stroke of the paddle prescribed in each procedure. The rate of the displacement increment is decreased as the desired wave height is increased. Furthermore, the maximum paddle velocity (which occurs at mid-stroke) is obtained for this range of wave heights. Again, the slope of the maximum paddle velocity is decreases as dimensionless amplitude increases.

The dimensionless phase speed (or Froude number, $F = C / \sqrt{g h}$ ) is plotted on Figure 13(a) for the selected wave height range. It shows that the Third order method phase speed is matches to the Byatt-Smith numerical estimation in the studied range of ε. The dimensionless outskirts decay coefficient $(2 α = β . h)$ versus the dimensionless amplitude is plotted on Figure 13(b). This outskirts decay coefficient describes the way free surface elevation tends towards the mean level at infinity. Stokes showed that $β$ is a solution of the following equation, also used by Byatt-Smith [36]:

$F^{2} = \frac{\tan β}{β} (25)$

As a matter of fact the Third order method for $ϵ < 0.4$ matches the Byatt-Smith reference but for greater $ϵ$ some discrepancy is observed.

Conclusion and Future Works

The ISPH numerical method is used to simulate the solitary wave generation and propagation along the channel at different wave amplitudes and water depths in an ISPH-based numerical wave flume. The numerical results are compared with analytical data in terms of free surface displacements, fluid particle velocity, phase speed, flow field counters and some other wave parameters.

In the first section, the free surface profile variations over time, position and through solitary wave amplitude ranges are assumed. The numerical free surface profiles are compared with analytical results at various times and it is proved that studied Grimshaw Third order method has acceptable relative variation for 5 seconds after the wave propagation, about 10%.

In the second section, solitary wave paddle motions, paddle velocities and accelerations, their displacements, phase speeds and outskirts decay coefficients are assumed. Results were in a good agreement with the analytical data. Then, maximum paddle velocity and displacement of the wave generation procedures are derived and as seen they coincide with the analytical data. Based on the obtained results, I am continuing this methodology in the two solitary waves runup, the C-wave in sloped beachesm ship motions in the regular waves and breaking criteria for two solitary waves investigations in the in-hand studoes.

References

Copyright

© 2016 Farhadi A, et al. This is an open-aaacess article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.