Home | About | Journals | Submit | Contact Us | Français |

**|**Bioengineering (Basel)**|**v.4(1); 2017 March**|**PMC5590438

Formats

Article sections

- Abstract
- 1. Introduction
- 2. Results from Finite Element Model
- 3. Results from Experiments
- 4. Study of Drug Particle Velocity and Angular Deviation
- 5. Conclusions
- References

Authors

Related links

Bioengineering (Basel). 2017 March; 4(1): 24.

Published online 2017 March 14. doi: 10.3390/bioengineering4010024

PMCID: PMC5590438

Ramana Pidaparti, Academic Editor and Hu Yang, Academic Editor

Department of Aerospace Engineering, Indian Institute of Science, Bangalore 560012, India; Email: moc.liamg@kevivunat

Received 2017 February 18; Accepted 2017 March 13.

Copyright © 2017 by the authors.

Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

This paper presents the design optimization of diaphragms for a micro-shock tube-based drug delivery device. The function of the diaphragm is to impart the required velocity and direction to the loosely held drug particles on the diaphragm through van der Waals interaction. The finite element model-based studies involved diaphragms made up of copper, brass and aluminium. The study of the influence of material and geometric parameters serves as a vital tool in optimizing the magnitude and direction of velocity distribution on the diaphragm surface. Experiments carried out using a micro-shock tube validate the final deformed shape of the diaphragms determined from the finite element simulation. The diaphragm yields a maximum velocity of 335 m/s for which the maximum deviation of the velocity vector is 0.62°. Drug particles that travel to the destination target tissue are simulated using the estimated velocity distribution and angular deviation. Further, a theoretical model of penetration helps in the prediction of the drug particle penetration in the skin tissue like a target, which is found to be 0.126 mm. The design and calibration procedure of a micro-shock tube device to alter drug particle penetration considering the skin thickness and property are presented.

Micro-scale drug delivery devices have recently drawn a great attention due to their advantages over hypodermic needles. Several interesting developments such as liquid jet injectors, powder injectors, micro-needles and thermal micro-ablation help in tackling pain and needle phobia [1,2,3,4]. Apart from being expensive, micro-needles can cause local inflammation and skin irritation [1]. Thermal micro-ablation procedure is slow and the liquid jet injector causes bleeding [1]. Powder injectors are non-invasive drug delivery systems that have received much attention because the drug delivery is painless and the chances of spreading of autoimmune and communicable diseases are none. Needleless drug delivery devices [5,6,7,8,9,10,11,12,13,14] use a shock tube to deliver painkillers, contraceptives, genetic material and insulin, which are frequently used treatments. Often, these devices involve a contoured shock tube with pressurized gas filled in first the compartment sealed by a diaphragm [5,6,7,8,9,10,11]. In such devices, the drug particles are placed on the surface of the diaphragm facing the pressurized gas. The drug delivery procedure ruptures the diaphragm, causing the pressurized gas to pass through the nozzle section, simultaneously accelerating the drug particle along the flow. The drug particle and the gas escape the nozzle (second compartment) and hit the skin target. Another category of shock tubes uses a ballistic gun that involves a moving piston, stopping suddenly to allow the particles catapulted towards the target [12,13,14]. However, this method involves bulky devices, which have parts subjected to wear and tear. Recently, researchers have focused on optimizing the needleless drug delivery systems for their size and cost. In this direction, Jagadeesh and Takayama [15] have reviewed applications of micro shock waves generated by piezoceramics, electro-hydraulics, micro-explosives and the pulsed laser. In addition, Jagadeesh et al. [16] have proposed a typical needleless drug delivery device consisting of a diaphragm and explosive driver to propel the liquid drug to the target. However, the study did not report the device calibration considering drug particle and skin properties.

Many researchers have reported the analysis of flow and drug particle penetration in the contoured shock tube-based particulate drug delivery device [15,16,17,18,19,20,21,22,23]. The velocity, density and size of the drug particle were the main parameters governing drug particle penetration. Kendal [10] studied particle flow in a contoured shock tube using particle image velocimetry, which captured 500 m/s velocity of gold particles with a diameter of 1 μm. A unified penetration model predicted the particle penetration of 20–30 μm, which matched with experiments. Mitchell et al. [11] used a semi-empirical relation to determine particle penetration depth in mucous tissue. They studied the penetration of particle with 0.8–3.7 μm size and with a density of 16,800 kg/m^{3}. A velocity of particles 550 m/s gave particle penetration of 20–16 μm. Numerical simulation of gas–particle interaction using a CFD approach predicted the possibility of achieving 570 m/s velocity with polystyrene particles of 39 μm [19]. Optimization of exit-to-throat area ratio using fluid mechanics yielded a velocity of 1050 m/s and 400 m/s in a contoured and conical nozzle respectively. The magnitude of the velocity and the direction of drug particle govern the penetration in the target. Tissue damage is greater if the direction of drug particle deviates significantly from the desired direction. Therefore, understanding the velocity distribution is also important for the efficient design of the diaphragm and the delivery mechanism. In addition to the velocity of drug particle, the density and yield stress of the skin influences the particle depth penetration. These properties of skin vary depending upon the body site and age of a person [24,25,26]. Thus, the drug delivery device has to be calibrated based on the condition of human skin before drug delivery. The calibration requirement needs an accurate estimate of the velocity profile and the drug particle penetration.

In the present work, we report the capability of explosive driven shock tube-based drug delivery device developed by Jagadeesh et al. [16] to impart velocity to the drug particles. The objective of this work is to optimize the diaphragm for efficient drug delivery via simulations since experiments are expensive. In this device, as shown in Figure 1a, the expendable entities are the diaphragm with a coat of drug particles and a polymer tube with an explosive coating. Drug particles gain velocity entirely by the motion of the diaphragm due to the incident shock wave produced by igniting explosive coated in a polymer tube [16]. The exploded particles and gases remain trapped within the shock tube after its operation, making it the cleanest and safest drug delivery device among the devices driven by explosives [5,6,7,8,9,10,11]. The entire duration of shock tube operation is a few milliseconds. Thus, the heat generation is negligible, by which the device does not heat up. For large-scale applications, it is necessary to design the diaphragm for optimal material usage, easy manufacturing, safety, maximizing drug particle velocity and optimal deviation of the velocity trajectories for specific target delivery. In this study, we consider the design of copper, aluminium and brass diaphragms for a single use in an explosive driven shock tube. The paper also presents a method to estimate drug particle penetration in human skin with various different biomechanical properties. We also propose a calibration procedure of the device to achieve the desired drug particle penetration. The next section presents the shock pressure history and diaphragm material models used in the finite element simulations.

