Numerical prediction of the piezoelectric transducer response in the acoustic nearfield using a one-dimensional electromechanical finite difference approach

TL;DR

This paper introduces a one-dimensional finite difference model to predict the response of a piezoelectric transducer to optoacoustic pressure waves in the acoustic nearfield, facilitating better interpretation of experimental signals.

Contribution

The study develops a simple electromechanical finite difference model that accurately predicts piezoelectric transducer responses in nearfield conditions, considering impedance variations and electromechanical coupling.

Findings

Simulated transducer signals match experimental data well.

The model effectively predicts the influence of impedance changes.

The approach simplifies nearfield transducer response analysis.

Abstract

We present a simple electromechanical finite difference model to study the response of a piezoelectric polyvinylidenflourid (PVDF) transducer to optoacoustic (OA) pressure waves in the acoustic nearfield prior to thermal relaxation of the OA source volume. The assumption of nearfield conditions, i.e. the absence of acoustic diffraction, allows to treat the problem using a one-dimensional numerical approach. Therein, the computational domain is modeled as an inhomogeneous elastic medium, characterized by its local wave velocities and densities, allowing to explore the effect of stepwise impedance changes on the stress wave propagation. The transducer is modeled as a thin piezoelectric sensing layer and the electromechanical coupling is accomplished by means of the respective linear constituting equations. Considering a low-pass characteristic of the full experimental setup, we obtain the resulting transducer signal. Complementing transducer signals measured in a controlled laboratory experiment with numerical simulations that result from a model of the experimental setup, we find that, bearing in mind the apparent limitations of the one-dimensional approach, the simulated transducer signals can be used very well to predict and interpret the experimental findings.

Tables

Table 1. Table 1: Configuration and simulation parameters of the computational domain. From left to right: layer material, [ z − , z + ] subscript 𝑧 subscript 𝑧 [z_{-},z_{+}] -range of layer (in mm mm {\rm{mm}} ), density ρ 𝜌 \rho (in mg / mm 3 mg superscript mm 3 {\rm{mg/mm^{3}}} ), wave velocity c 𝑐 c (in mm / ms mm ms {\rm{mm/ms}} ), absorption coefficient μ a subscript 𝜇 𝑎 \mu_{a} (in mm − 1 superscript mm 1 {\rm{mm}^{-1}} ) and layer impedance Z 𝑍 Z (in mg / ( ms ​ mm 2 ) = 10 2 ​ Rayl mg ms superscript mm 2 superscript 10 2 Rayl {\rm{mg/(ms\,mm^{2})}}=10^{2}{\rm{Rayl}} ).

Material	$[z_{-},z_{+}]$	$\rho$	$c$	${\mu }_a$	$Z$
PMMA	$[0.000,4.980]$	1.18	2.77	0.0	3.3
Glue	$[4.980,4.995]$	1.00	2.50	0.0	2.5
PVDF	$[4.955,5.005]$	1.78	2.25	0.0	4.0
Glue	$[5.005,5.020]$	1.00	2.50	0.0	2.5
PMMA	$[5.020,5.500]$	1.18	2.77	0.0	3.3
Ink	$[5.500,5.550]$	1.00	1.80	100.0	1.8
Glass	$[5.550,9.550]$	2.23	5.60	0.0	12.5

Equations

$$K(z){\partial }_{z}u(z,t)+{\partial }_{t}p(z,t)=0,$$

$$\rho (z){\partial }_{t}u(z,t)+{\partial }_{z}p(z,t)=0.$$

$$[K{D}_{z}u+D_{t}p{]}_i^n=0,$$

$$[\rho D_{t}u+D_{z}p{]}_i^n=0,$$

$$\displaystyle[D_{t}u]_{i}^{n}=\frac{1}{\Delta t}\Big{(}u_{i}^{n+\frac{1}{2}}-u_{i}^{n-\frac{1}{2}}\Big{)},$$

$$\displaystyle[D_{z}u]_{i}^{n}=\frac{1}{\Delta z}\Big{(}u_{i+\frac{1}{2}}^{n}-u_{i-\frac{1}{2}}^{n}\Big{)}.$$

