Matlab notes

Matlab workshops

Table of Contents

• 1. Examples
  • 1.1. find – matlab
  • 1.2. reshape of averaged 2D data
  • 1.3. average vector from multiple files–Matlab
  • 1.4. radial profiles of mean streamwise velocity at X/D=3
  • 1.5. X velocity of 2nd stokes wave
    • 1.5.1. Error
    • 1.5.2. XVelocityStokes.m
  • 1.6. solve an equation to find pitch angle for implicit optimization of blade
    • 1.6.1. bem.m
  • 1.7. Power coefficient in post processing of Cm history
  • 1.8. column vector to row vector
  • 1.9. row vector to column vector
  • 1.10. find the root of an single nonlinear equation
  • 1.11. calculate Cp history
  • 1.12. assign column data to a variable
  • 1.13. Extract a column data and assign to a variable and do math calculations

1 Examples

 

1.1 find – matlab

find the location of a number in a matrix in matlab

matrix

x =
 2     2     3
 4     3     2
 6     4     8

now I want to get the location of a number 4.

I want ans like this :

ans=(2,1) (3,2)

as these are the locations for 4 in matrix.

[i,j]=find(x==4)

1.2 reshape of averaged 2D data

ans=(2,1) (3,2)
X=X.'   % transpose
x_array=reshape(X, 499,[]);  
save -ascii x_array.dat x_array

B = reshape(A,sz) reshapes A using the size vector, sz, to define size(B). For example, reshape(A,[2,3]) reshapes A into a 2-by-3 matrix. sz must contain at least 2 elements, and prod(sz) must be the same as numel(A).

exampple:

reshape a 1-by-10 vector into a 5-by-2 matrix

>> A=1:10；
>> B=reshape(A,[5,2])
B = 5×2

     1     6
     2     7
     3     8
     4     9
     5    10

1.3 average vector from multiple files–Matlab

n=359;
a=[];
b=[];
c=[];
% for loop
for i=1:n  
    filename=sprintf('output_%d.dat',i);
    fileinfo = importdata(filename);
    ux=fileinfo.data(:,7);   % extract column data
    uy=fileinfo.data(:,8);
    uz=fileinfo.data(:,9);
    a=[a,ux];   % put all Ux in an array
    b=[b,uy];
    c=[c,uz];

end
aver_Ux=sum(a,2)/n;  % sum of array elements
aver_Uy=sum(b,2)/n;

aver_Uz=sum(c,2)/n;

• sum

1.4 radial profiles of mean streamwise velocity at X/D=3

load aver_ux_array.dat;
load z_array.dat;
r=z_array(:,1);
r=r.'
r_j=0.00125;
r_nor=r/d;
ux_x_3d=aver_ux_array[41,:];  % extract 41 arrow, which is X/D=3
ux_x_3d=ux_x_3d.';  % transpose
u0=1102.9;
ux_nor=ux_x_3d/u0;
ux_vs_r_3d=[r_nor ux_nor];
save -ascii ux_vs_r_3d_nor.dat ux_vs_r_3d;
plot(r_nor,ux_nor);
xlabel('R/R_{j}');
ylabel('U/U_{j}');
legend({'X/D_{j}=3'},'Location','northeast');
print('ux_vs_r_3d','-depsc');  %Encapsulated PostScript® file
print('ux_vs_r_3d','-dpng');  %png image

1.5 X velocity of 2nd stokes wave

1. define a UDF named 1.5.2
2. type the following commands:

t=[0:1:360]*pi/180;
plot (t,XVelocityStokes(t));

output

1. get the Y coordinates in the plot, and assign it to a new variable, 'Y'

>> points=plot (t,XVelocityStokes(t));
>> Y=get(points,'ydata');

1. save figures

The print command allows you to send plots to you printer and to save plots in a variety of formats. For example,

print -dpsc

prints the current figure to a color PostScript printer. And,

print -deps foo.eps

saves the current figure to an encapsulated PostScript file called foo.eps.

1. assign X and Y coordinates to a new variable, 'XY", and save to a txt file

XY=[t;Y];

save -ascii  X_velocity.txt XY

1.5.1 Error

XVelocityStokes at line 12 column 2 error: evaluating argument list element number 2

1.5.2 XVelocityStokes.m

function u=XVelocityStokes(t) % x velocity of 2nd stokes wave U=0.6; H=0.08; g=9.81; L=4.8; k=2.0*pi/L; d=1.6; T=1.78; z=0; omega=2.0*pi/T; u= H*g*k/(2*(omega-k*U))*exp(k*z)*(1+exp(-2*k*(d+z)))/(1+exp(-2*k*d))*cos(omega*t);

