MATLAB Examples

Example

The friction velocity (u_star) is defined as:

$u_{*} = \left( \overline{u^{\prime}w^{\prime}}^{2}+\overline{v^{\prime}w^{\prime}}^{2}\right)^{1/4}$

It is here computed from the three wind components, with unknown instrument tilt. The method proposed by Klipp [1] is compared to the standard method, which is based on a direct calculation of the covariance between the three wind components. This comparison is done using the function "frictionVelocity" with an optional parameter named 'method', which is either 'standard' (by default) or 'Klipp2018'.

[1] Klipp, C. (2018). Turbulent friction velocity calculated from the Reynolds stress tensor. Journal of the Atmospheric Sciences, (2018).

Contents

Generate tilted wind velocity records

clearvars;close all;clc;
% load simlated wind velocities
load('data.mat','u','v','w','t')
[M,N]=size(u);

% rng(2)

yaw  = 35; % yaw angle (horizontal plane, between u and v)
elev = 3; % elevation angle (vertical plane, between w and u)
incl = 5; % inclination angle (vertical plane, between w and v)


R1 = [cosd(yaw),-sind(yaw),0;sind(yaw),cosd(yaw),0;0,0,1]; % matrix rotation around axis z
R2 = [cosd(elev),0,-sind(elev);0,1,0;sind(elev),0,cosd(elev)]; % matrix rotation around axis zy
R3 = [1,0,0; 0,cosd(incl),-sind(incl);0,sind(incl),cosd(incl)]; % matrix rotation around axis x


A = R3*R2*R1;  % 3D rotation matrix

% Construction of tilted velocity component
u_tilted = zeros(size(u));
v_tilted = zeros(size(u));
w_tilted = zeros(size(u));
for ii=1:M,
    dummy = A\[u(ii,:);v(ii,:);w(ii,:)];
    u_tilted(ii,:) = dummy(1,:);
    v_tilted(ii,:) = dummy(2,:);
    w_tilted(ii,:) = dummy(3,:);
end

Computation of the friction velocity

u_star2 = nan(1,M);
u_star1 = zeros(1,M);
R = nan(3,3,M);
for ii = 1:M,
    [u_star2(ii)] = frictionVelocity(u_tilted(ii,:),v_tilted(ii,:),w_tilted(ii,:),'method','Klipp2018');
    [u_star1(ii)] = frictionVelocity(u(ii,:),v(ii,:),w(ii,:));
end


clf;close all;
fig5 = figure;
plot(u_star2,u_star1,'b*');
daspect([1 1 1])
hold on;box on;

p = polyfit(u_star2,u_star1,1);
x = linspace(nanmin(u_star2(:)),nanmax(u_star2(:)),2);
y = polyval(p,x);
plot(x,y,'k')
[r p] = corr(u_star2',u_star1');
title(['R^2 = ' num2str(r.^2,4)]) ; % Add the value of r to the plot.
plot([min([u_star1,u_star2]),max([u_star1,u_star2])],[min([u_star1,u_star2]),max([u_star1,u_star2])],'r--')
legend('Measured','1:1','best fit','location','southEast')
ylabel('target friction velocity (m/s)')
xlabel('Retrieved friction velocity (m/s)')
set(gcf,'color','w')
xlim([min([u_star1,u_star2]),max([u_star1,u_star2])])
ylim([min([u_star1,u_star2]),max([u_star1,u_star2])])