(**a**) Schematic of a hand-held micro-shock tube-based drug delivery device (**b**) Deformed shape of the diaphragms (diameter = 8 mm and thickness indicated) subjected to explosive driven shock.

A uniform pressure distribution serves as a close approximation of pressure on the diaphragm, which results due to the release of compressed gas in a shock tube. In such calculations, friction or wall shear effect due to the shock tube is neglected [27]. Spherical pressure distribution closely approximates the pressure on diaphragm generated by explosives [16]. Such a pressure distribution is expressed as:

$$\begin{array}{c}\mathrm{P}\left(\mathrm{r},\mathrm{t}\right)={\mathrm{P}}_{0}\frac{{\left({\mathrm{r}}_{0}^{2}-{\mathrm{r}}^{2}\right)}^{1/2}}{{\mathrm{r}}_{0}}\left(\frac{\mathrm{t}}{{\mathrm{t}}_{\mathrm{i}}}\right);\text{\hspace{1em}}\mathrm{t}<{\mathrm{t}}_{\mathrm{i}}\\ \mathrm{P}\left(\mathrm{r},\mathrm{t}\right)={\mathrm{P}}_{0}\frac{{\left({\mathrm{r}}_{0}^{2}-{\mathrm{r}}^{2}\right)}^{1/2}}{{\mathrm{r}}_{0}};\text{\hspace{1em}}{\mathrm{t}}_{\mathrm{i}}<\mathrm{t}<{\mathrm{t}}_{0}\end{array}\}$$

(1)

where ${\mathrm{P}}_{0}$ is the shock pressure, ${\mathrm{r}}_{0}$ is the radius of the diaphragm, $\mathrm{t}$ is the time, ${\mathrm{t}}_{\mathrm{i}}$ is the initial rise time of shock and ${\mathrm{t}}_{0}$ is the time duration of shock. The tube axis passes through the centre of the diaphragm. This type of ramp is typical to shock tube problem. However, the micro-shock tube considered in the present study produces an explosive-driven shock with chemical kinetics at the inner surface of the tube where the explosive material layer is present and combustion process initiates. This leads to several complications. In order to simplify the problem, we assume that at any cross-section of the tube, micro-shocks emanate from the surface of the tube and all of these micro-shocks from the periphery of the cross-section finally give rise to a spherical shock. This spherical shock front remains unaltered due to chemical kinetic fronts propagating on the inner surface along the tube axis with a constant speed and peak pressure. We further assume that such a spherical shock first impinges on the diaphragm covering the tube cross-section and then the pressure builds up in the entire diaphragm, which is because of a larger radius of the diaphragm compared to the micro-shock tube diameter. A more realistic simulation would involve modelling this pressure build up in the outer periphery of the diaphragm. It would require a coupled model involving combustion-driven shock and elasto-visco-plastic deformation of the diaphragm. This will be considered in a separate article. We consider elasto-plastic analysis since plastic deformation is not important in the present study. This is because the maximum velocity imparted to the particle is mainly when the diaphragm is in an elastic state. In the plastic state, strain hardening causes the diaphragm to retard, reducing the velocity. The main hypothesis behind the present modelling approach is drawn on the preliminary calculations that the pressure build-up mechanism in the central region as well as in the outer periphery of the diaphragm can be approximated as a spherical shock without any discontinuity at r = r_{0}. The objectives of finite element simulation here are twofold. (1) To identify the shock pressure for the assumed distribution, which produces the experimentally-obtained deformed shape of the diaphragm (See Figure 1b) and (2) to understand the elasto-plastic response from diaphragms of various materials and the thickness and the velocity imparted to the drug particles.

Several analytical models analyse the diaphragms or plates undergoing plastic deformation upon loading due to the explosion, blast or shock wave [28,29,30,31,32,33]. Using these models, measurement of the resultant velocity distribution on the diaphragm surface due to impinging explosive pressure is difficult. Simulations based on finite element method provide an easy, cost-effective and reliable way to carry out a parametric sensitivity study and design of the diaphragm. In the present work, we have carried out finite element simulations to predict the maximum velocity of the diaphragm. To capture the plastic deformation of the diaphragm, an elastoplastic material model was considered for detailed simulations. Bilinear, multi-linear and exponential kinematic hardening functions available in the literature for different alloys of basic materials approximate the elasto-plastic deformation. A schematic representation of the stress-strain behaviour of the linear hardening model (shown in Figure 2) described by a linear hardening function [34] in terms of the tangent modulus (E_{T}) and Young’s modulus (E) is written as:

$$\mathrm{H}=\frac{{\mathrm{E}}_{\mathrm{T}}}{1-\left({\mathrm{E}}_{\mathrm{T}}/\mathrm{E}\right)}.$$

(2)

The present work uses a finite element model to analyse the deformation of the diaphragm using the elasto-plastic material property. Although the diaphragm is conventionally modelled using the shell element, we consider the solid element to capture the effect of thinning of the diaphragm plastically, as is the case observed experimentally. Considering a three-dimensional problem, the strain-displacement relations in matrix form are given by [35],

$$\left\{\begin{array}{c}{\mathsf{\epsilon}}_{\mathrm{x}\mathrm{x}}\\ {\mathsf{\epsilon}}_{\mathrm{y}\mathrm{y}}\\ {\mathsf{\epsilon}}_{\mathrm{z}\mathrm{z}}\\ {\mathsf{\epsilon}}_{\mathrm{x}\mathrm{z}}\\ {\mathsf{\epsilon}}_{\mathrm{y}\mathrm{z}}\\ {\mathsf{\epsilon}}_{\mathrm{x}\mathrm{y}}\end{array}\right\}=\left[\begin{array}{cccccc}\frac{\partial}{\partial \mathrm{x}}& 0& 0& \frac{\partial}{\partial \mathrm{z}}& 0& \frac{\partial}{\partial \mathrm{y}}\\ 0& \frac{\partial}{\partial \mathrm{y}}& 0& 0& \frac{\partial}{\partial \mathrm{z}}& \frac{\partial}{\partial \mathrm{x}}\\ 0& 0& \frac{\partial}{\partial \mathrm{z}}& \frac{\partial}{\partial \mathrm{x}}& \frac{\partial}{\partial \mathrm{y}}& 0\end{array}\right]\left\{\begin{array}{c}{\mathrm{u}}_{\mathrm{x}}\\ {\mathrm{u}}_{\mathrm{y}}\\ {\mathrm{u}}_{\mathrm{z}}\end{array}\right\},$$

