1 Introduction

Rarefied microscale gas flows are present in a wide variety of applications. In biomedical science, they are used in the study of the transport and inhalation of small droplets which can potentially carry respiratory viruses when they are expelled from an infected person during coughing or sneezing [1]. They are also applied in the design of microelectromechanical system (MEMS) devices where the relevant gas flows are slow rarefied flows [2]. In astrophysics, there are certain bodies in the Solar System whose atmospheres are of extremely low density such that the gases in question must be considered as rarefied gases [3]. It is widely known that microscale dilute gas flows cannot be accurately described with the traditional fluid dynamics framework based on the Navier–Stokes–Fourier closure and that the Navier–Stokes–Fourier equations can fail to capture the relevant phenomena even in a qualitative manner [4]. The study of rarefied gas motion goes back at least as far as the work of Maxwell, who used kinetic theory to study gas flows induced by temperature gradients [5]. There have since been many studies of rarefied gas flows and in particular, there has been a great deal of interest in thermophoresis on a sphere in a steady-state rarefied gas flow (thermophoresis being the force experienced by solid particles or surfaces in a rarefied gas in the presence of a temperature gradient) [6, 7]. Interest in rarefied gas flows in idealised geometries apart from microflows around isolated spherical particles has recently renewed, with examples including the rarefied gas flow between two infinite parallel plates [8, 9] and the drag force on a sphere moving close to a planar boundary through a highly rarefied gas [10]. Modelling of rarefied gas flows is generally complicated even for steady-state flows in simple geometries because the usual continuum NSF equations lose their accuracy as the Knudsen number \(\text {Kn}\) (the ratio of the mean free path in the gas to the characteristic length scale L for the problem) becomes larger.

For Knudsen numbers bigger than 0.01, the NSF equations typically give a picture which is not even qualitatively correct, missing out on the existence of effects such as heat flow from colder regions of the gas to warmer regions [11]. Once the gas flow becomes moderately rarefied, the Boltzmann equation is the only model which gives accurate answers, and must in general be solved numerically. Numerical solution of the Boltzmann equation can be implemented via direct simulation Monte Carlo (DSMC) [12] or deterministic methods [13], but these are usually prohibitively expensive in computational cost for the flows which interest us. DSMC becomes computationally expensive when used to model low-speed microflows where the flow velocity is much smaller than the thermal velocity and is particularly expensive for Knudsen numbers in the so-called transition regime between 0.01 and 1 at low speeds. In the flows in MEMS which interest us in applications, the signal-to-noise ratio is too low to make solution by DSMC computationally viable [13, 14]. Our statements on computational cost depend on the assumption that the interaction between the particles is described by a simple inverse power law potential rather than the more realistic Lennard-Jones potential but this is a standard and widely used assumption in the literature (as is our assumption that the gas is monatomic) [15]. There are also methods which use the smallness of the Knudsen number to greatly increase the efficiency of numerical solution of the Boltzmann equation for moderately small Knudsen numbers [16].

One possibility for effective modelling of rarefied gas flows in the transition regime which also might have some physical insight is to create continuum equations using moment approximations to the Boltzmann equation. The original approximation of this kind was made by Grad who expanded the distribution function in Hermite polynomials and derived 13-moment equations known as the Grad or G13 equations [17, 18]. The Grad equations suffer from a number of deficiencies, including inability to capture the effect of Knudsen boundary layers. Struchtrup and Torrilhon proposed a regularisation of the Grad equations to overcome these issues which are known as the R13 equations [19]. We will not attempt to explain the method of moments here or justify the Grad equations, as there are various detailed reviews in the literature which the reader could consult [20]. Informally, one can think of the Grad equations as an approximation to the Boltzmann equation which is suitable for low-speed microflows at the transitional Knudsen numbers which we are going to study. Moment equations do not currently have tractable methods of numerical solution in three dimensions. Techniques do exist for certain two-dimensional flows, but these are difficult to apply to low-speed external flows [21]. However, given that the flows which interest us are slow creeping flows, one can make significant simplifications by considering the linearised forms of the Grad and R13 equations and in this paper we will only work with the linearised versions. Even in linearised form, solution of the R13 equations is non-trivial and typically involves exact symmetry assumptions [22]. This is the main reason for studying analytic solutions to the Grad equations rather than the R13 equations, and there are already some existing examples in the literature.

We will now summarise the previous work which has been done in this respect. In an important study of thermophoresis on spherical particles, Young examined several models for particle thermophoresis and solved the Grad equations with spherical symmetry to obtain an analytic solution for the drag force on a sphere (Young explicitly states that he has chosen to study the Grad equations because of the complexity of the R13 equations even in their linearised form) [23]. Young also showed that the G13 method generates a hierarchy of expressions for the thermophoretic force on a sphere at low Knudsen numbers which includes all the classic results and that the G13 solution can be used to derive an interpolation formula for the transition regime [23]. The relevance of the G13 method is further clarified in [20]. Before Young, Dwyer attempted to give a theory for the thermophoretic force on a spherical particle based on an analytic solution of the Grad equations [24]. Lockerby and Collyer derived Green functions for the G13 equations, where a Green function is taken to be the flow response to a Dirac delta forcing term. We will discuss these solutions in more detail in Sect. 2 [14]. These solutions are useful from the numerical point of view because one can exploit linearity of the equations and represent a flow field with a superposition of Green functions [25]. The number of analytic solutions to the Grad equations is still somewhat limited given the physical interest of the equations and that these solutions are mostly restricted to the study of thermophoresis on spherical particles or to solutions corresponding to points placed in an infinite dilute gas.

The aim of this work is to study the classical analytic solutions and techniques from potential theory for Stokes flow and investigate the extent to which these techniques can be carried over to slow, steady-state rarefied gas flows. In the process, we hope to derive new exact and integral solutions which can be employed in numerical schemes or used to model steady-state dilute gas flows in certain geometries which are idealised but which at least include some kind of boundary. Given that the equations of motion and the boundary conditions for the analogous boundary value problems are more complicated, it would be natural to expect that these techniques will not carry over entirely and that it might not be possible to derive the analogue for certain analytic solutions of Stokes flow. This turns out to be the case, as we will see shortly.

The structure of the paper is as follows. In Sect. 2, we will fix notation and outline the methods, equations, boundary conditions and solutions which we will use in the rest of the paper. In Sects. 3 and 4, we will consider some applications of the method of images. The basic idea of this method in the context of fluid dynamics is to consider the flow fields generated by a point forcing or a point heat source close to a boundary surface in terms of an ‘image system’ of singularities on the other side of the wall. The related fields can then be found by using Fourier analysis to obtain the image system which satisfies the required boundary conditions. This method has not before been applied to dilute microscale gas flows as far as we know. Using the method of images we will derive the velocity and temperature fields due to point forcings and point heat sources in the presence of and without velocity slip and temperature jump at the boundary. The purpose of this is two-fold: to obtain a better physical understanding of velocity (temperature) fields due to a point forcing (point heat source) in dilute microscale gas flows and heat flows close to a plane boundary, and to derive novel exact solutions which can be implemented in the study of such flows in simple geometries and problems. This could potentially enable us to better understand the movement of slender bodies or small aspherical particulates in rarefied microscale gas flows close to walls. The solutions could also be employed in certain numerical schemes [14, 25].

In Sect. 5, we derive a novel solution for the unsteady Grad equations given a time-dependent point heat source. The motivation for this is to try to extend our understanding of dilute microscale gas flows slightly beyond the usual steady-state case and to obtain a solution which complements the existing time-dependent fundamental solutions obtained for unsteady Stokes flow [26]. This solution could be implemented in a numerical method which uses time-dependent fundamental solutions and used to model an unsteady heat flow in a very simple geometry [27]. Our assumption of one dimension seems prohibitive, but is quite typical in the existing literature [28]. Many of the solutions which we derive refer to a flow with constant velocity and pressure. This is obviously quite a restrictive condition, but does correspond to some physically realistic scenarios of current interest (for example, Fourier flow, where one has a static gas between two stationary parallel plates and the flow structure as it exists is due only to the temperature difference between the plates) [29]. We will refer to these types of flow as heat flows to emphasise that the medium is not moving.

Finally, in Sect. 6, we derive the analogue of the ‘rotlet’ from studies of Stokes flow (the rotlet being a point torque rather than a point force). As well as extending the result from Stokes flows to rarefied microscale gas flows for the sake of interest, this result has potential applications in biology and biomedical science since micro-organisms can exhibit helical beating motions as they move close to flat walls, in which case angular momentum needs to be accounted for [30]. To use the rotlet which we derive in applications, one would need to place it close to a no-slip boundary and determine the associated velocity and pressure fields. One might wonder if this could be made more realistic by placing the rotlet close to a partial-slip boundary, but we show in Sect. 4 that obtaining the relevant fields even for the simpler case of a point forcing in a rarefied gas flow close to a partial-slip plane boundary is unfortunately intractable.

2 Preliminaries

For the sake of clarity, we will now fix all the notation and equations which are used in the rest of the paper so that the reader can refer back to them as necessary. The basic geometry which we will consider for all the results in Sects. 3 and 4 is shown in Fig. 1. The domain on which our PDEs are defined is the upper half-space given by the set of all triples \(S = \{ (x,y,z) : z > 0\}\) and the boundary of the domain (denoted in grey in Fig. 1) is the flat infinite plane given by all triples \(\{ (x,y,z) : z = 0\}\). At a height \(z=h\) above the plane, one then places either a point source of force \(\mathbf{f } \delta (\mathbf{r })\) or a point source of heat \(g \delta ( \mathbf{r })\) at the point \((x,y,z) = (0,0,h)\), where \(\mathbf{f }\) is the point force vector, g is the strength of the heat source, \(\delta \) denotes the Dirac delta function and \(\mathbf{r } = (x,y,z-h)\) is the vector from the source to the point of observation. On the opposite side of the wall outside of the domain on which the PDEs are defined, one imagines that there is a mirror image point source of force \(\overline{\mathbf{f }} \delta ( \overline{\mathbf{r }})\) or a point source of heat \(g \delta (\overline{\mathbf{r }})\), where \(\overline{\mathbf{f }}\) is the associated point force vector chosen to satisfy the no-flux boundary condition, and \(\overline{\mathbf{r }} = (x,y,z+h)\) is the position vector for the point at which the image source is observed. The idea of the method of images is to decompose the velocity field \(\mathbf{v }\) into three parts

$$\begin{aligned} \mathbf{v } = \mathbf{u } + \mathbf{u }' + \mathbf{u }'', \end{aligned}$$
(1)
Fig. 1
figure 1

The point force or heat source is located at a height h above the boundary at \(z=0\) and there is a mirror image point force or heat source whose strength is equal in magnitude at a height h below the boundary outside the domain in which the PDEs are defined. Note that the point source of force can be perpendicular to the plane at \(z=0\), or parallel to it

Full size image

where \(\mathbf{u }\) is the velocity field due to the point forcing placed in the flow, \(\mathbf{u }'\) is the velocity field due to the mirror image point forcing placed on the other side of the wall and \(\mathbf{u }''\) is the unknown velocity field due to some system of singularities on the other side of the wall. The basic reason for needing more singularities in the image system apart from the mirror image point force is that velocity is a vector quantity and not a scalar, hence an image point force by itself would not be able to cancel those components of the flow which are both tangential and normal to the boundary surface. Note that the prime does not denote differentiation and we hope that this is not misleading, as we feel that it is more confusing to denote the different parts of the total velocity field \(\mathbf{v }\) with different letters. \(\mathbf{v }\) as a whole is required to solve the boundary value problem, so one uses Fourier analysis to find the unknown velocity field \(\mathbf{u }''\) and finally adds the three parts together to obtain the total flow field due to a point force or heat source close to the wall. The analysis in question requires two-dimensional Fourier transforms. For reference, the Fourier transform and the inverse Fourier transform are defined via

$$\begin{aligned} \hat{{u}}_i ( k_1, k_2, z)= & {} \frac{1}{2 \pi } \int \int u_i (x,y,z) \mathrm{e}^{\mathrm{i} k_1 x + \mathrm{i} k_2 y} \, \text {d}x \, \text {d}y, \end{aligned}$$
(2a)
$$\begin{aligned} u_i ( x, y, z)= & {} \frac{1}{2 \pi } \int \int \hat{{u}}_i ( k_1, k_2, z) \mathrm{e}^{-\mathrm{i} k_1 x - \mathrm{i} k_2 y} \, \text {d}k_1 \, \text {d}k_2. \end{aligned}$$
(2b)

Useful tables of Fourier transforms can be found at [31, 32]. Quantities with hats indicate Fourier transformed quantities. In Fig. 1, we have mentioned that one must consider separately the cases where the point forcing is perpendicular or parallel to the wall: this is because the presence of the wall breaks isotropy of particles moving close to it. A superscript \(\perp \) for a quantity indicates that we are considering the case where the point forcing is perpendicular to the boundary (in the z-direction), and a superscript \(\parallel \) for a quantity indicates that we are considering the case where the point forcing is parallel to the boundary (in the x-direction). To be clear, this has nothing to do with the actual components of the velocity field itself. By the no-flux boundary condition, the z-component of the velocity field vanishes at the boundary whether the point forcing is in the x- or the z-direction.

2.1 Stokes equations with point forcing

We now outline the PDEs, boundary conditions and solutions which we will be referring to throughout the paper. To begin with, the Stokes equations are given by

$$\begin{aligned}&\nabla \cdot \mathbf{v } = 0 , \end{aligned}$$
(3a)
$$\begin{aligned}&\nabla p + \nabla \cdot \mathbf{S } = 0 , \end{aligned}$$
(3b)
$$\begin{aligned}&\mathbf{S } = - 2\, \text {Kn}' \, \overline{\nabla \mathbf{v }} , \end{aligned}$$
(3c)

where \(\text {Kn}'\) is a rescaled version of the Knudsen number which is relevant for the gas flows which we are studying

$$\begin{aligned} \text {Kn}' = \sqrt{\frac{2}{\pi }} \text {Kn} \end{aligned}$$
(4)

and the overline denotes the symmetric trace-free part of a tensor. The equations which we consider are dimensionless and linearised about an equilibrium state which is defined by a constant reference density \(\rho _0\) and a reference temperature \(\theta _0\). We only consider small deviations away from equilibrium in this work. The relations between the variables for dimensionless deviation away from the equilibrium state (denoted with tilde symbols) and variables with dimension are

$$\begin{aligned} \theta= & {} \theta _0 (1 + {\tilde{\theta }}), \,\,\,\, \,\,\,\, \rho = \rho _0 ( 1 + {\tilde{\rho }}), \,\,\,\, \,\,\,\, p = p_0 (1 + {\tilde{p}} ), \,\,\,\, \,\,\,\, \mathbf{r } = L \tilde{\mathbf{r }}, \end{aligned}$$
(5a)
$$\begin{aligned} \mathbf{v }= & {} \sqrt{R \theta _0} \tilde{\mathbf{v }}, \,\,\,\, \,\,\,\, \mathbf{q } = \rho _0 (R \theta _0)^{3/2} \tilde{\mathbf{q }} ,\,\,\,\, \,\,\,\, \mathbf{S } = \rho _0 R \theta _0 \tilde{\mathbf{S }}, \end{aligned}$$
(5b)

where L is a characteristic length scale, \(p_0\) is the equilibrium pressure and R is the ideal gas constant. In the rest of this work, tilde symbols will not be shown and all variables will be dimensionless unless stated otherwise.

The Green function for the Stokes equations with a point forcing term on the right-hand side of the momentum equation is extremely well studied in the literature (recall that the point forcing for the Stokes equations is called the Stokeslet) [33]. The reader should bear in mind that with a point forcing term applied, the momentum equation (3b) is actually

$$\begin{aligned} \nabla p + \nabla \cdot \mathbf{S } = \mathbf{f } \delta ( \mathbf{r } ). \end{aligned}$$
(6)

Assuming that the pressure and the magnitude of the velocity and stress tensors vanish as one goes to infinity, it is well known that one obtains [14, 32]

$$\begin{aligned} \mathbf{v } = \mathbf{f } \cdot \mathbf{G }, \,\, \, \, \,\, \, \, p = \frac{\mathbf{f } \cdot \mathbf{r }}{4 \pi | \mathbf{r }|^3} , \,\, \, \, \,\, \, \, \mathbf{S } = - \frac{\mathbf{f } \cdot \mathbf{r }}{4 \pi } \Bigg ( \frac{\mathbf{I }}{|\mathbf{r }|^3} - \frac{3 \mathbf{r } \mathbf{r }}{| \mathbf{r } |^5} \Bigg ), \end{aligned}$$
(7)

where \(\mathbf{I }\) is the identity matrix and the Green function \(\mathbf{G }\) is given by

$$\begin{aligned} \mathbf{G } = \frac{1}{8 \pi \text {Kn}'} \Bigg ( \frac{\mathbf{I }}{|\mathbf{r }|} + \frac{\mathbf{r } \mathbf{r }}{|\mathbf{r }|^3} \Bigg ) . \end{aligned}$$
(8)

Note that this is the free-space solution and that we have not yet introduced a boundary condition at \(z=0\). Using the parallel and perpendicular notation which we mentioned above, one has for parallel and perpendicular Stokeslets

$$\begin{aligned} \mathbf{v }^{\perp }= & {} f \mathbf{G }_{x} , \end{aligned}$$
(9a)
$$\begin{aligned} \mathbf{v }^{\parallel }= & {} f \mathbf{G }_{z} , \end{aligned}$$
(9b)

where f is the magnitude of the force, so one takes the component of the Green function which corresponds to the direction of the point forcing. For the velocity field due to the mirror image point forcing, the sign in front of the force is reversed if one takes the force in the perpendicular direction

$$\begin{aligned} \mathbf{v }'^{\perp }= & {} -f \mathbf{G }_{x} , \end{aligned}$$
(10a)
$$\begin{aligned} \mathbf{v }'^{\parallel }= & {} f \mathbf{G }_{z} , \end{aligned}$$
(10b)

which often leads to some cancellations in the computations. The sign of the force is chosen in this way to enforce the no-flux condition at the boundary. The no-slip boundary condition is simply the requirement that every component of the velocity vanishes at the boundary

$$\begin{aligned} \mathbf{v } = 0 \end{aligned}$$
(11)

and for a stationary wall the slip boundary condition in terms of the stress tensor \(\mathbf{S }\) is [14]

$$\begin{aligned} \mathbf{v } = - \sqrt{\frac{\pi }{2}} \mathbf{n } \cdot \mathbf{S } \cdot (\mathbf{I } - \mathbf{n } \mathbf{n } ), \end{aligned}$$
(12)

where \(\mathbf{n }\) is the unit normal vector to the boundary surface. In our case, the boundary surface is a flat plane at \(z=0\) and the unit normal vector is the unit vector in the z-direction, so the slip boundary condition in terms of each component of the velocity becomes

$$\begin{aligned} v_1= & {} - \sqrt{\frac{\pi }{2}} S_{31} , \end{aligned}$$
(13a)
$$\begin{aligned} v_2= & {} - \sqrt{\frac{\pi }{2}} S_{32} , \end{aligned}$$
(13b)
$$\begin{aligned} v_3= & {} 0 . \end{aligned}$$
(13c)

2.2 Stokes equations with point heat source

One can instead look for solutions when the conservation equations are subject to a steady-state point heat source rather than a point force source (in [14], this point heat source is called a thermal Stokeslet). In terms of the stress tensor and heat flux, the conservation equations are

$$\begin{aligned}&\nabla \cdot \mathbf{v } = 0, \end{aligned}$$
(14a)
$$\begin{aligned}&\nabla p + \nabla \cdot \mathbf{S } = 0 , \end{aligned}$$
(14b)
$$\begin{aligned}&\nabla \cdot \mathbf{q } = g \delta ( \mathbf{r }) , \end{aligned}$$
(14c)

where \(\mathbf{q }\) is the heat flux and g is the strength of the point heat source. The equations for the stress tensor and heat flux are given by the NSF closure:

$$\begin{aligned} \mathbf{S }= & {} - 2\, \text {Kn} ' \overline{\nabla \mathbf{v }} , \end{aligned}$$
(15a)
$$\begin{aligned} \mathbf{q }= & {} - \frac{15}{4} \text {Kn}' \nabla \theta . \end{aligned}$$
(15b)

In the literature, one then typically looks for isobaric zero flow solutions such that \(\mathbf{v } = p = 0\) which gives us the equations

$$\begin{aligned} \nabla \cdot \mathbf{q }= & {} g \delta ( \mathbf{r }) , \end{aligned}$$
(16a)
$$\begin{aligned} \mathbf{q }= & {} - \frac{15}{4} \text {Kn}' \nabla \theta . \end{aligned}$$
(16b)

Taking the boundary condition to be that the temperature \(\theta \) scaled by the mass goes to its equilibrium value at infinity, one then finds that for the thermal Stokeslet applied to an infinite rarefied gas, the temperature \(\theta \) and the heat flux \(\mathbf{q }\) are

$$\begin{aligned} \theta= & {} \frac{g}{15\, \text {Kn}' \pi |\mathbf{r }|}, \end{aligned}$$
(17a)
$$\begin{aligned} \mathbf{q }= & {} \frac{g \mathbf{r }}{4 \pi |\mathbf{r }|^3} . \end{aligned}$$
(17b)

The no-temperature jump boundary condition is the requirement that the temperature vanishes at the boundary

$$\begin{aligned} \theta = 0 \end{aligned}$$
(18)

and the temperature jump boundary condition is

$$\begin{aligned} \theta = - \frac{1}{2} \sqrt{\frac{\pi }{2}} \mathbf{n } \cdot \mathbf{q } . \end{aligned}$$
(19)

Since the unit normal vector is the unit vector in the z-direction, the temperature jump boundary condition becomes

$$\begin{aligned} \theta = - \frac{1}{2} \sqrt{\frac{\pi }{2}} q_3 . \end{aligned}$$
(20)

2.3 Grad equations with point forcing

In this work, we are primarily interested in the Grad equations for rarefied gas flows. For a point forcing applied to the momentum equation, these consist of the usual equations for conservation of mass, momentum and internal energy

$$\begin{aligned} \nabla \cdot \mathbf{v }= & {} 0 , \end{aligned}$$
(21a)
$$\begin{aligned} \nabla p + \nabla \cdot \mathbf{S }= & {} \mathbf{f } \delta (\mathbf{r }) , \end{aligned}$$
(21b)
$$\begin{aligned} \nabla \cdot \mathbf{q }= & {} 0, \end{aligned}$$
(21c)

combined with equations for the stress tensor and the heat flux given by Grad’s closure

$$\begin{aligned} \mathbf{S }= & {} - 2\, \text {Kn}' \overline{\nabla \mathbf{v }} - \frac{4}{5} \text {Kn}' \overline{\nabla {\mathbf{q }}} , \end{aligned}$$
(22a)
$$\begin{aligned} \mathbf{q }= & {} - \frac{15}{4} \text {Kn}' \nabla \theta - \frac{3}{2} \text {Kn}' \nabla \cdot \mathbf{S } . \end{aligned}$$
(22b)

The point forcing in the case of the Grad equations is known as the Gradlet [14]. Solving these equations for an infinite rarefied gas, one finds that the velocity and pressure fields associated to a Gradlet in an infinite domain far away from boundaries are

$$\begin{aligned} \mathbf{v }= & {} \frac{\mathbf{f }}{8 \pi \text {Kn}'} \Bigg ( \frac{\mathbf{I }}{| \mathbf{r }|} + \frac{ \mathbf{r } \mathbf{r }}{|\mathbf{r }|^3} \Bigg ) - \frac{3 \text {Kn}' \mathbf{f }}{20 \pi } \Bigg ( \frac{\mathbf{I }}{| \mathbf{r }|^3} - \frac{3 \mathbf{r } \mathbf{r }}{|\mathbf{r }|^5} \Bigg ), \end{aligned}$$
(23a)
$$\begin{aligned} p= & {} \frac{\mathbf{f } \cdot \mathbf{r }}{4 \pi |\mathbf{r }|^3} . \end{aligned}$$
(23b)

The slip boundary condition for the Grad equations is more complicated than the slip boundary condition for the Stokes equations and includes a term which depends on the heat flux [34]:

$$\begin{aligned} \mathbf{v } = - \sqrt{\frac{\pi }{2}} \mathbf{n } \cdot \mathbf{S } \cdot (\mathbf{I } - \mathbf{n } \mathbf{n }) - \frac{1}{5} \mathbf{q } \cdot (\mathbf{I } - \mathbf{n } \mathbf{n }). \end{aligned}$$
(24)

For a flat plane as the boundary, these become

$$\begin{aligned} v_1= & {} - \sqrt{\frac{\pi }{2}}S_{3 1} - \frac{1}{5}q_{1} , \end{aligned}$$
(25a)
$$\begin{aligned} v_2= & {} - \sqrt{\frac{\pi }{2}}S_{3 2} - \frac{1}{5}q_{2} , \end{aligned}$$
(25b)
$$\begin{aligned} v_3= & {} 0 . \end{aligned}$$
(25c)

2.4 Grad equations with point heat source

As with the Stokes equations, one can instead consider the Grad equations with a point heat source applied to the conservation equations:

$$\begin{aligned}&\nabla \cdot \mathbf{v } = 0 , \end{aligned}$$
(26a)
$$\begin{aligned}&\nabla p + \nabla \cdot \mathbf{S } = 0, \end{aligned}$$
(26b)
$$\begin{aligned}&\nabla \cdot \mathbf{q } = g \delta (\mathbf{r }), \end{aligned}$$
(26c)
$$\begin{aligned}&\mathbf{S } = - 2\, \text {Kn}' \overline{\nabla \mathbf{v }} - \frac{4}{5} \text {Kn}' \overline{\nabla {\mathbf{q }}} , \end{aligned}$$
(26d)
$$\begin{aligned}&\mathbf{q } = - \frac{15}{4} \text {Kn}' \nabla \theta - \frac{3}{2} \text {Kn}' \nabla \cdot \mathbf{S } . \end{aligned}$$
(26e)

In this case, the point heat source is known as a thermal Gradlet. As before, it is usual to simplify by taking \(\mathbf{v } =p =0\) which gives us

$$\begin{aligned} \nabla \cdot \mathbf{S }= & {} 0, \end{aligned}$$
(27a)
$$\begin{aligned} \nabla \cdot \mathbf{q }= & {} g \delta (\mathbf{r }), \end{aligned}$$
(27b)
$$\begin{aligned} \mathbf{S }= & {} -\frac{4}{5} \text {Kn}' \overline{\nabla \mathbf{q }}, \end{aligned}$$
(27c)
$$\begin{aligned} \mathbf{q }= & {} - \frac{15}{4} \text {Kn}' \nabla \theta . \end{aligned}$$
(27d)

One then can solve to find the heat flux \(\mathbf{q }\), the stress tensor \(\mathbf{S }\) and the temperature \(\theta \) rescaled by the mass for a point heat source in an infinite rarefied stationary gas [14]:

$$\begin{aligned} \theta= & {} \frac{g}{15\, \text {Kn}' \pi |\mathbf{r }|}, \end{aligned}$$
(28a)
$$\begin{aligned} \mathbf{q }= & {} \frac{g \mathbf{r }}{4 \pi |\mathbf{r }|^3} , \end{aligned}$$
(28b)
$$\begin{aligned} \mathbf{S }= & {} - \frac{\text {Kn}' g}{5 \pi } \Bigg ( \frac{\mathbf{I }}{| \mathbf{r }|^3} - \frac{3 \mathbf{r } \mathbf{r }}{| \mathbf{r }|^5} \Bigg ) . \end{aligned}$$
(28c)

Again, as before, there is some additional complexity in the temperature jump boundary condition for the Grad equations which now involves a term which depends on the stress tensor [34]:

$$\begin{aligned} \theta = - \frac{1}{2} \sqrt{\frac{\pi }{2}} \mathbf{n } \cdot \mathbf{q } - \frac{1}{4} \mathbf{n } \cdot \mathbf{S } \cdot \mathbf{n } . \end{aligned}$$
(29)

For a flat surface, this simplifies to

$$\begin{aligned} \theta = -\frac{1}{2} \sqrt{\frac{\pi }{2}} q_3 - \frac{1}{4} S_{33} . \end{aligned}$$
(30)

The coupling between the stress tensor and the heat flux is an effect which is present in rarefied gas flows but not in standard Stokes flows.

3 Exact solutions from the method of images

The basic question is to what extent the fields due to a point force or a point heat source placed in an infinite rarefied gas are modified in the presence of a flat planar boundary, and this is what we will now proceed to study in some detail via the method of images. Along the way, we will also derive some solutions for Stokes flows, some of which are new as far as we are aware.

3.1 Solutions for no-slip and no-temperature jump boundaries

We will begin by focussing on the no-slip and no-temperature jump boundary conditions given by (11) and (18). The situation is simplest for a boundary value problem where the temperature is prescribed to vanish at the boundary. In this setting, the image system required to solve the boundary value problem is just the mirror image of whatever was placed in the domain on the other side of the wall.

As an example, we start with the temperature and heat flux fields associated to a point heat source placed in isobaric stationary-state Stokes flow close to a no-temperature jump surface (Eqs. (14), solutions (17) and boundary condition (18)). As emphasised before, this type of flow would be better characterised as a heat flow, since the medium is not moving. We will abuse notation slightly and decompose the total temperature field as \(\theta + \theta ' + \theta ''\), where \(\theta \) denotes the temperature field for the point heat source placed in the domain close to the boundary, and does not denote the total temperature field unless explicitly stated otherwise (this should only be at the end of a section when we are ready to write down total fields). Hopefully this is not too confusing, as it saves on having to introduce extra letters in the notation. If the total temperature field satisfies the no-temperature jump boundary condition, then we must have

$$\begin{aligned} \theta + \theta ' + \theta '' = 0 . \end{aligned}$$
(31)

After taking Fourier transforms of each quantity, we have

$$\begin{aligned} {\hat{\theta }} +{\hat{\theta }}' + {\hat{\theta }}'' = 0. \end{aligned}$$
(32)

\(\theta \) is given by (17a), so taking the Fourier transform we have

$$\begin{aligned} {\hat{\theta }} = \frac{g}{15\, \text {Kn}' \pi } \frac{1}{k} \mathrm{e}^{-kz} , \end{aligned}$$
(33)

where \(k=\sqrt{k_1^2 + k_2^2}\). The next point also has some potential for confusion, so we will clarify. The unknown temperature field \(\theta ''\) and unknown heat flux \(\mathbf{q }''\) by themselves satisfy Eqs. (14), but without the point heat source term on the right-hand side of the conservation equation:

$$\begin{aligned} \nabla \cdot \mathbf{q }= & {} 0 , \end{aligned}$$
(34a)
$$\begin{aligned} \mathbf{q }= & {} - \frac{15}{4} \text {Kn}' \nabla \theta . \end{aligned}$$
(34b)

This needs to be made clear as it may be used implicitly in the rest of the text. It should be clear from a physical point of view because the unknown fields are associated with a fictitious system of singularities on the other side of the wall outside the physical domain. To find \({\hat{\theta }}''\), one takes two-dimensional Fourier transforms of the above equations and solves the resulting differential equations. The above equations simplify to

$$\begin{aligned} \Delta \theta '' = 0 . \end{aligned}$$
(35)

The two-dimensional Fourier transform of this is

$$\begin{aligned} \frac{\partial ^2 {\hat{\theta }}''}{\partial z^2} - k^2 \theta ''=0 , \end{aligned}$$
(36)

where \(k^2 = k_1^2 + k_2^2\). The solution to this differential equation is

$$\begin{aligned} {\hat{\theta }}'' = A \mathrm{e}^{-kz} , \end{aligned}$$
(37)

where A is a constant to be determined and the positive exponential term has been neglected because the temperature goes to zero at spatial infinity. Substituting all the hatted quantities into the Fourier transformed boundary condition at \(z=0\), we have

$$\begin{aligned} A = - \frac{2g}{15\, \text {Kn}' \pi } \frac{1}{k}, \end{aligned}$$
(38)

so

$$\begin{aligned} {\hat{\theta }} '' = - \frac{2g}{15 \text {Kn}' \pi } \frac{1}{k} \mathrm{e}^{-kz} . \end{aligned}$$
(39)

From Eq. (14c), we have

$$\begin{aligned} q_1'' = - \frac{15}{4} \text {Kn}' \frac{\partial \theta }{\partial x},\quad q_2'' = - \frac{15}{4} \text {Kn}' \frac{\partial \theta }{\partial y} ,\quad q_3'' = - \frac{15}{4} \text {Kn}'\frac{\partial \theta }{\partial z}. \end{aligned}$$
(40)

Taking Fourier transforms as usual, we have

$$\begin{aligned} {\hat{q}}''_1= & {} \frac{15 \mathrm{i} k_1 }{4} \text {Kn}' {\hat{\theta }} = - \frac{g}{2 \pi } \frac{\mathrm{i} k_1}{k} \mathrm{e}^{-kz}, \end{aligned}$$
(41a)
$$\begin{aligned} {\hat{q}}''_2= & {} \frac{15 \mathrm{i} k_2 }{4} \text {Kn}' {\hat{\theta }} = - \frac{g}{2 \pi } \frac{\mathrm{i} k_2}{k} \mathrm{e}^{-kz}, \end{aligned}$$
(41b)
$$\begin{aligned} {\hat{q}}''_3= & {} \frac{15}{4}k \,\text {Kn}' A \mathrm{e}^{-kz} = - \frac{g}{2 \pi } \mathrm{e}^{-kz} . \end{aligned}$$
(41c)

We finish by taking the inverse Fourier transforms to arrive at

$$\begin{aligned} \theta ''= & {} - \frac{2g}{15\,\, \text {Kn}' \pi } \frac{1}{|\mathbf{r }|} , \end{aligned}$$
(42a)
$$\begin{aligned} \mathbf{q }''= & {} -\frac{g \mathbf{r }}{2 \pi } \frac{1}{|\mathbf{r }|^3}, \end{aligned}$$
(42b)

and adding the unknown fields to the fields due to the point heat source in the domain and the mirror image heat source on the other side to obtain the expected total fields for a heat source close to a no-temperature jump surface

$$\begin{aligned} \theta= & {} \frac{g}{15\, \text {Kn}' \pi } \Bigg ( \frac{1}{|\overline{\mathbf{r }}|} - \frac{1}{|\mathbf{r }|} \Bigg ), \end{aligned}$$
(43a)
$$\begin{aligned} \mathbf{q }= & {} \frac{g \mathbf{r }}{4 \pi } \Bigg ( \frac{1}{|\overline{\mathbf{r }}|^3} - \frac{1}{|\mathbf{r }|^3} \Bigg ). \end{aligned}$$
(43b)

We have been somewhat explicit here, but hope that the method of calculation should be clear in the rest of the paper even if some steps are missed out for brevity. Similar computations confirm that the fields due to a point heat source placed in a isobaric stationary-state rarefied gas flow close to a no-temperature jump surface (Eqs. (27), solutions (28) and boundary condition (18)) are equivalent to the ones obtained in (43a) and (43b).

3.2 Velocity distribution for a gradlet close to a no-slip surface

The next problem is to determine the fields due to a point forcing (known in this context as a Gradlet) placed in a rarefied gas flow close to a no-slip boundary (Eqs. (22), solutions (23), and boundary conditions (24)). As mentioned earlier, we do have to consider as separate cases when the point forcing is perpendicular to the wall or parallel to it, because the wall breaks isotropy of nearby particle motions. The Grad equations are generally combined with a slip boundary condition. The reason for this is that gas flows in microscale or nanoscale devices whose dimensions are of the order of the mean free path of the gas molecules typically exhibit a significant amount of gas slip. However, it is possible to have certain rarefied gas flows where the effect of gas slip can be approximately neglected [35, 36] and in this case, the exact solution which we are about to derive would be very useful and applicable in mesh-free numerical methods [37].

Starting with Eqs. (21) and removing the point forcing term applied to the momentum equation, it can be shown with some substitutions [23] that

$$\begin{aligned}&\nabla ^2 p'' = 0, \end{aligned}$$
(44a)
$$\begin{aligned}&\nabla p'' = \text {Kn}' \nabla ^2 \mathbf{u }'', \end{aligned}$$
(44b)
$$\begin{aligned}&\nabla ^2 \mathbf{q }''=0 . \end{aligned}$$
(44c)

Taking Fourier transforms and solving the resulting differential equations, we get

$$\begin{aligned} {\hat{u}}_1''= & {} \Bigg ( \frac{h}{8 \pi \text {Kn}'} \Bigg ) (B_1 + \mathrm{i} k_1 A z ) \mathrm{e}^{-kz}, \end{aligned}$$
(45a)
$$\begin{aligned} {\hat{u}}_2''= & {} \Bigg ( \frac{h}{8 \pi \text {Kn}'} \Bigg ) (B_2 + \mathrm{i} k_2 A z ) \mathrm{e}^{-kz}, \end{aligned}$$
(45b)
$$\begin{aligned} {\hat{u}}_3''= & {} \Bigg ( \frac{h}{8 \pi \text {Kn}'} \Bigg ) (B_3 + k A z ) \mathrm{e}^{-kz}, \end{aligned}$$
(45c)
$$\begin{aligned} {\hat{p}}= & {} \Bigg (\frac{h}{4 \pi } \Bigg ) A \mathrm{e}^{-kz} , \end{aligned}$$
(45d)

where h is the height at which the Gradlet has been placed above the boundary. The goal is now to use the boundary conditions and the continuity equation (21a) to determine the constants A and \(B_i\). We will start with the case where the Gradlet is perpendicular to the surface. The no-slip boundary condition is

$$\begin{aligned} \mathbf{u }^{\perp } + \mathbf{u }^{\perp '} + \mathbf{u }^{\perp ''} = 0 . \end{aligned}$$
(46)

The meaning of the primes in the decomposition of the total velocity field \(\mathbf{v }\) into three parts is the same as in the previous section and was also explained in the preliminary section. Next, add the velocity fields (23) for the Gradlet and the mirror image Gradlet when the Gradlet is in the perpendicular direction (i.e. one takes the z-component of the force). If one took the point forcing to be in the parallel direction, the expressions would be slightly different as one would take the x-component of the force. This leaves us with

$$\begin{aligned} u_1^{\perp } + u_1'^{\perp }= & {} - \frac{1}{4 \pi \text {Kn}'} \frac{xh}{r_h^3} - \frac{9 \text {Kn}'}{10 \pi } \frac{x h}{r_h^5} , \end{aligned}$$
(47a)
$$\begin{aligned} u_2^{\perp } + u_2'^{\perp }= & {} - \frac{1}{4 \pi \text {Kn}'} \frac{yh}{r_h^3} - \frac{9 \text {Kn}'}{10 \pi } \frac{y h}{r_h^5} , \end{aligned}$$
(47b)
$$\begin{aligned} u_3^{\perp }= & {} 0 , \end{aligned}$$
(47c)

where \(r_h^2 = x^2 + y^2 + h^2\) and h is always the height at which the point force or point heat source is placed above the boundary. We Fourier transform this to get

$$\begin{aligned} {\hat{u}}_1^{\perp } + {\hat{u}}_1^{\prime \perp }= & {} - \mathrm{i} \frac{h}{4 \pi \text {Kn}'} \frac{k_1}{k} \mathrm{e}^{-kh} - \mathrm{i} \frac{3 \text {Kn}'}{10 \pi } k_1 \mathrm{e}^{-kh} , \end{aligned}$$
(48a)
$$\begin{aligned} {\hat{u}}_2^{\perp } + {\hat{u}}_2^{\prime \perp }= & {} - \mathrm{i} \frac{h}{4 \pi \text {Kn}'} \frac{k_2}{k} \mathrm{e}^{-kh} - \mathrm{i} \frac{3 \text {Kn}'}{10 \pi } k_2 \mathrm{e}^{-kh} , \end{aligned}$$
(48b)
$$\begin{aligned} {\hat{u}}_3^{\perp } + {\hat{u}}_3^{\prime \perp }= & {} 0 . \end{aligned}$$
(48c)

Finally, we add all the hatted quantities at the boundary at \(z=0\) to obtain expressions for the coefficients \(B_i\):

$$\begin{aligned} \frac{h}{8 \pi \text {Kn}'} B^{\perp }_1= & {} - \mathrm{i} \frac{h}{4 \pi \text {Kn}'} \frac{k_1}{k} - \mathrm{i} \frac{3 \text {Kn}'}{10 \pi } k_1 , \end{aligned}$$
(49a)
$$\begin{aligned} \frac{h}{8 \pi \text {Kn}'} B^{\perp }_2= & {} - \mathrm{i} \frac{h}{4 \pi \text {Kn}'} \frac{k_2}{k} - \mathrm{i} \frac{3 \text {Kn}'}{10 \pi } k_2 \end{aligned}$$
(49b)
$$\begin{aligned} \frac{h}{8 \pi \text {Kn}'} B^{\perp }_3= & {} 0 . \end{aligned}$$
(49c)

Equation (49c) implies that \(B_3^{\perp } =0\), so we only need a third equation to close the system and determine \(A^{\perp }\). This is provided by the continuity equation (21a). Fourier transforming the continuity equation and making some cancellations, we end up with

$$\begin{aligned} A^{\perp } = \frac{\mathrm{i}}{k}(k_1 B^{\perp }_1 + k_2 B^{\perp }_2 ). \end{aligned}$$
(50)

Substituting these constants back into Eqs. (49) and taking inverse Fourier transforms, we obtain the unknown pressure and velocity fields which can be added to the fields for the Gradlet and the mirror image Gradlet to obtain the total fields associated with a point forcing placed in a rarefied gas flow perpendicular to the boundary.

$$\begin{aligned} u_1^{\perp ''}= & {} -\frac{h}{4 \pi \text {Kn}'} \Bigg ( \frac{\partial }{\partial z} \Bigg ( \frac{h x}{|\mathbf{r }|^3} \Bigg ) - \frac{\partial }{\partial z} \Bigg ( \frac{ z x}{|\mathbf{r }|^3} \Bigg ) \Bigg ) + \frac{3\, \text {Kn}'}{10 \pi } \Bigg ( \frac{\partial }{\partial z} \Bigg ( \frac{3 z x}{|\mathbf{r }|^5} \Bigg ) - \frac{\partial }{\partial z} \Bigg ( \frac{ 3 h x}{|\mathbf{r }|^5} \Bigg ) \Bigg ) , \end{aligned}$$
(51a)
$$\begin{aligned} u_2^{\perp ''}= & {} -\frac{h}{4 \pi \text {Kn}'} \Bigg ( \frac{\partial }{\partial z} \Bigg ( \frac{h y}{|\mathbf{r }|^3} \Bigg ) - \frac{\partial }{\partial z} \Bigg ( \frac{ z y}{|\mathbf{r }|^3} \Bigg ) \Bigg ) + \frac{3\, \text {Kn}'}{10 \pi } \Bigg ( \frac{\partial }{\partial z} \Bigg ( \frac{3 z y}{|\mathbf{r }|^5} \Bigg ) - \frac{\partial }{\partial z} \Bigg ( \frac{ 3 h y}{|\mathbf{r }|^5} \Bigg ) \Bigg ) , \end{aligned}$$
(51b)
$$\begin{aligned} u_3^{\perp ''}= & {} -\frac{h}{4 \pi \text {Kn}'} \Bigg ( \frac{\partial }{\partial z} \Bigg ( \frac{hz}{|\mathbf{r }|^3} - \frac{1}{|\mathbf{r }|} - \frac{z^2}{|\mathbf{r }|^3} \Bigg ) \Bigg )+ \frac{3\, \text {Kn}'}{10 \pi }\Bigg ( \frac{\partial }{\partial z} \Bigg ( \frac{3z^2}{|\mathbf{r }|^5} - \frac{1}{|\mathbf{r }|^3} - \frac{3 h z}{|\mathbf{r }|^5} \Bigg ) \Bigg ) , \end{aligned}$$
(51c)
$$\begin{aligned} p^{\perp ''}= & {} -\frac{h}{4 \pi } \Bigg ( \frac{\partial }{\partial z} \Bigg ( \frac{2 z}{|\mathbf{r }|^3} \Bigg ) + \frac{12\, \text {Kn}'^2}{5 h} \frac{\partial }{\partial z} \Bigg ( \frac{3 z^2}{|\mathbf{r }|^5} - \frac{1}{|\mathbf{r }|^3} \Bigg ) \Bigg ) . \end{aligned}$$
(51d)

It is straightforward to repeat the analysis for a point forcing in the parallel direction, as this just changes some of the expressions slightly due to certain expressions no longer cancelling out. The corresponding unknown velocity and pressure fields in this case are

$$\begin{aligned} u_1^{\parallel ''}= & {} \frac{h}{4 \pi \text {Kn}'} \Bigg ( \frac{\partial }{\partial x} \Bigg ( \frac{h x}{|\mathbf{r }|^3} \Bigg ) - \frac{\partial }{\partial x} \Bigg ( \frac{ z x}{|\mathbf{r }|^3} \Bigg ) \Bigg ) - \frac{3\, \text {Kn}'}{10 \pi } \Bigg ( \frac{\partial }{\partial x} \Bigg ( \frac{3 z x}{|\mathbf{r }|^5} \Bigg ) - \frac{\partial }{\partial x} \Bigg ( \frac{ 3 h x}{|\mathbf{r }|^5} \Bigg ) \Bigg ), \end{aligned}$$
(52a)
$$\begin{aligned} u_2^{\parallel ''}= & {} \frac{h}{4 \pi \text {Kn}'} \Bigg ( \frac{\partial }{\partial x} \Bigg ( \frac{h y}{|\mathbf{r }|^3} \Bigg ) - \frac{\partial }{\partial x} \Bigg ( \frac{ z y}{|\mathbf{r }|^3} \Bigg ) \Bigg ) - \frac{3\, \text {Kn}'}{10 \pi } \Bigg ( \frac{\partial }{\partial x} \Bigg ( \frac{3 z y}{|\mathbf{r }|^5} \Bigg ) - \frac{\partial }{\partial x} \Bigg ( \frac{ 3 h y}{|\mathbf{r }|^5} \Bigg ) \Bigg ) , \end{aligned}$$
(52b)
$$\begin{aligned} u_3^{\parallel ''}= & {} \frac{h}{4 \pi \text {Kn}'} \Bigg ( \frac{\partial }{\partial x} \Bigg ( \frac{hz}{|\mathbf{r }|^3} - \frac{1}{|\mathbf{r }|} - \frac{z^2}{|\mathbf{r }|^3} \Bigg ) \Bigg ) - \frac{3 \, \text {Kn}'}{10 \pi }\Bigg ( \frac{\partial }{\partial x} \Bigg ( \frac{3z^2}{|\mathbf{r }|^5} - \frac{1}{|\mathbf{r }|^3} - \frac{3 h z}{|\mathbf{r }|^5} \Bigg ) \Bigg ) , \end{aligned}$$
(52c)
$$\begin{aligned} p^{\parallel ''}= & {} \frac{h}{4 \pi } \Bigg ( \frac{\partial }{\partial x} \Bigg ( \frac{2 z}{|\mathbf{r }|^3} \Bigg ) + \frac{12\, \text {Kn}'^2}{5 h} \frac{\partial }{\partial x} \Bigg ( \frac{3 z^2}{|\mathbf{r }|^5} - \frac{1}{|\mathbf{r }|^3} \Bigg ) \Bigg ) . \end{aligned}$$
(52d)

For the reader who is familiar with the literature on singularities of Stokes flows, what we have shown here is that the image system which is required to solve the boundary value problem is the analogue of the corresponding image system needed to solve the same boundary value problem when one instead has a Stokeslet close to a flat no-slip wall [38]. Recall that the singular point forcing we have defined is only the simplest type of singularity of a flow, and that one can take derivatives to obtain higher-order singularities [39]. In the terminology of the literature on singularities of viscous flows, if the force is in the z-direction the additional singularities in the image system on the other side of the wall apart from the mirror image Gradlet are just the analogous source dipole in the z-direction and a z-dipole of z-Gradlets [30].

3.3 Pressure distribution for a stokeslet perpendicular to a slip surface

We will now derive the total pressure field due to a point forcing in a Stokes flow which is close to a slip boundary and perpendicular to the boundary (Eqs. (3), solutions (7) and boundary conditions (12)). Taking Fourier transforms of Eq. (5c) for the unknown stress tensor field, we obtain

$$\begin{aligned} {\hat{S}}_{1 3}^{\prime \prime \perp }= & {} {\hat{S}}_{3 1}^{\prime \prime \perp } = -\text {Kn}' \Bigg ( \frac{\partial {\hat{u}}_1^{\perp }}{\partial z} - \mathrm{i} k_1 {\hat{u}}_3^{\perp } \Bigg ) , \end{aligned}$$
(53a)
$$\begin{aligned} {\hat{S}}_{2 3}^{\prime \prime \perp }= & {} {\hat{S}}_{3 2}^{\prime \prime \perp } = -\text {Kn}' \Bigg ( \frac{\partial {\hat{u}}_2^{\perp }}{\partial z} - \mathrm{i} k_2 {\hat{u}}_3^{\perp } \Bigg ) . \end{aligned}$$
(53b)

We have only written down the transforms of the components which we will need for our computations. Making some rearrangements with the Stokes equations [38], we see that the unknown pressure satisfies the Laplace equation

$$\begin{aligned} \Delta p'' =0 , \end{aligned}$$
(54)

so as before, after taking Fourier transforms and solving the differential equation one obtains

$$\begin{aligned} {\hat{p}}'' = A \mathrm{e}^{-kz}. \end{aligned}$$
(55)

Again, as with Sect. 3.2, the Fourier transformed components of the unknown velocity field are

$$\begin{aligned} {\hat{u}}_1''= & {} \Bigg ( \frac{h}{8 \pi \text {Kn}'} \Bigg ) (B_1 + \mathrm{i} k_1 A z ) \mathrm{e}^{-kz}, \end{aligned}$$
(56a)
$$\begin{aligned} {\hat{u}}_2''= & {} \Bigg ( \frac{h}{8 \pi \text {Kn}'} \Bigg ) (B_2 + \mathrm{i} k_2 A z ) \mathrm{e}^{-kz}, \end{aligned}$$
(56b)
$$\begin{aligned} {\hat{u}}_3''= & {} \Bigg ( \frac{h}{8 \pi \text {Kn}'} \Bigg ) (B_3 + k A z ) \mathrm{e}^{-kz}. \end{aligned}$$
(56c)

Adding the relevant components for the velocity and stress tensor fields due to a Stokeslet and a mirror image Stokeslet in the perpendicular direction, we obtain

$$\begin{aligned} u_1^{\perp } + u_1^{\perp '}= & {} - \frac{h}{4 \pi \text {Kn}'} \frac{x}{r_h^3} , \end{aligned}$$
(57a)
$$\begin{aligned} u_2^{\perp } + u_2^{\perp '}= & {} - \frac{h}{4 \pi \text {Kn}'} \frac{y}{r_h^3} , \end{aligned}$$
(57b)
$$\begin{aligned} u_2^{\perp } + u_2^{\perp '}= & {} S_{31}^{\perp } + S_{31}^{\perp '} = S_{32}^{\perp } + S_{32}^{\perp '} = 0. \end{aligned}$$
(57c)

The Fourier transformed boundary conditions are

$$\begin{aligned}&{\hat{u}}_1^{\perp } + {\hat{u}}_1^{\perp '} + {\hat{u}}_1^{\perp ''} = - \sqrt{\frac{\pi }{2}} ( {\hat{S}}_{31}^{\perp } + {\hat{S}}_{31}^{\perp '} + {\hat{S}}_{31}^{\perp ''}), \end{aligned}$$
(58a)
$$\begin{aligned}&{\hat{u}}_2^{\perp } + {\hat{u}}_2^{\perp '} + {\hat{u}}_2^{\perp ''} = - \sqrt{\frac{\pi }{2}} ( {\hat{S}}_{32}^{\perp } + {\hat{S}}_{32}^{\perp '} + {\hat{S}}_{32}^{\perp ''}), \end{aligned}$$
(58b)
$$\begin{aligned}&{\hat{u}}_3^{\perp } + {\hat{u}}_3^{\perp '} + {\hat{u}}_3^{\perp ''} = 0 . \end{aligned}$$
(58c)

Taking Fourier transforms of (57) and adding everything at the boundary at \(z=0\), we obtain

$$\begin{aligned}&{\hat{u}}_1^{\perp ''} - \text {Kn}\frac{\partial {\hat{u}}_1^{\perp ''}}{\partial z} = \mathrm{i} \frac{h}{4 \pi \text {Kn}'}\frac{k_1}{k}, \end{aligned}$$
(59a)
$$\begin{aligned}&{\hat{u}}_2^{\perp ''} - \text {Kn}\frac{\partial {\hat{u}}_2^{\perp ''}}{\partial z} = \mathrm{i} \frac{h}{4 \pi \text {Kn}'}\frac{k_2}{k}, \end{aligned}$$
(59b)
$$\begin{aligned}&{\hat{u}}_3^{\perp ''} = 0 . \end{aligned}$$
(59c)

Closing the system with the same Fourier transformed continuity equation as before, we obtain for the coefficients

$$\begin{aligned} B_1^{\perp }= & {} \frac{2 \mathrm{i} k_1}{k(1 + 2 k\, \text {Kn})} , \end{aligned}$$
(60a)
$$\begin{aligned} B_2^{\perp }= & {} \frac{2 \mathrm{i} k_2}{k(1 + 2 k\, \text {Kn})} , \end{aligned}$$
(60b)
$$\begin{aligned} B_3^{\perp }= & {} 0 , \end{aligned}$$
(60c)
$$\begin{aligned} A^{\perp }= & {} \frac{- 2}{1 + 2 k\, \text {Kn}} . \end{aligned}$$
(60d)

Note that the Knudsen number is now the standard Knudsen number because the scaling has cancelled out. After dimensionalisation to account for the slip length \(\lambda \), these are the same Fourier coefficients as the ones obtained in [39] i.e. when \(\text {Kn} \) is replaced by \(\lambda \). In fact, the Fourier coefficients are the same as those in the no-slip case (when the slip length \(\lambda \) is zero) divided through by \(1 + 2 \lambda k\), so by differentiating one can write down the relation

$$\begin{aligned} \Bigg ( 1 - 2 \lambda \frac{\partial }{\partial z} \Bigg ) {\hat{p}}^{\perp ''} (\mathbf{r }, \lambda ) = {\hat{p}}^{\perp ''} (\mathbf{r },0). \end{aligned}$$
(61)

This implies a differential equation for the inverted quantities:

$$\begin{aligned} \Bigg ( 1 - 2 \lambda \frac{\partial }{\partial z} \Bigg ) p^{\perp ''} (\mathbf{r }, \lambda ) = p^{\perp ''} (\mathbf{r },0). \end{aligned}$$
(62)

In [38], it is shown that

$$\begin{aligned} p^{\perp ''}(\mathbf{r },0) = \frac{h}{2 \pi } \frac{\partial }{\partial z} \Bigg (\frac{z}{|\mathbf{r }|^3} \Bigg ), \end{aligned}$$
(63)

so integrating by parts and adding the other contributions \(p^{\perp }\) and \(p^{\perp '}\) due to the Stokeslet and the mirror image Stokeslet, we find that the total pressure due to a point forcing in Stokes flow close to and perpendicular to a flat slip boundary is

$$\begin{aligned} p^{\perp } = \frac{1}{4 \pi } \Bigg ( \frac{z}{|\mathbf{r }|^3} - \frac{z}{|\overline{\mathbf{r }}|^3} + \frac{h}{\lambda } \int ^{\infty }_0 \text {d} s \, \Bigg [ \frac{\partial }{\partial z} \Bigg (\frac{z}{|\mathbf{r }|^3} \Bigg ) \Bigg ] (\mathbf{r } + (h+s) \mathbf{e }_z) \mathrm{e}^{-s/2 \lambda } \Bigg ) , \end{aligned}$$
(64)

where \(\mathbf{e }_z\) is the unit vector in the z-direction.

3.4 Temperature distribution for a thermal stokeslet close to a temperature jump surface

Finally, we consider the case of a point heat source (thermal Stokeslet) in heat flow close to a temperature jump surface (Eqs. (14), solutions (17) and boundary condition (19)). By the same analysis as before, one obtains for the Fourier transformed temperature

$$\begin{aligned} {\hat{\theta }} = A \mathrm{e}^{-kz}. \end{aligned}$$
(65)

The Fourier transformed boundary condition is

$$\begin{aligned} {\hat{\theta }} + {\hat{\theta }}' + {\hat{\theta }}'' = - \frac{1}{2} \sqrt{\frac{\pi }{2}}( {\hat{q}}_3 + {\hat{q}}_3' + {\hat{q}}_3'' ) . \end{aligned}$$
(66)

By the same or very similar analysis as Sect. 3.1, we substitute in the Fourier transformed quantities at the boundary where \(z=0\) and then rearrange for the coefficient A, which is found to be

$$\begin{aligned} A = -\frac{2g}{15\, \text {Kn}'\pi k (1 + \frac{15}{8} \text {Kn}\, k)} . \end{aligned}$$
(67)

Note that this is the same Fourier coefficient as the one which we obtained for the same problem in Sect. 3.1 with a no-jump boundary after dividing through by \(1 + 15\,\text {Kn}\,k/8\). Differentiating as before implies a differential equation for \({\hat{\theta }}''\), which is inverted to give

$$\begin{aligned} \Bigg ( 1 - \frac{15}{8} \text {Kn} \frac{\partial }{\partial z} \Bigg ) \theta '' = \theta _{nj}'', \end{aligned}$$
(68)

where \(\theta _{nj}''\) is the unknown temperature field for the heat source close to the no-temperature jump surface. Substituting in \(\theta ''_{nj}\) which was found in Sect. 3.1, integrating by parts and adding \(\theta ''\) to the temperature fields \(\theta \) and \(\theta '\) contributed by the heat source and the mirror image heat source, we obtain the total temperature field due to a point heat source placed in isobaric stationary Stokes flow close to a temperature jump boundary:

$$\begin{aligned} \theta = \frac{g}{15\, \text {Kn}' \pi } \Bigg ( \frac{1}{|\mathbf{r }|} - \frac{1}{|\overline{\mathbf{r }}|} + \frac{16}{15\, \text {Kn}} \int ^{\infty }_0 \text {d}s \, \Bigg [ \frac{1}{{|\mathbf{r }|}} \Bigg ] (\mathbf{r } + (h+s) \mathbf{e }_z ) \mathrm{e}^{-8 s/ 15\, \text {Kn}} \Bigg ). \end{aligned}$$
(69)

The components of the heat flux can be obtained in a similar way.

4 Solutions for rarefied gas flows close to boundaries

In this section, we will return to rarefied gas flows and heat flows modelled by the Grad equations. We will find in this section that the method of images can be pushed far enough to obtain the solution in the case of a point heat source in a heat flow close to a flat temperature jump surface, but that obtaining the solution in the case of a point forcing in a rarefied gas flow close to a flat velocity slip surface is intractable (even if we assume for simplicity that the point forcing is perpendicular rather than parallel to the boundary).

4.1 Temperature distribution for a thermal gradlet close to a temperature jump surface

We start with the case of a point heat source (thermal Gradlet) placed in a heat flow close to a temperature jump surface (Eqs. (26), solutions (28) and boundary condition (30)). The Fourier transformed temperature jump boundary condition is now quite complicated:

$$\begin{aligned} {\hat{\theta }} + {\hat{\theta }}' + {\hat{\theta }}'' =-\frac{1}{2} \sqrt{\frac{\pi }{2}} ( {\hat{q}}_3 + {\hat{q}}_3' + {\hat{q}}_3'' ) - \frac{1}{4} ({\hat{S}}_{33} + {\hat{S}}_{33}' + {\hat{S}}_{33}''). \end{aligned}$$
(70)

If we Fourier transform Eqs. (26) and solve the resulting differential equations, we obtain

$$\begin{aligned} {\hat{\theta }}''= & {} A \mathrm{e}^{-kz} , \end{aligned}$$
(71a)
$$\begin{aligned} {\hat{q}}_3''= & {} \frac{15}{4}k\, \text {Kn}' A \mathrm{e}^{-kz}, \end{aligned}$$
(71b)
$$\begin{aligned} {\hat{S}}_{33}''= & {} 3\, \text {Kn}'^2 k^2 A \mathrm{e}^{-kz} , \end{aligned}$$
(71c)

Using the solutions (28), we also have at the boundary

$$\begin{aligned} \theta + \theta '= & {} \frac{2 g}{15\, \text {Kn}' \pi r_h} , \end{aligned}$$
(72a)
$$\begin{aligned} q_3 + q_3'= & {} 0, \end{aligned}$$
(72b)
$$\begin{aligned} S_{33} + S_{33}'= & {} \frac{2\, \text {Kn}' g}{5 \pi } \Bigg ( \frac{3h^2}{r_h^5} - \frac{1}{r_h^3} \Bigg ). \end{aligned}$$
(72c)

Taking Fourier transforms of the above and adding everything at the boundary, we arrive after some rearrangements at

$$\begin{aligned} A = - \frac{2g + \frac{15}{10}\text {Kn}'^2 g k^2}{15\, \text {Kn}' \pi k ( 1 + \frac{15}{8}\text {Kn}\, k + \frac{3}{4} \text {Kn}'^2 k^2)}. \end{aligned}$$
(73)

Substituting this back into (71a) and differentiating, we obtain after some algebra the following relation:

$$\begin{aligned} \Bigg ( 1 - \frac{15}{8} \text {Kn} \frac{\partial }{\partial z} + \frac{3}{4} \text {Kn}'^2 \frac{\partial ^2}{\partial z^2} \Bigg ) {\hat{\theta }}'' = - \frac{2g}{15 \text {Kn}' \pi k}\mathrm{e}^{-kz} - \frac{ \text {Kn}' g k}{10 \pi } \mathrm{e}^{-kz} . \end{aligned}$$
(74)

This implies a differential equation for the unknown temperature field \(\theta ''\) after inverting:

$$\begin{aligned} \Bigg ( 1 - \frac{15}{8} \text {Kn} \frac{\partial }{\partial z} + \frac{3}{4} \text {Kn}'^2 \frac{\partial ^2}{\partial z^2} \Bigg ) \theta '' = -\frac{\text {Kn}' g}{10 \pi } \Bigg ( \frac{1}{|\mathbf{r }|^3} \Bigg ( \frac{3 z^2}{|\mathbf{r }|^2} - 1 \Bigg ) \Bigg ) -\frac{2g}{15 \text {Kn}' \pi } \frac{1}{|\mathbf{r }|}. \end{aligned}$$
(75)

Solving this linear PDE for \(\theta ''\), we obtain an explicit solution

$$\begin{aligned} \theta '' = \frac{16 g}{15 C \text {Kn}'^2 \pi }&\Bigg ( \int _0^{\infty } \text {d}s \, K \mathrm{e}^{-(15 \sqrt{2} \sqrt{1/\pi } + C)s/24 \text {Kn}' } (\mathbf{r } + (h+s) \mathbf{e }_z) \nonumber \\&\quad \quad - \int _0^{\infty } \text {d}s \, K \mathrm{e}^{(-15 \sqrt{2} \sqrt{1/\pi } + C)s/24 \text {Kn}' } (\mathbf{r } + (h+s) \mathbf{e }_z) , \end{aligned}$$
(76)

where

$$\begin{aligned} C= & {} \sqrt{450 \pi - 768}, \end{aligned}$$
(77a)
$$\begin{aligned} K= & {} \frac{ (3\, \text {Kn}'^2 - 4y^2 -4z^2)x^2 + (3\, \text {Kn}'^2 - 4z^2) y^2 - 6\, \text {Kn}'^2 z^2 - 2x^4 - 2y^4 - 2z^4 }{|\mathbf{r }|^5}. \end{aligned}$$
(77b)

The two arbitrary constants which appear in this solution vanish, because the temperature goes to zero at infinity and the constants both multiply exponential terms raised to a positive multiple of z. Adding \(\theta ''\) to \(\theta \) and \(\theta '\) finally results in the total temperature field due to a point heat source which is placed close to a flat temperature jump surface.

4.2 Gradlet perpendicular to a slip surface

We will end our discussion of the method of images by considering the case of a point forcing placed in a rarefied gas flow close to and perpendicular to a slip surface (Eqs. (21) and (22), solutions (23) and boundary conditions (25)). \({\hat{p}}^{\perp ''}\), \({\hat{u}}_1^{\perp ''}\), \({\hat{u}}_2^{\perp ''}\) and \({\hat{u}}_3^{\perp ''}\) are the same as in Sect. 3.2 where we studied the point forcing close to a no-slip boundary. The Fourier transformed boundary conditions are

$$\begin{aligned}&u_1^{\perp } + u_1^{\perp '} + u_1^{\perp ''} = - \sqrt{\frac{\pi }{2}}( S_{3 1}^{\perp } + S_{3 1}^{\perp '} + S_{3 1}^{\perp ''} ) - \frac{1}{5} (q_{1}^{\perp } + q_{1}^{\perp '} + q_{1}^{\perp ''} ), \end{aligned}$$
(78a)
$$\begin{aligned}&u_2^{\perp } + u_2^{\perp '} + u_2^{\perp ''} = - \sqrt{\frac{\pi }{2}}( S_{3 2}^{\perp } + S_{3 2}^{\perp '} + S_{3 2}^{\perp ''} ) - \frac{1}{5} (q_{2}^{\perp } + q_{2}^{\perp '} + q_{2}^{\perp ''} ) , \end{aligned}$$
(78b)
$$\begin{aligned}&u_3^{\perp } + u_3^{\perp '} +u_3^{\perp ''}= 0 . \end{aligned}$$
(78c)

Using Eqs. (21), we obtain

$$\begin{aligned} u_1 + u_1^{\perp '}= & {} - \frac{1}{4 \pi \text {Kn}'} \frac{xh}{r_h^3} - \frac{9\, \text {Kn}'}{10 \pi } \frac{xh}{r_h^5}, \end{aligned}$$
(79a)
$$\begin{aligned} u_2 + u_2^{\perp '}= & {} - \frac{1}{4 \pi \text {Kn}'} \frac{yh}{r_h^3} - \frac{9\, \text {Kn}'}{10 \pi } \frac{yh}{r_h^5}, \end{aligned}$$
(79b)
$$\begin{aligned} u_3 + u_3^{\perp '}= & {} 0 , \end{aligned}$$
(79c)
$$\begin{aligned} q_1^{\perp } + q_1^{\perp '}= & {} -\frac{9\, \text {Kn}'}{8 \pi } \frac{hx}{r_h^5} , \, \, \, \, q_2^{\perp } + q_2^{\perp '} = -\frac{9\, \text {Kn}'}{8 \pi } \frac{hy}{r_h^5}, \end{aligned}$$
(79d)
$$\begin{aligned} S_{31}^{\perp } + S_{31}^{\perp '}= & {} S_{32}^{\perp } + S_{32}^{\perp '} =0 . \end{aligned}$$
(79e)

Fourier transforming Eq. (22a) for the stress tensor \(\mathbf{S }\) we get

$$\begin{aligned} {\hat{S}}''_{ 3 1}= & {} - \text {Kn} \Bigg ( \frac{\partial {\hat{u}}''_1}{\partial z} - \mathrm{i} k_1 {\hat{u}}''_3 \Bigg ) - \frac{4}{5}\text {Kn}' \Bigg ( \frac{\partial {\hat{q}}''_1}{\partial z} - \mathrm{i} k_1 {\hat{q}}''_3 \Bigg ), \end{aligned}$$
(80a)
$$\begin{aligned} {\hat{S}}''_{ 3 2}= & {} - \text {Kn} \Bigg ( \frac{\partial {\hat{u}}''_2}{\partial z} - \mathrm{i} k_2 {\hat{u}}''_3 \Bigg ) - \frac{4}{5} \text {Kn}' \Bigg ( \frac{\partial {\hat{q}}''_2}{\partial z} - \mathrm{i} k_2 {\hat{q}}''_3 \Bigg ), \end{aligned}$$
(80b)

It is possible to Fourier transform equations (80) and add everything at the boundary as we have now done several times throughout this article, but after solving for the Fourier coefficients, the problem of taking inverse two-dimensional Fourier transforms is completely intractable.

5 Solution for unsteady thermal gradlet in one dimension

So far all the flows which we have studied have been time-independent. In dimensions two and three, it is also possible to study the solution to the unsteady Stokes equations with a time-dependent point forcing or point heat source applied to the conservation equations, but the analysis is somewhat involved [19, 23]. The next question is whether one can solve the unsteady Grad equations with a time-dependent point heat source. In order to simplify the equations far enough that they can be solved in this setting, we will have to assume one dimension: this is a standard assumption found elsewhere in the literature on the unsteady Grad and R13 equations [40]. We start with the linearised Grad equations from [28] in one dimension with a time-dependent point heat source applied to the conservation of energy equation. This gives us the linearised dimensionless conservation equations

$$\begin{aligned}&\frac{\partial \rho }{\partial t} + \frac{\partial v_1}{\partial x} = 0 , \end{aligned}$$
(81a)
$$\begin{aligned}&\frac{\partial v_1}{\partial t} + \frac{\partial S_{11}}{\partial x} + \frac{\partial p}{\partial x} = 0 , \end{aligned}$$
(81b)
$$\begin{aligned}&\frac{3}{2} \frac{\partial \theta }{\partial t} + \frac{\partial v_1}{\partial x} + \frac{\partial q_1}{\partial x} = g \delta (x) \delta (t) , \end{aligned}$$
(81c)

along with the linearised dimensional equations for the heat flux and the stress tensor

$$\begin{aligned}&\frac{\partial S_{11}}{\partial t} + \frac{8}{15} \text {Pr} \frac{{\overline{w}}_3}{{\overline{w}}_2} \frac{\partial q_1}{\partial x} = - \frac{2}{{\overline{w}}_2} \frac{1}{\text {Kn}} \Bigg ( S_{11} + \frac{4}{3} \text {Kn} \frac{\partial v_1}{\partial x} \Bigg ), \end{aligned}$$
(82a)
$$\begin{aligned}&\frac{\partial q_1}{\partial t} + \frac{5}{4\, \text {Pr}} \frac{\theta _4}{\theta _2} \frac{\partial S_{11}}{\partial x} = - \frac{1}{\theta _2} \frac{5}{2\, \text {Pr}} \frac{1}{\text {Kn}} \Bigg ( q_1 + \frac{5}{2\, \text {Pr}} \text {Kn} \frac{\partial \theta }{\partial x} \Bigg ) , \end{aligned}$$
(82b)

where \(\delta \) is the Dirac delta function and g is the strength of the point heat source. The other coefficients are defined in [28] and depend on the collision model which is being used. The Prandtl number is \(\text {Pr} = \mu c_p/k\), where \(\mu \) is the shear viscosity, \(c_p\) is the isobaric specific heat and k is the thermal conductivity. For simplicity, we will take the values for Maxwell molecules [28]:

$$\begin{aligned} \theta _2 = \frac{45}{8}, \,\,\,\,\, {\overline{w}}_2 = 2, \,\,\,\,\, {\overline{w}}_3 = \theta _4 = 3, \,\,\,\,\, \text {Pr} = \frac{2}{3} . \end{aligned}$$
(83)

As usual for the point heat source, we search for solutions with constant and uniform pressure and velocity. From the ideal gas law, this implies that \(\rho = - \theta \), so the equations become after some re-arranging

$$\begin{aligned}&\frac{\partial q_1}{\partial x} = g \delta (x) \delta (t) , \end{aligned}$$
(84a)
$$\begin{aligned}&\frac{\partial S_{11}}{\partial t} + \frac{8}{15} \frac{\partial q_1}{\partial x} = - \frac{1}{\text {Kn}} S_{11} , \end{aligned}$$
(84b)
$$\begin{aligned}&\frac{\partial q_1}{\partial t} = \frac{2}{3} \frac{1}{\text {Kn}} q_1 + \frac{5}{2} \frac{\partial \theta }{ \partial x} . \end{aligned}$$
(84c)

It is then straightforward to take separate Laplace transforms with respect to x and t, solve the resulting linear system for the Laplace transformed quantities \({\hat{\theta }}\), \({\hat{q}}_1\) and \({\hat{S}}_{11}\) and finish by taking inverse Laplace transforms to arrive at

$$\begin{aligned} q_1= & {} g \delta (t) , \end{aligned}$$
(85a)
$$\begin{aligned} S_{11}= & {} -\frac{8g}{15}\mathrm{e}^{-t/\text {Kn}} \delta (x) , \end{aligned}$$
(85b)
$$\begin{aligned} \theta= & {} gx \Bigg ( \frac{2}{5} \frac{\partial }{\partial t} \delta (t) - \frac{4}{15} \frac{1}{\text {Kn}} \delta (t) \Bigg ) . \end{aligned}$$
(85c)

The use of Laplace transforms for both variables makes sense as we only have first-order derivatives for space and time. In the case of Stokes flow, the preference for using the Laplace transform for the time domain and the Fourier transform for the spatial domain is slightly subtle, and not simply so that one can solve an ODE of at most second order [41]. This solution could be employed in a mesh-free numerical method to model an extremely simple time-dependent heat flow such as the one-dimensional heat flow within a rarefied vapour phase situated between two liquid bodies [28]. An example of such a numerical method (the method of fundamental solutions) is described in [27]. This method enables one to model a flow using fundamental solutions for the equations of motion. The boundary integral method also uses fundamental solutions but unlike MFS, this method requires the use of numerical integrations which can become computationally expensive. In formal terms, MFS avoids use of a computational mesh by approximating the solution via a series of radial basis functions [42, 43]. As far as we are aware, using the solution we have just derived in this method would not be too challenging, as the method has already been applied in three dimensions. We also attempted to solve the unsteady Grad equations in one dimension with a time-dependent point forcing applied to the momentum equation, but found obtaining non-trivial solutions to be intractable.

Fig. 2
figure 2

Streamlines for a Stokes Rotlet b Grad Rotlet

Full size image

6 Grad analogue of the rotlet

In the literature on singularities of viscous flow, the point forcing is only the simplest possible type of singularity [30]. One can also obtain a higher-order singularity of Stokes flow corresponding to rotational motion: in other words, one can study the response of the flow field to a point torque rather than a point forcing. This singularity is discussed in detail by Batchelor in the context of bulk stress in a suspension of non-spherical particles [44] and it should be simple to generalise to the Grad equations. If one introduces a Stokeslet in the flow, it is well known that one obtains higher-order singularities like the Stokes dipole (also called the doublet in the literature) and the Stokes quadrupole by taking the gradient [30]. Below, we do the same with the Gradlet and obtain the fields corresponding to a Grad dipole:

$$\begin{aligned} \mathbf{v }= & {} \Bigg ( \frac{\mathbf{D } - \mathbf{D }^\mathrm{T} - \text {Tr}\, \mathbf{D }}{|\mathbf{r }|^3} + \frac{3 \mathbf{r } \mathbf{r }}{|\mathbf{r }|^5} \Bigg ) \frac{\mathbf{r }}{8 \pi \text {Kn}'} - \Bigg ( \frac{\mathbf{D } - \mathbf{D }^\mathrm{T} - \text {Tr}\, \mathbf{D }}{|\mathbf{r }|^5} + \frac{5 \mathbf{r } \mathbf{r }}{|\mathbf{r }|^7} \Bigg ) \frac{9 \,\text {Kn}' \mathbf{r }}{20 \pi } , \end{aligned}$$
(86a)
$$\begin{aligned} p= & {} \frac{1}{4 \pi } \Bigg ( \frac{3 \mathbf{r } \cdot \mathbf{D } \mathbf{r }}{|\mathbf{r }|^5} - \frac{\mathbf{e } \cdot \mathbf{D } \mathbf{e }}{|\mathbf{r }|^3} \Bigg ) , \end{aligned}$$
(86b)
$$\begin{aligned} \mathbf{q }= & {} \Bigg ( \frac{\mathbf{D } - \mathbf{D }^\mathrm{T} - \text {Tr}\, \mathbf{D }}{|\mathbf{r }|^5} + \frac{5 \mathbf{r } \mathbf{r }}{|\mathbf{r }|^7} \Bigg ) \frac{9\, \text {Kn}' \mathbf{r }}{8 \pi } , \end{aligned}$$
(86c)

where \(\mathbf{D }\) is a rank two tensor of strengths and \(\mathbf{e }\) is a unit vector.

By analogy with the singularities of Stokes flow, we could refer to the symmetric part of the fields as the Grad stresslet (corresponding to straining motion) and the anti-symmetric part as the Grad rotlet (corresponding to rotational motion). The streamlines due to the Grad rotlet are circles around the line through the origin. In Fig. 2, we plot the streamlines for the Stokes and Grad rotlets for comparison when \(\text {Kn} =1\). The colourmap and colour scale are the same for both figures, so one can see that the gradient for the Grad rotlet is much larger. Note that the Grad rotlet has two terms, the first of which will dominate as \(\text {Kn}\) tends to zero, leading to the streamlines which are observed for the Stokes rotlet and described in [3]. The rotlet singularity could have applications in biomedical science because it corresponds to a point torque. If utilised in the correct way it could potentially capture angular momentum of cylindrical or aspherical dust particles close to boundaries in dilute gas flows (relevant for studying inhalation of particulates and droplets) [30]. The next step in that endeavour would be to derive the velocity fields when a Grad rotlet is placed close to a flat no-slip boundary in a rarefied gas flow (as mentioned earlier, a no-slip boundary condition might still be a good approximation for certain microscale rarefied gas flows). It would be desirable for realistic modelling purposes to derive the velocity fields associated to a Grad rotlet close to a slip surface in a rarefied gas flow, but we have already shown that this is an intractable problem even for the simpler case of a point forcing perpendicular to a slip surface.

7 Conclusions

In summary, we have studied some of the classical exact solutions and techniques from potential theory for Stokes flow (in particular, the method of images) and investigated the extent to which these carry over to rarefied microscale gas flows and heat flows modelled by the linearised Grad equations. As a result we arrive at a number of new exact solutions which could be applied to the study of these flows in idealised geometries or used in numerical schemes for modelling purposes [37]. In particular, the integral solution for a thermal Gradlet close to a temperature jump surface could be employed in boundary integral techniques for numerical solution of linear PDEs, where one computes flow and pressure fields by solving for distributions of singularities of the flow [25]. These numerical techniques are very familiar for Stokes flow calculations in physics and engineering and it might be interesting to use them for rarefied gas flow and heat flow calculations. We have also derived several other solutions which could be employed in boundary integral techniques for rarefied gas flow computations (and in fact, the derivation of the pressure field even for the Stokeslet perpendicular to a slip surface is new as far as we know). In some important respects, the techniques do not carry over completely and it is impossible to obtain an analogous solution for a point forcing placed in a rarefied gas flow close to a slip surface or a solution for the unsteady Grad equations in dimension higher than one. As mentioned, the time-dependent solution to the unsteady Grad equations with a heat source could be employed to study very simple one-dimensional heat flow problems [28] or it could be employed in a time-dependent mesh-free numerical method [27, 45]. It also seems possible that the Grad rotlet which we derived could find applications in biomedical science via the study of aspherical particles and droplets which move in dilute gas flows close to flat walls and have non-negligible angular momenta.