$$\displaystyle u_{i+\frac{1}{2}}^{n+\frac{1}{2}}~{}=~{}u_{i+\frac{1}{2}}^{n-\frac{1}{2}}~{}-~{}\rho_{i+\frac{1}{2}}^{-1}\frac{\Delta t}{\Delta z}\Big{(}p_{i+1}^{n}-p_{i}^{n}\Big{)},$$

$$\displaystyle p_{i}^{n+1}~{}=~{}p_{i}^{n}~{}-~{}K_{i}\frac{\Delta t}{\Delta z}\Big{(}u_{i+\frac{1}{2}}^{n+\frac{1}{2}}-u_{i-\frac{1}{2}}^{n+\frac{1}{2}}\Big{)},$$

$$S=s^{D}T+gD,$$

$$D=dT+{ϵ}^{T}E,$$

$${\partial }_{t}E(z,t)+h{\partial }_{z}u(z,t)=0$$

$$[D_{t}E+h{D}_{z}u{]}_i^n=0,$$

$$E_{i}^{n+1}~{}=~{}E_{i}^{n}~{}-~{}h\frac{\Delta t}{\Delta z}\Big{(}u_{i+\frac{1}{2}}^{n+\frac{1}{2}}-u_{i-\frac{1}{2}}^{n+\frac{1}{2}}\Big{)}.$$

$$\displaystyle U(t)=-\frac{d\ell}{\epsilon^{T}}\Big{(}\frac{1}{\ell}\int_{z_{-}}^{z_{+}}p(z,t)~{}\mathrm{d}z\Big{)}.$$

$$s_i^{\prime }=a{s}_i+(1-a)s_{i-1}^{\prime }.$$

Full text

1 Introduction

Optoacoustics (OAs) can be considered a compound phenomenon, consisting of two distinct processes that occur on different time-scales: fast optical absorption of laser energy inducing a photothermal heating of the absorbing media, and, subsequently, the emission of comparatively slow acoustic stress waves due to thermoelastic expansion and stress field relaxation [1, 2, 3, 4, 5]. Albeit thermoelastic expansion succeeds the absorption of laser light, it is by no means the only energy conversion process that supports the production of optoacoustic signals [6, 4, 3]. E.g., in case of laser generation of stress waves in liquid one might identify three relevant mechanisms, i.e. thermoelastic expansion, vaporization, and dielectric breakdown, occurring at increasing deposited power-densities [6]. Owing to the possible generation and propagation of transverse vibrations, the production of acoustic waves in solid is somewhat more intricate [2]. However, in the absence of any change in state of the underlying medium, thermoelastic expansion due to the absorption of laser energy can be considered the dominant conversion mechanism [6, 4]. Compared to the propagation of acoustic stress waves, which proceeds on a microsecond timescale, the optical absorption is assumed to occur instantaneously. Assuming the optical absorption to be instantaneous has consequences for the theoretical treatment of the problem inherent dynamics [7, 8, 5]. It not only allows to decouple the optical absorption problem from the acoustic propagation problem but also allows to simplify the latter as discussed in the remainder.

Here, we present a combined study, complementing measurements on a controlled experimental setup with custom numerical simulations in terms of a finite difference model of the underlying physical processes. In contrast to recent studies wherein we discussed measurement, simulation and approximate inversion of OA signals observed for layered PVA-H (polyvinyl alcohol hydrogel) phantoms in the acoustic farfield [9, 10], the specific object of the presented study is to model the observed PVDF (polyvinylidenfluorid) transducer response resulting from the subtleties of the source volume in the acoustic nearfield. The acoustic properties of the source volume, determining the propagation of the stress waves and their behavior upon crossing inter-layer boundaries, are considered to be piecewise homogeneous and not too rapidly varying (as, e.g., opposed to the scenario considered by Refs. [11, 12]). If the experimental setup satisfies nearfield conditions, i.e. under the assumption of plane acoustic waves, and for a translational symmetry of the region of interest in the plane perpendicular to the direction of the propagating stress waves, the evolution of the stress profile within the domain can be modeled in terms of the equations of one-dimensional ($1$D) linear elasticity [13, 4, 14]. The response of a piezoelectric sensing layer to trespassing stress waves is then included by the electromechanical coupling to the constituting equations of linear piezoelectricity [15]. The article presents an application for a simple and efficient effectively $1$D approach for the solution of a $3$D problem. Note that quite similar $1$D approaches where considered in the literature to study, e.g., complex transients in pulsed OA spectroscopy [16], surface heating and energy transfer in pulsed microwave catalytic systems [17], and, thin-film piezoelectric ultrasonic sensors [18].

