M Jimenez-Diaz
email: m.jimenez.diaz@tue.nl

J Van Dijk

M Van Der Mullen

J J A M Van Der Mullen

Effect of remote field electromagnetic boundary conditions on

published or not. The documents may come

Introduction

Microwave Induced Plasmas [START_REF] Moisan | Microwave excited plasmas volume 4 of Plasma Technology[END_REF] (MIPs) are interesting subjects for academic studies and offer many possibilities for industrial applications. The latter can be found in the fields of materials processing [START_REF] Boulos | Thermal plasmas: fundamentals and applications[END_REF][START_REF] Paquin | [END_REF][4], gas volume cleaning [5] and for the construction of sources of ions [6], radicals [START_REF] Montaser | Inductively coupled plasmas in analytical atomic spectrometry[END_REF] and photons [START_REF] Wharmby | [END_REF].

Depending on the type of microwave coupling, MIPs can be divided into plasmas inside resonant cavities, surface-wave sustained plasmas and free expanding plasma torches. This paper is devoted to the last category for which we mention as examples the microwave plasma torch (MPT) [9] and the axial injection torch (TIA) [10]. In both cases the plasma is created in the zone where the transition takes place from the electromagnetic (EM) structure of a coaxial waveguide to that of an expansion zone. The expansion can take place in open air or in a vessel. The plasma gas, in most cases argon, is generally provided via the hollow inner conductor of the coaxial waveguide.

The working power of the TIA can reach 2 kW whereas a microwave plasma torch (MPT) is limited to around 200 W. This limitation is due to the use of coaxial cables instead of rectangular waveguides [9] as the connection to transfer the EM energy from the power supply to the torch.

In the past several experiments have been carried out on both MPT and TIA plasmas. Applying absolute intensity measurements [11,[START_REF] Iordanova | Poly-diagnostic validation of spectroscopic methods: in-depth monitoring of microwave induced plasmas[END_REF] of atomic lines and the continuum it was possible to determine the electron density n e and temperature T e ; these values have been compared with the results of the Stark intersection method and Thomson scattering [START_REF] Jonkers | [END_REF][14][15][16][17]. Typical reported values for the TIA are n e = 3 • 10 21 m -3 and T e = 22000 K. For the MPT the typical values are n e = 3 • 10 20 m -3 and T e = 18000 K (see [START_REF] Timmermans | Atomic and molecular excitation processes in microwave induced plasmas a spectroscopic study[END_REF]).

These properties of the electron gas were found to depend on experimental settings like the power level, the flow and the chemical composition of the gas flow (e.g. Ar versus He [START_REF] Van Der Mullen | [END_REF]20]).

Especially the chemical composition was found to have an important impact and in line with this observation, it was understood that the internal plasma properties can also be changed by the entrainment of air. Oxygen and nitrogen molecules sucked in by the argon flow will enhance electron loss-mechanisms and will thus affect the properties of the electrons gas and the shape of the plasma flame. To study the effect of the environment, metallic plasma vessels were constructed in [11,[START_REF] Timmermans | Atomic and molecular excitation processes in microwave induced plasmas a spectroscopic study[END_REF]21,22], with the purpose to change the gas-environment of the plasma (e.g. by replacing air by argon). It was found that chemical features of the environment indeed have an impact on the plasma. As explained in [START_REF] Jonkers | [END_REF], when the plasma is created in a vessel filled with Ar, we get higher n e -and lower T e -values than when it is operated in air . Typical changes are an increase of n e from 6 • 10 20 m -3 in air to 1.2 • 10 21 m -3 in argon and a change in T e from 19000 K in air to 16000 K in argon. This was associated with an increase of the plasma radius from around 1.28 mm (air) to 1.79 mm (argon).

However the question rises in how far placing the plasma into a metallic gas-purification vessel will change the electromagnetic aspects of the plasma. A metal vessel placed around the plasma can cause back-reflection of the electromagnetic waves so that the electric field at the plasma location may change as well. More generally the question can be formulated as in how far the electromagnetic features of the environment will change the properties of microwave induced plasma. This question also addresses the impact of metallic grids that are sometimes placed around MIPs to avoid EM interferences with devices in the surroundings. It is not impossible that these grids will change the EM structure at the plasma location as well.

This study on the impact of the environment on the operation of MIP torches is based on numerical calculations. The modular toolkit PLASIMO [23,24] was employed to calculate the fields for three different type of configurations; open, semi-open and closed systems; they are shortly denoted by, O, S or C configurations. The role of the plasma is played by a cylindrical bar or cone of given permittivity and conductivity.

Realizing that the plasma, formed at a tip of a nozzle has a radial size of less than 1 mm, while the vessels used in practice have radii of about 6 cm, the aim of this paper can be formulated as the question in to what extend the remote ‡ field EM boundary conditions will effect the operation ‡ The distance of the EM boundaries to the plasma is far compared to the plasma radius but not large in comparison and the inner part of a MIP.

In the past several numerical studies were published on one of the configurations O, S or C. In [25] a numerical EM study of the MPT is given. It deals with an open configuration. In [26] the coaxial microwave torch (CMT), i.e. a semi-open configuration (metal vessel working as a circular waveguide) was studied. After giving the EM module, it is also coupled with a plasma fluid model where flow, temperature and species equations in local thermal equilibrium are included. In [27] a closed configuration is used. It consists of a circular waveguide with a metal wall at the top (shortcut).

In all three studies cited above the EM field is found by solving a time-harmonic wave equation for the azimuthal magnetic component (i.e. an equation similar to the Helmholtz equation). In our approach we apply a full vector formulation instead, meaning that apart from the azimuthal magnetic component we also solve the radial and axial electric field components in a coupled manner. With this formulation and the use of the staggered Yee grid [28], the boundary conditions can be fulfilled in a more natural way, which, as we will see, speeds up the convergence process. Together with the application of a direct matrix-vector solver this leads to an enhancement of the convergence speed of more than a factor hundred as compared to the work of [27]. Another difference with the previous studies is that our study is not confined to just one configuration but that instead all three types of configurations, O, S and C are modelled and compared with each other. In this way a better understanding of the EM interplay of the plasma and the environment can be obtained. It was found that the O configuration is the most robust. This configuration is free from resonances so that a feed-back of the plasma on resonances can be avoided. Nevertheless we will see that the plasma still plays a role in the relative importance of the absorption, transmission and absorption of EM energy. The plasma has an antenna function in converting the wave from a TEM mode of the coaxial waveguide into that of radiation into open air. In order to get insight in this conversion the role of various antenna elements is investigated by changing the plasma density and the geometry of the transition region (gap).

This paper is organized as follows. Section 2 contains the definition of the configurations, the plasma and geometry of the setup. The physical and numerical model are explained in section 3. The results and discussion given in section 4. We end with the conclusions in section 5.

Configuration, Plasma and Geometry

The calculations were performed on a model structure that simplifies the basic EM structure of the TIA and MPT or more generally MIP torches. As stated above, both the TIA and MPT plasmas can be seen as extensions of the inner conductor of a coaxial waveguide, meaning that in both cases the plasma forms part of the EM transition from the coaxial waveguide to the expansion region. In the model we can take for this expansion zone a metallic free environment or a circular wave-guide. For the latter we can take various structures that can be divided into closed and semi-open circular waveguides. This leads to the following configurations:

C (Closed)

The plasma expands into a closed circular waveguide, i.e. a circular wave guide with a shortcut at the end (see [START_REF] Timmermans | Atomic and molecular excitation processes in microwave induced plasmas a spectroscopic study[END_REF]27]).

S (Semi-Open)

The plasma expands into a open circular waveguide (see [26]).

O (Open)

The plasma expands into open air (see [25]).

