{ "cells": [ { "cell_type": "markdown", "id": "61c6f09c-d285-42a9-b3d3-9c166ab5271f", "metadata": {}, "source": [ "# Velocity data from a Superfluid wind Tunnel\n", "We are now going to apply the various methods we learned in the previous class to look at a velocity signal obtained in a superfluid wind tunnel, downstream an obstacle. The aim will mostly be to reproduce the results presented in the original article by Salort et *Al* : [Energy cascade and the four-fifths law in superfluid turbulence](https://hal.science/hal-00635578v2/file/2012-EPL-Cascade.pdf) The details of the experimental setup can be found in the paper, but briefly, longitudinal velocity measurements are obtained by recording the stagnation pressure at the tip of a mignature pressure probe facing the main flow.\n", "\n", "It is shown that, assuming that both the normal and superfluid components of He II have similar velocities ($v_s\\approx v_n$), the stagnation pressure $P_{stag}$ is related to the 'momentum' velocity $v$ (i.e. the velocity of the the center of mass for a fluid particle comprising both components) by the usual relation $P_{stag} \\approx P_{stat} + \\rho v^2/2$. \n", "\n", "The data are available for this class only, courtesy of P.E. Roche and J. Salort. They will eventually be made public with restrictions on their use, so **please do not use/distribute those data outside this class**.\n", " \n", "The original data are complex (acquisition at 12 kHz of in-phase and in-quadrature output of a lockin amplifier, see this [other article](https://hal.science/hal-00508872v2/file/Salort_TSF_PoF_2010.pdf)) and were pre-processed to simplify this class:\n", "the raw voltage signal, measured with a mean velocity $\\approx 1,15$ m/s at T=1,56K, has been cast to a velocity signal using a calibration obtained in-situ from mean velocity measurements. It is available as a Matlab 7 binary file `v_toupie_1.56K.mat`.\n", "\n", "Note that contrary to the Modane data, those where acquired downstream an obstacle. You will thus be able to see the signature of the periodic vortex shedding, the period of which will serve as the reference time scale of the flow. " ] }, { "cell_type": "code", "execution_count": 1, "id": "0bf0028f-c1e4-4b72-9c77-26057d43ad7a", "metadata": {}, "outputs": [], "source": [ "base_path = '/homelocal/pd218333/Documents/Toupie_Data/DataTD/'" ] }, { "cell_type": "markdown", "id": "51955606-8fc6-430c-b360-257bfd6d983b", "metadata": {}, "source": [ "## Load the data" ] }, { "cell_type": "code", "execution_count": null, "id": "bcf8d3be-6d68-4214-83ca-99065ed6382a", "metadata": {}, "outputs": [], "source": [ "from os.path import join\n", "from scipy.io import loadmat\n", "v_raw = loadmat (join (base_path, \"v_toupie_1.56K.mat\"))['v'].flatten()\n" ] }, { "cell_type": "code", "execution_count": null, "id": "c3a67350-f189-46c4-9f43-a2906d51543b", "metadata": {}, "outputs": [], "source": [ "## Show subsampled signal\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "_,axs = plt.subplots(1,2)\n", "\n", "npts = 1000\n", "stride = int(np.ceil (v_raw.size/npts))\n", "axs[0].plot (v_raw[::stride])\n", "\n", "axs[1].hist(v_raw, bins=50, density=True, histtype=\"step\")" ] }, { "cell_type": "markdown", "id": "f5291baa-3d3b-49c3-8d59-5b9f96972f4f", "metadata": {}, "source": [ "## Velocity spectrum\n", "* Compute and display the raw spectrum\n", "* Compensate the spectrum with $f^{5/3}$ to show that it can reasonably be compared with a Kolmogorov spectrum\n", "* Identify the the detachment frequency and compute the Strouhal number. Compare this number with the one in the paper.\n", "* Identify the energy peak at large frequency and filter it out to obtain the same base signal as in the paper.\n", "* Compute the mean velocity, the standard deviation (using either the Parceval theorem or the `np.std` method) and the turbulence intensity ratio. Also compute the turbulent Reynolds number, based on the vortex shedding length $L_0$, and based on the $\\lambda$. Compare those numbers with the paper." ] }, { "cell_type": "code", "execution_count": null, "id": "9d6a6d6f-6a99-4df5-b0a0-7e83c96080e3", "metadata": {}, "outputs": [], "source": [ "## Compute and show raw velocity spectrum\n" ] }, { "cell_type": "code", "execution_count": null, "id": "d00947a8-8293-4868-953b-f43e0845af1d", "metadata": {}, "outputs": [], "source": [ "## Compensate the spectrum and add a subplot to zoom in and show the vortex shedding peak" ] }, { "cell_type": "code", "execution_count": null, "id": "3d4a09fb-e3fd-4780-a3f3-324bfa692fd7", "metadata": {}, "outputs": [], "source": [ "## Filter out large frequency peak" ] }, { "cell_type": "code", "execution_count": null, "id": "dde70570-ff65-4adc-b90a-ccf086baf625", "metadata": {}, "outputs": [], "source": [ "## Compute standard deviation from Parceval theorem and basic np.std. Use that to obtain the turbulent intensity and the Reynolds number" ] }, { "cell_type": "markdown", "id": "f661ebbe-dabe-4eaa-8103-2a97e350a728", "metadata": {}, "source": [ "## Compare global vs local Taylor hypothesis t->x tranform\n", "In the article the authors have integrated the velocity time series to obtain the vector of equivalent spatial positions. Lets compare both!" ] }, { "cell_type": "code", "execution_count": null, "id": "fc114efa-aa00-4e71-9ed4-a53fd4af7b29", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "afb5fba5-2637-4074-a564-3f424ba27527", "metadata": {}, "outputs": [], "source": [ "## Look at the autocorrelation of the signal" ] }, { "cell_type": "code", "execution_count": null, "id": "33d90e86-c1d5-4842-b593-c0fbd8a68514", "metadata": {}, "outputs": [], "source": [ "from scipy.signal import correlate, correlation_lags,butter,filtfilt,detrend\n", "from scipy.interpolate import CubicSpline\n", "def autocorr (v):\n", " ## Compute the autocorrelation 'function' of the signal\n", " auto_corr = correlate (v, v, 'full')\n", " ## Normalize to obtain R\n", " norm = np.sum (v**2)\n", " R = (auto_corr/norm)\n", " lags = correlation_lags (v.size, v.size, 'full')\n", " ## Only retain positive times\n", " R = R[lags >= 0]\n", " lags = lags[lags >= 0]\n", " return R,lags\n", "\n", "## Several methods to remove the mean velocity and any eventual long term trends\n", "method = \"raw\" # raw, chunks or highpass\n", "vm = np.mean (v_mod)\n", "\n", "if method == \"raw\":\n", " v = detrend(v_mod)\n", " sigmav = np.std (v)\n", " R,lags = autocorr(v-vm)\n", "elif method == \"chunks\":\n", " # Cut the signal into smaller peaces so that eventual large time trends\n", " # can be removed peace-wise\n", " v = v_mod-vm\n", " ## Try to avoid the energy ray around 8 Hz\n", " tau_l = 1/9;\n", " stride = int (tau_l * fps);\n", " nchunks = int(np.floor(v.size / stride))\n", " for ii in range(nchunks):\n", " print (f\"{ii} / {nchunks-1}\", end=\"\\r\")\n", " vchunk = v_mod[(ii*stride):(ii*stride+stride)]\n", " if ii == 0:\n", " #R,lags = autocorr(vchunk - np.mean(vchunk))\n", " R,lags = autocorr(detrend (vchunk))\n", " sigmav = np.std(vchunk)\n", " else:\n", " R += autocorr(detrend (vchunk))[0]\n", " sigmav += np.std(vchunk)\n", " R /= nchunks\n", " sigmav /= nchunks\n", "elif method == \"highpass\":\n", " ## High pass filter at 2s (0.5Hz) ~ 80m\n", " fnyq = fps/2;\n", " tau_l = 560 * (0.18/41);\n", " #fhp = 1/tau_l;#Hz\n", " fhp = 2\n", " b, a = butter(3, fhp/fnyq, 'highpass')\n", " v = filtfilt(b, a, v_mod)\n", " R, lags = autocorr(v)\n", " sigmav = np.std (v)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "64c9902a-b974-4c3f-aacd-d41f8aa84712", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "8844d462-8281-4273-9034-b2b080432e29", "metadata": {}, "source": [ "## Dissipation rate from velocity increments structure function\n", "As explained in the original paper, *In the Kolmogorov 1941 framework, the scaling $S^{\\|}_n(r) \\propto (\\epsilon r)^{n/3}$ is expected for inertial range spacial increments* $r$, where $$S^{\\|}_n(r) = \\langle \\left[(\\mathbf{v}(x+r)-\\mathbf{v}(x)) \\cdot \\frac{\\mathbf{r}}{\\mathbf{r}}\\right]^n \\rangle $$ denotes the longitudinal structure function of order $n$ of the velocity increments across distances $r$ ($\\mathbf{r} = r\\mathbf{e_r}$).\n", "\n", "In this section we are going to use the 2nd and 3rd order structure functions to determine the dissipation rate $\\epsilon$ and then use classical HIT relation to derive the Taylor length scale and the Kolmogorov, dissipative, length scale.\n", "\n", "### Second order structure function\n", "For the second order structure function, the pre-factor of the scaling $S_2 = C_2 (\\epsilon r)^{2/3}$ has been determined empirically to be $C_2\\approx 2.1\\pm 0.1$.\n", "\n", "One can show that $S_2(r)$ can be obtained directly from the auto-correlation coefficient of the velocity:\n", "\n", "$$S_2(r) = 2\\langle v^2\\rangle\\left(1-R(r)\\right)$$\n" ] }, { "cell_type": "code", "execution_count": null, "id": "fda13276-3f4b-4b74-af5c-b6c1ac2d70d6", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [], "source": [ "## Compute S2 directly\n", "try:\n", " from numba import njit, prange\n", "except ModuleNotFoundError:\n", " from nonumba import njit\n", " prange = range\n", " \n", "@njit(parallel=True)\n", "def compute_moment_2 (v, n):\n", " S2 = np.zeros ((len(n)), dtype=np.float32).flatten()\n", " for i in prange(len(n)):\n", " nn = int (n[i])\n", " S2[i] = np.mean ((v[nn:]-v[:-nn])**2)\n", " return S2\n", "\n", "maxlen = 0.5\n", "npts = round (maxlen/vm*fps);\n", "\n", "\n" ] }, { "cell_type": "markdown", "id": "bc06596a-1ed1-48df-a608-bfcadc86e1a3", "metadata": {}, "source": [ "### Third order structure function\n", "The celebrated four-fith law is another way to devise the dissipation rate, with the advantage that it is an exact result from HIT theory that does not involve any empirical constant (like $C_2$ above).\n", "For sufficiently large Reynolds numbers this law writes\n", "\n", "$$S_3(r) = -\\frac{4}{5} \\epsilon r$$\n", "There is no simple method (other than loops) to obtain $S_3$ so we have to compute it manually:" ] }, { "cell_type": "code", "execution_count": 119, "id": "1cb2d9c9-9072-48c7-a13c-0f623399b90c", "metadata": {}, "outputs": [], "source": [ "try:\n", " from numba import njit, prange\n", "except ModuleNotFoundError:\n", " from nonumba import njit\n", " prange = range\n", " \n", "@njit(parallel=True)\n", "def compute_moment_3 (v, n):\n", " S3 = np.zeros ((len(n)), dtype=np.float32).flatten()\n", " for i in prange(len(n)):\n", " nn = int (n[i])\n", " S3[i] = np.mean ((v[nn:]-v[:-nn])**3)\n", " return S3\n", "\n", "maxlen = 0.5 #m\n", "npts = round (maxlen/vm*fps);\n" ] }, { "cell_type": "code", "execution_count": null, "id": "195c3063-8829-405c-95d6-5c99b24489f2", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "2da23b6d-b5c1-4a53-b389-70a81554f4cc", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.3" } }, "nbformat": 4, "nbformat_minor": 5 }