Design of Functionally Graded Alloys for Locks Highly Resistant to Ultrasonic Detector Attacks

1. Introduction

Mechanical locks have been continuously improved [1,2,3], while mechanical assemblies of electric locks still have problems to be solved. This is why mechanical locks are still an important link in the chain of security. Along with the development of mechanical locks, different attack techniques are also being improved [4,5,6,7,8]. The ultrasonic method is one of them.

The motivation for this paper is to increase mechanical locks’ resistance to active and passive ultrasonic attacks. This is not to be done by significant changes in their mechanical assemblies. Shapes and dimensions of all other mechanical components (including the driver pins) will be left the same. Only the chemical composition of the key pins will be changed. This way, a high-security lock’s protections against standard mechanical attacks [9,10,11] (e.g, bumping, picking or raking using Kopsch’s tools, mechanical decoding with "Sputnik" tool, drilling, etc.) in the form of pick-proof pins, trapping calottes (Fangkalotten) [2,3], double pins, pins in multiple planes, cascaded locks (Kaskadensystem) [3,9], anti-drilling rods, etc, will be kept, and extra protections will be added.

Active ultrasonic detectors [12] read the lock pins’ lengths by measuring the ultrasound echo delay time on the tip of a pin. Passive ultrasonic detectors [13], which are cheaper, register vibrations on a pin’s lateral surface (Figure 1) using a laser beam perpendicular to the pin. The laser microphone records pin vibrations at the pin’s impact on the plug/rotor, after the pin’s spring is compressed and abruptly released. Measured natural frequencies are then correlated to the pins’ lengths.

For both types of ultrasonic attacks to succeed, the Young’s modulus E and density $\rho$ of key pins must be known to the attacker, and pins of different lengths must have different natural frequencies. Both key pins and driver pins on all known high-security pin-tumbler locks are always made of homogeneous metal alloys, most often brass or steel. Driver pins are all the same length, while key pins from one manufacturer are usually made in 10 standard lengths. A code defining cuts in a key to match, e.g., a particular 5-pin lock is, for example, 2-4-6-5-8. This is why it’s easy to correlate a key pin’s natural frequency to its length, which makes every such lock vulnerable to ultrasonic attacks.

A lock resistant to ultrasonic attacks should have key pins with equal natural frequencies regardless of their lengths. Furthermore, the passive ultrasonic attack won’t succeed if the laser probe is pointed at an area on a pin close to a standing wave node. This is why the key pins should be made of alloys with variable $E(x)$ and density $\rho (x)$ along the pin’s longitudinal axis x.

Natural frequency f of longitudinal vibrations of a thin rod [14] of length L made of a homogeneous material is

$$f={\displaystyle \frac{1}{2L}}\sqrt{{\displaystyle \frac{E}{\rho }},}$$

(1)

where E is Young’s modulus and $\rho$ is density. This means that an equal frequency for all 10 pins of different lengths can be achieved with different homogeneous alloys. This is true, but this way, the positions of standing wave nodes $x_0$ (Figure 2) are fixed to a point $x_0<L/2$ and can’t be moved. If it’s possible to move $x_0$ to a point close to pin’s tapered tip, like in Figure 2, which is a spot where a laser probe will be pointing, it will be more difficult for an attacker to record a useful signal because of the small amplitude of vibrations.

Even for a homogeneous material, equation (1) is valid only if the rod can be considered "thin"(i.e, if the ratio $L/d>1.7$, for error < 0.1%), and if cross-section $A(x)$ is constant. Gradient $dA(x)/dx$ increases the speed of sound $v(x)$, thus increasing the frequency f and moving the node point $x_0$ to the left from $L/2$. Its preferred position is actually to the right, as shown in Figure 2.

In this paper, FGMs are proposed to be used for making key pins, instead of usual homogeneous alloys. FGMs have been researched since mid-1980s [15]. Many new FGMs along with new methods of manufacturing have been invented since then [16,17,18,19]. Their first significant use case was for space shuttle wings, to combine tensile strength of steel with superior high-temperature properties of ceramic materials. Simply laminating steel wings with a layer of ceramics makes the wing structure too weak, which is why a gradual transition from pure ceramics on the surface to metal in the bulk was devised to make the wing structure mechanically and thermally stronger.

FGM will improve a lock’s resistance to ultrasonic attacks, because by varying $E(x)$ and $\rho (x)$ along the x-axis, both pin’s natural frequency and position of standing-wave node can be manipulated. This way, it will be difficult to correlate a pin’s natural frequency to its length. The method to produce such FGM pins, i.e. 3D printing from metal powders is already well established - a so-called direct energy deposition (DED) method. It uses several different metal powders, mixed at different ratios on every point along the x-axis, to achieve the desired $E(x)$ and $\rho (x)$ values. The volume ratios of metal powders are calculated from $E(x)$ and $\rho (x)$ according to Voigt’s model.