(3)

The constitutive model relates stress to strain by:

$$\left\{\begin{array}{c}{\mathsf{\sigma}}_{\mathrm{x}\mathrm{x}}\\ {\mathsf{\sigma}}_{\mathrm{y}\mathrm{y}}\\ {\mathsf{\sigma}}_{\mathrm{z}\mathrm{z}}\\ {\mathsf{\sigma}}_{\mathrm{x}\mathrm{z}}\\ {\mathsf{\sigma}}_{\mathrm{y}\mathrm{z}}\\ {\mathsf{\sigma}}_{\mathrm{x}\mathrm{y}}\end{array}\right\}=\left[\begin{array}{cccccc}{\mathrm{C}}_{11}& {\mathrm{C}}_{12}& {\mathrm{C}}_{12}& 0& 0& 0\\ {\mathrm{C}}_{12}& {\mathrm{C}}_{11}& {\mathrm{C}}_{12}& 0& 0& 0\\ {\mathrm{C}}_{12}& {\mathrm{C}}_{12}& {\mathrm{C}}_{11}& 0& 0& 0\\ 0& 0& 0& \mathrm{G}& 0& 0\\ 0& 0& 0& 0& \mathrm{G}& 0\\ 0& 0& 0& 0& 0& \mathrm{G}\end{array}\right]\left\{\begin{array}{c}{\mathsf{\epsilon}}_{\mathrm{x}\mathrm{x}}\\ {\mathsf{\epsilon}}_{\mathrm{y}\mathrm{y}}\\ {\mathsf{\epsilon}}_{\mathrm{z}\mathrm{z}}\\ {\mathsf{\epsilon}}_{\mathrm{x}\mathrm{z}}\\ {\mathsf{\epsilon}}_{\mathrm{y}\mathrm{z}}\\ {\mathsf{\epsilon}}_{\mathrm{x}\mathrm{y}}\end{array}\right\},$$

(4)

where ${\mathrm{C}}_{11}=\mathrm{E}\left(1-\mathsf{\upsilon}\right)/\left(1-2\mathsf{\upsilon}\right)\left(1+\mathsf{\upsilon}\right);$ ${\mathrm{C}}_{12}=\mathrm{E}/\mathsf{\upsilon}\left(1-2\mathsf{\upsilon}\right)\left(1+\mathsf{\upsilon}\right);$ $\mathrm{G}=\mathrm{E}/2\left(1+\mathsf{\upsilon}\right)$. The equation for momentum balance of a three-dimensional elastic solid is expressed as:

$$\nabla \u2022\mathsf{\sigma}+\mathrm{f}=\mathsf{\rho}\ddot{\mathrm{u}},$$

(5)

By substituting Equation (4) in Equation (5), we get the standard Navier’s equation, given by:

$$\mathsf{\rho}\frac{{\partial}^{2}\mathrm{u}}{\partial {\mathrm{t}}^{2}}-\nabla .\left(\overline{\mathrm{C}}\nabla \mathrm{u}\right)=\mathrm{f},$$

(6)

where ${\overline{\mathrm{C}}}_{\mathrm{i}\mathrm{j}\mathrm{k}\mathrm{l}}=\left({\mathrm{C}}_{\mathrm{i}\mathrm{j}\mathrm{k}\mathrm{l}}+{\mathrm{C}}_{\mathrm{k}\mathrm{l}\mathrm{i}\mathrm{j}}\right)/2$ and ${\mathrm{C}}_{\mathrm{i}\mathrm{j}\mathrm{k}\mathrm{l}}$ is the elasticity tensor whose matrix form is given in Equation (4). Next, Equation (6) is transformed to a first-order state variable form [36]. Defining the velocity field as $\mathrm{v}={\left[{\mathrm{u}}_{\mathrm{x},\mathrm{t}}{\mathrm{u}}_{\mathrm{y},\mathrm{t}}{\mathrm{u}}_{\mathrm{z},\mathrm{t}}\right]}^{\mathrm{T}}$ and considering Rayleigh-type damping [37] of the form $\mathrm{d}=\mathsf{\alpha}\mathsf{\rho}+\mathsf{\beta}\overline{\mathrm{C}}$, Equation (6) is rewritten as

$$\left[\begin{array}{cc}\mathrm{I}& 0\\ 0& \mathsf{\rho}\end{array}\right]\frac{\partial}{\partial \mathrm{t}}\left[\begin{array}{c}\mathrm{u}\\ \mathrm{v}\end{array}\right]-\nabla \u2022\left[\left[\begin{array}{cc}0& 0\\ \overline{\mathrm{C}}& \overline{\mathrm{C}}\mathsf{\beta}\end{array}\right]\nabla \left[\begin{array}{c}\mathrm{u}\\ \mathrm{v}\end{array}\right]\right]+\left[\begin{array}{cc}0& -\mathrm{I}\\ 0& \mathsf{\rho}\mathsf{\alpha}\end{array}\right]\left[\begin{array}{c}\mathrm{u}\\ \mathrm{v}\end{array}\right]=\left[\begin{array}{c}0\\ \mathrm{f}\end{array}\right].$$

(7)

Elastic-plastic analysis procedure involves solving the system of equilibrium equations given by:

$$\left\{\mathrm{F}\right\}=\left[\mathrm{K}\right]\left\{\mathrm{U}\right\},$$

(8)

where, $\left\{\mathrm{F}\right\}$ is the external load matrix, $\left[\mathrm{K}\right]$ is the system stiffness matrix and $\left\{\mathrm{U}\right\}$ is the system displacement matrix. By loading the structure in an incremental way, the finite element model extends the elastic problems in the plastic range. Von Mises yield criteria is used in the present simulation to identify yielding in an element which is given in terms of principal stresses ${\mathsf{\sigma}}_{1}$, ${\mathsf{\sigma}}_{2}$ and ${\mathsf{\sigma}}_{3}$ (calculated at nodal points) and yield stress ${\mathsf{\sigma}}_{0}$ given by:

$${\left({\mathsf{\sigma}}_{1}-{\mathsf{\sigma}}_{2}\right)}^{2}+{\left({\mathsf{\sigma}}_{2}-{\mathsf{\sigma}}_{3}\right)}^{2}+{\left({\mathsf{\sigma}}_{3}-{\mathsf{\sigma}}_{1}\right)}^{2}=2{\mathsf{\sigma}}_{0}^{2}.$$

(9)

After identifying the yielding, following equilibrium equation in terms of incremental displacement $\left\{\Delta \mathrm{U}\right\}$, elastic load $\left\{\Delta \mathrm{F}\right\}$ and plastic load $\left\{\Delta {\mathrm{F}}^{\mathrm{p}}\right\}$ is solved to obtain $\left\{\Delta \mathrm{U}\right\}$, that is,

$$\left[\mathrm{K}\right]\left\{\Delta \mathrm{U}\right\}=\left\{\Delta \mathrm{F}\right\}+\left\{\Delta {\mathrm{F}}^{\mathrm{p}}\right\}.$$

(10)

The plastic part of the strain is used for increments in the iterative loop while the elastic part of the strain is used to update nodal displacements. The subsequent steps use strain displacement relation to compute strain increment $\left\{\Delta \mathsf{\epsilon}\right\}$. Now, the finite element procedure estimates elastic strain $\left\{\Delta {\mathsf{\epsilon}}^{\mathrm{e}}\right\}$ by subtracting plastic strain $\left\{\Delta {\mathsf{\epsilon}}^{\mathrm{p}}\right\}$ from strain increment $\left\{\Delta \mathsf{\epsilon}\right\}$. Through the stress–strain relations involving the elasto-plastic constitutive model, new stress $\left\{\Delta \mathsf{\sigma}\right\}$ is computed that gives stress redistribution. This stress strain relation uses a hardening function H, elastic stiffness matrix [D] and elasto-plastic stiffness matrix ${\left[\mathrm{D}\right]}_{\mathrm{e}\mathrm{p}}$ as given below [38],

$$\left\{\Delta \mathsf{\epsilon}\right\}={\left[\mathrm{D}\right]}_{\mathrm{e}\mathrm{p}}^{-1}\left\{\Delta \mathsf{\sigma}\right\},$$

(11)

$${\left[\mathrm{D}\right]}_{\mathrm{e}\mathrm{p}}=\left[\mathrm{D}\right]-\left[\mathrm{D}\right]\left\{\frac{\partial \mathrm{Q}}{\partial \mathsf{\sigma}}\right\}{\left\{\frac{\partial \mathrm{F}}{\partial \mathsf{\sigma}}\right\}}^{\mathrm{T}}\left[\mathrm{D}\right]{\left[\mathrm{H}+{\left\{\frac{\partial \mathrm{F}}{\partial \mathsf{\sigma}}\right\}}^{\mathrm{T}}\left[\mathrm{D}\right]\left\{\frac{\partial \mathrm{Q}}{\partial \mathsf{\sigma}}\right\}\right]}^{-1},$$

(12)

where Q is elastic potential and F is von Mises flow rule, given by:

$$\mathrm{Q}={\left({\mathsf{\sigma}}_{1}-{\mathsf{\sigma}}_{2}\right)}^{2}+{\left({\mathsf{\sigma}}_{2}-{\mathsf{\sigma}}_{3}\right)}^{2}+{\left({\mathsf{\sigma}}_{3}-{\mathsf{\sigma}}_{1}\right)}^{2},$$

(13)

$$\mathrm{F}=\sqrt{{\left({\mathsf{\sigma}}_{\mathrm{x}\mathrm{x}}-{\mathsf{\sigma}}_{\mathrm{y}\mathrm{y}}\right)}^{2}+{\left({\mathsf{\sigma}}_{\mathrm{y}\mathrm{y}}-{\mathsf{\sigma}}_{\mathrm{z}\mathrm{z}}\right)}^{2}+{\left({\mathsf{\sigma}}_{\mathrm{z}\mathrm{z}}-{\mathsf{\sigma}}_{\mathrm{x}\mathrm{x}}\right)}^{2}+3{\mathsf{\sigma}}_{\mathrm{x}\mathrm{y}}^{2}+3{\mathsf{\sigma}}_{\mathrm{y}\mathrm{z}}^{2}+3{\mathsf{\sigma}}_{\mathrm{z}\mathrm{x}}^{2}}-{\mathsf{\sigma}}_{0}$$

(14)

The diaphragm is circular and symmetric about its axis. Thus, the analysis considers only a quarter model. Due to a fixed condition of the diaphragm along its circumference in its holder, the boundary along the circumference is fully restrained in the finite element model. The finite element model uses a moving mesh, which captures large deformation. Finite element mesh uses Lagrangian quadratic elements with mid nodes on the edge of the solid elements. The iterative GMRES solver is employed to solve the finite element system.

Experiments were carried out using aluminium, copper and brass diaphragms by a research team at the Shock Wave Laboratory, Indian Institute of Science. The data were obtained for the purpose of validation of the present study. Dynamic response measurements with the diaphragm were carried out using a non-contact laser interferometer by the present authors. Finite element simulations (using elasto-plastic material properties of aluminium) were carried out to assess the magnitude and duration of loading with reference to the deformed shape of the diaphragm. The following material properties were used for aluminium: Young’s modulus E = 70 GPa, Poisson’s ratio $\mathsf{\upsilon}$ = 0.33, yield strength = 20 MPa, and tangent modulus E_{T} = 7 GPa. The following material properties are used for copper: Young’s modulus E = 110 GPa, Poisson’s ratio $\mathsf{\upsilon}$ = 0.35, yield strength = 33 MPa, and tangent modulus E_{T} = 11 GPa. Mass proportional damping coefficient and stiffness proportional damping coefficient are taken as 1 and 0.001, respectively. A moving mesh was employed since deformation is very large. A ramp pressure distribution is assumed with a pressure of 600 MPa with an initial rise time of 1 μs. The magnitude of pressure was identified as the one that gives the deformed shape and peak deflection at the centre of the diaphragm identical to the ones obtained experimentally.

