1 Introduction

Recently, terahertz (THz) radiation in the range between 0.1 and 3 THz has been used to analyze the internal structure of a large number of materials and elements. Such materials also include multilayer composite materials, which are opaque for visible and infrared light. Moreover, the analysis using THz radiation allows to determine thicknesses of internal layers and to define hidden defects and non-uniformities inside those materials [1–3]. It is worth noticing that THz radiation provides a non-destructive, non-ionizing, and contactless method for material evaluation. THz techniques can compete or at least be used complementary with thermographic [4], ultrasonic [5], acoustics [6], X-ray [7], and microwave [8] material examination methods.

The possibility to examine different composite materials in the THz range strongly depends on the attenuation caused by the respective materials. As a consequence, carbon fiber composites cannot be investigated with THz waves due to the high reflectivity for THz radiation which almost precludes penetration inside the internal structure [9]. On the other hand, materials exhibiting much better properties for the examination in the THz range are glass fiber-reinforced composites. These and many other materials are being investigated with THz radiation even though there are still problems with the transmission through the material and in addition with scattering in the internal structure [3, 10].

For our investigations, we have chosen ultra-high molecular weight polyethylene (UHMWPE) as a composite material because it has a low absorption coefficient varying between 0.3 and 1.5 cm−1 over the whole spectral range. Moreover, the refractive index of n = 1.53 [11] is almost constant over a very wide range of frequencies and, thus, the influence of dispersion is very small. Furthermore, such composites are used in manufacturing bulletproof jackets, helmets, and ballistic shields because of their high mechanical resistance. The determination of thicknesses of their layers, investigating delaminations and other defects, is crucial for reliable characterization of the manufacturing process and quality control. Up to now, the UHMWPE composites have not been intensely studied in the THz range [12, 13].

The non-destructive evaluation of composite materials [1–3, 14–16] was mostly carried out by the use of THz time-domain spectroscopy (TDS). Nevertheless, there are many other techniques based, for example, on FMCW radar [17], CW imaging using quantum cascade lasers [18], a gas laser [9], or a Gunn diode [19]. In this work, our THz TDS system is used in reflection configuration in order to obtain information about the structure of the multilayer sample including the distribution and thicknesses of particular internal layers. Additionally, we can get information about the type and the position of various defects. Transmission imaging is only used to determine composite material parameters (cumulative information about the sample [20]).

In reflection geometry, the THz pulses are partly reflected from the front surface and from the interfaces between the media having different refractive indices (like layers/plies). In addition, the THz pulses are scattered or reflected from the back surface. Because of the multilayered structure of composite materials, a sequence of pulses shifted in time is observed each pulse coming from a reflection from a particular interface. Generally, knowing the refractive index of the single layer, we can determine its thickness on the basis of the time difference between two pulses (their peaks) reflected from the front and back surfaces of this layer [13]. However, for thin layers, roughly about 50 μm, pulses reflected from adjacent layers overlap and cause errors in determining thicknesses of consecutive layers.

To solve this problem, van Mechelen et al. [21] assume that the investigated material that composed of a few layers is a stratified system, and for each layer, both its thickness and its refractive index are optimized. Krimi et al. have recently proposed a time-domain fitting procedure for a structure composed of few thin layers of different thickness, useful in inspection of automotive paints [22]. This approach is based on the transfer matrix method (TMM) and requires prior knowledge of the material parameters to determine the thicknesses of the layers.

Here, the algorithm presented in [22] was adapted to a composite sample composed of about 200 layers with a quasi-binary profile of the refractive index, which means that values of the refractive index for subsequent layers are as follows: n + Δn 1, n − Δn 2, n + Δn 3,…, respectively; n > > Δn i and Δn 1≈Δn 2≈Δn 3…. The optimization of both parameters, i.e., the refractive index and the thickness for each layer independently, would be computationally demanding. Therefore, Krimi’s algorithm [22] was improved in such a manner that calculations were extended to the optimization of thicknesses of all 200 layers. We assume Δn 1 = Δn 2 = Δn 3 = Δn and Δn was in general optimized for the whole sample. Due to the larger number of layers and thus more degrees of freedom for the optimization algorithm, it is very difficult to apply directly the optimization approach proposed in [22]. For each layer, the search interval was about 60 μm. For 200 layers, this corresponds to 20060 possible solutions if we have a resolution of 1 μm. The probability that the algorithm finds a better solution in this vast space is very low. For this reason, at first, we use a brute force method with a resolution of 10 μm to reduce the search space. With this hybrid approach, we reduce the search interval per layer from 60 to 10 μm and, thus, we significantly increase the optimization efficiency and accuracy. Using optimized values of thicknesses, we reconstruct the waveform. It is worth noticing that there was excellent agreement between measurements and simulations.

The paper is organized as follows: The internal structure and the properties of composites made from UHMWPE are described in Section 2. Next, the results of the refractive index measurements for two bulk samples and for the sample composed of four single plies (in two orthogonal orientations of the electric field vector) are presented in Section 3. Afterwards, the theory of the reflection and the description of the simulation method are presented in Section 4. Section 5 describes the experimental TDS setup and the results of our measurements. The fitting procedure was applied to the reflected waveforms and allowed to calculate thicknesses of all layers. To the best of our knowledge, it is the first time that thicknesses of all layers from a multilayer (~200 plies) structure were determined by means of THz TDS technology with such a precision. Moreover, the origin of the observed birefringence is discussed in Section 6.

2 The Characteristic of the Samples