Mechanical vibrations on FGM structures [20,21,22] have been modeled by different application-specific models, which is a further contribution of this paper, along with a method to optimize $E(x)$ and $\rho (x)$ profiles. This means that mathematical models of pins’ vibrations on a mechanical transmission line with variable $E(x)$, $\rho (x)$ and $A(x)$ have to be devised first; this is done in Section 2. Simplified models, appropriate for use in an optimization procedure (i.e, fast enough to calculate 500-1000 successive simulations in a reasonable time) have been derived, in Section 3. Initial calculations for pins with two distinct sections ($E(x)$ and $\rho (x)$ as step functions) as a starting point for the optimization procedure (to reach set referent values both for f and $x_0$) are done in Section 4, along with the transition to gradual $E(x)$ and $\rho (x)$ derived from Butterworth polynomials. An optimization procedure to calculate the optimal $E(x)$ and $\rho (x)$ profiles is defined and performed in Section 5. Final adjustments to ameliorate the resistance against ultrasonic attacks (by compensating for short stick effects and introducing some random variables) are performed in Section 6.

2. Mathematical Models of Pins Vibrations

Key pins excited by an axial impulse (from a piezo transducer probe of an active ultrasonic decoder, or an impact with a plug/rotor in front of a laser probe of a passive ultrasonic decoder) behave as a mechanical transmission line with free ends. Both ends of the key pins are not affixed, and the loading spring is very soft compared to a metal pin; hence, the termination impedance can be considered practically infinite at both ends. According to [23], the error in natural frequency of one-dimensional thin-rod axial model is less than -2.7% for the smallest key pins (i.e, for pins with ratio $L/d=1.4$, which is the lowest value among all standard pin-tumbler locks). The difference in natural frequencies between the two shortest standard pin lengths (pin-1 and pin-2, pin numbering goes up to pin-10 for the longest pin, for a lock with 10 pin lengths) is usually higher than 10%.

Partial differential equation (PDE) for one-dimensional mechanical line (shown in Figure 3 - K is stiffness in [N/m] of a line segment of width $\Delta x$, A is cross-section area, F is force on segment $\Delta x$, $\Delta m$ is its mass, D is a damping coefficient in [Ns/m], and $\xi$ is a dimensionless damping factor) [21,24] with variable parameters (2) can’t be solved analytically in the general case.

$${\displaystyle \frac{\partial }{\partial x}}\left(E\left(x\right)A\left(x\right){\displaystyle \frac{\partial u\left(x,t\right)}{\partial x}}\right)=\rho \left(x\right)A\left(x\right){\displaystyle \frac{{\partial }^{2}u\left(x,t\right)}{\partial t^2}}$$

(2)

Damping factor $\xi$ for metallic alloys is very low (in the order of ${10}^{-3}$), hence its term in (2) is neglected. It has no influence on the natural frequency and position of standing wave nodes. Equation (2) is then written in a more appropriate form:

$${\displaystyle \frac{\partial E_A\left(x\right)}{\partial x}}{\displaystyle \frac{\partial u(x,t)}{\partial x}}+E_A(x){\displaystyle \frac{{\partial }^{2}u(x,t)}{\partial x^2}}={\rho }_A\left(x\right){\displaystyle \frac{{\partial }^{2}u(x,t)}{\partial t^2}}$$

$$E_A\left(x\right)=E\left(x\right)A(x),{\rho }_A\left(x\right)=\rho \left(x\right)A(x)$$

(3)

2.1. FEM Model

FEM models are accurate, but slow. A FEM model is coded in MATLAB to be used as a reference for the evaluation of faster and simplified models designed in Section 3. Sampling time is $\Delta t=T_d$, number of finite elements is $N=500$, pin length is L, length of one FEM segment is $\Delta x=L/N=l$. Signal $u(x,t)$ discretized in space and time is $u(j,k)$. Using standard Euler’s approximations for derivatives, the following equations (4)-(6) derived from (3) will define the FEM model.

$$u\left(j,k+1\right)=w^2{v}^2\left(j\right)u\left(j+1,k\right)+\left[w^2{v}^2\left(j\right)·\right.$$

$$·\left(a\left(j\right)-2\right)+2\left]·u\left(j,k\right)+w^2{v}^2\left(j\right)\left[1-a\left(j\right)\right]·\right.$$

$$·u\left(j-1,k\right)-u(j,k-1)$$

(4)

$$w={\displaystyle \frac{T_d}{l}},a\left(j\right)=1-{\displaystyle \frac{E\left(j-1\right)A(j-1)}{E\left(j\right)A(j)}}$$

(5)

$$v(j)=\sqrt{{\displaystyle \frac{E(j)}{\rho (j)}}},x(j)={\displaystyle \frac{jL}{N}}=jl$$

(6)

The boundary condition for the left-hand ($j=1$) free end excited by impact force $F_1$ is described by (7), and analogously for the right-hand free end ($j=N$).

$$u\left(1,k+1\right)=\left(F_1\left(k\right)+\left(u\left(2,k\right)-u(1,k)\right){\displaystyle \frac{E_1{A}_1}{l}}\right){\displaystyle \frac{T_d^2}{{\rho }_1{A}_{1}l}}+$$

$$+2u\left(1,k\right)-u\left(1,k-1\right)$$

(7)

A condition (8) for a minimum sampling time $T_d$ is required for a stable FEM simulation.

$$T_d<\begin{array}{c}min\\ j\end{array}\sqrt{{\displaystyle \frac{\rho (j)}{E(j)}}}l$$

(8)

2.2. Some Exact Solutions of the PDE