In selecting an appropriate model-plasma we are guided by the experimental values given in [START_REF] Jonkers | Excitation and Transport in small scale Plasmas[END_REF] and the observation given in [START_REF] Timmermans | Atomic and molecular excitation processes in microwave induced plasmas a spectroscopic study[END_REF][25][26][27]. These are depicted in figure 1 where we see that, in practice, the plasma has two zones ( [START_REF] Timmermans | Atomic and molecular excitation processes in microwave induced plasmas a spectroscopic study[END_REF]): an ionizing and a recombining zone. In constructing the model we confine ourselves to the active (ionizing) part where most of the electrons are found and we will assume that this model-plasma has the simple structure of a cylindrical bar, with a constant relative permittivity r and conductivity σ (see section 3.1). The values of r and σ are obtained from a simple model for the electron momentum balance [27,[START_REF] Hagelaar | [END_REF] that requires input with the wavelength; We prefer use remote field rather than far-field to avoid confusion with the terms used in antenna theory where far-field is used for large distances compared to the wavelength. values for the electron density n e , the elastic collision frequency ν eh and the electromagnetic field frequency ω (see [START_REF] Van Den Donker | Modelling microwave plasmas for deposition purposes : exploring the freedom in space and chemistry[END_REF]). In the current study, we follow [START_REF] Álvarez | [END_REF] in taking ν eh = 100ω, with ω = 2πf while f = 2.45 GHz. The electron density is set to n e = 1 • 10 21 m -3 , except in section 4.2, where n e is varied between 10 20 m -3 and 10 22 m -3 in a parameter study on the effect of n e on the power response of the system.

Geometry

Figure 2(a) shows the geometry of the setup which exhibits a 2D axial-symmetric structure. Figure 2(b) depicts the structure of the corresponding computational grid. The coaxial waveguide has an inner radius R 2 = 5.5 mm, an outer radius R 3 = 14.5 mm and ends at Z 1 = 100 mm.

The nozzle starts at Z 1 , transforms into the tip at Z 2 = 106 mm and ends at Z 3 = 118 mm. The zone into which the plasma expands (in open air or into a circular waveguide) has a radius of R 4 and ends at Z 4 , both of which are parameter in this study. These parameters values are given in table 1. As stated above, the plasma is located at the transition of the coaxial waveguide to the expansion region. As shown in figure 2(c), the plasma region is defined by the following (r,z)values (R 1 , Z 3 ), (0, Z 3 ), (R 1 , Z p ) and (0, Z p ). Al the three configurations O, S, and C are described by the same geometry. The difference lies, as we will see, in the boundary conditions for the side and top of the circular waveguide.

Model

This section describes the electromagnetic model. The electromagnetic equations are given in section 3.1, the grid is given in section 3.2 while their respective boundary conditions are introduced in section 3.3. We will comment on the numerical implementation in section 3.4. A summary of the model is given in Table 2 and a comparison with other models in Table 3.

Main characteristics of the current electromagnetic model Electromagnetic Problem Deterministic (not the modes but the fields are determined) Harmonic fields (exp(iωt)) 2D (symmetry axis along z, [25], [26], and [27]. Here, TIA stands for Torch à Injection Axiale (Axial Injection Torch), MPT for Microwave Plasma Torch and CMT for Coaxial Microwave Torch.

∂/∂φ = 0) Full vector Maxwell Equations Transversal Magnetic ( E r , E z , H φ ) Projection Finite Differences Grid staggered (Yee grid) Ohm's law J = σ E Solver Direct Solver (superLU)

Maxwell Equations

For a dissipative and non-magnetic linear medium and assuming harmonic time dependency exp(iωt), the Maxwell equations are

∇ × E = -iωµ 0 H, (1)
∇ × H = J + iω r 0 E, (2)
where H, E and J are respectively the complex magnetic and electric fields and the current density, while r is the relative permittivity of the neutral and non-dispersive medium, 0 the vacuum permittivity and µ 0 the vacuum permeability. The angular frequency is ω = 2πf with f the EM wave frequency.

The axial symmetry is so that the electromagnetic quantities are azimuth-independent (∂/∂φ = 0) and the problem can be described in two dimensions. This axial-symmetry is also refered to as the m = 0 mode.

The dimensions of the coaxial waveguide, R 2 and R 3 , are selected so that only a transverse electromagnetic (TEM) mode can propagate. Due to the axial symmetry, m = 0, and because the expansion zone is being excited by a TEM mode where H z = 0, only TM 0 modes are supported in the expansion zone [33][34][START_REF] Harrington | Time-harmonic electromagnetic[END_REF]. This means that only the E r , E z and H φ field components of a TM 0 mode are needed in the model.

We assume that the current density is determined by the electrons solely, and that it linearly depends on the electric field as expressed by Ohms' law, J = σ E, where σ is the complex conductivity. The last, together with the harmonic assumptions, leads to the definition of the complex permittivity [START_REF] Van Den Donker | Modelling microwave plasmas for deposition purposes : exploring the freedom in space and chemistry[END_REF][START_REF] Bowers | High frequency electron resonances and surface waves in unmagnetized bounded plasmas[END_REF]

ˆ r = r -iσ( 0 ω) -1 . (3)
In the case of plasma r = 1 and

σ = 0 ω 2 p ν eh + iω (4)
with the plasma frequency ω p = e 2 n e /m e 0 and ν eh the collision frequency.

The complex conductivity [START_REF] Van Den Donker | Modelling microwave plasmas for deposition purposes : exploring the freedom in space and chemistry[END_REF][START_REF] Bowers | High frequency electron resonances and surface waves in unmagnetized bounded plasmas[END_REF] in the plasma is obtained from the current density J = -en e v e and the cold collisional electron approximation of the electron momentum balance equation for harmonic fields which reads

iω v e + ν eh v e = - e m e E, thus v e = -e E m e (ν eh + iω) . (5)
In this manner we show that J = σ E and obtain σ.

With the above assumptions, the Maxwell equations included in the model get the form

∂E z ∂r - ∂E r ∂z = k 0 Hφ , (6)
∂ Hφ ∂z = k 0 ˆ r E r , (7) 1 r
∂r

Hφ ∂r = -k 0 ˆ r E z (8)
where k 0 = ω √ µ 0 0 is the vacuum wavenumber, Hφ = iZ 0 H φ the normalized azimuthal magnetic field component, and Z 0 = µ 0 / 0 the impedance of free space. The time-averaged power density absorbed in the domain, Q (W/m -3 ), is computed as

Q = 1 2 Re J * E = 1 2 Re (σ) | E| 2 . ( 9
)
So that the total power absorbed in volume V equals

P abs = V QdV = 1 2 V Re (σ) | E| 2 dV. ( 10
)
The complex Poynting vector is obtained in this case as

S = 1 2 E × H * = 1 2 êz êr êφ E z E r 0 0 0 H * φ = 1 2 E r H * φ êz - 1 2 E z H * φ êr , (11)
and the time-average Poynting vector is ¯ S = Re( S).

The conservation of electromagnetic energy is given by the Poynting theorem

∇ • ¯ S + Q = 0, ( 13
)
which in integral form reads

B ¯ S • n db + P abs = 0. ( 14
)
This merely states that the net § power of the EM radiation entering the volume V enclosed by surface B (with outward normal vector n) equals the power absorbed in V.

Grid

Equations ( 6), ( 7) and ( 8), and the corresponding boundary conditions (see Section 3.3), are projected onto a staggered grid using finite differences following the method of Yee [28]. Figure 3(a). shows a cell for the Yee grid, where H φ is defined at the nodal points and E r and E z at the edges of the grid cell. The location of the discrete fields components makes it possible to satisfy the continuity of the field at the edges. The computational grid is composed of interior and boundary cells. The electromagnetic fields are only solved in the interior cells. Each nodal point of an interior cell has, depending on the type of material, certain values of the permittivity r and conductivity σ. In this study the materials used are air with r = 1 and σ = 0, and plasma for which r = 1 and σ is given by equation ( 4). The values at possible interfaces between two interior cells are determined by interpolation. In the boundary cells, the electromagnetic field components are not solved. These cells are only used to set the boundary conditions for the adjacent interior cell. Boundary conditions are applied to the electric field components, E r , E z as they are the only field component defined at the interface of boundary cells and interior cell. They are described in section 3.3. The computational grid is stretched, meaning that the density of grid points is a function of radial position (see [START_REF] Mobley | [END_REF], [23]). The density of points close to the axis is six times higher than at radial edge of expansion zone. .

Boundary Conditions

In general three different types of boundary conditions (BC) can be defined for a scalar field F at a surface B (i) Dirichlet:

F ( x) = p( x), (ii) Neumann: n • ∇F ( x) = q( x), (iii) Mixed: n • [∇F ( x) + h( x)F ( x)] = w( x)
. § The net power entering a region refers to the power associated to the ingoing minus that of the outgoing waves.

Here x is the position vector of a point on B, n the normal on B while p( x), q( x), h( x), and w( x) are explicitly known functions of x . If p( x), q( x) or w( x) = 0 for all x ∈ B, the BC is named homogeneous, otherwise the boundary is called inhomogeneous. As stated above we only need to consider the BCs of the field components E r and E z . So they can play the role of F given above. Below we will sketch the various types of boundaries of the configuration under study and the corresponding conditions imposed to E r and/or E z .

Homogeneous

• Symmetry-BC (r = 0). The axial symmetry requires at the axis a homogeneous Neumann BC i.e. type (ii) with F = E z , applied at positions x ∈ {Z p ≤ z ≤ Z 4 , r = 0} and q( x) = 0; i.e.

∂E z ∂r x = 0. ( 15
)
• Perfect Electrical Conductors (PEC)-BC. For metallic boundaries we assume that the E field component parallel to the surface is zero which demands for a homogeneous Dirichlet boundary conditions i.e. type (i) with q( x) = 0. For example at the inner conductor of the coaxial waveguide the condition reads

E z | x = 0 (16)
∂E r ∂z + iβE r = 0. ( 17
)
Here β is the wavenumber of the wave. With this axial PBC, we assume that the EM energy leaves the domain in axial direction with a wave in the T M 01 circular waveguide mode (see [START_REF] Pozar | Microwave engineering[END_REF]) with wavenumber β = k 2 -(p 01 /R) 2 , in which R is the radius of the circular waveguide, k = ω √ µ 0 0 r the wavenumber of the medium, while p 01 = 2.405 is the first zero of the Bessel function J 0 of first kind [START_REF] Abramowitz | Mathematical handbook of functions[END_REF]. Note that this boundary condition assumes a plane wave propagating in the positive z-direction E r ∝ exp(-iβz). This boundary condition is described in detail in [38,[START_REF] Mur | [END_REF]42]. • PBC at side-boundary for the radial propagation. A homogeneous mixed boundary (type (iii))

is used at the side-boundary

x ∈ {Z 1 ≤ z ≤ Z 4 , r = R 4 , }, for case O, with F = E z , q( x) = k H {2} 1 (kr) H {2} 0 (kr)
, and w( x) = 0. Thus

∂E z ∂r + k H {2} 1 (kr) H {2} 0 (kr) E z = 0 ( 18
)
with H a Hankel function [START_REF] Abramowitz | Mathematical handbook of functions[END_REF]. This BC, described in detail in [42], assumes an outward propagating radial wave with E z ∝ H {2} 1 (kr). • Excitation-BC. This inhomogeneous mixed BC, used at the bottom boundary x ∈ {z = 0, R 2 ≤ r ≤ R 3 }, has two components, one for the incident and one for the reflected wave. We

This BC is often referred to as the absorbing BC (see [38]), however since we use absorption in another meaning and since this BC deals with the propagation of the wave outside the domain we prefer to call this the Propagation BC assume that incident power driven into the plasma has a TEM mode so that the corresponding E-field component reads

E r | + (r, z) = C TEM r exp(-iβz) ( 1 9 )
with

C TEM = Z 0 |P i |/π/ ln(R 3 /R 2 )
, and |P i |, the (imposed) incident power. Due to the discontinuity present in the domain the incident wave will generate a wave that is reflected back to the entrance. For the corresponding component E - r we apply a PBC such that

∂E - r ∂z -ikE - r = 0. ( 20
)
In this case the wavenumber for the TEM mode coincides with the medium wavenumber, k. Note that this BC imposes at the entrance a plane wave propagating in the negative zdirection i.e. E - r ∝ exp(+ikz). This is the opposite to the top boundary where exp(-iβz) is used instead. As the BC should be applied to the total E r field, we substitute the relation 20), which, using [START_REF] Van Der Mullen | [END_REF], gives

E - r = E r -E + r into equation (
∂E r ∂z -ikE r = -2ik C T EM r exp(-ikz). ( 21
)
As we can see the excitation-BC is of type (iii) with F = E r , q( x) = ik, and w( x) = 2ikC T EM exp(-ikz)/r.

To summarize, we give an overview of the BC of the three types of configurations C, S, and O (figures 4(a), 4(b) and 4(c), respectively) The configurations have in common the excitation-BC at the entrance of the coaxial waveguide, the PEC-BC at the inner and outer metallic walls of the coaxial waveguide and the symmetry-BC at the central axis of the expansion region. The differences are (see figure 4):

C. Closed the PEC-BC are applied at the side-and top-boundaries, S. Semi-open the PEC-BC is applied at the radial side while an PBC acts at the top, O. Open PBCs are applied to both the side (radial PBC) and top boundaries (axial PBC).

Numerics

The discretisation of equations ( 6), ( 7) and ( 8) into the Yee grid [28,[START_REF] Van Den Donker | Modelling microwave plasmas for deposition purposes : exploring the freedom in space and chemistry[END_REF] generates a system of linear equations

A • v = b ( 22
)
where A is a matrix that represents the operator for the Maxwell equation and the boundary conditions, v a vector formed by the discrete electromagnetic field components located at the grid positions, and b the source vector related to the excitation BC.

The creation and management of the corresponding system of linear equations is performed using the LinSys library (see [23]). A direct solver (SuperLU, see [43]) is used to compute the solution of equation ( 22), i.e. the electromagnetic fields at the grid positions. From the distribution of the electromagnetic fields, the power density Q is computed using a discrete version of equation ( 9) and the total power is computed integrating it over the volume. The Poynting vector is computed using equation (11) for each grid cell.

The Poynting vector is integrated over the surface area at the boundaries of bottom (excitation) entrance, top, and lateral sides. This gives the net power through these boundaries which is used in section 4 to define the reflected and transmission power coefficients.

The energy conservation equation ( 14) is applied to each grid cell. For a grid with N ∼ 10 5 , the computation time is t c ∼ 10s, with a numerical precision of 10 -12 , more specifically for N = 75932 the problem is solved within t c = 4.94 s in a desktop computer using a processor of type AMD PHENOM II X4 955 at 3.2 GHZ and 530 MB out of 8 GB of memory of type DDR3 at 1333 MHz. The operating system is Ubuntu 10.04. This makes the EM model adequate for a future application in a self-consistent description of the fluid and electromagnetic aspects of MIPs.

Results and Discussion

In order to study the power response of the various torch configurations, we define at the excitation boundary the incident P i and reflected power P r . Together they form the net power, P ns = P i +P r , supplied to the MIP. At the boundaries (of the domain) we have P t , that is the power transmitted to the surroundings. This can be divided into P t = P t,z + P t,r thus the transmission at the top (axially emitted) and lateral side (radially emitted) (see figure 4). A positive sign refers to power leaving the domain (e.g. P t,z and P t,r ), whereas a negative sign refers to power entering the domain (e.g. P i ). Also P abs is positive; the energy used for heating the electrons can be regarded as leaving the EM domain. With this sign-convention equation ( 14 We can see in figure 5, that for the Closed configuration, the values for the transmission are absent. This is evident as the presence of both the radial and axial metallic limitation obstructs propagation. So only the values of the A-and R-coefficients remain. They strongly depend on the value of R 4 . Strong absorption peaks ('resonances') are seen for R 4 = 49 mm, 78 mm and 98 mm. The positions of these peaks are comparable to those found with the analytical solution of a circular waveguides [27,[START_REF] Pozar | Microwave engineering[END_REF]. However the presence of a plasma can change the positions and shapes of these peaks drastically. We can see this in figure 12, where curves A 1 and A 2 correspond to n e = 10 21 and 5 • 10 21 m -3 , respectively, both with ν/ω = 100. Note how the increment in n e modifies the position and strength of the absorption peaks.

Below a critical radius (R c = 47 mm for f = 2.45 GHz) the character of the propagation changes and both C-and S-configurations share the same power response. The reason is that for R < R c , the dispersion relation [START_REF] Álvarez | [END_REF] provides an imaginary wavenumber and therefore only non-propagative solutions exist in the circular waveguide.

Another trend is found in the S-configuration as we can see in figure 6. There is no radial transmission across the metal wall, so that T r = 0. Again the absorption and reflection depend on the radial position but this dependency is not so strong as in the C-case. The absorption peaks at R 4 = 49 mm and 78 mm found in the C case are no longer present; only a the broad absorption peak at R 4 = 98 mm is left over. Therefore, the operation of the torch around R 4 = 98 mm looks reasonable.

As we said before for R < R c , both the S-and C-configurations behave the same due to the presence of non-propagating solutions.

For the O-configuration we see in figure 7, that the power emitted in axial direction can be neglected and that all the power response coefficients are more or less independent of the radius R 4 for R 4 > 40 mm. The latter can be seen as a validation of the propagation boundary condition as placed at R 4 . However, for R 4 < 40 mm, the PBC is so close to the plasma that the electromagnetic wave is different from that assumed at the boundary and the reflections increases. Therefore the PBC formulated in equation ( 18) is not adequate for these low distances. Moreover, for R 4 = 14.5 mm the radial transmission becomes negative i.e. power enters from the side instead of leaving. Of course, this cannot happen in reality and so we have to conclude that it can not be used for so small values of R 4 .

The relative change is defined as ∆f/f = max(f )-min(f )] max(f )

, with f = A,T, R. For R 4 > 40 mm, we find slight changes with a relative change ∆A/A = 3%, which is even lower for the range R 4 > 80 mm (where ∆A/A = 1% is found) in T, A and R that can be attributed to the presence of numerical dispersion i.e. the wavenumber of the wave changes at it propagates through the numerical grid ( see for instance section 3.4.4 in [START_REF] Rumpf | Design and optimization of nano-optical elements by coupling fabrication to optical behaviour[END_REF]). This can be improved with an optimum grid refinement (for a detail explanation on numerical dispersion for different spacial and time discretization schemes see [START_REF] Liu | [END_REF]).

The low value of T z for the O-configuration is better understood when looking to the radiation pattern of a monopole antenna close to a metal ground. As we can see in figure 8, the monopole antenna emits much more power in the radial than in the axial direction.

The figures 9, 10 and 11 show the influence of the axial extension of the expansion zone. In Figure 9, giving A and R as a function of L = Z 4 -Z 1 for the C-case, we see, just as in figure 5, sharp changes in A and thus R. The maximum in A takes place at L = 150 mm, which is the value used for Z 4 = 250 mm in figures 5, 6 and 7.

For the O-and S-configuration, figures 7 and 6, we find only a slight dependence on Z 4 , which again can be seen as a validation for propagation boundary condition performed at Z 4 .

In the O-configuration, for L > 124 mm the relative change for the absorption is ∆A/A < 1%. Between 124 mm and 84 mm, ∆A/A ≤ 3%. However, for L = 64 mm, the relative change becomes ∆A/A > 9%. This mean that the PBC is not adequate when placed at L < 84 mm i.e. a distance to the plasma of 54 mm for a plasma length of 20 mm.

In the S-Configuration, for L between 64 mm and 210 mm the relative change is less than 10%. This high relative change can be attributed to the presence of numerical dispersion ( [START_REF] Rumpf | Design and optimization of nano-optical elements by coupling fabrication to optical behaviour[END_REF]).

The above makes clear that especially the closed configuration shows absorption peaks as functions of the sizes of the enclosure. In previous studies it was found that the dependency of the resonances on the dimensions of the enclosure, strongly depends on plasma conditions (see [27]). If for instance the cavity is tuned such that the power absorption is at maximum for a plasma driven in argon, it will be brought out of resonance if a small amount of H 2 or other molecular species is introduced. After such a change in plasma chemistry the cavity must be retuned. This implies that the C-configuration is not so suitable for plasma applications. It is illustrative to refer to the Evenson cavity [46] and the Beenakker cavity [47], which are not longer popular in spectrochemical analysis due to the resonances.

From now on we confine ourselves to the open case which indeed nowadays is the most often used configuration.

Effect of the electron density value

As stated before the plasma is located at the transition region where the electromagnetic structure of the coaxial waveguide converts into that of the expansion region. In order to investigate in how far the plasma plays a role in this conversion we studied the power response of the open configuration as a function of the electron density. In all cases the plasma shape is fixed to that of a cylindrical bar with a radius and length of 1 mm and 20 mm respectively. The results, given in figure 13, show that, just as in figures 7 and 11, the transmission in axial direction can be neglected. Thus the dominant power channels are the radial transmission T r , the absorption A

A, R, T z , T r [%] R 4 [mm] Open R c A R T r T z Figure 7.
Power response of the Open configuration with the same height and model plasma as that used for figures 5 and 6; the absorption, reflection coefficients and transmission are given as functions of the radius R 4 . Note that there are no absorption peaks for R 4 ≥ 40 mm where the relative change is around 3%. For R 4 < 40 mm the propagation boundary condition is apparently not adequate. For R 4 = 14.5 mm the radial transmission even becomes negative i.e. power enters from the radial side. The Open-Configuration behaves as a monopole antenna near a ground metal (see figure 8) and this explains that Tr is much larger than Tz.

and reflection R. It is seen that, starting at n e = 10 21 m -3 the A-value increases from 20 to almost 50% for increasing n e and that this increase takes place at the expense of the R-value.

Above n e = 10 21 m -3 the radial transmission increases sharply thereby decreasing the A -value. Apparently the plasma serves as an antenna and the higher n e , the better the antenna is. However, the better the antenna function the less energy is left over for absorption. This sets an optimum value of the plasma density. At lower n e -values more energy is absorbed thereby increasing the plasma density. However, due to this n e -increase the plasma improves its antenna function, so that more energy is emitted and thus less is used for the absorption in the plasma. This tends to decrease n e which in turn leads to an increase in absorption. This negative feed-back might be (one of) the reason(s) for the fact that a MIP torch show an instable behavior.

Effect of the shape of the transition region (the 'gap')

In order to understand the antenna function we played with other elements of the transition region, hereafter shortly denoted as the gap. The possible structures are given in figure 14 where the model-plasma remains unaltered; it is bar-shaped with a electron density of 10 21 m -3 while the radius and length are 1 mm and 20 mm, respectively. Figure 14(a) shows the influence of the extension of the inner conductor given by the insertion depth d = Z 2 -Z 1 , figure 14(b) shows a perpendicular partial coverage of the gap, while figure 14(c) deals with an oblique partial coverage parallel to the tip of the inner conductor. The perpendicular coverage of figure 14(b) is inspired by [10] (see figure 1(a)). The oblique coverage is inspired by [48]. The results are schematically given in figure 15.

It is found that lowering the inner conductor from d = 6 mm to -6 mm leads to an enhancement of the reflection from about 20% to 31% but also that the absorption increases The Open-Configuration behaves as a monopole antenna near a ground metal (see figure 8) and this explains that Tr is much larger than Tz. Note that the absorption peaks are not longer present and that for L > 124 mm the relative change for the absorption is ∆A/A < 1%. Between 124 mm and 84 mm, ∆A/A ≤ 3%. However, for L = 64 mm or lower, the relative change becomes ∆A/A > 9%.

from about 46% to 55%. In figure 16, the extension of the nozzle is further reduced and even withdrawn into the outer cylinder of the coaxial waveguide. As we can see for d < -25 mm the complete plasma is hidden in the outer conductor wall and the power response is that of the C-and S-configurations, figures 5 and 6, for R 4 = R 3 = 14.5 mm. This operation seems lucrative since large absorption-values are obtained. However such a plasma in a cylinder is not very practical as it might cause ablation of the cylinder wall. Moreover, since the main plasma part will not operate in an open region, it will not be applicable for spectrochemical analysis or for the treatment of surfaces and volumes. The coverage of the gap has a substantial effect on the reflection which changes from 20% (no coverage) to about 32% (oblique coverage) and 52% (horizontal coverage). So the negative effect Power response of the Open configuration with the same expansion zone radius as that of figures 9 and same height of that of figures 5; the absorption, reflection and transmission coefficients are given as functions of the electron density ne. At lower ne-values more energy is absorbed thereby increasing the plasma density. However, due to this ne-increase the plasma improves its antenna function, so that more energy is emitted and thus less is used for the absorption in the plasma. This tends to decrease ne which in turn leads to an increase in absorption.

of coverage is that more power is send back to the magnetron power generator; the positive side is that from the net power admitted to the torch a larger fraction is absorbed by the plasma.

Conclusions

A flexible and versatile EM model is constructed using the PLASIMO platform [23,24]. It is employed to explore the EM features of axi-symmetric microwave induced plasma (MIP) torches.

The basis of the model is formed by the Maxwell equations used in a full vectorial form, meaning that for the TM mode the three vector components E r , E z and H φ are solved in coupled manner. This is an important difference with most of the others methods found in literature for which the EM problem is described as a Helmholtz-like wave equation for the H φ component solely. to 31% but also that the absorption increases from about 46% to 55%. For d < -25 mm the complete plasma is hidden in the outer conductor wall and the power response resembles that of the C-and S-configurations, figures 5 and 6, for R 4 = R 3 = 14.5 mm Note that Tz = 0, because the region below the outer conductor for z between Zp and Z 3 acts as a circular waveguide which radius R 3 is below the critical radius Rc = 47 mm and therefore waves cannot propagate.

The full-vector approach has the advantage that the boundary conditions can be implemented in a natural way. This, together with the employment of a direct matrix solver, leads to a speed-up in convergence by a factor 100 when compared to [27]. For a grid of N ∼ 10 5 cells the solution is reached in t c ∼ 10 s whereas the iterative process depicted in [27], demands t c ∼ 10 3 s for N ∼ 10 4 .

Three types of configurations were investigated; they are classified as Closed C-, Semi-open Sand Open -configurations. Changing from one configuration-type to an other is realized by simply changing the boundary condition of the expansion zone.

It is found that the C-configuration, for which the expansion takes place in a closed cylindrical cavity, strongly depends on the radial R and axial extension L. For various R-and L-values sharp resonances (strong absorption peaks) where found. Especially the radial sizes for which resonance are obtained are extremely sensitive to plasma features. This makes the C configuration less suitable for applications in the field for which (sudden) changes in the plasma conditions can be expected such as found in the field of spectrochemistry and the treatment of surfaces and volumes.

The features of the S-configurations, for which the wave expansion takes place in an open circular wave guide, are much less dependent of the sizes and therefore not so sensitive to plasma chemistry. These S-configuration can be seen as a model for the practice often used in laboratories to place MIP torches in a open cylindrical grid with the purpose to protect the environment against microwave EM radiation. Calculations show that this is good practice provided the radius of this grid is around 90 mm.

Most attention is paid to torches of the open configuration as these are most often used in practice. It is found that the power balance as expressed by the values of the absorption, transmission and reflection depends on the electron density of the plasma. The reason is that plasma is an essential element in the conversion of the EM waves of the co-axial structure to the EM-structure of the expansion region. This antenna function is better for increasing n e -values. However, this implies that for increasing n e the absorption decreases so that the plasma does not get enough energy to keep n e high. This negative feed-back might be the reason for the instabilities often seen in the operation of open configured MIP torches.

Apart form the electron density we also investigated other aspects of the antenna-say gapregion. So for instance it is found that the plasma absorption strongly depends on the presence of a metallic obstruction placed at the end of the coaxial waveguide and the degree at which the inner conductor of this coaxial waveguide extends into the expansion zone. The influence of two different gap obstructions were investigated namely a horizontal and an oblique one. It is found that these both enhance the reflection and absorption.

Due the flexibility of this EM model it can be used in future for further explorations of the impact of the various construction elements on the performance of torches and other types of MIPs. Moreover, it is ready to be used in combination with a set of transport equations to get a self-consistent description of the plasma-EM interaction.

Figure 1 .

 1 Figure 1. Left: a schematic view of a TIA. From the EM point of view the TIA basically consists of a coaxial structure placed perpendicular to a rectangular waveguide (WR-340). The plasma gas flows inside the inner conductor of the coaxial waveguide from the bottom via the nozzle to the expansion region. The EM energy generates the plasma just above the nozzle, thus in the region where the EM-transition takes place from the coaxial waveguide to the expansion region. Right: a schematic view of the two main plasma regions: the ionizing and recombining zone.

Figure 2 .

 2 Figure 2. (a) the schematic geometry of the electromagnetic structure for the TIA showing 2D axial-symmetry; (b) the computational grid (see Table.1 for values of geometry parameters); (c) detailed view of the model plasma, in this case a bar.

Figure 3 .

 3 Figure 3. (a): general structure of a (Yee) cell used in the numerical method; (b): example of interfaces between different types of material cells.

  at boundary positions x ∈ {0 ≤ z ≤ Z 1 , r = R 2 }. • Propagation boundary condition (PBC) for axial propagation . A homogeneous mixed boundary, i.e. type (iii), is used at the top-boundary positions x ∈ {z = Z 4 , 0 ≤ r ≤ R 4 }, of the S and O configuration with F = E r , q( x) = iβ, and w( x) = 0. Thus

Figure 4 .

 4 Figure 4. Different types of expansion regions: (a) Closed: a circular waveguide with short-cut at top i.e. a PEC-BC at the side and at the top, (b) Semi-Open: a circular waveguide open at the top side i.e a PEC-BC at the side and axial PBC at the top and (c) Open: Open air i.e. radial PBC at the side and axial PBC at the top.

  ) transforms into P = P i + P r + P t + P abs = 0(23) Substituting P i = -|P i | and dividing all terms by |P i |, equation (23) can be written as1 = R + T + A(24)whereR = |P r |/|P i |, T = |P t |/|P i |and A = |P abs |/|P i | are the reflection, transmission and absorption (power) coefficients of the configuration. The transmission coefficient can be further specified as T = T z +T r . The results given below are all based on simulations taking |P i | = 100 W. The results of R, T , and A are obtained after normalizing the corresponding powers by |P i |, thus dividing by 100 W. The relative importance of R,T and A will be denoted as the power response of the configuration.

4. 1 .

 1 Effect of the Remote Field Boundary ConditionsThe power responses as functions of the vessel radius R 4 of the three different configurations C, S, and O are given in figures 5, 6 and 7. They have in common that the model-plasma has the same electron density n e = 10 21 m -3 and shape: a cylindrical bar with a radius of 1 mm and an axial extension of 20 mm. The height of the expansion zone, L = Z 4 -Z 1 , in all three cases equals 150 mm which corresponds to Z 4 = 250 mm. It should be noted, however, that only in the C-case the axial extend is limited by a physical (i.e. metallic) boundary. For the O-and S-cases the axial extend represents the position of the axial propagation BC. Mutatis mutandis the same applies to the radial extend as given by the R 4 -value. In the C-and S-case this is a physical boundary, in the O-configuration the R 4 -value represents the position of the radial propagation BC.

Figure 5 .

 5 Figure 5. Power response of the Closed configuration with a bar-shaped plasma with of radius of 1 mm, axial extension 20 mm and an ne-value of 10 21 m -3 . The height of the expansion zone equals 150 mm which corresponds with Z 4 = 250mm; the absorption and reflection coefficients are given as functions of the radius R 4 . As the configuration is closed the transmission is zero. Clear absorption peaks (resonances) can be seen at R 4 = 49 mm, 78 mm and 98 mm. For R 4 below the critical radius Rc = 47 mm, the behaviour is the same as in the Semi-open case because there is no propagating solution through the circular waveguide and the results are independent of the presence of a top metal wall.

Figure 6 .

 6 Figure 6. Power response of the Semi-open configuration with the same height and model plasma as that of figure5; the absorption, reflection and transmission coefficients are given as functions of the radius R 4 . As the configuration is closed in radial direction Tr is zero. Note that the absorption peaks at R 4 = 49 mm and 78 mm found in the C case are no longer present; only a more broad absorption peak at R 4 = 98 mm is left over. For R 4 below the critical radius Rc = 47 mm, the behaviour is the same as in the closed case because there is no propagating solution through the circular waveguide and the results are independent of the presence of a top metal wall.

Figure 8 .Figure 9 .

 89 Figure 8. Radiation (transmission) pattern of a monopole antenna showing that an antenna close to a metallic ground transmits much more power in the radial than in the axial direction. Obtained using EZNEC (http://www.eznec.com/).

Figure 10 .Figure 11 .

 1011 Figure 10. Power response of the Semi-open configuration with the same expansion zone radius and model plasma as that of figure 9; the absorption, reflection and transmission coefficients are given as functions of the height of the expansion zone L = Z 4 -Z 1 . As the configuration is closed in radial direction Tr is zero. Note that the absorption peaks are not longer present and that the relative change is around 10%.

2 Figure 12 .Figure 13 .

 21213 Figure 12. Absorbed power response as function of R 4 in the closed configuration with the same height as that of figure 5. The model plasma is the same with exception that ne equals 10 21 m -3 in A 1 and 5 • 10 21 m -3 in A 2 . Note how the increment in ne modifies the position and strength of the absorption peaks.

Figure 14 .Figure 15 .Figure 16 .

 141516 Figure 14. Different types of transition or 'gap-' regions: (a.1) Open gap with nozzle insertion of d = Z 2 -Z 1 = 6 mm, (a.2) Open gap with nozzle insertion of d = -6 mm, (b) ObliqueObstruction (similar to[48]) and (c) Obstruction perpendicular to the nozzle; similar to the gap in the TIA design[10] . The extension of the nozzle in (b) and (c) is that of (a.2). All of them are in O-configuration.

Table 1 .

 1 Table.1 for values of geometry parameters); (c) detailed view of the model plasma, in this case a bar. Parameter distances as used in the model-setup (see figure 2(b))

		Parameters of the set-up
	Parameter	Value	Parameter	Value
	R 1	1 mm	Z 1	100 mm
	R 2	5.5 mm	Z 2	106 mm
	R 3	14.5 mm	Z 3	118 mm
	R 4	variable	Z 4	variable

Table 2 .

 2 Summary of the main characteristics of the electromagnetic model.

	Main differences with other electromagnetic models
		current	Bilgic	Gritsinin	Alves
			[25]	[26]	[27]
	Fields (TM)	H φ , E r , E z
	Equation	Full Vector	Azimuthal Magnetic Component
	(Comps. )	(H φ , E r , E z )		(H φ )
	Grid	Staggered		Centered
		(Yee)
	BC for	E r , E z		H φ
	Torch	TIA, MPT	MPT	CMT	TIA
		and CMT
	Different 'Gaps'	3	1	1	1
	Configurations	C, S, O	O	S	C
	Open Boundary Plane Wave	Plane Wave Plane Wave	No
	at Exp. Reg.	Radial Wave
	Solver	Direct	Direct	Iterative	Iterative
	Grid Points, N	∼ 10 5	∼ 10 4	Not given	∼ 10 4
	Comp. Time t c	∼ 6 s	Not given	Not given	∼ 600 s

Table 3 .

 3 Comparison of our electromagnetic model with those of

Acknowledgments

This work is supported and funded by the Dutch Technology Foundation STW, the Eindhoven University of Technology and Draka Communications. The authors would like to thank A. Sola and A. Gamero from the University of Cordoba (Spain) for the fruitful discussions.


W W F Chong

M Teodorescu
email: m.s.teodorescu@cranfield.ac.uk

H Rahnejat

Effect of lubricant molecular rheology on formation and shear of ultra-thin surface films

Keywords: Molecular level tribology, solvation, viscous shear, cavitation, nano-scale non-Newtonian friction

Physics of molecularly thin fluid films, formed between surface features at close range is investigated. It is found that the interplay between discrete lubricant drainage from such contacts and localised contact deflection plays an important role both on the load carrying capacity of these asperity level conjunctions as well as on friction. Small spherical molecules tend to solvate near assumed smooth surfaces of asperities at nano-scale. Their discrete drainage at steadily decreasing gaps adds to the viscous friction of any bulk lubricant film. However, at the same time the generated solvation pressures increase the load carrying capacity. Conversely, long chain molecules tend to inhibit solvation, thus show a decrease in the load carrying capacity, whilst through their wetting action reduce friction. Consequently, real lubricants should comprise molecular species which promote desired contact characteristics, as indeed is the case for most base lubricants with surmised properties of certain additives. The methodology presented underpins the rather empirical implied action of surface adhered films. This is an initial approach which must be expanded to fluids with more complex mix of species. If applicable, this could also be an alternative (potentially time saving) approach to Monte-Carlo simulations for molecular dynamics.

Introduction

There is a progressive demand for accurate models to predict boundary friction between real (rough) surfaces. The main goal is to reduce frictional losses, and therefore, improve mechanical efficiency. Most frictional losses are as the result of any direct interaction of rough surface asperities. Nevertheless, a significant part of friction is also due to the shear of a thin lubricant film, usually trapped between opposing asperities of the contacting surfaces. In some applications an ultra-thin surface film is only possible, where molecular behaviour of liquid lubricant near surface asperities is of significance. These are usually lightly loaded contacts. However, in the micro-scale conjunctions of light to moderate load intensity a thin lubricant film adsorbed or entrained in between asperities can also behave in the same manner. This is also true of any lubricated conjunction at start-up or during inlet boundary reversal, where a coherent lubricant film cannot be assured.

New manufacturing techniques can now provide better control of surface topography (e.g. roughness, texturing). Molecular rheology of lubricant can also be improved by choosing an appropriate composition. This means that additives, having desired molecular interactions with solid bounding surfaces, are added to the base lubricants. It is hoped that the bulk rheology of the mixture would improve the load carrying capacity of the overall conjunction and possess low shear resistance. At the same time molecular behaviour near the smooth solid boundaries or around surface asperities would reduce friction when a coherent film is not formed. Most tribological conjunctions are subject to transient conditions, such as in start-up, inlet boundary reversals or because of application of cyclic loads. Therefore, interactions occur between the various phenomena at different physical scales; between bulk micro-scale lubrication and nano-scale localised interactions (e.g. around asperity tip pair contacts).

At micro-scale, viscous shear is the dominant mechanism underlying friction, as in hydrodynamic and elastohydrodynamic conditions. These mechanisms are better understood. At nano-scale, boundary interactions are affected by other kinetic laws such as intermolecular and surface energy effects. There is, of course, the critical interaction between these phenomena, particularly in the intervening scale. This interaction is between surface topography and lubricant rheology at short range (thin surface films of dimensions of several molecular diameters of the base lubricant or its additives). There is no deterministic model representing the formation of such ultra-thin films, particularly with respect to their interactions with the bulk tribological conditions. This study is an initial attempt to address this shortcoming.

At one end of the interaction range one can consider a single hemisphere. This represents a surface feature (such as an introduced dimple) or an idealised asperity tip, sliding against a perfectly smooth substrate, furnished with a uniform film of a basic lubricant (of a single molecular composition). Such a lubricant may represent a thin protective layer or a sealant, which when subjected to shear forms a thin surface adsorbed film. Octamethycyclotetrasiloxane (OMCTS) is used as a low friction protective layer or a sealant in some applications, for example in siliconised rubber seals. These are often used in inhalation devices. OMCTS is a non-polar fluid with molecular diameter of around 1nm. Al-Samieh and Rahnejat [START_REF] Al-Samieh | Ultra-thin lubricating films under transient conditions[END_REF] have shown that near molecularly smooth surfaces the behaviour of such a fluid is dominated by the solvation effect for films of up to several molecular diameters deep. Chan and Horn [START_REF] Chan | The drainage of thin liquid films between solid surfaces[END_REF] have shown that there is discretised drainage of OMCTS from such conjunctions, indicating its near surface solvation effect. It is postulated that non-polar small molecular diameter spherical molecules (such as those of OMCTS) reduce localised viscous friction near asperity tips at nano-scale, which one may assume to be ideally smooth at this range. This is because solvation leads to dewetting of contacting surfaces.

At the other limit of this investigation, there are tribological conjunctions in microscale, where the load carrying capacity of the contact must be assured. Therefore, dewetting of the contiguous surfaces is not considered as a viable option. Consequently, long chain molecules are used to discourage density perturbation of the lubricant near the bounding contacting solid surfaces. Israelachvili [START_REF] Israelachvili | Intermolecular and surface forces[END_REF] suggests that long chain molecules are best suited for this purpose and indeed are used regularly as additives to base lubricants as boundary friction modifiers. This is another postulate which should be investigated. In the current study the same sliding contact model is used, but with an idealised lubricant film, based on long chain molecules of Hexadecane or Tetradecane. These are used in diesel fuel, where they act as the lubricant in off-road vehicle fuel pumps. They intervene between the sliding teeth surfaces of pump gears.

The mechanism of delivery of the lubricant into the contact in both the examples investigated is assumed to be through hydrodynamics. The method of solution, therefore, combines hydrodynamics of the conjunction with solvation of the lubricant near the smooth surface of an assumed hemispherical feature. The effects of starvation and cavitation beyond the contact exit are also taken into account. Starvation affects the supply of lubricant into the contact, thus the thickness of film which influences solvation. Inclusion of cavitation and lubricant reformation influences the load carrying capacity of the contact, which determines the instantaneous equilibrium load shares of viscous and surface force contributions. Hitherto, such a study has not been reported in open literature.

Mathematical Model

Figure 1 is a schematic representation of contact cross-section between a rigid hemispherical feature, sliding relative to a smooth flat semi-infinite elastic half space. Under the usual conditions encountered in lubricated contacts, a number of distinct regions would occur in the conjunction. There is the inlet wedge region, where the lubricant is entrained into the contact (in the opposite direction to the sliding motion), from a curved meniscus. This is followed by a region of parallel film (corresponding to high pressures) with any small strain localised deformation in accord with Hertzian assumptions (elastohydrodynamic region). At the exit region from the Hertzian zone, there is film rupture boundary, beyond which falling pressures below that of lubricant vaporisation leads to streamers; formation of vapour fingers in the lubricant flow. This is the cavitation region. Somewhere downstream of this region, lubricant film is assumed to reform. This rather orderly and now well understood lubricant flow characteristics is based on the principle of continuity of flow, which is fundamental to the bulk lubricant behaviour. Abd Al-Samieh and Rahnejat [START_REF] Al-Samieh | Nano-lubricant film formation due to combined elastohydrodynamic and surface force action under isothermal conditions[END_REF] showed that any near surface effect plays an insignificant role in continuity of flow, until the film thickness is gradually reduced to several molecular diameters of the lubricant. Then, discrete lubricant molecular behaviour pursues, followed by a stepwise drainage from the contact, roughly in accord with the molecular diameter of lubricant, which are assumed to be spherical. This finding conformed to the experimental observations of Chan and Horn [START_REF] Chan | The drainage of thin liquid films between solid surfaces[END_REF] and predictions of Matsuoka and Kato [START_REF] Matsuoka | An ultra-thin liquid film lubrication theory calculation method of solvation pressure and its applications to EHL problem[END_REF].

Motion Figure 1. Lubricated conjunction (Corrected)
Therefore, for fairly thin films, the pressures generated in the conjunction are due to a number of mechanisms, whose individual contributions vary according to the gap size, h(x, y) :

h(x, y) = h 0 + h s (x, y) + δ(x, y) ( 1 )
where h s (x, y) is the local gap at any location within the conjunction and h 0 is the minimum clearance.

For inert surfaces and non-polar lubricants, assuming no direct contact between the surfaces, the total conjunctional pressure is due to hydrodynamics, van der Waals pressure and solvation:

p T = p h + p w + p s (2)
where: p h is the hydrodynamic component of the pressure, p w is the van der Waals component and p s is the solvation contribution.

It should be noted that the contribution due to van der Waals interaction can become significant at almost vanishing separations.

The load carrying capacity of the contact is obtained as:

W = p T dA (3)
where p T > 0 means that negative (attractive) pressures do not contribute to the load carrying capacity.

Hydrodynamic Pressure

Al-Samieh and Rahnejat [START_REF] Al-Samieh | Ultra-thin lubricating films under transient conditions[END_REF] and Matsuoka and Kato [START_REF] Matsuoka | An ultra-thin liquid film lubrication theory calculation method of solvation pressure and its applications to EHL problem[END_REF] solved the Reynolds equation to obtain the hydrodynamic contribution. For the point contact geometry, the two dimensional Reynolds equation is:

∂ ∂x ρh 3 η . ∂p h ∂x + ∂ ∂y ρh 3 η . ∂p h ∂y = 12 U ∂ρh ∂x + V ∂ρh ∂y + d dt (ρh) (4)
where, U is the speed of entraining motion of the lubricant into the conjunction, being the average speed of the contacting surfaces (or half the speed of sliding in the case studied here) and V is the speed of any side-leakage of lubricant out of the contact area.

The hydrodynamic load carrying capacity of the contact is confined to the regions of full film, prior to a region of cavitation beyond a film rupture point. Swift-Stieber exit boundary conditions is used: P = ∂ P /∂x| cavitation = 0 red [START_REF] Swift | The stability of lubricating in journal bearings[END_REF][START_REF] Stieber | Hydrodynamische Theorie des Gleitlagers, Das Schwimmlager[END_REF]. Cavitation, though, plays a very significant role in the disruption of free surface film uniformity as desired for thin protective surface films. Additionally, Reynolds equation is not suited to the prediction of the cavitation region, where the interface between the lubricant and air causes loss of laminar flow. This is an underlying assumption for Reynolds equation.

To overcome this shortcoming the cavitation model proposed by Elrod [START_REF] Elrod | A cavitation algorithm[END_REF] is used here. This model modifies the Reynolds equation to take into account the boundary between a full film region and that representing a cavitation region. The approach approximates the fully demarcated vaporization cavitation boundary condition proposed by Jakobsson and Floberg [START_REF] Jakobsson | The finite journal bearing considering vaporization[END_REF] and Olsson [START_REF] Olsson | Cavitation in dynamically loaded bearings[END_REF] (now referred to as the JFO boundary condition). In the simplified, but computationally more efficient Elrod method a fraction film content ratio, 0 < θ ≤ 1 is defined. The hydrodynamic pressure can now be defined as:

p h = gβlnθ + p c ( 5
)
where cavitation is assumed to occur at the constant pressure, p c under the current isothermal analysis, β is the lubricant bulk modulus and g is a switching function:

g = 1 ⇒ Full film, if θ ≥ 1 0 ⇒ Cavitation, if 0 < θ < 1 (6)
Now for the full film region, Reynolds equation is modified to:

∂ ∂x ρ c h 3 η gβ ∂θ ∂x + ∂ ∂y ρ c h 3 η gβ ∂θ ∂y = 12 U ∂θρ c h ∂x + V ∂θρ c h ∂y + d dt (θρ c h) (7)
Referring to [START_REF] Swift | The stability of lubricating in journal bearings[END_REF], it is clear that in the cavitation region the Couette flow only operates, thus:

u ∂ ∂x + v ∂ ∂y + ∂ ∂t (θρ c h) = 0 (8)
It is reasonable to assume negligible side leakage from the contact with thin films and at relatively low sliding speeds, thus: V = 0 . With this assumption taken into account, it is clear that the cavitated region is governed by a balance between the squeeze film effect (the ultimate term in the last equation) and the drainage of the fluid from the contact on the account of hydrodynamics alone. Separation of the surfaces on account of squeeze; ∂θh/∂t > 0 increases the chance of cavitation, which is the same as an increase in the speed of entraining motion (Al-Samieh and Rahnejat [START_REF] Al-Samieh | Physics of lubricated impact of a sphere on a plate in a narrow continuum to gaps of molecular dimensions[END_REF]). One should then expect reduced cavitation at lower sliding speeds as an outcome of the analysis, which would encourage better lubricant distribution on the free surface post film reformation. This is a desired outcome in any reciprocating contact.

Lubricant rheology, therefore, is critical to the maintenance of a continuum. Yet, density alters as a direct result of the different kinetic laws that co-exist in the nanoscale. In bulk it alters with hydrodynamic pressure, and near the solid barrier it has a higher density profile due to the constraining effect of solvation pressure (section 2.2 ). The bulk density variation with the contact pressure for mineral oils is given by Dowson and Higginson [START_REF] Dowson | Elasto-hydrodynamic lubrication: The fundamentals of roller and gear lubrication[END_REF]:

ρ = ρ o 1 + 0.6 × 10 -9 × p h 1 + 1.7 × 10 -9 × p h (9)
As a first approximation this is taken to also represent the density-pressure dependency in this analysis. Matsuoka and Kato [START_REF] Matsuoka | An ultra-thin liquid film lubrication theory calculation method of solvation pressure and its applications to EHL problem[END_REF], Abd Al-Samieh and Rahnejat [START_REF] Al-Samieh | Nano-lubricant film formation due to combined elastohydrodynamic and surface force action under isothermal conditions[END_REF] and Al-Samieh and Rahnejat [START_REF] Al-Samieh | Ultra-thin lubricating films under transient conditions[END_REF][START_REF] Al-Samieh | Physics of lubricated impact of a sphere on a plate in a narrow continuum to gaps of molecular dimensions[END_REF] also assumed the same for the bulk hydrodynamic behaviour of the lubricants such as those described above and used in the current analysis. Clearly, for the cavitated region: ρ = ρ c when p = p c .

Lubricant viscosity-pressure variation can be predicted using the Roelands' equation [START_REF] Roelands | Correlational aspects of the viscosity-temperature-pressure relationships of lubricating oils[END_REF]:

η = η o e αp h ( 10
)
where α =

Solvation Pressure

For very smooth surfaces, the solvation effect dominates in gaps of several molecular diameters of the intervening fluid. The phenomenon is as the result of density variation of liquids near solid boundaries [START_REF] Israelachvili | Intermolecular and surface forces[END_REF]. The effect is significant in nano-scale conjunctions, inducing a structureless conjunctional environment which disrupts the comfort of a hydrodynamic continuum and thus proper wetting of the contiguous solids. Paradoxically, the dewetting action of solvation can guard against the formation of a meniscus, which otherwise would contribute to the work of adhesion.

It is postulated that the effect of solvation is more pronounced for small fluid molecules such as those of perfluoropolyether and Octamethylcyclotetrasiloxane, which have nominally spherical molecules of 1-1.5nm diameter. This effect is almost negligible for long chain molecules [START_REF] Israelachvili | Intermolecular and surface forces[END_REF]. The contact density variation near the smooth solid surfaces takes the form of (Henderson and Lozada-Cassou [START_REF] Henderson | A simple theory for the force between spheres immersed in a fluid[END_REF] and Evans and Perry [START_REF] Evans | Liquids at interfaces: What can a theorist contribute?[END_REF]) :

ρ = KT [ρ(h) -ρ(h → ∞)] ( 11
)
where ρ(h → ∞) is the corresponding value for a single surface, implying that h >> a. Therefore, it tends to the bulk density of the fluid, whose variation may be assumed to follow equation [START_REF] Jakobsson | The finite journal bearing considering vaporization[END_REF]. The density variation near the constraining barriers gives rise to a monotonic exponential pressure variation, which is alternately attractive-repulsive. This variation is explained in some detail by Israelachvili [START_REF] Israelachvili | Intermolecular and surface forces[END_REF], Teodorescu et al [START_REF] Teodorescu | Physics of ultra-thin surface films on molecularly smooth surfaces[END_REF] and, Gohar and Rahnejat [START_REF] Gohar | Fundamentals of tribology[END_REF], resulting in the discrete fluid drainage from the contact conjunction. Israelachvili [START_REF] Israelachvili | Intermolecular and surface forces[END_REF] proposed a solvation pressure of the following form:

p s = -Ce -h/a cos 2π h a ( 12
)
where: C = KT ρ(h → ∞) , which for the cases considered here: for OMCTS C = 172MP a, Hexadecane C = 62MP a and Tetradecane C = 49MP a (Chan and Horn [START_REF] Chan | The drainage of thin liquid films between solid surfaces[END_REF]). Note that the value of C is reduced as the size of the molecule increases. Size is not represented by molecular diameter for long chain molecules, which are rather similar to strings. In fact the diameter of Hexadecane is only 0.4nm. It is clear that the solvation pressure is only significant for small values of h/a .

Van der Waals Pressure

In nano-scale conjunctions other forms of kinetic interaction may also contribute, depending on the molecular structure of the fluid and the atomic structure of the bounding solids. If the fluid is considered as non-polar and the surfaces as inert, then long range van der Waals interactions between the fluid molecules and atoms of the surfaces and between the surfaces themselves play a role. These are weak attractive forces in the nano-scale, but can become very significant at closer range.

The van der Waals forces arising from these interactions can be obtained using Lifshitz theory [START_REF] Lifshitz | The theory of molecular attractive forces between solids[END_REF]. For nano-scale conjunctions, the procedure is detailed by Abd Al-Samieh and Rahnejat [START_REF] Al-Samieh | Nano-lubricant film formation due to combined elastohydrodynamic and surface force action under isothermal conditions[END_REF]. The current study addresses the nano-scale range, where some marginal contributions due to van der Waals interactions may be expected. Depending on the contact geometry, Israelachvili [START_REF] Israelachvili | Intermolecular and surface forces[END_REF] provides a simple expression. For flat surfaces (which may be assumed here because of the localized elastic deformation in the central small contact area in figure 1):

p w = - A h 6πh 3 (13)
where A h is the Hamaker constant, its value depends on the atomic structure of contacting bounding solids and the molecular composition of the intervening fluid. Its value is usually in the range: 10 -21 -10 -19 J , in this case: A h = 1.4 × 10 -20 J.

Contact Deflection

Contact deflection, required in equation ( 1) can now be calculated, using the elasticity potential equation (Johnson [START_REF] Johnson | Contact Mechanics[END_REF] and Gohar and Rahnejat [START_REF] Gohar | Fundamentals of tribology[END_REF]):

δ(x, y) = 1 πE p T (x , y ) (x -x ) 2 + (y -y ) 2
dx dy [START_REF] Henderson | A simple theory for the force between spheres immersed in a fluid[END_REF] where the deflection at a point (x, y) is determined due to all the generated pressures at points (x , y ). When a computational grid is made with the pressure distribution, p k,l , this equation can be stated as:

δ i,j = k l D k,l i,j p k,l (15)
where D k,l i,j are the influence coefficients (Johnson [19]). The reduced elastic modulus is

E = 1 -ν 2 1 E 1 + 1 -ν 2 2 E 2 -1 (16)

Method of Solution

A solution can now be obtained for θ, h, p, δ, ρ and η, by representing equation ( 7) in finite differences and obtaining a solution through low relaxation effective influence Newton-Raphson method with Gauss-Seidel iterations (Al-Samieh and Rahnejat [START_REF] Al-Samieh | Ultra-thin lubricating films under transient conditions[END_REF], Teodorescu et al [START_REF] Teodorescu | Physics of ultra-thin surface films on molecularly smooth surfaces[END_REF] and Chong et al red [START_REF] Chong | Cavitation induced starvation for pistonring/liner tribological conjunction[END_REF]). The procedure places the hemispherical feature at consecutive reducing separations h 0 from the surface of the flat semi-infinite plane. With an assumed speed of entraining motion (U ) and time step, this corresponds to the imposition of a rate of rigid body squeeze (∂h/∂t). Therefore, the conditions at each quasi-static step at various separations h s (x, y) may be regarded as rheonomic nonholonomic constraints. The simulation results in section 5 explore contact conditions with different rates of approach. The convergence criterion used in each quasi-static step is:

i j θ n i,j -θ n-1 i,j i j θ n i,j ≤ θ ( 17
)
where n is the iteration counter. If the above criterion is not satisfied, then:

(gθ

) n i,j = (gθ) n-1 i,j + ωΔθ n i,j (18)
The process is then repeated at the same separation. When the convergence criterion is met, the separation h 0 is reduced by a molecular diameter of the fluid and the entire iterative procedure is repeated again. In this manner the conjunctional behaviour is noted for different quasi-static steps at progressively reduced separations.

Conjunctional Friction

Friction generated by the ultra-thin adsorbed films on nominal molecularly smooth surfaces follow non-Newtonian shear due to chemical reactions which are based on thermal activation. Eyring [START_REF] Eyring | Viscosity, plasticity, and diffusion as examples of absolute reaction rates[END_REF] developed a model which described the viscosity of fluids under such circumstances. Eyring's model assumes that the motion of a volume of fluid molecules takes place in the presence of a cage-like potential barrier by their mere closely-packed arrangement in the bulk. In order to overcome this potential and escape the cage, the fluid molecules have to surmount the activation barrier potential. This is the thermal activation model which Eyring modified to include the effects of prevailing pressure and shear (Figure 2) as:

E y = Q -pΩ -τ φ ( 19
)

Energy

Reaction Path

Figure 2. The potential barrier in thermal activation

The shear stress acts on a stressed activation volume, φ , which may be conceived as a process coherence volume. This can be interpreted as the size of a moving segment in unit shear, whether it would be a molecule or a dislocation line. Therefore, during sliding motion the potential barrier is repeatedly overcome. The average barrier-hopping time for a single molecule is predicted using the Boltzmann distribution with a regular series of barriers and a high stress limit (i.e. τ φ/KT > 1 ). According to Briscoe and Evans [START_REF] Briscoe | The shear properties of langmuir-blodgett layers[END_REF], the shear stress is a function of velocity, pressure and temperature. For constant velocity and under isothermal conditions, it can be expressed as:

τ = τ 0 + ξp ( 20
)
where :

ξ = Ω φ
The Eyring shear stress is determined as :

τ 0 = 1 φ KT ln( U v 0 ) + Q .
where v 0 is the characteristic velocity related to the frequency of the process.

The sliding velocity, however, can change according to the operating speed of the contact. To account for this, the value of τ 0 is adjusted at a constant pressure, p v as:

τ 0 = KT φ lnU + τ 1 -ξp v ( 21
)
where

τ 1 = 1 φ [Q -p v Ω + KT lnv 0 ]
, and for the fluids used in this analysis: KT /φ ≈ 3.4/A , τ 1 ≈ 11.3/A , p v ≈ 100/A all in nN/m 2 and ξ ≈ 0.16 red [START_REF] He | Effect of interfacial liquid structuring on the coherence length in nanolubrication[END_REF] . Note, the contact area is taken as the area of the computational grid: A = dxdy and thus friction is obtained as:

F = τ dxdy; τ > 0 (22)

Results and Discussion

Series of simulations have shown that while the conjunction is separated by a molecularly thin film, the combined effects of solvation and hydrodynamics dominates, while the long range van der Waals interactions play an almost insignificant role. When the gap is greater than several molecular diameters of the intervening fluid, hydrodynamics dominate. As the gap is reduced the effect of solvation becomes increasingly important. For the fluid used (OMCTS) in the first part of this analysis, with molecular diameter a ≈ 1nm , an incipient solvation effect begins around gaps of the order of 8nm (figure 3). Nevertheless, hydrodynamics always plays the key role of entraining a thin film of lubricant into the conjunction through inlet wedge effect. 4) reformation zone (as postulated in description of figure 1). Each region plays a key role. The inlet meniscus acts as the reservoir supplying the contact with a film of lubricant. Lubricant ahead of the contact provides the meniscus liquid bridge through wetting action of the lubricant on the contiguous surfaces, and is indicated by the inward curved isobar ahead of the contact. Under these circumstances, for simplicity it is assumed that a meniscus is formed at a pressure close to that of lubricant's vaporization pressure p c . With entraining motion of the lubricant through the inlet wedge the generated pressures rise. The contact proper comprises a series of almost concentric isobars with increasing pressures towards the rear lubricant exit constriction from the contact. This is a characteristic of hydrodynamic pressures, except that solvation effect in such diminutive gaps causes pressure fluctuations superimposed upon the hydrodynamic pressure profile. Figure 4(b) is a cut through the centre line of the isobaric plot, which shows these pressure fluctuations.

The region of high pressures is immediately followed by oil film rupture, not unlike the behaviour of any hydrodynamic conjunction. With lubricant pressures falling below its vaporization pressure beyond the lubricant rupture point, cavitation occurs. The depth of cavitation is governed by the speed of entraining motion in the contact, as well as the gap size. The former is contributed by hydrodynamics, whilst the latter determines the extent of solvation. The fluid in the cavitation region is a mix of lubricant and vapour. Beyond this, in the trailing edge of the contact, the pressures are assumed to return to the atmospheric value. This implies the existence of a film reformation region, where the vapour escapes through the liquid. Figure 4(b) shows the trend towards reformation. The larger the extent of cavitation and reformation regions become, the greater is the deviation from uniformity of lubricant height upon the free surface of the semi-infinite plane after the exit constriction. It is natural to assume that clumps of lubricant may be formed on any surface beyond the film rupture region if many such hemispherical features or other textures slide along the flat substrate. The analysis here shows that even with assumed steady state condition the transient nature of solvation and its interplay with hydrodynamics can result in non-uniform distribution of lubricant on surfaces. (a) Isobaric pressure plot (MPa) Now returning to figure 3, as the gap is reduced, fluid film discretization occurs (as noted by Matsuoka and Kato [START_REF] Matsuoka | An ultra-thin liquid film lubrication theory calculation method of solvation pressure and its applications to EHL problem[END_REF], and Al-Samieh and Rahnejat [START_REF] Al-Samieh | Ultra-thin lubricating films under transient conditions[END_REF]). Layers of lubricant molecules are drained in a step-wise fashion from the conjunction, in this case at gap intervals of 1nm; the diameter of spherical molecules of OMCTS. In each step a further increase in load or a corresponding squeeze action is needed before another layer of lubricant is ejected out of the contact. This means that at a given sliding speed, whilst the hydrodynamic inlet flow follows Newtonian slow viscous action, the drainage from the contact does not conform to the continuity of flow condition. Solvation in effect accounts increasingly for load carrying capacity of the contact. It disrupts the structured nature of fluid flow and essentially promotes its dewetting. There is, of course, a limit to this height, determined by the meniscus pressure at the inlet (Abd Al-Samieh and Rahnejat [START_REF] Al-Samieh | Nano-lubricant film formation due to combined elastohydrodynamic and surface force action under isothermal conditions[END_REF]). The dotted lines in figure 3 indicate the oscillatory behaviour of solvation, in alternate attractive-repulsive action. The loss of load carrying capacity, indicated by the dotted line is purely theoretical as the hydrodynamic pressure at the nib of the contact guards against this effect, unless no entraining motion is to take place (such as in cessation of sliding).

Figure 5 shows a series of pressure distributions with the corresponding film shapes and fraction film ratios for steadily decreasing gaps in line with the steps indicated in figure 3. As the gap is reduced (e.g. from A to A') at the same sliding speed, the contribution due to solvation pressure is increased. This reorders the molecular content of the conjunction within the contact proper, followed by the drainage of a row of molecules. The film shapes show this discrete molecular drainage. In the intervening intervals (e.g. A' to B) an increase in the contact load (i.e. integrated pressure distribution) is due to molecular disposition to regiment rows of molecules within the conjunction. With increased pressures and thin films the depth of cavitation is increased directly as the result of the dewetting action of solvation. With larger cavitation and assuming a fairly uniform surface topography, a longer time, hence sliding distance, is required before film reformation is achieved through fluid surface tension. Therefore, under dynamic conditions with varying height of the surface features, it is clear that uneven distribution of film would result on the surfaces. The results indicate that this non-uniformity would be exacerbated at lower gap sizes. Hydrodynamic entrainment of the lubricant also has its own better known influence upon film rupture and cavitation. The height of lubricant at the meniscus and its distance to the leading edge of contact is critical to its wetting action (Hamrock and Dowson [24], Wedeven et al [START_REF] Wedeven | Optical analysis of a ball bearing[END_REF]). These are affected by the sliding speed. A faster sliding speed requires a longer inlet distance (i.e. a larger meniscus) in order to avoid starvation of the contact. However, this can also lead to thicker contact films and greater depth of cavitation. These effects lead to conditions that are contrary to the purpose of maintaining a better distribution of liquid lubricant in many applications. Figure 6 shows that the discretisation effect is slightly delayed with increasing speed of entraining motion. This delayed effect is caused by a rise in film thickness, thus more pronounced hydrodynamic action. Figure 7 shows a series of pressure distributions with their fraction film ratios, corresponding to different speeds of entraining motion in figure 6. Note that the rise in pressures is due to hydrodynamics in the inlet nib to the contact (this diminishes with smaller gaps). The solvation contribution remains the same for all the speeds in the same gap size. Interestingly, the depth of cavitation is increased as normally would with enhanced contact pressures, directly as the result of increased speed at the same effective contact load. Hence, reduced gap and increased sliding speed cause non-uniformity in free surface oil film, the former due to solvation and the latter due to hydrodynamic action.

Another key factor is friction. Assuming very smooth surfaces, at the molecular level separations here, friction is mainly due to viscous action of a lubricant film or introduced by a layer of molecules in drainage in accord with the activation energy required to dislocate them. There is clearly insignificant friction in the cavitation region. Figure 8 shows a plot of friction versus the conjunctional load (at a steadily reducing gap). The quasi-linear relationship between any subsequent instances of drainage is due to shear of a lubricant film, with a slope, which may be interpreted as a coefficient of friction. This alters slightly with both the sliding speed and contact pressures. The linear slope is as the result of changes in lubricant viscosity and film thickness, both affected by increasing load between steps of lubricant discretisation. The slope follows viscous friction. Seth et al red [START_REF] Seth | Influence of short range forces on wall-slip in microgel pastes[END_REF] studied slip of microgel paste particles against a solid boundary in the presence of water molecules. They showed that at each step of separation the slip of microgel particles along the wall follow viscous shear.

Referring back to figure 8, the change with pressure is almost insignificant at a given gap size as this is dominated by almost unchanged solvation pressure. However, during drainage a rise in friction occurs with no additional contact load. This is the effort expended to eject a row of molecules. In such instances interpretation of friction based on viscous action (or an effective coefficient of friction) is inappropriate. Here one may regard friction as breaking of cohesive intermolecular bonds during discrete drainage from the contact region or a volume of fluid (i.e. a row of molecules) overcoming an activation barrier to drain from the contact. With increased sliding speed thicker films cause delays in the same characteristics.

Friction characteristics of thin films are quite complex, because they are affected by their actions near solid barriers. This is not only true for OMCTS studied so far, but also for low shear strength films formed by all additives in micro-scale conjunctions in a large variety of engineering applications. One case already mentioned is the action of additive molecules such as Hexadecane in diesel, where the latter also acts as a lubricant, for example in off-road vehicles' fuel pumps. In these cases the surface effects are likely to be concentrated around smooth tips of asperities, which are larger in dimension than the molecular size of the surface active additives. Figure 9(a) shows a similar friction characteristic for Hexadecane, if it were to be used instead of OMCTS in the current study under otherwise exact simulation conditions. However, note that in the case of Hexadecane a small rise in friction is only discerned during lubricant drainage. This is because Hexadecane has long and narrow molecular structure which inhibits solvation (figure 9(b)). Thus, the discrete nature of drainage is reduced. This improves friction, at the expense of better wetting of contact. Hence, a quest to reduce wetness of the contact can have significant drawbacks in terms of friction. However, in certain cases this is a desired function, such as in seals of inhalation devices containing volatile mixtures.

The current analysis is isothermal, thus the effect of friction is understated. It is clear that the heat generated as the result of friction would deplete the film thickness, the load carrying capacity of the contact and exacerbate the effect of cavitation. The analysis shows that for some conjunctions one paradoxically needs the dewetting action of solvation, but uniformity, which is afforded by hydrodynamics. Hence, control of molecular composition of lubricants may be fundamental to future research in ultra-thin film devices.

Conclusion

There is a strong correlation between lubricant molecular structure and tribological behaviour of a conjunction. This becomes even more obvious when the conjunction is very small, at asperity level. The current paper investigates the interaction between local asperity or small feature deflection and molecular-level lubricant drainage as well as their impact on friction and local load carrying capacity of the lubricant. It is found that while long chain molecules inhibit discrete drainage of the lubricant from the contact leading to lower friction, they also lower asperity level load carrying capacity (see figure 9(a)) because of the absence of any significant solvation pressure. This can potentially lead to increased metal-to-metal interactions. Therefore, an optimal additive blend must consist of a variety of molecular species. The optimal proportion in the blend, should consider surface topography (e.g. roughness, texturing) together with molecular composition of the lubricant. For real lubricants, comprising an assortment of molecular species, the current model needs to be extended to take into account molecular mix as the result of competition near the solid boundaries. This would be affected not only by solvation of the molecules near the boundary solids, but also by the interaction between the various molecules themselves. More work will be required in the future to take these issues into account.

Figure 3 .

 3 Figure 3. Discretized contact characteristics in ultra-thin smooth conjunctions

Figure 4 (

 4 Figure4(a) shows an isobaric plot of conjunctional pressures. This comprises 4 regions of interest: (1) the inlet meniscus/wedge, (2) the contact proper, (3) a cavitation region and (4) reformation zone (as postulated in description of figure1). Each region plays a key role. The inlet meniscus acts as the reservoir supplying the contact with a film of lubricant. Lubricant ahead of the contact provides the meniscus liquid bridge through wetting action of the lubricant on the contiguous surfaces, and is indicated by the inward curved isobar ahead of the contact. Under these circumstances, for simplicity it is assumed that a meniscus is formed at a pressure close to that of lubricant's vaporization pressure p c . With entraining motion of the lubricant through the inlet wedge the generated pressures rise. The contact proper comprises a series of almost concentric isobars with increasing pressures towards the rear lubricant exit constriction from the contact. This is a characteristic of hydrodynamic pressures, except that solvation effect in such diminutive gaps causes pressure fluctuations superimposed upon the hydrodynamic pressure profile. Figure4(b) is a cut through the centre line of the isobaric plot, which shows these pressure fluctuations.The region of high pressures is immediately followed by oil film rupture, not unlike the behaviour of any hydrodynamic conjunction. With lubricant pressures falling below its vaporization pressure beyond the lubricant rupture point, cavitation occurs. The depth of cavitation is governed by the speed of entraining motion in the contact, as well as the gap size. The former is contributed by hydrodynamics, whilst the latter determines the extent of solvation. The fluid in the cavitation region is a mix of lubricant and vapour. Beyond this, in the trailing edge of the contact, the pressures are assumed to return to the atmospheric value. This implies the existence of a film reformation region, where the vapour escapes through the liquid. Figure 4(b) shows the trend towards reformation. The larger the extent of cavitation and reformation regions become, the

Figure 4 .

 4 Figure 4. Pressure variation in an ultra-thin gap of 8 nm

Figure 5 .

 5 Figure 5. Conjunctional conditions at a steadily reduced flying height (Corrected)

Figure 6 .

 6 Figure 6. Conjunctional characteristics with increasing slider speed

Figure 7 .

 7 Figure 7. Liquid film discretization at increasing speed of entraining motion

Figure 8 .Figure 9 .

 89 Figure 8. Friction characteristics of OMCTS

T 1

 1 Temperature ( o C) U Speed of entraining motion (m/s) V Speed of side-leakage (m/s) W Load carrying capacity (N ) Z Pressure-viscosity index (-) a Molecular diameter of lubricant (m) g Switch function (-) h Gap size (m) h 0 Minimum rigid clearance (m) h s Local gap (m) i, k, x Coordinates along x-axis (-) j, l, y Coordianates along y-axis (-) n Iteration counter (-) p T Total pressure (P a) p c Cavitation pressure (P a) p h Hydrodynamic pressure (P a) p s Solvation pressure (P a) p v Constant pressure for varying velocities (P a) p w van der Waals pressure (P a) r Contact curvature radius (m) t Time step (s) v 0 Characteristic velocity (m/s) Ω Lubricant volume where pressure, p T is acting on (m 3 ) α Pressure-viscosity coefficient (P a -1 ) α o Viscosity coefficient at p = 0 and 60 o C (P a -1 ) β Lubricant bulk modulus (P a) δ Elastic deformation (m) θ Convergence criteria for θ (-) η Lubricant viscosity (P a.s) η 0 Lubricant viscosity at p = 0 and 60 o C (P a.s) ν Poisson's ratio (-) ω Relaxation factor for θ (-) φ Activation volume (m 3 ) ρ Lubricant density (kg/m 3 ) ρ c Lubricant density during cavitation (kg/m 3 ) ρ 0 Lubricant density at p = 0 and constant temperature (kg/m 3 ) τ Shear stress (P a) τ Eyring shear stress (Velocity dependent) (P a) τ 0 Eyring shear stress (Pressure dependent) (P a) θ Non-dimensional density or Fractional film content (-) ξ Eyring model constant (-)

Acknowledgement

The authors acknowledge the technical support from partners and sponsorship provided by the EPSRC through the ENCYCLOPAEDIC program grant.

Barrier height for Eyring model (J) E Modulus of elasticity (P a) E Reduced modulus of elasticity (P a) F Friction force (N ) K Boltzman constant (JK -1 ) Q Process activation energy (J)