UHMWPE consists of extremely long chains of polyethylene (Fig. 1a) with a simple repeated atomic structure and a molecular mass of about 2 million u [23]. The intermolecular strength of this material is high as it is caused by longitudinal polyethylene chains bonded by many Van der Waals bonds. UHMWPE is a very tough material, and it is resistant to many external influences and has high impact strength [24, 25]. This composite is available on the market in the form of fibers or bulk material, mostly as sheets or rods.

Fig. 1
figure 1

Sketch of the long chains of UHMWPE (a). SEM photograph of fibers of BT10 composite (b). Photographs of the BT10 (c) and HB50 (d) composites

Full size image

This paper concentrates only on fiber-based composites. Fibers are extruded from the heated gel of UHMWPE using a spinneret, and more than 95 % of polymer chains have parallel orientation [25, 26]. Fibers having different diameters are frequently used for manufacturing ropes, nets, and protective textiles. Moreover, flexible tapes made of perpendicular layers (often called plies) consisting of fibers are tailored for helmets, vests, or ballistic shields [12].

Based on two kinds of tapes from Dyneema® company (BT10 and HB50), two kinds of composites were manufactured in the Polish company, Lubawa. Therefore, both composites have their names derived from the names of the tapes. The BT10 is the laminated resin-free tape which consists of a few woven and interleaving plies of PE fibers in two perpendicular orientations (0°/90°). The diameter of a single fiber is about 15 μm (Fig. 1b), and the interleaving period is about 100 mm. The thickness of a single ply was about 40–50 μm [13] as measured with a micrometer screw as well as with a microscope. The arrangement of plies inside an HB10 tape in this process is random, and due to its woven nature, there are regions with regular (0°–90°–0°–90°) and irregular (90°–0°–0°–90°) sequence of layers (Fig. 2). More information about this tape was not available to the authors.

Fig. 2
figure 2

The arrangement of UHMWPE tapes and composites made out of BT10 and HB50 tapes

Full size image

The non-woven HB50 tape consists of four alternately perpendicular plies of PE fibers with diameters of about 17 μm each [24]. The thickness of a single ply is about 55–70 μm, and a plastic rubber (17 % of styrene-isoprene-styrene triblock copolymer (SISTC) [24]) is used as a matrix. In this case, the orientation of plies is regular (0°–90°–0°–90°) due to the fact that the non-woven HB50 tape consists of a regular sequence of plies having mutually perpendicular orientations (Fig. 2).

Long tape reels of both kinds were cut into smaller sheets which were then layered and hot pressed to form a multilayer composite material (Fig. 2). The photographs of BT10 and HB50 composite samples having the dimensions of 50 × 50 mm2 and thicknesses 10.40 and 9.45 mm, respectively, are shown in Fig. 1c, d. The thickness of each sample varies in the range ±0.1 mm. The exact number of plies in both samples is not known to the authors.

3 Characterization of UHMWPE in Transmission

We investigated the multilayer BT10 sample having a smaller thickness of q = 3.25 mm using a standard time-domain spectroscopy system (Spectra 300) from TeraView. A relatively thin sample was used, because the measurement of a thicker sample may be distorted due to the relatively short Rayleigh range of the beam in the system used for the measurements. The determination of the refractive index (n) and the extinction coefficient (κ) is very important for further analysis of experimental results and simulation parameters. Both values (n and κ) were determined using THz TDS as described in [27].

Figure 3a presents three experimental waveforms: the reference and two waveforms after transmission through the sample for two orientations (0° and 90°) of the electric field with respect to the arbitrary chosen edge of the sample. The spectral characteristics of the refractive index and the extinction coefficient are presented in Fig. 3(b, c). Due to many weak internal reflections inside the multilayer sample, the measured waveform shows weak oscillations after the main pulse, what explains the characteristic spectral features around 1.0 and 1.4 THz (Fig. 3(b, c)). The absorption peak of polyethylene at about 2.2 THz can also be clearly identified.

Fig. 3
figure 3

The characterization of the 3.25-mm-thick BT10 sample: one reference and two transmitted pulses for 0° and 90° orientations of the electric field (a), refractive index characteristics (b), and extinction coefficient characteristics (c). The characterization of the 180-μm-thick BT10 four-ply sample: zoom on the pulses (d) and refractive index characteristics for 0° and 90° orientations of the electric field (e)

Full size image

One can see (Fig. 3(b)) that the refractive index (n) varies in the range of 1.5300–1.5360. For the sake of simplicity for further analysis, we assumed that n is constant throughout the range of interest, what implies the group index is nominally equal to the phase index: the material is dispersionless. Although this assumption is not absolutely true, it is a sufficiently good approximation for samples having a thickness of about 10 mm which was proven experimentally and described in Section 5.2. The average refractive index (n av ) of the sample, equal to 1.5326, was calculated from the time difference (Δt) between the reference and sample pulse using the following formula: n av  = Δt · c / q, with Δt = 5.767 ps. For the same reason, the absorption of all the samples is considered to be constant and the average extinction coefficient (κ av ) needed for further simulations was determined to be 2 × 10−3.

In order to measure the optical parameters of a single layer, a single 45-μm-thick ply was mechanically detached from the sample. Because of the small thickness of a single layer, the pulses transmitted through the sample for both orientations are almost indistinguishable, showing no detectable time difference. In addition, instabilities of the laser used in the TDS setup did not allow for the determination of a time delay between the reference and probe pulse. Therefore, a sample composed of four single plies was created by folding four parts of the same layer and pressing them together, resulting in a total thickness of about 180 μm. A clear time difference was observed for the propagation time of two orthogonal orientations of the electric field, although the thickness of the sample can be assumed to be constant in both directions. Using the measured time delay (Δt), we determined the refractive indices for parallel and perpendicular orientation of the electric field to be n 90 = 1.510 and n 0 = 1.545, respectively. This indicates that the unidirectional fibers exhibit a birefringence of Δn = n 0–n 90 = 0.035.