The article is organized as follows. In section 2 we discuss the computational model for the numerical simulation of stress wave sensing in piecewise homogeneous elastic media in detail. Therefore, we discuss the finite difference stencils used to approximate the underlying differential equations of the continuous models (i.e. linear elasticity and linear piezoelectricity) and illustrate their implementation in terms of python modules111For conciseness, the presented code listings are not documented thoroughly. In general, this has to be considered bad programming style. However, note that the code in the supplementary material under Ref. [19] is documented well. that are later used for the numerical experiments discussed in section 4. In section 3 we point out limitations of the presented approach and conclude with a summary in section 5.

2 Modeling stress wave sensing in a piecewise homogeneous elastic medium

In our one-dimensional approach, the computational domain represents a heterogeneous elastic medium for which the acoustic displacement $r(z,t)$ at the field point $z$ defines the material velocity field via $u(z,t)=\dot{r}(z,t)$. Material velocity $u(z,t)$ and excess pressure $p(z,t)$ are related via the set of coupled first-order differential equations of linear elasticity [13, 14]

[TABLE]

Therein $\rho(z)$ denotes the local density and $K(z)\equiv\rho(z)c(z)^{2}$ signifies the bulk modulus of elasticity wherein $c(z)$ refers to the local speed of sound.

Subsequently we discuss the finite-difference approach used to model the piezoelectric “sensing” of acoustic stress waves. In subsection 2.1 we introduce the data structure used to represent the piecewise homogeneous medium. In subsection 2.2 we elaborate on the mechanism of optoacoustic signal generation, responsible for the generation of initial stress profiles within our model. The finite difference approach for the propagation of the acoustic stress waves in piecewise homogeneous elastic media is illustrated in subsection 2.3. The coupling to the equations of state of linear piezoelectricity is discussed in subsection 2.4. Finally, a simple postprocessing strategy to account for the disturbance of the transducer signal by the experimental setup is discussed in 2.5.

2.1 Discretization of the computational domain

In principle, the above set of equations describes a continuous model. However, so as to be able to numerically solve these equations for reasonable boundary conditions and initial values, we need to consider a discrete realization of the computational domain that holds local material properties such as $\rho(z)$ and $c(z)$. Later, we will consider the propagation of laser generated acoustic stress waves within the medium. Therefore, the optical properties of the medium, responsible for the absorption and scattering of laser light are of importance. In our simplified $1$D model we will consider the limiting case of purely absorbing media. Therefore we here introduce a further local material property, namely the absorption coefficient $\mu_{a}(z)$. In order to set up a data structure that holds all these local material properties, consider a discrete $z$-grid with constant mesh width, i.e. $z_{i}=z_{\rm{min}}+i\Delta z$ for $i=0\ldots N-1$ and $\Delta z=(z_{\rm{max}}-z_{\rm{min}})/N$, and let, e.g., $\rho_{i}$ refer to $\rho(z_{i})$. Then, a proper domain data structure might be implemented by the python class Domain in code listing 1. Subsequently we assume that the respective code is available in the python module file domain.py.

2.2 Optical absorption of laser energy