Variables x and t can always be separated, because standing waves will be formed on a line with two free ends. Furthermore, the amplitude $U(x)$ is important for optimization, not the instantaneous value of the longitudinal vibrations $u(x,t)=U(x)exp(i\omega t)$. Boundary conditions will thus always be $dU(0)/dx=0$ and $dU(L)/dx=0$. For the simplest case, i.e, a thin pin with constant E, $\rho$, and A, the natural frequency is (1).

2.2.1. Constant $E_0$ and ${\rho }_0$, Exponential $A(x)$

For a pin with cross-section described with an exponential function $A(x)=A_{0}exp(-hx)$, and constant E and $\rho$, the exact natural frequency can be calculated by analytical solution of (3):

$$f={\displaystyle \frac{1}{2L}}\sqrt{{\displaystyle \frac{E_0}{{\rho }_0}}}\sqrt{1+{\displaystyle \frac{h^2{L}^2}{4{\pi }^2}}}$$

(9)

This natural frequency is always higher than for the constant cross-section (1). This indicates that the gradient $dA(x)/dx$ increases the propagation speed $v(x)$ in the same material.

2.2.2. Constant $A_0$, Exponential $E(x)$ and $\rho (x)$

For the Young modulus and density described as $E(x)=E_{0}exp(ax)$, $\rho (x)={\rho }_{0}exp(cx)$, the propagation speed and natural frequency can also be calculated by analytical solution of equation (3):

$$v(x)=\sqrt{{\displaystyle \frac{E_0}{{\rho }_0}}}·exp\left({\displaystyle \frac{a-c}{2}}x\right)$$

(10)

$$f={\displaystyle \frac{c-a}{4}}\sqrt{{\displaystyle \frac{E_0}{{\rho }_0}}}·{\displaystyle \frac{1}{exp\left(\frac{c-a}{2}L\right)-1}}$$

(11)

These equations can also be analytically calculated by integration of time ${\displaystyle \frac{T}{2}}={\int }_0^L{\displaystyle \frac{dx}{v(x)}}$, which will be a basis for devising fast and simple models of FGM key pins vibrations in next section.

3. Fast Mathematical Models for Optimization of FGM Key Pins Profiles

Starting from PDE (3) and the standing wave equation $u(x,t)=U(x)exp(i\omega t)$, the ordinary differential equation (ODE) (12) is derived:

$${\displaystyle \frac{d^{2}U\left(x\right)}{{dx}^2}}=-{\displaystyle \frac{d{E}_A\left(x\right)}{dx}}{\displaystyle \frac{1}{E_A\left(x\right)}}{\displaystyle \frac{dU\left(x\right)}{dx}}-{\omega }^2{\displaystyle \frac{{\rho }_A\left(x\right)}{E_A\left(x\right)}}U\left(x\right)$$

(12)

This is not analytically solvable for every $E_A(x)$ and ${\rho }_A(x)$ function. However, for a known $E_A(x)$ and ${\rho }_A(x)$, the solution (f and $x_0$) can be found with sufficient precision with less than 10 successive simulations, using a Simulink model in Figure 4. Variable x is treated as time. Integrators work with 4^th order Runge-Kutta with constant time step (actually the FEM element length $l=L/N$). Initial value on integrator $U(0)$ is set to any max value (1 $\mu m$ in this case), and on the other integrator to $U^{\prime }(0)=0$ (this is one boundary condition). Hence, if $\omega =2\pi f$ is set to the pin’s correct natural frequency, the second boundary condition $U^{\prime }(L)=0$ will be met at the end of the simulation (Figure 5b) for f=476 kHz. The standing wave node position $x_0=2.92$ mm is then read directly from $U(x)$ in Figure 5a.

One Simulink simulation takes 0.3 seconds, which means 3 seconds will be needed to establish $x_0$, f and $U(x)$ profile. This is 10 times faster than one FEM model simulation in Matlab (cca. 30 seconds). This model is more useful for checking f and $x_0$ calculated by some other method, because it is then 100 times faster than FEM, and its precision (i.e, f, $x_0$ and $U(x)$) is the same as FEM’s. However, a faster model will be needed for FGM pins’ profile optimization procedure, since 500-1000 successive simulations are needed in simplex optimization of one pin to reach an optimum profile.

3.1. Speed Profile $v(x)$ Approximation

If $A(x)=A_0$ is constant, the ultrasound speed at point x can be calculated as $v\left(x\right)=\sqrt{E(x)/\rho (x)}$. Natural frequency $f=1/T$, and standing wave node position $x_0$, similarly to (11), can then be calculated by numerical integration of time (13):

$${\displaystyle \frac{T}{2}}={\int }_0^L{\displaystyle \frac{dx}{v(x)}},{\displaystyle \frac{T}{4}}={\int }_0^{x_0}{\displaystyle \frac{dx}{v(x)}}={\int }_{x_0}^L{\displaystyle \frac{dx}{v\left(x\right)}}$$

(13)