Refractive indices calculated for both orientations (Fig. 3e) exhibit big differences over the whole measured range of frequencies in comparison with Fig. 3(b), where the two curves almost overlap with each other. This phenomenon is resulting from the fact that the pulse traveling through successive plies is reflected several times from particular layers giving rise to an additional contribution to the pulse propagating directly through the sample (etalon effect). In our case, layers are thin (~45 μm) and, therefore, such reflections add to the main pulse, which causes a change not related to spectral features of the medium. As a result, both values of the refractive indices have a relatively large decrease around 1 THz. The difference for all other frequencies is quite uniform and is equal to about 0.03–0.04. In Fig. 3e, two additional calculated values of refractive indices (n 0 and n 90) are shown. They correspond to the average values of both indices in the frequency range under investigation and will be used in further calculations.

It must be pointed out that the measurement of thin layers is very demanding also due to non-perfect stability of the TDS system and problems with the precise determination of the sample thickness. Therefore, the obtained value of the birefringence (Δn) must be considered as an approximation. Moreover, since the calculated average refractive index (n av ) can be treated as a mean of about 70 perpendicular plies, the three indices fulfill the following condition: n 0 > n av  > n 90, which is in accordance with our expectations.

The same measurements were repeated for the 2.95-mm-thick multilayer HB50 sample. The complex refractive index (N) is equal to N = 1.521 + i · 0.002. Due to the fact that the plies of HB50 are glued to each other by a resin, it was impossible to separate a single ply. This fact prevented us from determining the birefringence of a single layer. The origin of this birefringence is discussed in Section 6.

4 Theory of Reflection and Simulation Method

This section contains the simplified theory for the determination of thickness from reflection measurements in case of a multilayer structure. Moreover, a description of the TMM used in our simulations is presented together with the fitting procedure to determine the thickness of each layer.

4.1 Fresnel Equations

In general, if a short pulse with an amplitude equal to 1 is incident on a multilayer structure with relatively thick layers, the reflected waveform consists of a series of pulses reflected from the interfaces between layers with different refractive indices. Both amplitudes and phases of the reflected pulses are governed by Fresnel equations. For perpendicular incidence on an interface between two media with refractive indices n i and n j (j = i + 1), the amplitude transmission (t ij ) and reflection (r ij ) coefficients can be written as

$$ {t}_{ij}=\frac{2{n}_i}{n_i+{n}_j}, $$
(1)
$$ {r}_{ij}=\frac{n_i-{n}_j}{n_i+{n}_j}. $$
(2)

If n j  > n i , then r ij  < 0, leading to the 180° phase change of the electric field during reflection. Let us assume a structure with k-alternating layers with a binary profile of the refractive index, namely n 1 = n low = 1.510, n 2 = n high = 1.545, n 1, n 2,…. The structure is surrounded from both sides by air with the refractive index (n 0 = n k + 1 = 1). For simplicity, the extinction coefficient of the sample was assumed to be small and was neglected in the considerations. All individual amplitudes (a q ) in the waveform can be calculated from

$$ {a}_1={r}_{01}, $$
(3)
$$ {a}_q={\left(-1\right)}^2\cdot {r}_{\left(q-1\right)q}\cdot \left({t}_{01}\cdot {t}_{10}\right)\cdot {\left({t}_{12}\cdot {t}_{21}\right)}^{q-1};\kern1em q=2,\kern0.5em 3,\kern0.5em \dots \kern0.5em k. $$
(4)

For a 20-ply structure, we can calculate the particular amplitudes which are equal to a 1 = −0.203, a 2–19 = ±0.009, and a 20 = 0.019.

4.2 Reflection from UHMWPE Composites

Although Fresnel equations are always valid, the above calculations do not take into account the fact that in UHMWPE composites, particular layers are no longer thick compared to the terahertz pulse duration. Figure 4a shows the electric pulses (green) registered using the detector coming from the reflection from the multilayer structure of the pulse generated by a source with the electric field vector (E) perpendicular to the XZ plane. Due to the fact that the diameters of UHMWPE fibers are much smaller than the wavelength of the incident THz radiation, the plies of fibers can be substituted with homogenous plies with the constant refractive index, n low or n high. This binary refractive index modulation is caused by plies with alternating orientation of fibers relative to the E vector.

Fig. 4
figure 4

Reflections of the terahertz pulse from the multilayer structure (a) and the terahertz pulse illuminating this structure (b)

Full size image

All thicknesses of plies are in the range of about 35–100 μm. The terahertz pulse generated in the TDS system, described in Section 5.1, consists of the main peak with a full width at half maximum (FWHM) equal to 0.3 ps and two side lobes at distances ±0.4 ps with amplitudes of about −0.44 and −0.17, respectively (Fig. 4b). However, for a 40-μm-thick layer with n = 1.5, the transit time of a pulse reflected from the backside is equal to 0.4 ps and thus similar to the time difference between the peak and the lobes. This leads to a temporal overlap of pulses reflected from consecutive interfaces. In the simplest case presented in Fig. 4a, we can assume that pulse II, which is reflected from the second interface (II), is partially superimposed with pulses I and III reflected from the first interface and the third one, respectively. This superposition generally depends on the thicknesses and the refractive indices of the layers. The phase flip is observed for pulses I and II (n 0 < n 1 for the first interface and n 1 < n 2 for the second).