In our model, the propagation of acoustic stress waves is triggered by an initial distribution of acoustic stress $p(z,0)=p_{0}(z)$ for a medium at rest $u(z,0)=0$. The nonzero initial pressure profile results from a thermoelastic conversion of deposited laser energy to mechanical stress within the medium. In comparison to the typical timescale of mechanical response, the deposition of laser energy can be considered instantaneous. Below, the material parameter $\mu_{a}(z)$ accounts for the absorption of photons within the medium. As a consequence, considering a light flux in $z$ direction, a laser beam will experience a decrease of its incident fluence $f_{0}$ with increasing depth $z$ following the Beer-Lambert decay law $f(z)=f_{0}\exp\{-\int_{z_{\rm{min}}}^{z}\mu_{a}(z^{\prime})~{}\mathrm{d}z^{\prime}\}$ [20]. For such a purely absorbing medium, the amount of locally absorbed laser energy $W(z)$ is related to the respective fluence decay through $W(z)=-\mathrm{d}f(z)/\mathrm{d}z=\mu_{a}(z)f(z)$ [21]. Finally, the efficiency of the conversion of deposited laser energy to acoustic stress is governed by the Grüneisen parameter $\Gamma$, by means of which $p_{0}(z)=\Gamma\,W(z)$. Note that the initial pressure pulses obtained in this manner exhibit abrupt changes along the $z$ axis, signaling a sudden increase or decrease of the absorption coefficient. Using these “shockwave” initial conditions in the acoustic propagation algorithm implemented below might cause numerical artefacts in the observables. In such a situation, if the absorption coefficient is nonzero in a small range $z\in[z_{p},z_{p}+\delta]$ only, a remedy might be to use a simple Gaussian function with peak intensity $f_{0}$, $1/e$ extension of $\delta$, and centered at $z=z_{p}+\delta/2$, instead. Albeit this does not properly model the exponential attenuation of laser fluence, it allows to study the principal distortion of pressure profiles upon propagation. Another possibility to circumvent the challenges such a shockwave might cause for a finite-difference scheme is to adopt a high-precision finite-volume procedure, see the discussion in sec. 5.

Once the computational domain is initialized and the optical properties declared, a virtual laser beam can be propagated through the medium by means of the function absorbLaserBeam in code listing 2. On input it expects three arguments, i.e. the details of the computational domain Glob, the Grüneisen parameter Gamma and the initial laser fluence f0. On output it yields a tuple consisting of the initial pressure profile as a function of $z$ and the fluence of the transmitted part of the beam. Further, an implementation of a simple Gaussian pressure profile as discussed above is provided by the function deltaPulse.

2.3 Acoustic propagation via finite differences

The above hyperbolic set of partial differential equations might be discretized via a staggered-grid leapfrog scheme [22]. Therefore, consider an additional $t$-grid with constant mesh width $\Delta t$, i.e. let $t_{n}=t_{\rm{min}}+n\Delta t$ for $n=0\ldots M-1$ and $\Delta t=(t_{\rm{max}}-t_{\rm{min}})/M$, and let, e.g., $p_{i}^{n}$ refer to $p(z_{i},t_{n})$. Then, discretization yields the numerical approximations

[TABLE]

of Eqs. (1a) and (1b), respectively, wherein the centered half-step grid-derivatives for the exemplary variable $u$ read

[TABLE]

For convenience, shifting the velocity field to the intermediate coordinates, i.e. letting $n\to n+1/2$ in Eq. (2a) and $i\to i+1/2$ in Eq. (2b), and using the harmonic mean $\rho_{i+\frac{1}{2}}^{-1}=(\rho_{i}^{-1}+\rho_{i+1}^{-1})/2$, we obtain the evolution equations

[TABLE]

for the acoustic velocity and excess pressure fields. On either end of the computational domain, the boundary conditions implement pressure-release boundaries, i.e. we impose $p_{0}^{n}=0$ and $p_{N-1}^{n}$ at each timestep, defining a free surface. A vectorized python implementation of the above formulae is shown in code listing 3. The first $4$-tuple of parameters specify the initial values for the excess pressure (p0) and the details of the computational domain, i.e. sonic velocity (v), density (rho), and mesh width (dz). The parameter measure refers to a function that facilitates the monitoring of relevant observables. This callee will be defined later on. However, note that it receives four arguments that specify the current time step (n), the time increment (dt), the full velocity field (u) and excess pressure (p). Finally, Nt signifies the overall number of simulation time steps until termination of the propagation process.

2.4 Realizing the piezoelectric coupling

The implementation of a piezoelectric sensing layer requires the electromechanical coupling of the equations of linear elasticity to the state equations of $1$D piezoelectricity. Here we consider a piezoelectric layer as a stress sensing device, relying on the direct piezoelectric effect in which a mechanical load, i.e. stress within the material, is converted to an electric field. If electrodes are connected to the opposing layer surfaces, the resulting potential difference across the layer might be measured. The theoretical framework in which the piezoelectric response might be described depends on mechanical and electrical boundary conditions (BCs). Depending on the subtleties of the experimental setup that needs to be modeled, one distinguishes a mechanically free setup (i.e. at constant stress), labelled “T”, and a mechanically clamped setup (i.e. at constant strain), labelled “S”, as well as an electrical short-circuit setup (i.e. at constant electrical field), labelled “E”, and an electrical open-circuit setup (i.e. at constant electrical displacement), labelled “D”. Here, we consider the idealized case of a mechanically free setup with an open-circuit. The respective $1$D constituting relations read