The centre of the diaphragm attains its maximum velocity within the initial rise-time period of 1 μs as shown in Figure 3. The simulation results for all types of materials and diaphragm thickness shows such velocity profile. Simulations using copper give a lower velocity as compared to aluminium. This is primarily due to a relatively higher modulus of elasticity of copper. Figure 4 shows the deformation history at the centre of the diaphragm in the form of out of plane displacement. The deformation stabilizes with an exponentially decaying velocity. This corresponds to the accumulation of residual deformation left out corresponding to plastic deformation. Copper has higher strength and settles with comparatively smaller plastic deformation.

Out-of-plane velocity history at the centre of diaphragm for various different materials and diaphragm thickness.

Figure 1b shows a set of diaphragms exposed to shock in a micro-shock tube. Jagadeesh et al. [16] have documented the expendable materials and the method followed to operate a micro-shock tube device. In the present study, the experiments were performed with the same setup with different diaphragm materials. Figure 5a shows the deformed shape obtained from finite element simulation for the given pressure loading and a moving mesh. Figure 5b shows the shape of the experimental specimen after the shock event. Table 1 compares the deformation at the centre of the diaphragm obtained by simulations to that of experiments. The deformed shape varies along radial coordinate (See Figure 6). The material properties used are the properties at annealed state. Pressure history was assumed based on the explosive driven shock mechanisms in the micro-shock tube as discussed in Section 1.1. Pressure history estimated from the analysis represents the incident pressure on the diaphragm by the explosive used with the shock tube. The obtained pressure history by matching the deformation of the diaphragm analytically and experimentally serves as an input to study the drug particle travel, velocity and deviation. The next section presents the simulation of drug particle travel, velocity and deviation followed by the estimation of drug particle penetration in the skin.

Final deformation of a 0.1-mm thick copper diaphragm, (**a**) deformed shape obtained from finite element simulation and (**b**) deformed shape of the experimental specimen.

Comparison between finite element and experimental deformations of the diaphragm with various thicknesses. Material properties of aluminium were used in the simulations.

The distribution of drug particle deviation and velocity needs to be studied for effective and accurate drug delivery. The analytical displacement and velocity history obtained from Section 2 was post processed to simulate the path of travel and velocity of the drug particle deposited on the diaphragm. The drug particles are usually deposited circularly at the centre, thus the particles deposited from the centre to the radial location of about r = R were considered. Neglecting diffusion, the travel path of drug particles from the diaphragm to the target located at a distance of 4 mm from it was studied (See Figure 7a). Due to the symmetry of diaphragm and drug deposition, the path traced by drug particles was symmetric about the centre and thus the deviation of the path was same radially. Thus, few deposited particles on diaphragm along the cross section taken through the centre of the diaphragm were considered for the travel path simulation as shown in Figure 7a.

Simulation of deformation of aluminium diaphragm and drug particle travel (**a**) at t = 0 s; (**b**) at t = 0.8 μs; (**c**) at t = 0.9 μs; (**d**) at t = 0.95 μs; (**e**) at t = 2 μs and (**f**) at t = 1 ms.

The initial location of drug particles is on the outer surface of the diaphragm. On the bottom surface, the pressure history P(r,t) was applied (see Figure 7b). The velocity increases up to 0.8 μs and the particles remain attached to the surface due to the gain in the momentum. At 0.8 μs, the diaphragm experiences strain hardening due to high strain rate and suffer a decrease in velocity. A critical cohesive force reaches that is required to break out the particles from the surface. Particles detach from the surface of the diaphragm and travel towards the target with a certain deviation (deviation is evaluated in next section). The loss in energy of the particle due to detachment was considered negligible and deviation due to diffusion was neglected. The time at which particle experience retardation occurred was assumed as the time when the particle detached from the diaphragm. The breaking of cohesive bonding was not simulated in the present study. The particles after the detachment were traced further with their velocity and deviation at this instance. With the velocity distribution and initial position of the particles, the final positions at various time intervals were evaluated. Figure 7c,d shows the deformation of the diaphragm and particle positions while reaching the target. At 0.95 μs, entire particles reach the target, but the diaphragm still deforms plastically due to pressure still acting on it. The diaphragm further deforms as shown in Figure 7e and reaches a final deformed position as shown in Figure 7f.

The angle of deviation of the diaphragm surface velocity vector with respect to the central axis was estimated from the simulations. The angle of deviation ($\mathsf{\delta}$) as shown in Figure 8 is computed by:

$$\mathsf{\delta}={\mathrm{tan}}^{-1}\left(\frac{{\mathrm{V}}_{\mathrm{x}}}{{\mathrm{V}}_{\mathrm{y}}}\right)\frac{180}{\mathsf{\pi}},$$

(15)

where ${\mathrm{V}}_{\mathrm{x}}$ is the surface velocity component in the x-direction and ${\mathrm{V}}_{\mathrm{y}}$ is the surface velocity component in the y-direction. The angle $\mathsf{\delta}$ gives an estimate of how much the drug particles would deviate from straight path over a distance of travel before they hit the target. Greater deviation causes the drug particle to travel comparatively more distance, resulting in a loss of momentum. An oblique penetration can become inefficient in terms of the penetration depth in the target and depend on the target (tissue/cell structure). This extends the velocity requirement at points having a higher deviation, which is difficult to achieve. Deviation also causes greater damage to the tissue. Thus, the objective of the analysis is primarily to minimize the angle of deviation $\mathsf{\delta}$ and at the same time maximize the velocity for a particular shock-loading.