The influence of $dA(x)/dx$ along with $E(x)$ and $\rho (x)$ on speed $v(x)$ is difficult to calculate analytically. The main idea is to approximate the increase of speed $v\left(x\right)=\sqrt{E(x)/\rho (x)}$ at the pin’s tip (where $dA(x)/dx\ne 0$) accurately enough, then (13) can be used as a fast model for optimization procedure. This way, f and $x_0$ can be quickly calculated using (13), and then model (12) can be used to calculate $U(x)$ and to confirm the values of f and $x_0$. The simplest approach is to use an exponential function (14), where $v_a\left(x\right)=\sqrt{E(x)/\rho (x)}$, as a speed correction factor at the tip ($x>L_1$).

$$v(x)=v_a(x)exp(b(x-L_1))$$

(14)

Even if $v(x)$ thus doesn’t match exactly at every point on the x-axis, but f and $x_0$ match the ODE model (12), this will be precise enough for the pin optimization procedure. For a pin of certain dimensions (length L, and a tip of constant dimensions for all 10 standard lengths of a particular lock manufacturer), it will be shown that coefficient $b=b(L)$ can be calculated to achieve sufficient precision. Standard key pins of the locksmithing company Schlage will be further used as a case study.

3.2. Calculation of b Coefficient

Two straight lines and one circular arc approximate the tip of a standard Schlage key pin, Figure 6. The coefficient b was calculated for 10 homogeneous pins (made of brass, $v(x)=v_0=3870$ m/s). Then it will be shown that this exponential correction of speed is also accurate for arbitrary variable $E(x)$ and $\rho (x)$. Pulse propagation time from $x=0$ to $x=L$ is:

$$t\left(L\right)={\int }_0^L{\displaystyle \frac{dx}{v(x)}}={\int }_0^{L_1}{\displaystyle \frac{dx}{v_0}}+{\int }_{L_1}^L{\displaystyle \frac{dx}{v_{0}exp(b(x-L_1))}}=$$

$$={\displaystyle \frac{L_1}{v_0}}+{\displaystyle \frac{1}{b{v}_0}}\left(1-exp(b(L-L_1))\right)=T/2,$$

(15)

where natural frequency is $f=1/T$, propagation time from one end to standing wave node is $t(x_0)=t(L)/2$, and $x_0=v_0·t(x_0)$ for $x_0<L_1$.

The following substitutions will be introduced (16). The length of a tip $d_{sch}=1.4$ mm is constant for all 10 pins; $L_C$ is dependent on L.

$$d_{sch}=L-L_1,L_C={\displaystyle \frac{1}{b}}\left(1-exp(b(L-L_1))\right)$$

(16)

From (15) and (16), the implicit equation (18) is derived.

$$f={\displaystyle \frac{v_0}{2}}·{\displaystyle \frac{1}{L_1+L_C}},x_0={\displaystyle \frac{1}{2}}\left(L_1+L_C\right)$$

(17)

$$b={\displaystyle \frac{1}{v_0/(2f)-L_1}}\left(1-exp(b{d}_{sch})\right)$$

(18)

Equation (18) will be used to iteratively calculate the value of b according to the following procedure:

1.) For a standard pin length L, and $A(x)=d^2(x)\pi /4$ with $d(x)$ according to Figure 6, run a Simulink simulation according to (12) and Figure 4. This yields f and $x_0$.

2.) Using f, iteratively calculate $b=b(L)$ from (18)

3.) Calculate $L_C=L_C(b)$ from (16)

4.) Take the next standard value for L and go to step 1.

This way, for 10 standard Schlage pin lengths,

L[mm]=[4.19 4.57 4.95 5.33 5.72 6.10 6.48 6.86 7.24 7.62]

The following values of b and $L_C$ are obtained:

b[$m^{-1}$]=[516 520 52 4 528 532 537 541 545 549 553]

$L_C$[mm]=[1.00 1.00 0.99 0.99 0.98 0.98 0.98 0.97 0.97 0.97]

3.3. Comparison of Two ODE Models

Approximation of initial ODE model (12) with variable $A(x)$, $E(x)$ and $\rho (x)$ (called "A" model) is compared with the ODE model with constant cross-section $A(x)=A_0$, and variable $E(x)$ and $\rho (x)$ (called "v" model). In the second model, the $E(x)$ will be multiplied by the factor $exp(b(x-L_1))$ for $x>L_1$, and $\rho (x)$ will be divided by the same factor, which will effectively multiply the speed $v\left(x\right)=\sqrt{E(x)/\rho (x)}$ by that factor along the pin’s tip.

3.3.1. Natural Frequencies of Brass Schlage Pins

Pins are homogeneous, made of brass ($E_0$=120 GPa, ${\rho }_0$=8000 kg/m³). The results of simulations of ODE models "A" and "v" in Figure 7 show that there are very small differences between models (maximum frequency error is 0.88 kHz, and maximum standing wave node position error is 0.015 mm) for homogeneous pins. For comparison, ODE model simulations of cylindrical pins (such as driver pins, without tapered tips) are also performed, showing that the frequency error increases up to 10% (for pin-1) if the tapered tips are neglected.

3.3.2. Comparison of "A" and "v" Models of FGM Pins

ODE models simulations are performed for an FGM pin of length $L=$6.86 mm, as in the previous subsection. The natural frequency of the pin is $f=$364 kHz, the node position is $x_0$=3.98 mm, and the accuracy remains the same. Variable $E(x)$ and $\rho (x)$ in Figure 8 are shown multiplied by factor $A(x)/A_0$ for the "A" model, and with the exponential factor for "v" model. Both ODE models also show equally high accuracy compared with the FEM model.