[TABLE]

wherein $E$, $D$, $S$ and $T$ refer to the electric field, field displacement, strain and stress, respectively. Further, $s^{D}$ denotes the mechanical compliance, $d$ refers to the piezoelectric strain constant and $\epsilon^{T}$ signifies the dielectric coefficient at constant stress. Note that under the assumption of a vanishing charge density in the unstressed state of the layer, the field displacement $D$ is zero. Further, considering the displacement gradient formulation of the strain, i.e. $S(z,t)=\partial_{z}r(z,t)$ and thus $\partial_{t}S(z,t)=\partial_{z}u(z,t)$, allows to cast Eqs. (5a) and (5b) into the computationally convenient form

[TABLE]

of the direct piezoelectric effect, wherein $h=d/(\epsilon^{T}s^{D})$. Again, using centered half-step grid derivatives as above, this can be approximated by the finite difference equation

[TABLE]

providing an evolution equation for the electric field in the piezoelectric layer in the form

[TABLE]

Finally, the potential difference $U(t)=\phi(z_{+})-\phi(z_{-})$ between the opposing layer surfaces (located at $z_{-}$ and $z_{+}$) can be obtained from $E(z,t)=-\partial_{z}\phi(z,t)$ by numerical quadrature, e.g. using a trapezoidal rule. Note that, considering the constituting equation of the direct piezoelectric effect, i.e. Eq. (5b) together with the $1$D stress-pressure relation $T(z,t)=-p(z,t)$ of hydrostatic compressions [13], the potential difference $U(t)$ between the opposing layer surfaces (located at a distance $\ell=z_{+}-z_{-}$) can be found to be

[TABLE]

Thus, the potential difference is simply proportional to the average pressure within the piezoelectric sensing layer, i.e. $U(t)\propto\bar{p}(t)$. As pointed out earlier, the above finite difference scheme was derived for the idealized setup considering the mechanically free BCs. Simulations for a mechanically clamped setup can easily be done by setting $h=e/\epsilon^{S}$ in Eq. (8), wherein $e$ refers to the piezoelectric stress constant and $\epsilon^{S}$ signifies the dielectric constant at fixed stress. This follows from considering the piezoelectric constituting equations for BCs “S” and “D”.

A proper transducer data structure that implements the finite difference approximation to the linear piezoelectric constituting equations is given by the python class PiezoTransducer shown in code listing 4. An instance of the transducer class expects three arguments on input: an instance Glob of the discrete computational domain, a tuple (zMin, zMax) specifying the surface locations of the sensing layer and an optional argument h that combines the piezoelectric and mechanical system parameters as discussed above. Note that the class provides a method called measure that can be passed as callback function to the propagateStressWaveSLDE routine, assuming the role of the function measure in code listing 3.

2.5 Signal postprocessing in the time domain

Note that the experimental setup used to amplify and postprocess the piezoelectric response might have an impact on the shape of the signal itself. In order to mimic such a disturbance one might pursue one out of several postprocessing strategies, based on digitally filtering the detected signal in the time- or frequency domain [22]. Such a filter takes an input signal $s(t)$, i.e. a time-series of input points, and yields a modified output signal $s^{\prime}(t)$ subject to several physical constraints. E.g., a time-domain filter might work online, processing input points as they are recorded, or offline, as a convenient postprocessing device. While the former working mode naturally ensures causality since it prohibits the filter to access input points that are out of time, yet, the latter working mode allows also for a more general behavior wherein a certain input point might depend on earlier as well as later input points. Here, for convenience, we opt for an offline filter that works in the time domain and obeys causality. Without elaborating on the subtleties of the experimental setup used in our laboratory experiments, we here assume that the experimental setup has the characteristics of a first order low pass filter, suppressing high frequency features of the input signal. This choice is purely phenomenologic since the effect of such a filter is consistent with our observations.