Figure 7b with a 0.8-us snapshot shows the particles remaining attached to the diaphragm surface. However, the velocity vector orientations (angular deviation) at that time are not perpendicular for all particles. The angular deviations are shown in Figure 9. The angular values are due to the deformed curvature of the diaphragm. The particles located at the centre have a least deviation and highest velocity. Figure 9 shows the variation of deviation of particles radially at the time of detachment. Maximum deviation is associated with the particles located at the periphery of the deposition. The maximum deviation for various materials of the diaphragm and thickness are tabulated in Table 2. The deviation is mainly due to large deformation as seen for the copper material. In the case of copper, the deformation and the deviation is comparatively less. Thinner specimens have comparatively more deformation and deviation as seen in Table 2. Maximum deviation observed is 0.62° for an aluminium specimen of 0.1-mm thickness. For the drug particle with velocity 310 m/s, the deviation of 0.62° causes a shift of 0.043 mm in the *x*-direction for a travel of 4 mm, which is very useful for focusing drug delivery in micro-scale. Thus, the particle deposition can be extended to the periphery of diaphragm until limiting deviation is reached. Thus, the area of deposition available for a particular diaphragm can be estimated effectively. This study forms a basis to predict the effective penetration of drug particle in the target material, which is discussed next.

Having analysed the travel path of drug particles using finite element analysis, drug particle penetration in the skin or tissue is simulated. Upon reaching the skin, the drug particles impact the skin surface via their gained momentum and reach the epidermis by puncturing micron-sized holes through the stratum corneum [1]. The micro-pore closes immediately after penetration by the surrounding tissue due to pressure. Such, drug-particle penetration studies are performed using an analytical model derived by Dehn [39]. The penetration model is represented as:

$${\mathrm{d}}_{\mathrm{P}}=\frac{4{\mathsf{\rho}}_{\mathrm{P}}{\mathrm{r}}_{\mathrm{P}}}{3{\mathsf{\rho}}_{\mathrm{T}}}\left[\mathrm{ln}\left({\mathsf{\rho}}_{\mathrm{T}}{\mathrm{v}}_{\mathrm{P}}^{2}+3{\mathsf{\sigma}}_{\mathrm{y}}^{\mathrm{T}}\right)-\mathrm{ln}\left(3{\mathsf{\sigma}}_{\mathrm{y}}^{\mathrm{T}}\right)\right],$$

(16)

where ${\mathrm{d}}_{\mathrm{P}}$ is the particle penetration depth, ${\mathsf{\rho}}_{\mathrm{P}}$ and ${\mathsf{\rho}}_{\mathrm{T}}$ are the densities of the drug particle and tissue respectively, ${\mathrm{v}}_{\mathrm{P}}$ is the velocity, ${\mathrm{r}}_{\mathrm{P}}$ is the radius of the particle and ${\mathsf{\sigma}}_{\mathrm{y}}^{\mathrm{T}}$ is the yield stress of the tissue. The penetration obtained is the penetration along the deviated direction. Effective depth of penetration along the direction of shock tube axis as a function of deviation is given by:

$${\mathrm{d}}_{\mathrm{P}}^{\mathrm{E}\mathrm{f}\mathrm{f}}={\mathrm{d}}_{\mathrm{P}}\mathrm{cos}\left(\mathsf{\delta}\right).$$

(17)

The momentum and size of the particles affect the penetration of the skin as given by Equation (16). Since the penetration is instantaneous due to impact, the physiochemical properties do not affect the penetration. Using the penetration model (Equation (17)) and the velocity obtained from simulations, the drug particle penetration was estimated for different properties of the skin. The thickness and material property of the skin vary considerably depending upon the site of drug delivery and age of the person [40]. Wildnauer et al. investigated the yield stress of the stratum corneum of the skin and found the yield stress to vary between 4.9–22 MPa. Studies related to the density of different components of the skin report the density of epidermis to vary between 1110–1190 kg/m^{3} [41]. For the size of drug particle ${\mathrm{r}}_{\mathrm{P}}$ = 15 μm and density ${\mathsf{\rho}}_{\mathrm{P}}$ = 19,300 kg/m^{3}, the penetration depth along the radial coordinate of the diaphragm is shown in Figure 10. The penetration depth decreases with radial location due to the reduction of the velocity at the periphery. Most of the central region has a negligible variation in the penetration depth. At a radial location of 2 mm, the reduction of the penetration observed is 24.4%. Figure 10a shows a negligible variation of drug particle penetration with the density of the skin. However, the yield stress significantly influences drug particle penetration as seen in Figure 10b. At a lower yield stress, drug particle penetration in the skin is significantly higher. Lower yield stress levels of 5 MPa indicate soft skin zones on the body where drug particle penetration is 67% higher as compared to the hard zone having a yield stress of 25 MPa. Thus, to achieve desired drug particle penetration, the calibration of micro-shock tube device considering the thickness and yield stress of the skin is required. The required design and calibration procedure is explained in the next section.

The micro-shock tube based drug delivery device has a potential to deliver a wide variety of drugs like traditional, protein delivery, gene therapy and DNA vaccination [1] by depositing them on the diaphragm surface. To control the drug dosage, the density of drug particles deposited is controlled, keeping the surface area of deposition limited by the deviation of particles as discussed in Section 4.2. Each drug particle (frequently made up of gold or tungsten) has a coating of the drug, which gets directly absorbed in the viable epidermis after penetration. The finite element computational model, particle travel simulation and penetration model can further provide an easy optimization procedure of micro-shock tube-based drug delivery devices with less dependency on expensive experiments as shown in Figure 11. The procedure is briefly explained in the following steps:

Design optimization procedure of micro-shock tube based drug delivery device using simulation and penetration models.

Step 1: The first step of the procedure involves the determination of shock strength form the experiments by comparing the deformation of diaphragms obtained from experiments and computational model reported in Section 1.3.

Step 2: Next, the pressure of the shock is taken as an input for the finite element simulation model and velocity profile of drug particles at the time of detachment from the diaphragm is estimated using the computation model. Further, the velocity profile of the drug particle reaching the skin target is estimated using the simulation procedure elaborated in Section 4.1 and Section 4.2.

Step 3: In the next step, drug particle penetration is estimated using the penetration model as explained in Section 4.3.

Step 4: In this step, drug particle penetration is compared with the desirable range of penetration for the given site of the skin.

Step 5: In this step, the shock tube design is finalized if the drug particle penetration is desirable.

Step 6: In cases where particle penetration is undesirable, the explosive coating length in the explosive tube is altered to vary the shock strength followed by Step 1–5. These, in turn, alters the velocity profile of the diaphragm, thereby altering the penetration depth of the drug particles.