4.3 Simulations

Due to the complexity resulting from this partial superposition of pulses, the waveform reflected from the multilayer structure with the refractive index of ith layer (N i  = n i  + iκ (i = 1…k)) was simulated using the TMM [22]. The structure is surrounded from both sides by air with the complex refraction index (N 0 = N k + 1 = 1). In TMM, each layer is modeled with the following two matrices [22]:

  • D ij describing the behavior of the terahertz pulse at the interface between ith and jth (j = i + 1) layers containing the corresponding Fresnel coefficients

$$ {D}_{i,j}=\frac{1}{t_{ij}}\left[\begin{array}{cc}\hfill 1\hfill & \hfill {r}_{ij}\hfill \\ {}\hfill {r}_{ij}\hfill & \hfill 1\hfill \end{array}\right], $$
(5)

  • P i describing the propagation of the pulse through the ith layer

$$ {P}_i\left(\omega \right)=\left[\begin{array}{cc}\hfill \exp \left(\frac{i\omega {N}_i{d}_i}{c}\right)\hfill & \hfill 0\hfill \\ {}\hfill 0\hfill & \hfill \exp \left(-\frac{i\omega {N}_i{d}_i}{c}\right)\hfill \end{array}\right] $$
(6)

where ω is the angular frequency and c is the speed of light in vacuum.

The multilayer structure can be described as the multiplication of all layers’ matrices

$$ {M}_{\mathrm{total}}\left(\omega \right)={\displaystyle {\prod}_{i=0}^k{P}_i\left(\omega \right)\cdot {D}_{i,i+1}}=\left[\begin{array}{cc}\hfill {M}_{11}\left(\omega \right)\hfill & \hfill {M}_{12}\left(\omega \right)\hfill \\ {}\hfill {M}_{21}\left(\omega \right)\hfill & \hfill {M}_{22}\left(\omega \right)\hfill \end{array}\right]. $$
(7)

The transfer function for the reflection geometry (R(ω)) can be calculated from the following equation:

$$ R\left(\omega \right)=\frac{M_{21}\left(\omega \right)}{M_{11}\left(\omega \right)}. $$
(8)

Finally, the simulated signal (E r (t)) reflected from the multilayer sample can be calculated using the inverse numerical Fourier transform (F −1) according to the following formula:

$$ {E}_r(t)={F}^{-1}\left[R\left(\omega \right)\cdot F\left[{E}_0(t)\right]\right] $$
(9)

where F[E 0(t)] is the numerical Fourier transform of the incident THz pulse. The measured reference THz signal (see Figs. 4b and 6b) from the TDS setup was used as E 0(t) in the simulations presented here.

Figure 5 represents the simulation of the waveform reflected from a 20-ply sample as described above. Thicknesses of the layers (d 1–d 20) were chosen randomly in the range from 35 to 100 μm and are presented in Fig. 5 together with the binary refractive index profile (white and yellow colors). Each layer has the extinction coefficient (κ = 2 × 10−3). The x-axis was scaled from the time domain into the distance domain Z (μm), knowing that the total thickness of the sample (d = 1250 μm) corresponds to the time difference (Δt = 12.738 ps) between the front and back surface echoes. The average refractive index can be calculated using the following formula:

Fig. 5
figure 5

The simulation of the waveform reflected from the 20-ply structure (transfer matrix method (TMM))

Full size image
$$ {n}_{av}=\frac{0.5\cdot \varDelta t\cdot c}{d} $$
(10)

which is equal to 1.529 for the simulated structure.

One can notice that the waveform starts with the front surface echo with the negative amplitude (b 1 = −0.201). Next, we observe reflections from individual interfaces with positive and negative amplitudes in the range (b 2–19 = ±0.014) and the high positive pulse with the amplitude (b 20 = 0.018) connected with the back surface of the sample. Comparing a i calculated in Section 4.1 with b i simulated here, we can assess the influence of the side lobes on the individual amplitudes. Amplitudes of front (a 1 vs. b 1) and back (a 20 vs. b 20) surface reflections are very similar, while the amplitudes 2–19 differ significantly as caused by the superposition of the side lobes.

Moreover, due to the two-sided lobe shape of the terahertz pulse, the positions of maxima and minima are shifted with respect to the real positions of the interfaces: significantly at the beginning (18 μm) and less (±3 μm) in the middle of the waveform. Similar result can be obtained if we reverse the refractive index profile (n 1 = n high = 1.545, n 2 = n low = 1.510,…, n 20 = n low).

4.4 Determination of Thicknesses of All Layers

It can be concluded from Fig. 5 that the determination of layers’ thicknesses based on positions of maxima and minima will lead to significant errors. Therefore, the method described in [22] was adopted for precise calculation of these thicknesses. The algorithm requires the exact number of layers in the structure and the refractive index profile of the sample. The first value can be obtained by counting the maxima and minima in the waveform. The average refractive index of the whole sample can be calculated using Eq. 10. It is assumed that the refractive index profile of the structure is binary, which means that the refractive index of the ith layer can be described as follows: n i  = n av  ± 0.5 · Δn. In the case considered here, Δn = n high − n low = 0.035. However, for real samples, Δn depends on characteristic features of the sample (e.g., caused by the manufacturing process) and can be different. Therefore, it was necessary to repeat the procedure a few times for various values of Δn to find the optimal value.

The algorithm seeks for global minima of the function QERR described as