On the other hand, as it’s shown in Figure 8a, the ratio of the maximum and minimum value of $E_A(x)$ of the "A" model can be much higher than the ratio of respective values of the "v" model. This can increase numeric errors at ODE model integrators on the "A" model. This is why, for certain $E(x)$ profiles, the accuracy of the "A" model decreases compared to FEM and "v" models.

3.4. Testing of Fast Time Integrator Model

Since the accuracy of the ODE "v" model is confirmed, and hence the accuracy of speed correction with the exponential factor, the fast time integrator model (13) can now be tested for accuracy, since it uses that same approximation. For a pin with $E_v(x)$ and ${\rho }_v(x)$ profiles as in Figure 8, the fast time integrator model results are f=403 kHz, and $x_0$=4.07 mm. Therefore, the frequency error is >10%.

On the other hand, for a pin with $E_v(x)$ and ${\rho }_v(x)$ profiles as in Figure 9 (transition points of both E and $\rho$ are now close to $x_0$), accuracy is much better. ODE "v" model gives f=476.7 kHz, $x_0$=4.43 mm, while the fast time integrator gives f=476.2 kHz, $x_0$=4.43 mm.

The accuracy of this model increases as transition points approach the standing wave node. However, it will be shown (in Section 5) that it’s possible to set the optimization objective and penalty functions so that node position $x_0$ falls close to both transition points at the end of the optimization procedure.

This model is still accurate enough for the optimization procedure and the most appropriate for it. It can complete 1000 simulations in MATLAB in 3 seconds. When the optimum point is reached, the final result can be double-checked and confirmed by one simulation using the ODE "v" model.

4. Initial Calculations for FGM Pins

There is more than one way to achieve the optimization objective, but the simplex method’s convergence is not always guaranteed. This is why its initial vectors must be properly set. The following procedure will be described, which achieves reliable convergence.

4.1. Setting the Initial Parameters for Optimization

Setting the initial parameters for optimization is performed according to the following algorithm:

1.) A desired distance of the standing wave node from the pin cone $d_0=L_1-x_0$ is first set. Let $d_0=$1 mm.

2.) Calculate (according to (16)-(18)), for a homogeneous pin (e.g, made of brass, like in subsection 3.2), length $L_{ref}=L_C+2d_0+d_{sch}$=4.395 mm (since $L_C$ depends on L, a few iterations may be needed).

3.) For this $L_{ref}$, also according to (18), calculate $f=v_2/2/(L_{ref}-d_{sch}+L_C)$ = 476 kHz, for $v_2$ defined by initial homogeneuous material (brass), i.e, 3800 m/s.

This is not the length of a standard Schlage pin (it falls between pin-1 and pin-2). It will be used to set a reference for natural frequency (476 kHz) and node position ($d_0$=1 mm left from the tapered tip). The idea is to keep the right-hand side the same (made of brass), and change the composition of the homogeneous alloy on the left-hand side (e.g, like Figure 10), so its speed of sound $v_1$ satisfies the equation:

$${\displaystyle \frac{v_1}{L-d_0-d_{sch}}}={\displaystyle \frac{v_2}{L_{ref}-d_0-d_{sch}}}$$

(19)

4.) This way, the wave propagation time on the left-hand side remains the same, keeping the frequency f and node position $d_0=L_1-x_0$ unchanged, for every standard pin of length L:

L[mm]=

[ 4.19 4.57 4.95 5.33 5.72 6.10 6.48 6.86 7.24 7.62 ]

$v_1$[km/s]=

[ 3.41 4.13 4.86 5.58 6.32 7.05 7.77 8.50 9.22 9.95 ]

FEM simulations confirm this approach; natural frequency and node positions remain the same for all 10 standard lengths.

4.2. Using Butterworth Polynomials to Improve the Structure

A structure like that in Figure 10, with an abrupt transition, may be difficult to implement in practice. The reason is a problem with achieving sufficient adhesion along the separation plane ($x=L_2$), along with a possible problem of fitting their different crystal grids. Furthermore, when attacked with an active ultrasonic detector (ultrasonic pulse injected at the pin’s tip $x=L$), there will be a reflection from the separation plane ($x=L_2$), and then another from the flat, left-hand end ($x=0$). Their delay times will be equal for any pin length, but the ratio of amplitudes and phases of these two reflected pulses will be different for different materials on left-hand side. This is another reason why the transition must be made gradually (a gradual transition won’t create reflections), so it can be realized using methods like 3D printing/sintering from metal powders.

Butterworth polynomials magnitudes (20) will be used to implement gradual transitions of $E(x)$ (in Figure 11), and analogously, $\rho (x)$ as well. The transitions are thus located in the bulk, not close to surface like in usual planar FGM structures described in [15].

$$E\left(x\right)=E_1{\displaystyle \frac{1+{\left(\frac{x}{x_{1E}}\right)}^{2n_E}}{1+{\left(\frac{x}{x_{2E}}\right)}^{2n_E}}},\rho \left(x\right)={\rho }_1{\displaystyle \frac{1+{\left(\frac{x}{x_{1\rho }}\right)}^{2n_{\rho }}}{1+{\left(\frac{x}{x_{2\rho }}\right)}^{2n_{\rho }}}}$$