Step 7: This the final step, which involves the testing of many diaphragms for rupture during delivery. Rupture of the diaphragm may cause injury to the skin due to explosive material driving the shock wave. A minimum of 50 diaphragms are tested to after deciding the thickness and explosive coating length. In the case of diaphragm rupture, the thickness in increased and Steps 1–7 are repeated until the diaphragm is stable.

An alternate way to alter drug particle penetration is to change the diaphragm radius, which is difficult since it involves the redesign of diaphragm holder or the diameter of shock tube itself. Following the above procedure, diaphragm thickness and explosive coating length can be determined for possible types of skin types and sites on the body. We now present an example of the drug particle penetration of 0.12 mm in human skin with an aluminium diaphragm, which most of the times is the transdermal zone of the skin. The material properties of the skin usually considered for analysis by researchers [10,11,41] are ${\mathsf{\rho}}_{\mathrm{T}}$ = 1150 kg/m^{3}, ${\mathsf{\sigma}}_{\mathrm{y}}^{\mathrm{T}}$ = 45 MPa, ${\mathrm{r}}_{\mathrm{P}}$ = 15 μm and ${\mathsf{\rho}}_{\mathrm{P}}$ = 19,300 kg/m^{3}. Following the above design and calibration procedure, the average penetration depth for the skin was estimated to be 0.126 mm, which is close to the desired value. To avoid the deviation of the drug particle the peripheral deposition can be reduced and density can be increased. To maintain uniformity of dose, the deposition is to be uniform with minimal deviation. The procedure yields a diaphragm thickness of 0.1 mm and explosive coating length of 55 mm. The designed diaphragm also does not rupture during drug delivery. Thus, the present study demonstrates efficient ways to design shock tube-based drug delivery devices with any actuation material like explosive or pressurized gas.

The paper covers a brief review of various needleless drug-delivery devices. The paper describes the modelling of a diaphragm using a linear hardening model and finite element method for the transient analysis of diaphragm to determine the velocity induced on the drug particles. The final deformation obtained from finite element transient analysis and experimental specimens matched well. The paths of drug particles deposited at various locations on the diaphragm were simulated using the deformation and velocity history. The time histories of deformation and velocities of diaphragm gave the velocity of 335 m/s with a maximum deviation of 0.62° of the velocity vector from the shock tube axis. The range of velocities and deviations obtained for a various material of diaphragm suggests aluminium as the best choice for diaphragm material. Variation of drug particle penetration for different properties of human skin like yield stress and density was studied using an analytical drug particle penetration model. An estimated average particle penetration of 0.126 mm was seen in the human skin. Further, the design and calibration procedure of micro-shock tube is proposed, which enables optimal drug delivery at any site on the human skin. The study presented in the paper provides an easiest way to design and optimize the diaphragm for the required penetration of drug particles. The procedure optimizes the diaphragm and settings of the shock tube device for drug delivery considering the skin site and age of a person.

The authors acknowledge the initial input regarding shock tube design from G. Jagadeesh from Aerospace Engineering, IISc and visiting scientists S.R. Nagaraj and S.G. Rakesh for their important technical input and ideas, which served as background and motivation for this work.

Author Contributions

V.T.R. and D.R.M. conceived and designed the study; V.T.R. performed the experiments and simulations; V.T.R. and D.R.M. analyzed the data and wrote the paper.

Conflicts of Interest

The authors declare no conflict of interest.

1. Arora A., Prausnitz M.R., Mitragotri S. Micro-scale devices for transdermal drug delivery. Int. J. Pharm. 2008;364:227–236. doi: 10.1016/j.ijpharm.2008.08.032. [PMC free article] [PubMed] [Cross Ref]

2. Brown M.B., Martin G.P., Jones S.A., Akomeah F.K. Dermal and transdermal drug delivery systems: Current and future prospects. Drug Deliv. 2006;13:175–187. doi: 10.1080/10717540500455975. [PubMed] [Cross Ref]

3. Jain K.K. Drug delivery systems—An overview. Methods Mol. Biol. 2008;437:1–50. doi: 10.1007/978-1-59745-210-6_1. [PubMed] [Cross Ref]

4. Schramm J., Mitragotri S. Transdermal drug delivery by jet injectors: Energetics of jet formation and penetration. Pharm. Res. 2002;19:1673–1679. doi: 10.1023/A:1020753329492. [PubMed] [Cross Ref]

5. Austin G.D., Salo T.J., Colburn T.J. Single-Use Medicine Delivery Unit for Needleless Hypodermic Injector. 5938637. U.S. Patent. 1999 Aug 17;

6. Bellhouse B.J., Sarphie D.F., Greenford J.C. Needleless Syringe Using Supersonic Gas Flow for Particle Delivery. US 5899880 A. U.S. Patent. 1999 May 4;

7. Bellhouse B.J., Sarphie D.F., Greenwood J.C. Method of Delivering Powder Transdermally with Needless Injector. 5630796. U.S. Patent. 1997 May 20;

8. Weston T.E. Needle-Less Injector. 5891086. U.S. Patent. 1999 Apr 6;

9. Lee S., McAuliffe D.J., Kodama T., Doukas A.G. In vivo transdermal delivery using a shock tube. Shock Waves. 2000;10:307–311. doi: 10.1007/s001930000059. [Cross Ref]

10. Kendall M.A.F. The delivery of particulate vaccines and drugs to human skin with a practical, hand-held shock tube-based system. Shock Waves. 2002;12:23–30. doi: 10.1007/s001930200126. [Cross Ref]

11. Mitchell T.J., Kendall M.A.F., Bellhouse B.J. A ballistic study of micro-particle penetration to the oral mucosa. Int. J. Impact Eng. 2003;28:581–599. doi: 10.1016/S0734-743X(02)00150-1. [Cross Ref]

12. Emilio D.M. Needleless Prefilled Disposable Hypodermic Injector. US 3115133 A. U.S. Patent. 1963 Dec 24;

13. McKinnon C.M., Nakagawa T., Wilcox C.E. Needleless Hypodermic Injection Device. US 5064413 A. U.S. Patent. 1991 Nov 12;

14. Peterson S.F., McKinnon C.N., Jr., Smith P.E., Nakagawa T., Bartholomew V.L. Needleless Hypodermic Injection Methods and Device. US 5520639 A. U.S. Patent. 1996 May 28;