$$ \mathrm{QERR}\left({d}_1\kern0.5em \dots \kern0.5em {d}_k\right)={\displaystyle {\sum}_t\left[{\mathrm{signal}}_{\mathrm{measured}}(t)-{\mathrm{signal}}_{\mathrm{simulated}}\left({d}_1\kern0.5em \dots \kern0.5em {d}_k,t\right)\right].} $$
(11)

QERR(d 1…d k) describes the mean square deviation between the measurement and simulation. Our approach utilizes this formulation of the objective function which is based on the amplitude and the phase of the reflected terahertz subpulse at each boundary. Thicknesses of individual layers are frequently varied until an optimal correlation between the measured and the simulated pulse is achieved. The minimization of the value of QERR corresponds to the increase of the correlation factor and thus leads to the correct thicknesses of the layers. Due to the existence of several local minima, it is difficult for a deterministic optimization algorithm to find the global minimum. Therefore, we utilize a genetic optimization algorithm to estimate the thickness for each layer by minimizing the error function QERR. The applied stochastic algorithm is also called differential evolution [28, 29].

The operation of this evolutionary algorithm is very similar to other stochastic optimization algorithms. Firstly, an initial vector population (parent population) is randomly generated which consists of random thickness values for all layers. We assume a uniform probability distribution of these values over the entire given range

$$ \mathrm{thickness}(i)={\mathrm{thickness}}_{\mathrm{lower}}+\mathrm{rand}\left(0,\kern0.5em 1\right)\cdot \left({\mathrm{thickness}}_{\mathrm{upper}}-{\mathrm{thickness}}_{\mathrm{lower}}\right) $$
(12)

where i is the layer index, rand(0,1) is the random value from the (0,1) range, and thicknesslower and thicknessupper are the lower and upper boundaries of the searched thickness, respectively.

With a random mutation in the parent population based on Eq. (13), a new population is calculated

$$ {\mathrm{thickness}}_{\mathrm{new}}(i)={\mathrm{thickness}}_{\mathrm{old}}\left({r}_1\right)+\mathrm{rand}\left(-1.5,\kern0.22em 1.5\right)\cdot \left({\mathrm{thickness}}_{\mathrm{old}}\left({r}_2\right)-{\mathrm{thickness}}_{\mathrm{old}}\left({r}_3\right)\right) $$
(13)

where r 1–r 3 are random integer and mutually different indices.

In the next step, newly calculated thicknesses are compared to the old thicknesses using QERR(d 1…d k). If the new thicknesses’ values yield a smaller cost function value QERR than the old thicknesses, they are set to the next generation; otherwise, the old values are retained. This iteration is repeated until a given stopping criterion is reached. The last calculated population contains corresponds with a high probability to the global minimum for all layers.

The accuracy of the thickness calculation of layers depends on the accuracy of determination of material parameters (n and κ). In the above simulation, we know these parameters accurately, which leads to the precise reconstruction of the initial thicknesses with an error lower than 1 μm. In Section 5, we also apply this algorithm to the experimental data and we prove its correctness.

5 Experimental Setup and Results

This section contains results of experimental verification of the fitting procedure by means of a TDS setup in reflection configuration.

5.1 Experimental Setup

A THz TDS setup in reflection configuration (precisely described in [2]) was used for the investigation of UHMWPE samples. The TDS setup uses a Ti:sapphire laser, an InAs surface emitter, a GaAs photoconductive detector, and a mechanical chopper for generation and detection of THz pules (Fig. 6a). In the configuration shown in Fig. 6, the THz radiation hits the sample under normal incidence. The focal length of the second parabolic mirror is 105 mm, resulting in a spot size at the surface of the sample of about 2 mm (measured at 1 THz). A typical THz signal reflected from a mirror replacing the sample (reference) is presented in Fig. 6b. A long delay line (speed 0.33 ps/s) combined with a lock-in amplifier (time constant 30 ms, sensitivity 200 mV) provides a time delay of about 140 ps with 13,441 measuring points with spacing 0.0104 ps. Such precision is required for further analysis. The acquisition time of a single 140-ps-long waveform lasts about 7 min. After the main pulse, we observe the small secondary pulse at 15 ps with an amplitude of 2 connected with a back reflection from the inside of the setup. Moreover, the inset in Fig. 6b shows a zoom on the tail of the reference signal (20–125 ps) without any significant features with maximum amplitude equal to 0.6 and standard deviation 0.16. Assuming that the amplitude reflection coefficient (a 1) from the front surface of the sample is equal to 0.203, we can assume that the influence of this tail on the signal reflected from the sample is negligible. A zoom on the main pulse with FWHM = 0.3 ps and side lobes is also presented in Fig. 4b.

Fig. 6
figure 6

The experimental setup (a) and the reference signal (b). The upper and lower insets show the power spectrum and a zoom of the tail of the reference signal, respectively

Full size image

Due to the fact that in our experiment, we analyzed the time-domain data, the performance of the system was determined also in time domain, according to the procedure described in [30]. The signal-to-noise ratio and the dynamic range are equal to 90 and 1500, respectively. Moreover, a power spectrum of the reference pulse plotted in the inset of Fig. 6b shows a bandwidth up to 3 THz and a dynamic range of about 40 dB. Since the setup is situated in a box purged with dry air, the humidity is low (<1 %) and water features are not visible both in the time-domain signal and in the power spectrum.

The considered samples were placed on a rotating platform, and the THz radiation was precisely focused on their front surfaces using a z-stage. The BT10 sample was carefully aligned on the platform with respect to the parallel orientation of the electric field vector in relation to the direction of fibers of the front layer (0°). Under such conditions, the system measured the reflected waveform presented in Figs. 7 and 9. Next, the sample was rotated by 90° (to obtain the perpendicular orientation) and the waveform was also acquired (Fig. 9). HB50 was investigated in the same way, and the result for 90° is presented in Fig. 10. The rotation axis coincides with the measurement point in order to measure the same spot under different orientations.