In general, a low-pass filter is characterized by time constant $\tau$, defining its “cutoff frequency” $\omega_{\rm c}=\tau^{-1}$. Signal features with a frequency $\omega>\omega_{\rm c}$ appear suppressed. If the input signal is sampled at constant time increment $\Delta t$, one might also define the filter parameter $a=\Delta t/(\tau+\Delta t)$, with $0\leq a\leq 1$, to characterize its performance. Denoting input and output signals as $s$ and $s^{\prime}$, respectively, the effect of a simple first order low-pass filter might be cast into the recurrence relation

[TABLE]

A data structure that implements such a low-pass filter is shown in code listing 5. Therein an instance of the time domain filter class Filter expects two arguments on input: the ordered sequence of time samples t and the input signal s0. The low-pass is provided as method lowPass. As an argument, the method takes the characteristic smoothing parameter a of the filter and yields an output signal that contains the distortion of the input signal by the low pass. The filter described above can be used to postprocess the detector signals recorded via the piezoelectric transducer implemented in code listing 4.

3 Limitations of the one-dimensional approach

Since the theoretical approach implemented in sect. 2 is tailored to a particular physical scenario, there are limitations that restrict its range of application. Below, without any claim on completeness, we address some of these limitations.

In principle, considering a nonscattering $3$D setup wherein $\vec{r}=(x,y,z)$, the deposition of energy within the medium is governed by laser beam parameters such as the incident laser fluence at maximal intensity $f_{0}$, the transverse beam profile $f(x,y)$ and the temporal intensity profile $I(t)$ of the laser pulse, as well as the absorption coefficient $\mu_{a}(x,y,z)$ within the medium. Two prerequisites allow for an effectively $1$D treatment of the problem of optical absorption: (i) the absorption coefficient depends only on the depth coordinate $z$, i.e. $\mu_{a}(x,y,z)\equiv\mu_{a}(z)$, thus realizing a layered medium with translational invariance in the $(x,y)$-plane, and, (ii) the transverse beam profile is wide enough to ensure initially plane acoustic wavefronts. This can be accomplished by requiring the penetration depth $\ell$ of the laser into the medium, defined via $1\equiv\int_{z_{\rm{min}}}^{\ell}\mu_{a}(z^{\prime})~{}\mathrm{d}z^{\prime}$, to be smaller than the characteristic lengthscale, say, e.g., the beam diameter $d$, that defines the extension of the laser spot on the medium (assuming a constant laser intensity within the spot). Then, the validity of the $1$D approach is limited by the onset of diffraction during the propagation of the initial pressure profile along the $z$ axis.

In our numerical approach, the propagation of stress waves is accomplished by a finite-difference approximation of the equations of linear elasticity, i.e. Eqs. (1a) and (1b), under the assumption of initial plane acoustic waves. While this approach seems sufficient for our particular application, see the results reported in sect. 4, there are several options to extend the numerical procedure. E.g., if radial inertia following the optical absorption of laser energy might not be neglected, Love’s modified wave equation might be used to include the effect of dispersive waves during the propagation process [23]. If elastic media with losses, as, e.g. tissue, are considered, a Kelvin-Voigt model wave equation might be used to account for viscoelastic effects [15]. Further, note that if the assumption of initial plane acoustic waves and stress wave detection in the acoustic nearfield are not satisfied, a more complete $3$D description of the problem that naturally accounts for the effect of acoustic diffraction is needed. Also note that since we study the propagation of stress waves prior to thermal equilibration of the source volume, which occurs on larger timescales, we neglect diffusion of heat within the computational domain.

For our numerical experiments we consider a PVDF (polyvinylidene fluoride) polymer as piezoelectric sensing layer. Note that PVDF exhibits not only piezoelectric but also pyroelectric properties [24]. Thus, if the temperature of the PVDF layer is expected to increase notably within the duration of the measurement process, pyroelectric effects might be expected in addition to the direct piezoelectric effect discussed in subsect. 2.4. However, in the particular source volume configuration studied in the presented article, the initial stress pulse is separated from the PVDF layer by an approx. $0.5\,{\rm{mm}}$ wide PMMA layer. Since the thermal diffusivity of PMMA amounts to $\alpha=0.1054\,{\rm mm^{2}/s}$ (at a temperature of $T=25^{\circ}$, see Ref. [25]), the approximate thermal relaxation time for thermal equilibration over the former distance is $t_{\rm{r}}=9.49\,{\rm{s}}$, exceeding the observation time by several orders of magnitude. Hence, as pointed out above, we neglect heat diffusion and thus also pyroelectric effects in the PVDF layer.