(20)

$$E_2=E_1{\left({\displaystyle \frac{x_{2E}}{x_{1E}}}\right)}^{2n_E},{\rho }_2={\rho }_1{\left({\displaystyle \frac{x_{2\rho }}{x_{1\rho }}}\right)}^{2n_{\rho }}$$

(21)

Equations (20) are magnitudes of Butterworth filter transfer functions, with ±3 dB points at $x_1$ and $x_2$. Their exponents $n_E$ and $n_{\rho }$, however, don’t have to be integers. To estimate the width of the transition region $\Delta x_{trans}$ for profile $E(x)$ (to set maximum acceptable values of exponent $n_E$), a tangent in the inflection point $x_{inf}$ is used, in a span between $E_1$ and $E_2$ (22), and analogously for profile $\rho (x)$.

$$E^{\prime }(x_{inf})={\displaystyle \frac{E_1}{x_{inf}}}\left({\left({\displaystyle \frac{x_{2E}}{x_{1E}}}\right)}^{2n_E}-1\right){\displaystyle \frac{4n_E^2-1}{8n_E}}$$

(22)

$$x_{inf}=x_{2E}{\left({\displaystyle \frac{2n_E-1}{2n_E+1}}\right)}^{{\displaystyle \frac{1}{2n_E}}},\Delta x_{trans}={\displaystyle \frac{E_2-E_1}{E^{\prime }(x_{inf})}}$$

(23)

From (22) and (23) it can be estimated that $\Delta x_{trans}$ will be $\Delta x_{trans}$ >0.4 mm if $n_E,n_{\rho }<20$, hence it will be a bound for one of the penalty functions for simplex optimization. Standard metal powders with a particle diameter of 10 $\mu m$ are easily available to achieve a sufficiently gradual transition. Furthermore, the lower bound will be $n_E,n_{\rho }>2$, because if the exponent is too low, then the equations (21) for $E_2$ and ${\rho }_2$ are not valid.

4.3. Setting the Initial Vector for Simplex Optimization

The optimization vector h will be defined as follows:

$$h=[\Delta x_E\Delta x_{\rho }{x}_{mid}{\rho }_1{\rho }_2{E}_1{E}_2]$$

(24)

$x_{2E}=x_{mid}-\Delta x_E,x_{1E}=x_{mid}+\Delta x_E$

$x_{2\rho }=x_{mid}-\Delta x_{\rho },x_{1\rho }=x_{mid}+\Delta x_{\rho }$

The initial values of the h vector, for every pin length L and corresponding speed $v_1$ calculated by (19) are (25). Value $g_{inf}=|E^{\prime }(x_{inf})|/E_1$ is an acceptable gradient at the inflection point.

$v_2=3800$ m/s, ${\rho }_2=8000$ kg/m³, $E_2={\rho }_2{v}_2^2$, $n=10$

${\rho }_1={\rho }_2{v}_2/v_1$, $E_1=E_2{v}_1/v_2$, $L_2=L-d_0-d_{sch}$, $g_{inf}=1500m^{-1}$

$$\Delta x_E=\left|L_2(1-k)\right|,k=exp\left({\displaystyle \frac{ln(1+2g_{inf}{L}_2/n)}{4n}}\right),\Delta x_{\rho }=\Delta x_E,x_{mid}=L_2.$$

(25)

4.4. Approximation of E(x) and E’(x) for ODE Model

ODE "v" model requires gradient $E^{\prime }(x)$ for simulation, to confirm the optimization result for each pin length L. To simplify $E(x)$, and especially $E^{\prime }(x)$ for Simulink simulation, (20) will be approximated by $E_p(x)=E_1{\sum }_{i=0}^{10}{p}_i{x}^i$, fitted by the least-squares method.

5. Optimization of E(x) and ρ(x) Profiles

The model to be used in the simplex optimization procedure was defined in Section 3, followed by the optimization vector and its initial values in Section 4. The objective function and penalty functions are next to be defined. Since the goal is to get f and $x_0$ as close as possible to references $f_{ref}$=476 kHz and $x_{0ref}=L-d_0-d_{sch}$ for each pin length L, the main objective function is (26), with weight factors $k_x={10}^{10}$ and $k_f={10}^{-6}$.

$$F_{obj}=k_x{(x_0-x_{0ref})}^2+k_f{(f-f_{ref})}^2$$

(26)

First penalty function (27) is calculated from exponents $n_E$ and $n_{\rho }$. These exponents are calculated in each optimization step from values in the optimization vector h, according to (21). Penalty increases when $n_E,n_{\rho }<2$ or when $n_E,n_{\rho }>20$, as explained in subsection 4.2. Factors $k_{E1}$ or $k_{\rho 1}$ have a value of 100 if $n_E$<2 or $n_{\rho }$<2, otherwise zero. Similarly, $k_{E2}$ or $k_{\rho 2}$ have a value of 10 if $n_E$>20 or $n_{\rho }$>20, otherwise zero.

$$F_{pn}=k_{E1}{(n_E-2)}^2+k_{E2}{(n_E-20)}^2+$$