Fig. 7
figure 7

The waveform reflected from the BT10 sample (orientation 0°). The inset shows a zoom of the waveform with the binary profile of the refractive index

Full size image

The algorithm proposed in Section 4.4 is required to know the exact number of layers (k) in the structure. The condition which guarantees the reliable distinction of all interfaces is derived from [31]: the longitudinal (axial) resolution of the TDS system (g), understood as the ability to resolve two closely spaced reflections associated with adjacent layers, should be lower than the smallest thickness of a layer in the sample. g can be determined using the formula [31]

$$ g=\frac{c\cdot \tau }{2n}. $$
(14)

For the refractive index (n = 1.53) and the temporal width of the THz pulse (τ = 0.3 ps) (Fig. 4b), the longitudinal resolution is estimated to be about 30 μm, which is less than the smallest thickness of the layer both measured using a microscope and calculated using our algorithm. Therefore, the applied TDS system is capable to clearly distinguish all interface-derived maxima and minima in the reflected signal and to count the number of layers in the sample.

5.2 Thicknesses of Layers in BT10 Sample

Figure 7 shows the measured waveform reflected from the 10.40-mm-thick BT10 composite material presented in Fig. 1c. The electric field is parallel to the direction of fibers of the first ply (orientation 0°), which means that the refractive index of this ply is equal to n high. One can easily notice the front and back surface echoes and in between many reflections from all individual layers of the sample. The total propagation time of the pulse through the sample and back is 106.209 ps. The waveform was smoothed by means of a Chebyshev filter. Based on the thickness of the sample, the x-axis was scaled from the time domain into a distance (Z (mm)). The FWHM of the echo reflected from the beginning of the sample is τ 1 = 0.30 ps, and that from the back of the sample is τ 2 = 0.33 ps, respectively, which confirms that the dispersion of the sample having a thickness of around 10 mm is relatively small. This fact allows the assumption in our calculations of n = const.

The inset in Fig. 7 shows a zoom of a part of the waveform with few maxima and minima. We counted k = 189 plies in the waveform. Green dotted lines determine interfaces between the layers. Most of the layers have a calculated thickness of about 40–50 μm, while the third and the tenth ply are about two times thicker than the others. Taking into account the woven nature of the BT10 tape, this 90–100 μm thickness corresponds to two layers having the same orientation of fibers so that the refractive index step between them is missing.

In accordance with Fresnel equations, if the wave propagates from the medium with lower (higher) refractive index into the second medium with higher (lower) refractive index, the signal reflected at the interface has an inverted (non-inverted) phase. On this basis, the binary refractive index profile—the regions with higher and lower refractive indices—was determined (see inset in Fig. 7). The average refractive index of the sample was calculated using Eq. 10 (Δt = 106.209 ps, d = 10.40 mm) and is equal to n av  = 1.5317, which is in good agreement with the transmission calculation from Section 3. The same measurement and analysis were carried out for the 90° orientation, where the electric field is perpendicular to the fibers of the first layer, and the obtained waveform was also used for further analysis.

The TMM-based algorithm described in Section 4.3 requires the knowledge of k, n av , κ av , and Δn for the calculation of layers’ thicknesses. k and n av were determined precisely above, the extinction coefficient (κ av  = 0.002) was adopted from the transmission measurement (Section 3), while Δn was calculated using the trial and error method with a best fit of Δn = 0.05. We assume that the single ply of UHMWPE fibers slightly expanded during the process of separation from the sample which would explain the difference between Δn = 0.035 calculated in Section 3 and the value obtained for the plies in the sample.

Thicknesses of all layers were optimized by means of the TMM algorithm and the function QERR for both orientations of the electric field, and they are presented in Fig. 8(a). The simulation lasted about 5 h using a tabletop PC. We obtained the following results for both orientations 0° and 90°; 138 plies have an average thickness of 43 ± 2 μm, and 51 plies have an average thickness of 88 ± 5 μm. The maximum and minimum thicknesses are equal to 103 and 37 μm, respectively. Taking into account the occurrence of double plies (with thicknesses above 50 μm), we obtained 240 plies with the average thickness of 43 μm.

Fig. 8
figure 8

Thicknesses of layers (a) and differences between thicknesses calculated for both orientations of the electric field (0° and 90°) (b)

Full size image

Moreover, Fig. 8(b) shows differences between the thicknesses calculated for both orientations; 132 layers (70 %) have exactly the same thickness, 53 layers (28 %) differ by ±1 μm, and only 4 layers (2 %) differ by 2 or 3 μm.

In order to compare measurement results with simulations, we also simulated the waveform reflected from the structure with thicknesses calculated above. Figure 9 presents four waveforms: two measured reflected waveforms from the structure for both orientations of the electric field (0° and 90°) and two simulated reflected waveforms from a virtual sample with parameters corresponding to the experimental one also for two orientations of E field. Moreover, the inset illustrates the zoom of the region 4.0–4.6 mm (the same as in the inset of Fig. 7).

Fig. 9
figure 9

Four reflected waveforms—measured and simulated—for both orientations of electric field (0° and 90°). The inset shows the zoom of the waveforms in the region between 4.0 and 4.6 mm

Full size image