4 Results and Discussion

The numerical approach detailed in sect. 2 allows to implement a piezoelectric sensing sensing layer with or without acoustic backing layer. Both setups allow for notably different signals that might be expected. Here, in order to perform custom numerical simulations for an experimental setup we seek to match, we consider a PVDF layer glued (using Norland Optical Adhesive “NOA85V”) in between two layers of PMMA (polymethylmethacrylat), serving as backing- and frontlayer. As illustrated in Fig. 1, we assume that the deposition of laser energy in the computational domain occurs in a thin layer of black ink in between the PMMA frontlayer and a closing plate of glass. The configuration of the domain, i.e. the thickness of all layers along with their assumed optic and acoustic properties, is listed in Tab. 1. Further details regarding the experimental setup are discussed elsewhere [26].

Analysis of the transducer response

Starting from an initial stress profile, obtained following the OA signal generation procedure outlined in subsect. 2.2 and evolved using the staggered-grid leapfrog scheme presented in subsect. 2.3, Fig. 2(a) shows the time sequence of the potential difference $U(t)$ between the boundaries of the PVDF layer. As argued in subsect. 2.4, the resulting Voltage signal is proportional to the average pressure in the piezoelectric sensing layer. Assuming that the intrinsic lengthscale $\lambda_{\rm{ac}}=2\ell$ (wherein $\ell$ is the penetration depth of the laser defined in sect. 3) that might be used to define the acoustical wavelength in the given setup exceeds the thickness of the piezoelectric layer, the Voltage $U_{\rm{s}}(t)$ caused by the trespassing stress wave is simply proportional to the local pressure at, say, the surface facing the PMMA frontlayer, i.e. $U_{\rm{s}}(t)\propto p(z^{\rm{PVDF}}_{+},t)$ [27, 28] (in our numerical experiment $z^{\rm{PVDF}}_{+}=5.005\,{\rm{mm}}$). If the width of the PVDF layer exceeds $\lambda_{\rm{ac}}$ notably, both Voltage signals are expected to differ significantly [27]. We verified this behavior by means of further numerical experiments based on a simpler design of the computational domain (not shown). In this particular application we assume the width of the PVDF layer $\Delta z_{\rm{PVDF}}=10\,{\rm{\mu m}}$ and the acoustical wavelength $\lambda_{\rm{ac}}=20\,{\rm{\mu m}}$ to be approximately of the same extend. As evident from Fig. 2(a), both signals are (still) in reasonable agreement as expected for a “thin” sensing layer and a “wide” trespassing stress wave.

An implementation of the computational domain using the simulation parameters listed in Tab. 1 is shown in code listing 6. Provided that the imported modules are in the local searchpath, the necessary statements to produce a plot similar to Fig. 2(a) requires merely $48$ lines of python code.

Note that the time sequence of the transducer response $U(t)$ features several groups of features: the first group reaches the PVDF layer in the interval $t\in(0.15,0.45)\,{\rm{\mu s}}$, resulting from the initial stress wave passing the Ink-PMMA (reflection coefficient $C_{r}({\rm{Ink}},{\rm{PMMA}})\approx 0.29$) and PMMA-PVDF inter-layer boundaries also after multiple reflections within the Ink layer, enclosed bye PMMA and Glass ($C_{r}({\rm{Ink}},{\rm{Glass}})\approx 0.74$). Since both relevant reflection coefficients are positive, the first group of features is a train of compression peaks. The first reflection with a negative reflection coefficient occurs at the PMMA-Glue inter-layer crossing ($C_{r}({\rm{PMMA}},{\rm{Glue}})\approx-0.13$), triggering a rarefaction wave that travels back towards the PMMA-Ink interface where its sign changes again due to the apparent coefficient $C_{r}({\rm{PMMA}},{\rm{Ink}})\approx-0.29$, approaching the sensing layer, once again, as compression wave. However, due to multiple inter-layer reflections that involve coefficients $C_{r}<0$, rarefaction waves eventually trespass the PVDF layer, see, e.g. the second group of features reaching the PVDF layer in the interval $t\in(0.50,0.75)\,{\rm{\mu s}}$. The sequence of sign changes in the $1$D pressure profile $p(z,t)$ can be disentangled best using the space-time plot shown in Fig. 2(b).