$$+k_{\rho 1}{(n_{\rho }-2)}^2+k_{\rho 2}{(n_{\rho }-20)}^2$$

(27)

The following penalty function (28) is introduced because the fast time integrator model (13) used in the optimization procedure loses precision if $x_0$ is located outside of ($x_{1E}$,$x_{2E}$) or ($x_{1\rho }$,$x_{2\rho }$) range.

$$F_{pm}=k_m{\left({\displaystyle \frac{x_{1E}+x_{2E}}{2}}-x_0\right)}^2+k_m{\left({\displaystyle \frac{x_{1\rho }+x_{2\rho }}{2}}-x_0\right)}^2$$

(28)

The optimization algorithm will hence be steered to move both arithmetic mean values towards $x_0$. The weight factor is set to $k_m=4·{10}^8$.

Boundary values $E_1$, $E_2$, ${\rho }_1$ and ${\rho }_2$ must stay within some reasonable range. This is why brass (the most common material for pins in pin-tumbler locks) was chosen as the initial point, and then the ultrasound speeds $v_1$ for the left-hand side of each pin were calculated in Section 4 according to (19) and Figure 10. Although a much wider span of speeds of sounds is possible, ranging from 1200 m/s (common lead) to 20000 m/s (diamond), the goal is to stay within a span of common, easily obtainable metal alloys, with a sufficient minimum shear strength (in the order of a common brass, i.e, 200 GPa).

This is why a blue pentagon (Al-bronze-brass-Fe-Be) in Figure 12 is chosen as the permitted area. The following penalty function (29) will indicate if points $M_1$=($E_1$, ${\rho }_1$) or $M_2$=($E_2$, ${\rho }_2$) fall outside of the permitted area.

$$F_{pE\rho }=k_E\rho ·{\displaystyle \frac{{\sum }_{j=1}^2{\sum }_{i=1}^5{P}_{triangle,i}(E_j,{\rho }_j)-{2P}_{pent}}{{2P}_{pent}}}$$

(29)

The weight factor is $k_E\rho =2·{10}^3$. The constant $P_{pent}$ is the area of the pentagon. A sum of areas of five triangles is calculated for boundary points $M_1$ and $M_2$. For point $M_1$ its five triangles are: ($M_1$-bronze-Al), ($M_1$-Al-Be), ($M_1$-Be-Fe), ($M_1$-Fe-brass), ($M_1$-brass-bronze), and likewise for point $M_2$. If both boundary points are located within the pentagon, the sum of the areas of the five triangles equals the pentagon area, and the penalty function is zero. If a boundary point is located outside the pentagon, the sum of the areas of the five triangles is greater than the pentagon area, and the penalty function value (29) increases accordingly.

The cumulative penalty function is then calculated as a sum $F_p=F_{pn}+F_{pm}+F_{pE\rho }$.

According to [15] and [25], both Young’s modulus and density of a metal alloy (and a metal-ceramic composite as well) can be calculated (30) based on a simple rule-of-mixture with volume fractions ${\varphi }_i$ as weight factors, assuming additivity of volume (Voigt model). The Voigt model is reasonably accurate for metal alloys and metal-ceramic composites, but not for porous materials [26], for which a more precise model (Mori-Tanaka [15]) must be used. This means that any (E, $\rho$) combination within any triangle in Figure 12. can be achieved by a certain mixture of three or more metals.

$$E(x)=\sum _{i=1}^N{\varphi }_i(x)·E_i,\rho (x)=\sum _{i=1}^N{\varphi }_i(x)·{\rho }_i$$

(30)

Boundaries of the pentagon in Figure 12 are set considering the following criteria:

1.) A sufficient shear strength for soft metals with E<70 GPa is difficult to achieve, hence this is one of the limits. Duralumin with 95% aluminium will have a much higher shear strength than pure aluminium, but practically the same Young’s modulus and density, according to (30).

2.) The speed of sound in beryllium is very high (v=12500 m/s), and it’s not too expensive compared to the rest of the high-security lock, since it will be used to make only key pins. If it’s acceptable to use metal-ceramic composites (which are also possible to implement as FGMs, according to [15]), possible replacements for beryllium are alumina (Al₂O₃, E=400 GPa, $\rho$=4000 kg/m³, v=10000 m/s), silicon carbide (SiC, E=400 GPa, $\rho$=3150 kg/m³, v=11250 m/s) or boron carbide (B₄C, E=450 GPa, $\rho$=2500 kg/m³, v=13400 m/s). According to [27] and [28], alloys of beryllium, aluminium and iron are feasible in any ratio.

3.) The purpose of this paper is to prove that it’s possible to calculate profiles of $E(x)$ and $\rho (x)$ for standard key pin lengths to achieve desired f and $x_0$, within a reasonable permitted range of $E(x)$ and $\rho (x)$, using reasonably priced and obtainable materials, considering the advances in production of various FGMs. The exact composition of pins’ alloys (to achieve a minimum shear strength and other important properties) and the proper method of production (e.g, 3D printing and sintering using metal powders) are out of the scope of this paper.

5.1. Results of Optimization Procedure Using Simplex Method