Due to the almost perfect longitudinal arrangement of the UHMWPE fibers in the plies, the layers with higher refractive index for 0° orientation of E, have lower refractive index for 90° orientation and vice versa. Therefore, the reflected waveforms for those two orientations have opposite phase values, which is clearly visible in the inset of Fig. 9.

Simulated waveforms agree well with the experimental data when we take into account shape, amplitudes, and positions of echoes with the correlation coefficient equal to 0.985 and 0.979, for 0° and 90°, respectively. Some layers exhibit almost perfect agreement between simulation and experiment. The typical difference between peaks of measured and simulated pulses is about 3–5 μm, while the maximum difference is about 10 μm. Small differences in the amplitudes of reflected pulses can be related to the non-perfect microstructure state of the interface between fibers of adjacent layers.

It is obvious that all thicknesses determined for two orientations independently are almost identical, and therefore, further investigations can be limited to one orientation only. The calculation of all thicknesses was carried out, assuming a binary profile of the refractive index of the sample, which is not completely true in the real samples due to technological limitations. It is possible that individual plies have refractive indices slightly different from the assumed binary ratio n low/n high. This would explain the maximum discrepancies of about ±10 μm between the measured and the simulated waveform. However, these variations of individual refractive indices must be small, because thicknesses calculated for both orientations are almost equal. Until now, the software developed here is not capable of optimization of both the thickness and the refractive index for each layer in a 200-ply material.

5.3 Thicknesses of Layers in HB50 Sample

Measurements of the BT10 sample (Section 5.2) showed that calculations for both orientations of the electric field give very similar results. Therefore, in the case of the HB50 sample, we focus only on the 90° orientation. Figure 10 represents the waveform reflected from the 9.45-mm-thick HB50 composite material as shown in Fig. 1d. The waveform was smoothed by means of the Chebyshev filter, and the x-axis was scaled from the time domain into a distance (Z (mm)) based on the thickness of the sample. The inset (on the right) shows a zoom of part of the waveform with maxima and minima coming from reflections from consecutive layers. The electric field is perpendicular to the fibers of the first ply (orientation 90°), which means that the refractive index of this ply is equal to n low. The propagation time of the pulse through the sample and back is 95.85 ps, corresponding to n av  = 1.5171. This is a slightly lower value of the refractive index than in the case of BT10; nevertheless, this result is consistent with the results from transmission measurements (see Section 3).

Fig. 10
figure 10

Waveforms reflected from the HB50 sample—measured and simulated. The inset on the right shows the zoom of the marked part of waveforms. The inset on the left illustrates calculated thicknesses of particular layers

Full size image

Based on the algorithm presented above, we counted k = 154 plies in the waveform and determined the refractive index profile. Most of the layers have a thickness of about 60 μm, and double layers are not observed in this case. Using the trial and error method, we determined the parameters Δn = 0.04 and κ av  = 0.002 which give the best fit between the measured and the simulated waveform.

Finally, it should be pointed out that the average thickness of the layers is equal to 61 ± 7 μm, while the maximum and minimum thicknesses are equal to 82 and 45 μm, respectively. All calculated values of thicknesses are depicted in Fig. 10 (inset on the left). Almost 90 % of calculated thickness values are within a range of 50–70 μm, what is in good agreement with results presented in the literature [24].

Moreover, the simulation of the reflected signal was carried out for earlier calculated thicknesses of layers and for determined material parameters (n av , Δn, κ av , and k), and its result is also shown in Fig. 10. The similarity between the measured and the simulated waveforms is very good with the correlation coefficient equal to 0.969.

For the HB50 sample, differences between amplitude values and the distribution of successive maxima and minima are larger than for the BT10 sample. The typical difference between peaks of measured and simulated pulses is about 5–7 μm, with a maximum deviation of 17 μm. Such difference may result from the microstructural inhomogeneity of interfaces between successive layers, small bubbles of air existing between layers, or real thickness fluctuations of layers.

6 Discussion of the Birefringence

To describe the interaction between the polarized THz radiation and the UHMWPE fiber, we utilized the Polder and van Santen approach [32], adopted to the THz range by Jordens et al. [33]. This theory is valid, if the diameter of the fibers is much smaller than the wavelength of the radiation, which is satisfied in the considered case, because the diameter of the fibers is about 15 μm and the radiation having the frequency of 1 THz corresponds to the wavelength equal to 300 μm. According to this theory, the effective permittivity of a dielectric mixture composed of fibers embedded in a homogenous matrix for both orientations of the electric field (ε 0 and ε 90) can be calculated using the following formulas (Fig. 11):

Fig. 11
figure 11

Sketch of fibers in a single ply with the relevant permittivities

Full size image
$$ \delta =\frac{\varepsilon_0-{\varepsilon}_{\mathrm{Matrix}}}{\varepsilon_{\mathrm{Fiber}}-{\varepsilon}_{\mathrm{Matrix}}}, $$
(15)
$$ \delta =\frac{\left(1-\frac{\varepsilon_{\mathrm{Matrix}}}{\varepsilon_{90}}\right)\left({\varepsilon}_{90}+{\varepsilon}_{\mathrm{Matrix}}\right)}{2\left({\varepsilon}_{\mathrm{Fiber}}-{\varepsilon}_{\mathrm{Matrix}}\right)}, $$
(16)

where δ is the volumetric content of the fiber inside the matrix. For the resin-free BT10 ply, we assumed that δ = 0.95 and the effective permittivity of fibers ε fiber = (n av )2 = (1.5326)2. The matrix is with ε matrix = 1. With these values, the following permittivities were obtained: ε 0 = (1.508)2 and ε 90 = (1.499)2, resulting in a birefringence of Δn = 0.009. Although it was assumed that the composite consists of 5 % of air, which is an already rather high value for these samples, the calculated birefringence (0.009) is too low, and it must be concluded that this theory is insufficient to explain the experimentally observed value of Δn = 0.05. In the case of an HB50 ply, the refractive index of the matrix (plastic rubber) is not known for the authors, which prevents us from calculating the effective permittivity.