Assuming that the impact of the experimental setup on the transducer signal can be described by means of a low-pass filter with time constant $\tau\approx 0.05\,{\rm{\mu s}}$ effectively smoothes the input transducer signal and yields an output signal as shown in Fig. 3(a). Therein, the grouped features of the input signal are partially lost due to the attenuation of frequencies above the cutoff $\omega\approx 20\times 10^{6}\,{\rm{rad/s}}$ (i.e. $f=\omega/2\pi\approx 3.2\,{\rm{MHz}}$) in the ultrasonic frequency range.

This completes the interpretation and analysis of transducer signals obtained for the computational domain shown in Fig. 1. The comparison of the numerical results to experimental data follows below.

Comparison to experimental results

In Fig. 3 we compare the numerical simulations to experimental data, obtained from measurements of OA stress waves via a PVDF transducer using a configuration similar to Fig. 1. The figure allows to compare the predicted transducer response (including the assumed low-pass distortion due to the experimental setup), to the measured response, cf. the black dashed lines in Figs. 3(a),(c). The inferred transducer response, obtained by correcting for the transducer transfer-function in the experiment [26], also compares well to the calculated response of the PVDF layer, cf. the solid lines in Figs. 3(a),(c). As evident from the figure, the sequence of compression and rarefaction pulses observed in the numerical and laboratory experiments are in excellent agreement. Note that the group of rarefaction pulse signal features in the range $t\in(1.6,1.9)\,{\rm{\mu s}}$, resulting from a reflection of the acoustic wave on the far end of the glass layer, has a notably lower amplitude in the experiment. This might be due to the circumstance that the stress wave needs to traverse the liquid-solid ink-glass and ink-PMMA inter-layer boundaries possibly several times. Thus, loss effects that are relevant in crossing such a boundary, not accounted for in our approach, might be necessary to explain this observation.

Verification using a finite-volume approach

Modeling wave propagation in nonconservative hyperbolic systems, such as linear elastic wave propagation in varying heterogeneous media governed by Eqs. (1a)-(1b), render a challenge for finite-difference methdos. Alternatively one might consider high-resolution finite-volume methods originally developed for nonlinear problems [12, 29]. Albeit computationally more expensive than the finite difference approach illustrated in sect. 2, they allow for a solution of problem instances that are not amenable to finite-difference methods, such as, e.g. shock-wave propagation in highly discontinuous nonlinear media [30]. For numerical redundancy we verified the results, i.e. location of signal features and ordering of compression and rarefaction pulses, obtained using our (simple) finite-difference approach via an independent (elaborate) finite-volume implementation in terms of an acoustic $1$D Riemann solver, provided by the PyClaw CLAWPACK tool [31, 32], cf. Figs. 3(a) and (b).

5 Conclusions

In the presented article we considered the problem of optoacoustic generation and propagation of stress waves within a piecewise homogeneous material, with focus on their detection using a piezoelectric sensing layer. The specific objective of this study was to implement a simple numerical model that facilitates the simulation of the transducer response in the acoustic nearfield in terms of an effectively $1$D finite difference approach, and to complement experimental results via numerical simulations. Comparing numerical simulations and experimental data obtained for a certain layered setup of the source volume, we found that the portrayed numerical approach accurately predict its elastic and piezoelectric response. In turn, the modeling of the piezoelectric transducer response in the presented approach proved useful for interpreting measured transducer signals and for verifying the assumed transfer-function of the employed transducer.

Since the non-availability of code impedes transparency and reproducability of results in scientific publications [33, 34, 35] we considered it useful to make the concise research-code for the presented study, including but not limited to the code listings 1 through 6 along with all scripts needed to reproduce all figures, publicly available [19].

Acknowledgments

This research work received funding from the VolkswagenStiftung within the “Niedersächsisches Vorab” program in the framework of the project “Hybrid Numerical Optics” (HYMNOS; Grant ZN 3061). The software was developed and tested under OS X Yosemite (Version: 10.10.3) on a MacBook Air featuring a 1.7GHz Itel Core i5 processor and 4 GB DDR3 using python [36] version $2.7.6$ and numpy version $1.8.0{\rm{rc}}1$ [37].