The optimization procedure for all ten standard pin lengths, using the fast time integrator model (13), the objective function defined by (26), and penalty functions defined by (27)-(29), with initial values for each pin length set by (25), was performed. The resulting $E(x)$ and $\rho (x)$ profiles are shown in Figure 12 and Figure 13. The time needed for optimization procedure in MATLAB, for all ten pins, is cca. 30 seconds.

The initial span of speeds of sound (19) calculated for 2-section pins (Figure 10), as the initial values for simplex optimization, was 3410-9950 m/s. The speed span after the completed optimization (Figure 12) is 3350-10450 m/s, for pins with $E(x)$ and $\rho (x)$ modeled using Butterworth polynomials (20).

Red ’X’ symbols on Figure 12 indicate right-hand, tapered pin tips, i.e. ($E_2$, ${\rho }_2$) points, with properties close to brass. Black ’X’ symbols indicate left-hand ends, or ($E_1$, ${\rho }_1$) points. FGM pins profiles between red ’X’ and black ’X’ are described by (20) and (21).

FEM simulation tests for all ten optimized pins indicate a maximum ±0.8% deviation in natural frequency from reference $f_{ref}$=476 kHz, and a maximum ±0.05 mm deviation in positions of standing wave nodes $x_0$.

6. Final Adjustments

The initial 1-dimensional FEM model, along with derived ODE and fast time integrator models, is valid for slender sticks, i.e. for $L>2d$. It can be used even in the range $d<L<2d$, but the natural frequency calculated from simulation results will thus be higher (because of short stick effects) than the actual natural frequency, according to [23]. The pin’s natural frequency, for pin-1, pin-2 and pin-3 must therefore be corrected by increasing $E_1$ (and consequently $E_2$ and every value of $E(x)$ calculated by simplex optimization) and decreasing ${\rho }_1$ (and therefore ${\rho }_2$ and $\rho (x)$ profile) by a certain factor $k_c$

$$E_c(x)=E(x)·k_c,{\rho }_c(x)=\rho (x)/k_c$$

(31)

The speed of ultrasound is thus increased by the $k_c$ factor, effectively increasing the actual pin’s natural frequency by the $k_c$ factor. The actual pins will thus all have the same natural frequency $f_{ref}$=476 kHz, when short stick effects are compensated. According to [23], the $k_c$ factor value for Schlage pin-1 is therefore $k_c$=1.023, for pin-2 is $k_c$=1.011 and for pin-3 is $k_c$=1.002. It can be neglected for all longer pins, from pin-4 to pin-10.

6.1. Introducing Random Fluctuations

Producing pins with exact $E(x)$ and $\rho (x)$ profiles as calculated by the simplex method (as on Figure 12) will be an additional challenge. A certain method of production (e.g, metal powder 3D printing) may produce pins with slightly different natural frequencies f, which could still be correlated to their lengths.

This can be preventively remedied by additionally correcting every pin’s $E_c(x)$ and ${\rho }_c(x)$ profiles with a certain random factor $k_{rnd}$ calculated before production of each pin.

$$E_{rnd}(x)=E_c(x)·k_{rnd},{\rho }_{rnd}(x)={\rho }_c(x)/k_{rnd}$$

(32)

If a random number generator sets $k_{rnd}$ at 0.95<$k_{rnd}$<1.05 range, thus introducing a random $\pm 5\%$ frequency fluctuation, this may be (depending on production process tolerances) enough to fully obfuscate any possible remaining correlation between any pin’s length and its natural frequency.

7. Conclusion

Mechanical locks are an important link in the chain of security. They won’t be fully supplanted by electric locks. Simple methods for quickly defeating fingerprint ID or RFID readers, or simply bypassing electric locks have already been found. Picking mechanical locks requires more training and manual dexterity. Mechanical locks are hence still being researched and improved, along with methods of attack.

Instead of using usual homogeneous metal alloys, it is proposed to make a lock with key pins made of FGM alloys, to make the lock more resistant to ultrasonic attacks. The key pins made of FGM alloys will thus have equal natural frequencies. Furthermore, the standing wave nodes on pins will be placed in desired, the most convenient positions.

The proposed method for calculating the profiles of FGM pins using the simplex optimization procedure with the fast time integrator model (to quickly calculate f and $x_0$) derived from the ODE "v" model needs less than one minute to complete the calculations for 10 standard key pin lengths of any manufacturer. Although the original 1-dimensional PDE and FEM models can have up to cca. +3% error for the shortest pins, it can be corrected by factors for certain $L/d$ ratios. The additional randomness (for production procedures like 3D printing) can easily be introduced to increase security.

The security of a pin-tumbler lock is increased, making it less susceptible both to active and passive ultrasonic attacks at an additional price of more expensive FGM key pins, leaving all other mechanical lock parts unchanged.

Possible types of locks for further research are wafer locks (like Japanese MIWA) or disc-tumbler locks (like DOM Diamant). Unlike pins in pin-tumbler locks, wafers and discs are 2-dimensional elements. The following challenge for research is slider locks (for example, EVVA 4KS), since sliders are 3-dimensional elements. Being 2-D and 3-D structures, mathematical models for wafers, discs and sliders will be increasingly more complicated. Methods for precise production (like 3D printing) of pins, wafers, discs or sliders made of FGM alloys will also have to be further researched and improved in order to meet the new demands.