1.6 solve an equation to find pitch angle for implicit optimization of blade

description: Solve the following equation: variable: φ

TSR*X={sin(\phi)(2cos \phi -1)}/{(1+2cos\phi)(1-cos\phi)}

• X, radial distance, (0.1-1)

Use the “bem.m” function and “fzero” function in the command window: >> z=fzero(@bem,0.6)

http://uk.mathworks.com/help/matlab/ref/fzero.html?searchHighlight=fzero&s_tid=doc_srchtitle

1.6.1 bem.m

  y=TSR*X
x = \phi

function y = bem( x )
%Equation (1)
%   TSR*X=1.2,X=0.4 in this example
y=sin(x)*(2*cos(x)-1)/((1+2*cos(x))*(1-cos(x)))-1.2;
        
end

1.7 Power coefficient in post processing of Cm history

In Fluent, you can only get unscaled moment and drag To get power/thrust coefficient of a tidal turbine, you need to setup a user defined function in Excel or matlab

The expression of Cp:

Cp=moment*\omega /(0.5 \rho v^3  A)

where the unit of ω is in rad/s

1.8 column vector to row vector

.' #transposes

To get any vector to be a row vector simply first force it to be a column vector, then transpose it:

A = rand(3,1);
B = A(:).'; % (:) forces it to unwrap along columns, .' transposes

A row vector gets unwrapped to a column vector and then transposed to a row vector again, and a column vector gets unwrapped to a column vector (thus doesn't change) and then gets transposed to a row vector.

Note that any N-dimensional matrix will be forced into a row vector this way due to column unwrapping:

A = rand(3,3,3,3);
B = A(:).'; % 81 element row vector

https://stackoverflow.com/questions/53173615/how-to-force-a-vector-to-be-row-vector

1.9 row vector to column vector

a=a(:) reshapes all elements of A into a single column vector. This has no effect if A is already a column vector. https://stackoverflow.com/questions/13413930/change-row-vector-to-column-vector

1.10 find the root of an single nonlinear equation

Problem: find the root of the equation:

\begin{equation} TSR*x=\frac{sin \phi (2cos\phi -1)}{(1+2cos\phi)(1-cos\phi)} \label{eq:inflow angle} \end{equation}

soluton: define a function named "bem.m"

function y = bem( phi )
% find the inflow angle
%   X=r/R
x=1;
TSR=4.25;
y=sin(phi)*(2*cos(phi)-1)/((1+2*cos(phi))*(1-cos(phi)))-TSR*x;
end

Type the following command in the matlab command window

>> clear
>> fun=@bem;
>> z=fzero(@bem,0.6)

then the output: z =

0.7793

• use fzero function: root of nonlinear function

syntax:

>> x = fzero (fun,x0)

• fzero solves fun(x)=0

1.11 calculate Cp history

input data : unscaled time history of moment/thrust

example of input data, "moment.out"

4283 6.983 -0.1910873139538953
4284 6.984 -0.191019809738711
4285 6.985 -0.1909838904131738
4286 6.986 -0.190943968230172
4287 6.987 -0.1908886443401208
4288 6.988 -0.1908541205872921

1. To load input data:

> load moment.out

1. extract time and unscaled moment and assign it to a variable named "m"

> m=moment(:, 3); # assign the 3rd column data of "moment.out" to variable "m"

> time = moment(:,2);

1. calculating power coefficient, Cp

cp=coeff(m); # "coeff " is a function to calculate the Cp of 3 rotor based on input moment (one rotor )

1. assign "time" and "cp" variables to a new variable "cp_his"

> cp_his = [time cp];

1. save the workspace variable, "cp_his", to a txt file

> save -ascii cp_his.txt cp_his

## syntax: save ; file format; export file name; variable

1. plotting the time history of Cp

> plot (time, cp);

1.12 assign column data to a variable

> load moment.out # load file

> m=moment(:, 3); # assign the 3rd column data of "moment.out" to variable "m"

1.13 Extract a column data and assign to a variable and do math calculations

1. To load data file:

>> load data.txt

1. To assign the second column data to a variable named "m" :

> m = data(:,2) ;

## (:, 2) : row, all, column, 2nd

## ";" is used to supress the output

1. To export the data file to a text file

https://www.youtube.com/watch?v=E56egH10RJA

keywords:

Automatic arithmetic operation of column data 一列数据自动代数运算

use Readtable function

Author: ka

Created: 2020-01-11 六 16:12

Emacs 24.5.1 (Org mode 8.2.10)

Validate