15. Jagadeesh G., Takayama K. Novel applications of micro-shock waves in biological sciences. J. Indian Inst. Sci. 2002;82:1–10.

16. Jagadeesh G., Prakash G.D., Rakesh S.G., Allam U.S., Krishna M.G., Eswarappa S.M., Chakravortty D. Needleless vaccine delivery using micro-shock waves. Clin. Vaccine Immunol. 2011;18:539–545. doi: 10.1128/CVI.00494-10. [PMC free article] [PubMed] [Cross Ref]

17. Kendall M.A.F., Quinlan N.J., Thorpe S.J., Ainsworth R.W., Bellhouse B.J. Measurements of the gas and particle flow within a converging-diverging nozzle for high speed powdered vaccine and drug delivery. Exp. Fluid. 2004;37:128–136. doi: 10.1007/s00348-004-0792-4. [Cross Ref]

18. Arun K.R., Kim H.D. Computational study of the unsteady flow characteristics of a micro shock tube. J. Mech. Sci. Technol. 2013;27:451–459. doi: 10.1007/s12206-012-1259-9. [Cross Ref]

19. Liu Y., Kendall M.A.F. Numerical analysis of gas and microparticle interactions in a hand-held shock-tube device. Biomed. Microdevices. 2006;8:341–351. doi: 10.1007/s10544-006-9596-z. [PubMed] [Cross Ref]

20. Liu Y. Performance studies of particle acceleration for transdermal drug delivery. Med. Biol. Eng. Comput. 2006;44:551–559. doi: 10.1007/s11517-006-0050-4. [PubMed] [Cross Ref]

21. Quinlan N.J., Kendall M.A.F., Bellhouse B.J., Ainsworth R.W. Investigations of gas and particle dynamics in first generation needle-free drug delivery devices. Shock Waves. 2001;10:395–404. doi: 10.1007/PL00004052. [Cross Ref]

22. Rasel M.A.I., Taher M.A., Kim H.D. A study on the gas-solid particle flows in a needle-free drug delivery device. J. Therm. Sci. 2013;22:340–344. doi: 10.1007/s11630-013-0633-y. [Cross Ref]

23. Truong N.K., Liu Y., Kendall M.A.F. Gas and particle dynamics of a contoured shock tube for pre-clinical microparticle drug delivery. Shock Waves. 2006;15:149–164. doi: 10.1007/s00193-006-0034-1. [Cross Ref]

24. Robertson K., Rees J.L. Variation in epidermal morphology in human skin at different body sites as measured by reflectance confocal microscopy. Acta Derm. Venereol. 2010;90:368–373. doi: 10.2340/00015555-0875. [PubMed] [Cross Ref]

25. Liang X., Boppart S.A. Biomechanical properties of in vivo human skin from dynamic optical coherence elastograph. IEEE Trans. Biomed. Eng. 2010;57:953–959. doi: 10.1109/TBME.2009.2033464. [PMC free article] [PubMed] [Cross Ref]

26. Alexander H., Cook T. Variations with age in the mechanical properties of human skin in vivo. In: Kenedi R.M., Cowden J.M., editors. Bed Sore Biomechanics. Macmillan Education; London, UK: 1976. pp. 109–117.

27. Anderson J.D. Modern Compressible Flows. 2nd ed. McGraw-Hill Inc.; New York, NY, USA: 1990.

28. Krajcinovic D. Clamped circular rigid-plastic plates subjected to central blast loading. Comput. Struct. 1972;2:487–496. doi: 10.1016/0045-7949(72)90003-X. [Cross Ref]

29. Li Q.M., Jones N. Blast loading of fully clamped circular plates with transverse shear effects. Int. J. Solids Struct. 1994;31:1861–1876. doi: 10.1016/0020-7683(94)90196-1. [Cross Ref]

30. Neuberger A., Peles S., Rittel D. Spring back of circular clamped armor steel plates subjected to spherical air-blast loading. Int. J. Impact Eng. 2009;36:53–60. doi: 10.1016/j.ijimpeng.2008.04.008. [Cross Ref]

31. Nurick G.N., Gelman M.E., Marshall N.S. Tearing of blast loaded plates with clamped boundary conditions. Int. J. Impact Eng. 1996;18:803–827. doi: 10.1016/S0734-743X(96)00026-7. [Cross Ref]

32. Rajendran R., Lee J.M. Blast loaded plates. Mar. Struct. 2009;22:99–127. doi: 10.1016/j.marstruc.2008.04.001. [Cross Ref]

33. Stoffel M., Schmidt R., Weichert D. Shock wave-loaded plates. Int. J. Solid. Struct. 2001;38:7659–7680. doi: 10.1016/S0020-7683(01)00038-5. [Cross Ref]

34. Ahoranta M., Bukva P., Kováč P., Mikkonen R., Tarhasaari T. Estimation of the stress state of axially tensioned Bi-2223/Ag tapes. Physica C. 2005;432:239–249. doi: 10.1016/j.physc.2005.08.011. [Cross Ref]

35. Reddy J.N. An Introduction to Finite Element Method. 3rd ed. TATA McGraw-Hill; New York, NY, USA: 2005.

36. Chowdhury I., Dasgupta S.P. Computation of Rayleigh damping coefficients for large systems. Electron. J. Geotech. Eng. 2003;8:1–11.

37. Hall J.F. Problems encountered from the use (or misuse) of Rayleigh damping. Earthq. Eng. Struct. Dyn. 2006;35:525–545. doi: 10.1002/eqe.541. [Cross Ref]

38. Singh S. Theory of Plasticity and Metal Forming Processes. 3rd ed. Khanna Publishers; New Delhi, India: 2010.

39. Dehn J.A. Unified theory of penetration. Int. J. Impact Eng. 1987;5:239–248. doi: 10.1016/0734-743X(87)90041-8. [Cross Ref]

40. Wildnauer R.H., Bothwell J.W., Douglas A.B. Stratum corneum properties I. Influence of relative humidity on normal and extracted stratum corneum. J. Investig. Dermatol. 1971;56:72–78. doi: 10.1111/1523-1747.ep12292018. [PubMed] [Cross Ref]

41. Duck F.A. Physical Properties of Tissue. Academic Press; London, UK: 1990.

Articles from Bioengineering are provided here courtesy of **Multidisciplinary Digital Publishing Institute (MDPI)**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |