Dynamic Measurements with the Bicone Interfacial Shear Rheometer: The Effects of the Numerical Implementation of the Interfacial Boundary Condition

Abstract

The increasing popularity of flow field-based data analysis (FFBDA) techniques has a paradigmatic example in the routines already developed for the rotational oscillating bicone bob interfacial shear rheometer. Such routines use a second order centered finite difference (SOCFD) discretization scheme, in both the vertical and radial coordinates, for the velocity field in the bulk fluid subphase and a first order forward finite difference (FOFFD) scheme in the vertical coordinate for the velocity field at the air/water interface. Such a mixture of schemes causes non-smooth flow fields at the interface that can be tackled by appropriately devising a SOCFD scheme for the vertical coordinate at the interface using a line of “phantom” nodes that merely serve to adequately merge the Navier–Stokes equations and the Boussinesq–Scriven boundary condition at the interface. Here we report on a detailed analysis of the quantitative improvements of such a scheme over the previous one by comparing the structure of the flow fields at and close to the interface, the differences in the interfacial and bulk drag torques on the bicone bob, and the differences in the torque/displacement complex amplitude ratio.

1. Introduction

Interfacial shear rheology aims at characterizing the shear mechanical properties of thin layers supported by a Newtonian liquid subphase. The need for a careful consideration of the role of the liquid subphase was recognized [1,2,3] since the first interfacial viscometers were designed [4,5]. Hence, in order to achieve a proper characterization of the shear rheological properties of a sample, one needs to properly separate the effects of the interface and the subphase. Therefore, the formulation of manageable physical models including a realistic representation of the coupling between the flow fields at the subphase and the interface are required.

Moreover, the availability of increasingly fast hardware (multi-kernel microprocessors and graphic cards) and powerful software platforms (MATLAB^®, GNU Octave, Python^®, Mathematica^®, etc.) has dramatically decreased the computational time to the level that, nowadays, including such flow field-based data analysis (FFBDA) routines for real-time processing of the experimental data in modern interfacial shear rheometers is conceivable [6].

The first FFBDA scheme Reynaert et al. [7] was designed for the magnetic needle ISR in the Helmholtz coil configuration, that was soon followed by another scheme designed for the double wall-ring (DWR) rotational interfacial shear rheometer [8]. Some variants were proposed for the Helmholtz coil magnetic needle ISR [9,10] and the magnetic tweezers configuration too [11]. Very recently, a FFBDA scheme adequate for the micro-button rotational interfacial shear rheometer has been proposed too [6]. Let us mention here that data analysis is not the only source of inaccuracy in interfacial rheometry; for a thorough discussion of experimental problems that affect measurement accuracy see, for instance, reference [12].

The bicone bob rotational rheometer is one of the oldest (see, for instance [13,14,15,16] and references therein) and, together with the DWR, one of the most extended interfacial rheometer configurations, probably because of the simplicity of the elements needed to set it up if a modern rotational rheometer is available [17]. A FFBDA software package for the rotational bicone bob interfacial shear rheometer under oscillatory forcing [18,19,20] was recently described and made publicly available. The program was developed by adapting the bicone geometry the ideas already used in references [11], for the magnetic trap ISR, and [17], for the oscillating rotational bicone rheometer. Briefly, by using the flow field configurations, obtained from numerical solutions of the Navier–Stokes equation with adequate boundary conditions, in an iterative process, the program yielded a converged value of the complex Boussinesq number, $B{o}^{*}$, that was later used to obtain [21] the complex interfacial dynamic moduli, $G^{*}$, or interfacial shear viscosity, ${\eta }^{*}$.

In that implementation [19] (code at [20]), different discretization schemes were used in the horizontal (SOCFD) and vertical (FOFFD) spatial directions at the interface line, which cause a sub-optimal evaluation of the velocity profiles right at the bicone rim [22]. Such errors might affect the calculation of the interfacial drag on the bicone probe and, consequently, the converged values of the complex dynamic moduli. Moreover, imprecision in the evaluation of the velocity fields close to the bicone rim may demand high mesh resolution for a proper representation of the velocity field when non-linear velocity profiles appear at the interface or the subphase, which is precisely the regime where FFBDA schemes are mandatory. Such undesirable effects jeopardize the possibility of including FFBDA schemes in the standard data analysis routines of commercial interfacial rheometers.

Here we report on an improvement of the numerical scheme that now incorporates a fully reworked discretization method that is now second order centered finite differences (SOCFD) in both spatial directions at both the subphase and, more importantly, the interface. In order to achieve the formulation of the SOCFD at the interface, a line of phantom nodes has been defined above the interface. However, an adequate combination of the discretized Boussinesq–Scriven boundary condition and the Navier–Stokes equations yields expressions for the matrix elements corresponding to the interface nodes in which the values of the velocity field at the phantom nodes are not present and, consequently, the values of the velocity field at the phantom nodes do not appear explicitly in the numerical formulation of the problem. A new version of the program (PV-II, hereafter; the previous program version will be labeled as PV-I)) has been made publicly available already [22] (code at [23,24]). No detailed account of the consequential improvement in the results of the data analysis has been provided whatsoever. Here, we provide a comparative analysis of the performance of both program versions in a wide range of physical parameters.

The paper is organized as follows. In Section 2, we make a brief introduction to the mathematical representation of the physical model of the bicone interfacial rheometer and a brief description of the numerical implementation of the new disctretization close to the interface. In Section 3, we illustrate the improvements in the performance of the program by comparing results obtained with both versions of the software, in aspects such as the velocity field and profiles, the interfacial drag torque, and the values obtained for the dynamic moduli both for numerical and experimental data. Conclusions are stated in Section 4.

2. Methods

2.1. Brief Description of the Mathematical Model

For the sake of completeness, we briefly remind here the physical model behind our approach and its mathematical formulation. In interfacial rotational shear rheology with bicone probes, measurements are performed on samples contained in a cylindrical cup, the bicone is placed with its rim exactly at the interface level, and the interface spans between the bicone rim and the cup wall.

The velocity field is obtained by numerical integration of the Navier–Stokes equations for an non-compressible fluid imposing the following approximations: (i) the axisymmetric geometry imposes an azimuthal velocity field, $v(r,z,t)=v_{\theta }(r,z,t)$, (ii) the interface is flat and horizontal, (iii) the conical bob is represented by a null thickness disk, and (iv) Marangoni flows can be neglected. Several circumstances that might take place in the experiment, in particular the formation of a meniscus at the interface, that deviates its shape from being perfectly flat, or non-zero radial and vertical fluid velocities, can not be addressed with this model. The rotor is supposed to perform an oscillatory motion. Then, the azimuthal flow velocity can be assumed to be proportional to the disk edge velocity, $v_{\theta ,b}\left(t\right)$, by separating the time and space dependencies as follows

$$v_{\theta }(r,z,t)=g^{*}(r,z)v_{\theta ,b}\left(t\right)=i{R}_b\omega g^{*}(r,z){\theta }_0{e}^{i\omega t},$$

(1)

where $R_b$ is the conical bob radius, $\omega$ is the angular frequency, ${\theta }_0$ is the oscillatory displacement amplitude, being ${\theta }^{*}\left(t\right)={\theta }_0{e}^{i\omega t}$, and $g^{*}(r,z)$ is the complex non-dimensional spatial velocity amplitude, whose real and imaginary parts are in-phase and out-phase with respect to the bicone velocity.

Under these approximations, the non-dimensional Navier–Stokes equations read as follows

$$iR{e}^{*}{g}^{*}(\overline{r},\overline{z})=\frac{{\partial }^2{g}^{*}(\overline{r},\overline{z})}{\partial {\overline{r}}^2}+\frac{{\partial }^2{g}^{*}(\overline{r},\overline{z})}{\partial {\overline{z}}^2}+\frac{1}{\overline{r}}\frac{\partial g^{*}(\overline{r},\overline{z})}{\partial \overline{r}}-\frac{g^{*}(\overline{r},\overline{z})}{{\overline{r}}^2},$$

(2)

where $R{e}^{*}=\frac{\rho \omega R_c^2}{{\eta }^{*}}$ is the Reynolds number. The spatial coordinates, $\overline{r}$ and $\overline{z}$, have been made non-dimensional with the cup radius, $R_c$, while $1/\omega$ has been chosen as the characteristic time scale. We remark that the subphase viscosity, ${\eta }^{*}$, can be a complex number if the subphase is viscoelastic.

No-slip boundary conditions are imposed at the cup walls and floor, and at the bicone-fluid contact line, together with zero velocity along the cylindrical symmetry axis within the subphase. These boundary conditions are written as

$$\begin{array}{c}{g}^{*}(\overline{r},0)=g^{*}(1,\overline{z})=0,\hfill \\ g^{*}(0,\overline{z})=0,\hfill \\ g^{*}(\overline{r}\le {\overline{R}}_b,\overline{h})=\frac{\overline{r}}{{\overline{R}}_b}.\hfill \end{array}$$

(3)

At the air–water interface (${\overline{R}}_b<\overline{r}<1$ and $\overline{z}=\overline{h}$), we use the Boussinesq–Scriven boundary condition [7,25], which, in cylindrical coordinates, reads as

$$\frac{\partial g^{*}}{\partial \overline{z}}=B{o}^{*}\frac{\partial }{\partial \overline{r}}\left(\frac{1}{\overline{r}}\frac{\partial }{\partial \overline{r}}\left(\overline{r}{g}^{*}\right)\right),\mathrm{at}{\overline{R}}_b<\overline{r}<1,\overline{z}=\overline{h},$$

(4)

where the complex Boussinesq number is defined as $B{o}^{*}=\frac{{\eta }_s^{*}}{R_c{\eta }^{*}}$. Note that the interfacial viscosity, ${\eta }_s^{*}$, can be a complex number too.

The rotational rheometer provides the angular position of the rotor-bicone ensemble and the torque the rheometer is imposing on it. In order to relate this parameters with the hydrodynamic field solution we use the torque balance equation

$$T^{*}\left(t\right)+T_{\mathrm{sub}}^{*}\left(t\right)+T_{\mathrm{surf}}^{*}\left(t\right)=I\frac{{\partial }^2\theta \left(t\right)}{\partial t^2},$$

(5)

where $T^{*}$ is the forcing torque, $T_{\mathrm{sub}}^{*}$ and $T_{\mathrm{surf}}^{*}$ are, respectively, the subphase and surface drag torques and I is the moment of inertia of the rotor-bicone ensemble. The surface and subphase drag torques can be computed from the velocity profiles obtained after solving the hydrodynamic problem as:

$$\begin{array}{cc}\hfill T_{\mathrm{sub}}^{*}& =-i\omega 2\pi R_b{\eta }^{*}{\theta }_0{e}^{i\omega t}{\int }_0^{R_b}{r}^2{\left.\left(\frac{\partial g^{*}}{\partial z}\right)\right|}_{z=h}dr,\hfill \\ \hfill T_{\mathrm{surf}}^{*}& =i\omega 2\pi R_b^2{R}_{c}B{o}^{*}{\eta }^{*}{\theta }_0{e}^{i\omega t}\left(R_b{\left.\left(\frac{\partial g^{*}}{\partial r}\right)\right|}_{r=R_b,z=h}-1\right).\hfill \end{array}$$

(6)

Then, we make the ansatz that the forcing torque can be represented as

$$T^{*}\left(t\right)=T_0{e}^{i\omega t-\delta }=T_0{e}^{-i\delta }{e}^{i\omega t}=T_0^{*}{e}^{i\omega t}.$$

(7)

where $\delta$ is the phase difference between ${\theta }^{*}\left(t\right)$ and $T^{*}\left(t\right)$. Hence, one can define a complex amplitude ratio between the torque and the angular displacement, $A{R}^{*}$, as

$$\begin{array}{cc}\hfill A{R}^{*}& =\frac{T_0^{*}}{{\theta }_0},\hfill \\ \hfill \delta & =\arg \left(A{R}^{*}\right).\hfill \end{array}$$

(8)

Taking into account Equations (1), (5), (6), and (8), the amplitude ratio, the azimuthal velocity amplitude function $g^{*}$, and the complex Boussinesq number, $B{o}^{*}$, are related by the following expression

$$\begin{array}{cc}\hfill A{R}^{*}& =i\omega 2\pi R_b{\eta }^{*}\left[{\int }_0^{R_b}{r}^2\left(\frac{\partial g^{*}}{\partial z}\right){|}_{z=h}dr\right.\hfill \\ \hfill & \left.-R_b{R}_{c}B{o}^{*}\left(R_b\left(\frac{\partial g^{*}}{\partial r}\right){|}_{r=R_b,z=h}-1\right)\right]-I{\omega }^2.\hfill \end{array}$$

(9)

Solving for $B{o}^{*}$ in Equation (9) allows us to define an iterative procedure to obtain a converged value of the complex Boussinesq number from the experimental value of the complex amplitude ratio, $A{R}_{exp}^{*}$ [17,19]. The converged value of $B{o}^{*}$ can then be used to compute the interfacial dynamic shear moduli, $G^{*}$, or viscosity, ${\eta }^{*}$.

2.2. Numerical Methods: Implementation of the New Discretization at the Interface

The difference between PV-I and PV-II resides in the way the equations for the interface nodes are defined. In the previous version [19], we discretized the equation representing the Boussinesq–Scriven boundary condition (Equation (4)) by taking first order finite differences in the vertical coordinate, $\overline{z}$, and SOCFD in the radial direction, $\overline{r}$. In PV-II, on the contrary, we use SOCFD in both directions. This is done by introducing phantom nodes at the positions with indexes $(j,0)$, with ($\lfloor x\rfloor$ indicates the highest integer smaller than or equal to x.) $\lfloor N\overline{R_b}\rfloor +2\le j\le N$. Hence, we have included a line of phantom nodes right above the interface, and we now have ($N+1-(\lfloor N\overline{R_b}\rfloor +2)$) extra equations for the velocity field, $g_{j,0}^{*}$, at the extra nodes in the mesh.

Following [18], we can take advantage of the problem’s symmetry and solve for the flow field in the rectangle $0\le \overline{r}\le 1$ and $0\le \overline{z}\le \overline{h}$. A cartoon version of the mesh is illustrated in Figure 1.

The discrete version of the Boussinesq–Scriven boundary condition (Equation (4)), in the new scheme of SOCFD in both the horizontal and vertical directions, is

$$\begin{array}{cc}\hfill \frac{1}{2}\frac{M}{\overline{h}}\left(g_{j,0}^{*}-g_{j,2}^{*}\right)=B{o}^{*}{N}^2& \left(\frac{g_{j+1,1}^{*}-g_{j-1,1}^{*}}{2(j-1)}-\frac{g_{j,1}^{*}}{{(j-1)}^2}+g_{j-1,1}^{*}+g_{j+1,1}^{*}-2g_{j,1}^{*}\right),\hfill \\ \hfill \forall j& \in \mathbb{Z}/\lfloor N\overline{R_b}\rfloor +2\le j\le N,\hfill \end{array}$$

(10)

while the discrete version of the Navier–Stokes equations at the central nodes is

$$\begin{array}{cc}\hfill iRe{g}_{j,1}^{*}& =N^2\left(g_{j+1,1}^{*}+g_{j-1,1}^{*}-2g_{j,1}^{*}+\frac{g_{j+1,1}^{*}-g_{j-1,1}^{*}}{2(j-1)}-\frac{g_{j,1}^{*}}{{(j-1)}^2}\right)\hfill \\ \hfill & +{\left(\frac{M}{\overline{h}}\right)}^2\left(g_{j,0}^{*}+g_{j,2}^{*}-2g_{j,1}^{*}\right),\forall j\in \mathbb{Z}/\lfloor N\overline{R_b}\rfloor +2\le j\le N.\hfill \end{array}$$

(11)

Solving for $g_{j,0}^{*}$ in Equation (10) and substituting into Equation (11), the value of the velocity field at the phantom nodes, $g_{j,0}^{*}$, may be eliminated, and we are left with

$$\begin{array}{cc}\hfill iRe{g}_{j,1}^{*}& =N^2\left(1+2\frac{M}{\overline{h}}B{o}^{*}\right)\left(g_{j+1,1}^{*}+g_{j-1,1}^{*}-2g_{j,1}^{*}+\frac{g_{j+1,1}^{*}-g_{j-1,1}^{*}}{2(j-1)}-\frac{g_{j,1}^{*}}{{(j-1)}^2}\right)\hfill \\ \hfill & +2{\left(\frac{M}{\overline{h}}\right)}^2\left(g_{j,2}^{*}-g_{j,1}^{*}\right),\forall j\in \mathbb{Z}/\lfloor N\overline{R_b}\rfloor +2\le j\le N.\hfill \end{array}$$

(12)

Consequently there is no need to take the phantom nodes into account when solving the linear system corresponding to the discrete representation of the Navier–Stokes equation. We remark that the limit of negligible interfacial viscoelasticity, ($B{o}^{*}=0$), is well behaved because in such a case the above expression reduces to the one corresponding to a free interface condition ($\frac{\partial g^{*}}{\partial \overline{z}}=0$, at the interface).

Finally, by implementing the corresponding four node formula (see Figure 2) for the interface nodes in Equation (12), factorizing and rearranging the terms in $g_{j,k}^{*}$, we arrive at

$$\begin{array}{c}{g}_{j-1,1}^{*}\left[N^2\left(1+2\frac{M}{\overline{h}}B{o}^{*}\right)\left(1-\frac{1}{2(j-1)}\right)\right]\hfill \\ -g_{j,1}^{*}\left[iRe+N^2\left(1+2\frac{M}{\overline{h}}B{o}^{*}\right)\left(2+\frac{1}{{(j-1)}^2}\right)+2{\left(\frac{M}{\overline{h}}\right)}^2\right]\hfill \\ +g_{j+1,1}^{*}\left[N^2\left(1+2\frac{M}{\overline{h}}B{o}^{*}\right)\left(1+\frac{1}{2(j-1)}\right)\right]+g_{j,2}^{*}\left[2{\left(\frac{M}{\overline{h}}\right)}^2\right]=0,\hfill \\ \forall j\in \mathbb{Z}/\lfloor N\overline{R_b}\rfloor +2\le j\le N.\hfill \end{array}$$

(13)

Hence, the corresponding matrix elements are:

$$\begin{array}{c}A(j,j-1)=N^2\left(1+2\frac{M}{\overline{h}}B{o}^{*}\right)\left(1-\frac{1}{2(j-1)}\right),\hfill \\ A(j,j)=-iRe-N^2\left(1+2\frac{M}{\overline{h}}B{o}^{*}\right)\left(2+\frac{1}{{(j-1)}^2}\right)-2{\left(\frac{M}{\overline{h}}\right)}^2,\hfill \\ A(j,j+1)=N^2\left(1+2\frac{M}{\overline{h}}B{o}^{*}\right)\left(1+\frac{1}{2(j-1)}\right),\hfill \\ A(j,j+(N+1))=2{\left(\frac{M}{\overline{h}}\right)}^2,\hfill \\ \forall j\in \mathbb{Z}/\lfloor N\overline{R_b}\rfloor +2\le j\le N.\hfill \end{array}$$

(14)

3. Results: Comparative Performance

In this section we briefly sketch the most salient features concerning the improvements in the software performance by comparing the results of the new program version (PV-II, hereafter) against those of the previous one (PV-II, hereafter). The comparison will be performed on a typical experimental realization [17,18,19] of the bicone system with a bicone with radius $R_b=34$ mm, in a cup, with radius $R_c=40$ mm, the water bulk subphase being $H=22$ mm deep. The frequency of the oscillation has been set at $\omega =\pi$ rad/s. In the following, variables obtained by means of program versions PV-I and PV-II will be labeled with subscripts I and II, respectively. Regarding practical applications of such data analysis schemes in the lab, the balance between computational cost and computational errors is a key aspect of the problem. Hence, the effect of the mesh size on the numerical results is investigated in detail.

3.1. Velocity Field Configurations

In Figure 3, we show color coded views of the differences between the velocity fields obtained with both versions of the software. The left panel in Figure 3 shows the differences between the real parts of the velocity field, $\left|\Re \left[g_{II}^{*}(r,z)\right]-\Re \left[g_I^{*}(r,z)\right]\right|$, while the right panel shows the differences between the imaginary parts, $\left|\Im \left[g_{II}^{*}(r,z)\right]-\Im \left[g_I^{*}(r,z)\right]\right|$. It is clear that the most significant differences appear close to the interface and, particularly, at the bicone rim. Nevertheless, the observed values, for both the differences between real and imaginary parts are always smaller than $5\times {10}^{-3}$ (remember that the highest value of the velocity field amplitude is 1 just at the bicone rim).

According to Figure 3, relevant differences appear at the interface and just below it. In Figure 4, we show a comparison between the profiles of the real and imaginary parts of the non-dimensional velocity at the interface obtained with both versions of the software, $g_{s,I}^{*}\left(r\right)$ (panels to the left) and $g_{s,II}^{*}\left(r\right)$ (panel to the right), at two different mesh sizes, namely, $440\times 220$ (top row) and $2520\times 1260$ (bottom row).

In each graph, we consider three different dynamical situations where continuous and dashed lines represent the real and imaginary parts, respectively. Case A (black lines) represents a viscoelastic interface (${\eta }_s^{*}=(1-i)\times {10}^{-3}\mathrm{N}·\mathrm{s}/\mathrm{m}$) at high frequency ($f=5$ Hz). Case B (red lines) represents a purely viscous interface (${\eta }_s^{*}={10}^{-5}\mathrm{N}·\mathrm{s}/\mathrm{m}$) at an intermediate frequency ($f=0.5$ Hz). Finally, case C (blue lines) refers to a clean air–water interface at low frequency ($f=0.05$ Hz). It is clear that the slight piece-wise linear aspect of $\Re \left[g_{s,I}^{*}\left(r\right)\right]$ and $\Im \left[g_{s,I}^{*}\left(r\right)\right]$ at the bicone rim is not present anymore in $\Re \left[g_{s,II}^{*}\left(r\right)\right]$ and $\Im \left[g_{s,II}^{*}\left(r\right)\right]$. This effect should be more appreciable for interfacial velocity profiles far from linear, i.e., for low interfacial viscosity and high oscillation frequency.

3.2. Interfacial and Subphase Torques

The differences found in the flow field profile just at the rim of the bicone might bring noticeable changes in the interfacial and subphase drag torques and, consequently, in the amplitude ratio. In the following, we will illustrate the influence of mesh spacing and interfacial viscosity on the drag torques just for the case of purely viscous interfaces on top of a water subphase and at a fixed frequency of 0.5 Hz.

In Figure 5, we show the relative differences in the real, ${\Delta }_r(\Re \left(T_{\mathrm{sub}}^{*}\right))$ (solid symbols), and imaginary parts, ${\Delta }_r(\Im \left(T_{\mathrm{sub}}^{*}\right))$ (open symbols), of the subphase drag torque, $T_{\mathrm{sub}}$, between the solutions obtained with different mesh sizes taking as reference the solution obtained for the mesh with the highest resolution used ($2520\times 1260$ mesh). Black, red, and blue symbols correspond to meshes with $200\times 100$, $440\times 220$, and $1000\times 500$ nodes, respectively. The left and right panels display the results for PV-I and PV-II, respectively.

The results obtained for subphase drag with versions I and II of the program show no apparent differences. This fact can be easily explained, because the subphase drag occurs mostly at the lower surface of the bicone, where the difference between the flow field computed with PV-I and PV-II are negligible. Interestingly, the relative difference in $T_{\mathrm{sub}}$ is practically unaffected by the value of the interfacial viscosity ${\eta }_s^{*}$.

On the contrary, the interfacial drag torque is determined by the velocity field gradient right at the bicone rim and, consequently, the corresponding figure for the interfacial drag torque is expected to show some differences. Such differences can be appreciated in Figure 6, where we show the relative differences real part, ${\Delta }_r(\Re \left(T_{\mathrm{surf}}^{*}\right))$ (solid symbols), and imaginary part, ${\Delta }_r(\Im \left(T_{\mathrm{surf}}^{*}\right))$ (open symbols) of the interfacial drag torque, $T_{\mathrm{surf}}$, between the solutions obtained with different mesh sizes taking as a reference the $2520\times 1260$ mesh solution. Black, red, and blue symbols correspond to meshes with $200\times 100$, $440\times 220$, and $1000\times 500$ nodes, respectively. The left and right panels of Figure 6 correspond, respectively, to the results obtained with PV-I and PV-II.

In both versions of the program, ${\Delta }_r(\Re \left(T_{\mathrm{surf}}^{*}\right))$ and ${\Delta }_r(\Im \left(T_{\mathrm{surf}}^{*}\right))$ decrease with increasing ${\eta }_s$ until they eventually level at ${\eta }_s$ values between ${10}^{-5}$ and ${10}^{-4}\mathrm{N}·\mathrm{s}/\mathrm{m}$. However, ${\Delta }_r(\Re \left(T_{\mathrm{surf}}^{*}\right))$ is considerably larger in PV-I than in PV-II for all the mesh sizes analyzed. Considering ${\Delta }_r(\Im \left(T_{\mathrm{surf}}^{*}\right))$, version II again shows smaller values, but the difference between version I and II is in this case much smaller than for ${\Delta }_r(\Re \left(T_{\mathrm{surf}}^{*}\right))$.

This fact can be explained by the effect of the box size length scale [26]. There are several characteristic length scales that affect the flow behavior in the bicone geometry. Following Fitzgibbon et al. [26], we can distinguish the subphase viscous (Stokes) length scale, ${\ell }_{\omega }={\left(\frac{\nu }{\omega }\right)}^{1/2}$, where $\nu$ is the kinematic viscosity, the interfacial viscous length scale, ${\ell }_{\omega }^s={\left|B{o}^{*}\right|}^{1/2}{\ell }_{\omega }$, the box size, $\tilde{b}$, and finally the bob radius, $R_b$. In the terms of Reference [26], the box size is the interface radial extension, $\tilde{b}=R_c-R_b$. Then, the interfacial length scale for velocity attenuation becomes larger than the box size when

$${\ell }_{\omega }^s>\tilde{b}\Rightarrow \sqrt{|B{o}_{\omega }^{*}|}{\ell }_{\omega }>\tilde{b}\Rightarrow \sqrt{\frac{|{\eta }_s^{*}|}{\sqrt{\eta \rho \omega }}}>\tilde{b}.$$

(15)

In the configuration here considered, the cross over value for ${\eta }_s^{*}$ then should be

$$|{\eta }_s^{*}|\sim {\tilde{b}}^2\sqrt{\rho \omega \eta }\sim 6.4\times {10}^{-5}\mathrm{N}·\mathrm{s}/\mathrm{m},$$

(16)

that is in very good agreement with the cross over value found numerically. Below that value, the interfacial viscous length scale, which depends on the value of ${\eta }_s$, dominates the dynamics, while above that value the box size, which is constant, is the relevant length scale. In the bulk subphase the Stokes length is ${\ell }_{\omega }\sim 0.56$ mm; hence, in this study ${\ell }_{\omega }<H$ and, consequently, the relevant viscous length scale at the subphase is ${\ell }_{\omega }$.

The differences observed in the calculations of $T_{sub}^{*}$ and $T_{surf}^{*}$ certainly have consequences on the total drag torque too. Such differences can be appreciated in Figure 7, where we show the relative differences in the real part, ${\Delta }_r(\Re \left(T_{\mathrm{tot}}^{*}\right))$ (solid symbols), and imaginary part, ${\Delta }_r(\Im \left(T_{\mathrm{tot}}^{*}\right))$ (open symbols) of the total drag over the bicone + rotor assembly, between the solutions obtained with different mesh sizes taking as a reference the $2520\times 1260$ mesh solution. Black, red, and blue symbols correspond to meshes with $200\times 100$, $440\times 220$, and $1000\times 500$ nodes, respectively. The left and right panels of Figure 7 correspond, respectively, to the results obtained with versions I and II of the software.

Generally speaking, both the real and imaginary parts of the relative differences in the total torque are smaller in PV-II. Moreover, for all the mesh sizes analyzed and for both versions of the numerical approach, ${\Delta }_r(\Im \left(T_{\mathrm{tot}}^{*}\right))$ is considerably smaller than ${\Delta }_r(\Re \left(T_{\mathrm{tot}}^{*}\right))$ (see the different scale of the left and right vertical axes). Interestingly, the results for PV-I show a practically constant value of the relative differences in the real part of the total torque, ${\Delta }_r(\Re \left(T_{\mathrm{tot}}^{*}\right))$, while the imaginary part shows trend changes at two different values of ${\eta }_s$. Indeed, the imaginary part levels out below ${\eta }_s\sim {10}^{-5}$ N·s/m and above ${\eta }_s\sim {10}^{-3}$ N·s/m, while showing a steep decrease in between those values upon increasing ${\eta }_s$. On the contrary, for PV-II, the real part of the total torque is flat above ${\eta }_s\sim {10}^{-4}$ N·s/m and increases upon decreasing ${\eta }_s$ below that value, while the imaginary part shows a similar behavior although in this case the cross-over value is ${\eta }_s\sim {10}^{-3}$ N·s/m.

3.3. Complex Amplitude Ratio

The differences in the drag terms must show up in the complex amplitude ratio, $A{R}^{*}$, too. In Figure 8, we show the relative differences in the modulus and argument of the complex amplitude ratio for the two program versions (left and right panels correspond to PV-I and PV-II, respectively). The dependence with the mesh size has been studied and the relative differences are referred to the results of the highest resolution mesh used ($2520\times 1260$ nodes). Black symbol curves correspond to a $200\times 100$ mesh, red symbols to a $440\times 220$ mesh, and blue symbols to a $1000\times 500$ mesh. Filled and empty symbols correspond to relative differences in the modulus and argument, respectively.

The relative differences in the modulus and argument of the complex amplitude ratio show a more complex behavior as a function of ${\eta }_s$, because the contributions of the differences in the interfacial and subphase drags are combined together with the inertia term (see Equation (9)) and nonlinearly transformed into modulus and argument.

In Figure 8, the aforementioned three distinct regions appear too. In this case the cross-over values are located approximately at ${\eta }_s\sim {10}^{-4}$ and ${\eta }_s\sim {10}^{-2}$ N·s/m. Generally speaking, the values of the relative differences for PV-II are always smaller than those corresponding to PV-I, particularly regarding the argument. Moreover, the modulus relative differences shows a monotonic variation with ${\eta }_s$ for PV-II, while this is not true for PV-I. Overall, the convergence of the results as a function of mesh size is quite good because the relative differences always decrease upon increasing mesh spatial resolution. Indeed, the relative differences between the best resolution mesh ($2520\times 1260$ nodes) and the second best resolution mesh ($1000\times 500$ nodes) are about $0.1\%$ in the modulus and $0.005\%$ in the argument. Hence, if in a given application the computing time needs to be shortened, spatial meshes of $1000\times 500$ nodes or even $440\times 220$ nodes might be safely used.

3.4. Local Strain at the Bicone Rim

It is not surprising that physical quantities that involve the response of the whole interface, such as the total torque, the complex amplitude ratio, or the complex interfacial viscosity, show small differences between the values obtained with both versions of the program. The opposite is true concerning physical variables that depend mostly on the local values at the bicone rim, such as the interfacial torque or the interfacial strain close to the bicone rim. The changes in the interfacial torque have already been illustrated in Figure 6. Indeed, the interfacial strain close to the bicon rim is, typically, the reference variable to study whether the measurements are made in the linear regime or not.

In the mathematical representation of the problem used here, the expression for the interfacial strain at the bicone rim is

$${\gamma }^{*}\left(R_b\right)={\left.\left(\frac{\partial u_{\theta }}{\partial r}-\frac{u_{\theta }}{r}\right)\right|}_{r=R_b}={\left.R_b{\theta }_0{e}^{i\omega t}\left(\frac{\partial g^{*}}{\partial r}-\frac{g^{*}}{r}\right)\right|}_{r=R_b}.$$

(17)

Hence

$$\frac{{\gamma }^{*}\left(R_b\right)}{{\theta }_0{e}^{i\omega t}}=R_b{\left.\left(\frac{\partial g^{*}}{\partial r}\right)\right|}_{r=R_b}-g^{*}\left(R_b\right)=R_b{\left.\left(\frac{\partial g^{*}}{\partial r}\right)\right|}_{r=R_b}-1,$$

(18)

and, taking into account Equation (6), the expression for the relationship between the strain at the bicone rim and the interfacial torque is

$$\frac{T_{surf}}{{\theta }_0{e}^{i\omega t}}=i\omega 2\pi R_b^2{R}_{c}B{o}^{*}\eta \left(R_b{\left.\left(\frac{\partial g^{*}}{\partial r}\right)\right|}_{r=R_b}-1\right)=i\omega 2\pi R_b^2{R}_{c}B{o}^{*}\eta \frac{{\gamma }^{*}\left(R_b\right)}{{\theta }_0{e}^{i\omega t}}.$$

(19)

In Figure 9, we show the relative differences in the modulus and argument of the complex interfacial strain at the bicone rim for the two program versions (left and right panels correspond to PV-I and PV-II, respectively). The dependence with the mesh size has been studied and the relative differences are referred to the results of the highest resolution mesh used ($2520\times 1260$ nodes). Black symbol curves correspond to a $200\times 100$ mesh, red symbols to a $440\times 220$ mesh, and blue symbols to a $1000\times 500$ mesh. Filled and empty symbols correspond to relative differences in the modulus and argument, respectively.

On the one hand, the relative differences in the modulus of the complex interfacial strain at the bicone rim, ${\Delta }_r\left(\right|{\gamma }^{*}\left(R_b\right)\left|\right)$, do not differ significantly for both program versions. On the other hand, the relative differences in the argument of the complex interfacial strain at the bicone rim, ${\Delta }_r(arg\left({\gamma }^{*}\left(R_b\right)\right))$, show large differences for ${\eta }_s\le {10}^{-4}$ N·s/m. The relative differences in ${\Delta }_r(arg\left({\gamma }^{*}\left(R_b\right)\right))$ between program versions I and II are particularly important in the range ${10}^{-5}\lesssim {\eta }_s\lesssim {10}^{-4}$ N·s/m, where the relative differences obtained with version PV-I are higher by a factor of 10 than those corresponding to version PV-II. This variation, in favor of version PV-II, is important because the range of ${\eta }_s$ between ${10}^{-5}$ and ${10}^{-4}$ N·s/m is, precisely, where more gain should be expected from using FFBDA schemes in the bicone geometry, provided the rheometer’s resolution in torque and angular displacement is good enough.

3.5. Consistency Tests

Consistency tests consist of programming a given interfacial viscosity, solving the flow field configuration, obtaining the corresponding value of the complex amplitude ratio, and feeding the data analysis scheme with those data to check if the FFBDA scheme yields the same values that were programmed initially. The test is performed for a set of 30 values of the interfacial viscosity with logarithmic spacing in the interval ${10}^{-6}\le {\eta }_s\le {10}^0$ N·s/m. All of the results shown here have been obtained with a high resolution mesh ($1000\times 500$ nodes).

The consistency tests performed on the iterative scheme using PV-II, with the same protocol employed in [19], are shown in the left panels of Figure 10. Different rows correspond to the three aforementioned representative cases. The top row corresponds to a perfectly viscous interface, so that the real part of the complex interfacial viscosity should coincide with the programmed value of ${\eta }_s$ and the value of the imaginary part should be negligible, as it is shown in the top left panel of Figure 10, where the real part of the calculated values of the interfacial viscosity is always larger by at least three orders of magnitude than the corresponding values of the imaginary part. The middle row corresponds to a viscoelastic interface with equal real and imaginary components of the perfectly viscous interface, so that both the real and imaginary parts of the complex interfacial viscosity should coincide with the programmed value of ${\eta }_s$. This is precisely what is shown in the middle left panel of Figure 10 where the real and imaginary parts superpose nicely in a single symbol line. The bottom row corresponds to a perfectly elastic interface, so that the real part of the complex interfacial viscosity should be negligible, while the imaginary part should coincide with the programmed value of ${\eta }_s$. This happens again, as shown in the bottom left panel of Figure 10, but in a small region slightly above ${\eta }_s\sim {10}^{-5}$ N·s/m; such a behavior appeared also in PV-I [19].

The tests yielded very similar results to those reported for version I in [19]. This is not strange since in the consistency tests of each version of the software, the initial flow configuration is obtained using the same version that is used to process the data afterwards. The curves in the right side panels of Figure 10 show the number of iterations needed for convergence at each ${\eta }_s$ value for PV-II. They show very similar trends to those corresponding to PV-I [19] but the numbers of iterations needed for convergence with PV-II are slightly larger (2 to 3 more iteration steps, typically) than those corresponding to PV-I, probably a consequence of the more involved mathematical implementation of the boundary condition at the interface.

3.6. Comparative Results on Experimental Data

In the following, we illustrate the differences in the performance of both program versions on experimental data obtained in, first, thin silicone oil films of known interfacial viscosity and depth, and, second, on a pentadecanoic acid (PDA) Langmuir monolayer under an isothermal compression. All of the results shown here have been obtained with a high resolution mesh ($1000\times 500$ nodes); hence, the differences in the results obtained with program versions PV-I and PV-II have been minimized.

In Figure 11, we show the results of processing experimental data obtained at an air/water interface covered with three different silicone oil films. Red symbols: $\eta =33$ Pa·s, $d=150$ m; ${\eta }_s=5\times {10}^{-3}$ N·s/m. Black symbols: $\eta =33$ Pa·s, $d=75$ m; ${\eta }_s=2.5\times {10}^{-3}$ N·s/m. Blue symbols: $\eta =1$ Pa·s, $d=100$ m; ${\eta }_s={10}^{-4}$ N·s/m. In the left panel, circles and triangles correspond to the interfacial viscosity, ${\eta }_s\left(\omega \right)$, and interfacial loss modulus, $G_s^{″}\left(\omega \right)$, respectively. The values obtained with both program versions are indistinguishable in the logarithmic vertical scale chosen. Moreover, the values obtained for the interfacial viscosity agree very well with the expected value in the three cases, and the loss modulus shows a linear dependence on frequency (straight line with unity slope in the doubly logarithmic plot). Such data convergence was achieved with 5 iteration steps or less.

The corresponding relative differences in the values of the loss modulus, ${\Delta }_r\left({\left[G_s^{″}\right]}_{I,II}\right)=\frac{{\left[G_s^{″}\right]}_I-{\left[G_s^{″}\right]}_{II}}{{\left[G_s^{″}\right]}_{II}}$, and the interfacial viscosity, ${\Delta }_r\left({\left[{\eta }_s^{\prime }\right]}_{I,II}\right)=\frac{{\left[{\eta }_s^{\prime }\right]}_I-{\left[{\eta }_s^{\prime }\right]}_{II}}{{\left[{\eta }_s^{\prime }\right]}_{II}}$, processed with both versions of the program are shown in the right panel of Figure 11. As expected, the relative differences coincide for the loss modulus and the interfacial viscosity exception made of the points at which both moduli take values comparable to the resolution limit of the instrument ∼${10}^{-4}$ N/m [18], that according to the data shown in the left panel corresponds, approximately, to $\omega <2$ rad/s. Hence, the circles and triangles appear superposed in the right panel plot but at $\omega <2$ rad/s. Not surprisingly, the relative differences are much bigger for the case of the lower value of the interfacial viscosity because the case with ${\eta }_s={10}^{-4}$ N·s/m is close to the resolution limit of the instrument [18] and, consequently, small errors in the flow field solution yield larger errors in the determination of the values of the rheological properties. Moreover, the relative differences grow with frequency, which is not surprising because the interfacial velocity profile becomes increasingly nonlinear as the frequency is increased, and the difference between the profiles determined with both versions of the program should increase accordingly.

In the left panel of Figure 12, we show plots of the values obtained for the interfacial loss modulus, $G_s^{″}\left(\gamma \right)$, as a function of the strain, $\gamma$, for the silicone oil films with ${\eta }_s=5\times {10}^{-3}$ N·s/m (red triangles; expected loss modulus $G_s^{″}\sim 1.5\times {10}^{-2}$ N/m) and ${\eta }_s={10}^{-4}$ N·s/m (blue circles; expected loss modulus $G_s^{″}\sim 3\times {10}^{-4}$ N/m)). In these measurements the frequency was set at $f=0.5$ Hz ($\omega =\pi$ rad/s), and the oscillation amplitude was varied. The left panel shows that the measurements are well within the linear regime, and good agreement with the expected values of the loss modulus is obtained. The corresponding relative differences in the values of the loss modulus processed with both versions of the program are shown in the right panel of Figure 12. Again, the relative differences are larger for the case of the lower value of the interfacial viscosity.

In Figure 13, we show (red triangles) the values of ${\left[G_s^{″}\right]}_{II}$ (red triangles) as a function of the interfacial pressure, $\Pi$, of a pentadecanoic acid Langmuir monolayer [19] obtained with PV-II, together with the relative difference, ${\Delta }_r\left({\left[G_s^{″}\right]}_{I,II}\right)$, between the values obtained when analyzed with both versions of the program calculated with two different mesh resolutions, namely, $1000\times 500$ nodes (black dots) and $200\times 100$ nodes (blue dots). From the results obtained with the high spatial resolution mesh (black dots), it is clear that the values of the difference, due to higher errors in the computations made with PV-I, are very small when the interfacial viscosity is larger than ${10}^{-3}$ N·s/m, roughly corresponding to interfacial pressure values above 25 mN/m, while they grow up to values about $0.1\%$ for lower values of the interfacial viscosity. The comparison with the results obtained using the low spatial resolution mesh (blue dots) shows that the differences increase by, roughly, a tenfold factor.

4. Conclusions

In this report, we show a detailed study of the influence of the discretization scheme close to an air/water interface in a FFBDA scheme for the bicone bob interfacial rheometer under oscillatory forcing. We have compared the performance of two program versions, using first order forward (FOFFD) and second order centered (SOCFD) finite differences discretization schemes at the interface, respectively. The SOCFD scheme needs the introduction of a line of phantom nodes above the interface. Such nodes, however, do not appear explicitly in the linear system of equations that has to be solved numerically because they are eliminated when combining the Navier–Stokes Equation (11) and the Boussinesq–Scriven boundary condition at the interface (10).

The velocity profiles at the interface show significant differences at the bicone rim, where the solutions provided by the SOCFD scheme are much smoother than those pertaining to the FOFFD one (Figure 4). Such an advantage is more noticeable the larger the nonlinearity of the profile, i.e., the lower the interfacial viscosity and the higher the oscillation frequency and the lower the spatial resolution of the mesh used in solving the hidrodynamic problem. Interestingly, the advantage in the accuracy of the flow profiles, at the highest mesh resolution used here, shows up as a few percent differences on the interface and subphase drag torques (Figure 5, Figure 6 and Figure 7), on the complex amplitude ratio (Figure 8), and on the final values obtained for the complex viscosity or dynamic moduli (Figure 10). Hence, using high resolution meshes, the values of the dynamic moduli or the complex viscosity obtained with PV-I show only marginal differences with respect to the more correct ones obtained with PV-II. Consequently, the operating windows previously reported [18] for a bicone bob ISR supplemented with BiconeDrag PV-I should be valid too when using PV-II. Regarding the computation time, PV-II requires, typically, only about $20\%$ more iteration steps than PV-I to achieve convergence. Such an increase in time cost is not relevant for the typical applications of the program.