The main drawback of this theory is the fact that it assumes that fibers have uniform form and it does not take into account their microstructural nature. The internal structure of the fibers is connected with polyethylene molecules, and it should be underlined that the THz range of radiation corresponds to the molecular interaction of hydrogen bonding, van der Waals interactions, and lattice interactions [34, 35]. Another factor that should be taken into account is the influence of the strain and deformations that arise during the whole composite production process: fiber drawing, tape forming, and finally, hot pressing [36].

To sum up, the origin of the birefringence is not clear to the authors. In contrary to a bulk polyethylene, which is normally homogeneous, a single ply of UHMWPE fibers exhibits birefringence. This phenomenon can be the result of both the longitudinal shape of almost perfectly arranged fibers (Fig. 1b), the microscopic structure of many long PE chains with a simple repeated atomic structure (Fig. 1a), and the influence of the strain originating from the production process.

7 Summary, Conclusions, and Outlook

We reported on terahertz TDS investigations of two kinds of UHMWPE composite materials manufactured from BT10 and HB50 tapes from Dyneema®. At the beginning, the complex refractive indices were determined in transmission configuration of the TDS setup for both multilayered structures (BT10 and HB50) and for a structure composed of four plies (of BT10) having the same (parallel) orientation of PE fibers. For the HB50 sample, a birefringence was observed, which results in the dependence of the value of the refractive index on the orientation of the electric field vector in relation to the direction of fibers. This fact became the basis of the structure analysis in the experiment in the reflection configuration of TDS setup.

The mechanism of the formation of reflections inside the sample from subsequent fiber layers with alternate orientation (one perpendicular to the next one) was described based on the periodic modulation of the refractive index and was simulated using a transfer matrix method.

Using Fresnel’s equations, it was shown that reflections from particular interfaces between layers results in a temporal overlap of consecutive pulses. A simple analysis using only the peak position of the reflected pulses leads to an inaccurate determination of the thickness of particular layers.

Therefore, in order to determine thicknesses of layers of both UHMWPE composites, a time-domain fitting procedure was used, looking for a maximum correlation between the measured pulse and the simulated one. The reconstructed waveforms based on the optimized thicknesses show excellent agreement with the measurements with correlation coefficients above 0.97, resulting in typical differences in thickness between measurements and simulations in the order of 3–5 μm for BT10 and 5–7 μm for HB50, respectively.

It should be pointed out that fluctuations of the thickness of each layer may also be the reason for small differences between the amplitudes of the measured and the simulated pulse and the distribution of their peak values (minima and maxima). Moreover, these differences can also be caused by the microstructural inhomogeneity of the interfaces or by a small amount of air bubbles existing along the interface.

Another reason for these differences may be the fact that for the simulations, we assume a binary profile of the refractive index inside the investigated sample. This criterion may not be always satisfied due to small structural inhomogeneities. Our analysis shows that the transfer matrix method (TMM) is a very precise tool even if we use only average values of the refractive index (n av ) and extinction coefficient (κ av ) for the calculations.

As a result, thicknesses of multilayer structures were determined by means of the reflection TDS technology with high precision. In the case of a BT10 sample, we performed independent calculations for both orientations of the electric field in relation to the arbitrary chosen edge of the sample. Here, 240 layers were calculated having the average thickness of about 43 ± 2 μm (in this case, double layers existed—two consecutive layers with the same fiber orientation). Almost identical results were obtained for both orientations of the electric field. Ninety-eight percent of thickness values differ by less than 1 μm, which confirms the correctness of the applied method. For the HB50 sample, we counted 165 layers with the average thickness of about 61 ± 7 μm, which corresponds to data known in the literature for such samples. Here, double layers were not present due to the different structures of the composite material (non-woven).

The main restriction of the presented method is connected with the assumption of a binary profile of the complex refractive index of the samples, which is not completely true in real samples due to technological limitations. Moreover, the complex refractive index of the investigated samples was approximated to be spectrally constant, which seems to be a sufficiently good approximation for the 10-mm-thick samples. The accuracy of the calculation of layers’ thickness by means of the proposed TMM-based method can be improved by optimization of both thickness and the spectrally dependent complex refractive index independently for each layer. This computationally demanding task for a 200-ply material could be accelerated by using a graphical processor unit to analyze a multilayer system within a few minutes. With the hybrid approach, which combines the advantages of the brute force and the stochastic optimization approaches, we can significantly increase the efficiency of our approach. Finally, the proposed algorithm does not consider secondary reflections from the interfaces, so-called etalon effect. Although this phenomenon seems to be negligible for the considered structures due to the relatively small value of the reflection coefficients, this fact should be also taken into account in rigorous computations.

To sum up, the presented THz TDS-based method allows for non-ionizing and contactless examination of multilayer UHMWPE composite samples with accurate determination of number and individual thicknesses of layers. The growing number of threats of the contemporary world, including wars and terrorism hazards, leads to the need of finding the materials which can improve the people protection against different weapons. Therefore, UHMWPE composites, due to their high ballistic protection, seem to be appropriate response to these requirements. Moreover, new nondestructive evaluation methods, which can be applied for quality control of different products manufactured from UHMWPE (like personal armors, bulletproof vests, helmets), are really needed.