Simple Propeller example

In this example a simple propeller submitted to a constant uniform flow is considered. The geometry is given by three blades with circular cross section.

if ~(length(getenv('TESTS_RUN')) > 0 && strcmp( getenv('TESTS_RUN'), 'yes')), close all, clear all, end
%
addpath( genpath( [ pwd '/../../src'] ) );

Problem definition

The blades are considered considerable stiff and only lift is considered, thus a rigid rotation analytic solution can be used to verify the numerical solution.

The wind velocity is given by the windVel function

va = feval('windVel', 0,0) ;

and the blades are considered to have only lift, with a lift coefficient given by the function liftCoef

c_l = feval('liftCoef', 0) ;

the density and kinematic viscosity are given by

rhoA = 1.225 ; nuA = 1.6e-5 ;

The material parameters of the blades is given by steel

E = 210e9 ;  nu = 0.3 ; rho = 6000 ; 

and the geometric parameters length and diameter are

l = 3 ; d = 0.1;

Analytical solution

compute solution using the second cardinal, first the wind parameters are loaded lift load per unit of length:

fl = 1 / 2 * c_l * rhoA * norm(va) ^ 2 * d ;

the total moment induced in node 1 in x direction for is the sum for three blades:

moment1x = 3 * fl * l * l / 2 ;

and then the angular moment is:

bladeMass = rho * l * pi * d ^2 /4 ; 
Jrho =  3 * 1/3 * bladeMass  * l ^ 2 ; 
angleXnode1 = @(t)  moment1x / Jrho / 2 * t .^ 2 ;

Numerical solution

The material parameters are

materials                 = struct() ;
materials.hyperElasModel  = '1DrotEngStrain' ;
materials.hyperElasParams = [ E nu ]        ;
materials.density         = rho             ;

elements

The elements are

% nodes
elements             = struct() ;
elements(1).elemType = 'node'  ;
% blades
elements(2).elemType = 'frame' ;
elements(2).elemCrossSecParams = {'circle' ; d };
elements(2).massMatType =  'consistent'        ;
elements(2).dragCoefFunction = @(beta,Re) 0;
elements(2).liftCoefFunction ='liftCoef';
elements(2).momentCoefFunction = @(beta,Re) 0;
% auxiliar element
elements(3).elemType = 'truss' ;
elements(3).elemCrossSecParams = {'circle' ; 1.5*d };
elements(3).massMatType =  'lumped'        ;

boundary Conditions

boundaryConditions             = struct() ;
boundaryConds(1).imposDispDofs = [ 1 3 4 5 6 ] ;
boundaryConds(1).imposDispVals = [ 0 0 0 0 0 ] ;

mesh

mesh             = struct() ;
mesh.nodesCoords = [ 0        0              0            ; ...
                     0  l*sin( pi )        l*cos( pi )    ; ...
                     0  l*sin( pi/3  )     l*cos( pi/3 )  ; ... 
                     0  l*sin( 4*pi/3 )   -l*cos( 4*pi/3 ); ...
                     -d*.75 0 d ; ...
                     -d*.75 0 -l*1.5 ] ;

mesh.conecCell         = { } ;
mesh.conecCell{ 1, 1 } = [ 0 1 1   1   ] ;
mesh.conecCell{ 2, 1 } = [ 1 2 0   1 2 ] ;
mesh.conecCell{ 3, 1 } = [ 1 2 0   1 3 ] ;
mesh.conecCell{ 4, 1 } = [ 1 2 0   1 4 ] ;
% auxiliar elements
mesh.conecCell{ 5, 1 } = [ 0 1 1   5 ] ;
mesh.conecCell{ 6, 1 } = [ 0 1 1   6 ] ;
mesh.conecCell{ 7, 1 } = [ 1 3 0   5 6 ] ;

initial Conditions

% homogeneous initial conditions are considered, then an empty struct is set:
initialConds = struct() ;

analysisSettings

analysisSettings                        = struct() ;
analysisSettings.finalTime              =   400     ;
analysisSettings.deltaT                 =   5       ;
analysisSettings.methodName             = 'alphaHHT';
analysisSettings.stopTolIts             =   50      ;
analysisSettings.geometricNonLinearAero = true      ;
%analysisSettings.booleanSelfWeight      = false     ;
analysisSettings.stopTolDeltau          =   0       ;
analysisSettings.stopTolForces          =   1e-5    ;
analysisSettings.fluidProps = { rhoA ; nuA ; 'windVel' } ;

otherParams

otherParams.problemName =  'simplePropeller' ;
otherParams.plots_format = 'vtk' ;

Run ONSAS

[ matUs, ~ ] = ONSAS( materials, elements, boundaryConds, initialConds, mesh, analysisSettings, otherParams ) ; 

Verification

numerical time vector is given by:

timeVec = linspace(0, analysisSettings.finalTime, size(matUs, 2) ) ;

numerical rotation angle for nodal moment case is:

dofAngleXnode1 = 2 ;
angleXnode1Numeric = -matUs(dofAngleXnode1,:) ;

analytical rotation angle is:

angleXnode1Analytic = angleXnode1(timeVec) ;

verifBoolean = norm( angleXnode1Numeric - angleXnode1Analytic )  ...
                    < ( norm( angleXnode1Analytic ) * 5e-2 ) ;

Plots

lw = 2.0 ; ms = 10; plotfontsize = 22 ; spanPlotTime = 2 ;
fig1 = figure(1) ;
plot( timeVec(1:spanPlotTime:end), angleXnode1Analytic(1:spanPlotTime:end) ,'b-x' , 'linewidth', lw,'markersize',ms )
hold on, grid on
plot( timeVec(1:spanPlotTime:end), angleXnode1Numeric(1:spanPlotTime:end), 'ko' , 'linewidth', lw,'markersize',ms )
labx = xlabel('time(s)');   laby = ylabel('$\theta_x node 1$') ;
legend('analytic','numeric', 'location','North')
set(gca, 'linewidth', 1.2, 'fontsize', plotfontsize )
set(labx, 'FontSize', plotfontsize); set(laby, 'FontSize', plotfontsize) ;
title('simple propeller test')
if length(getenv('TESTS_RUN')) > 0 && strcmp( getenv('TESTS_RUN'), 'yes')
  fprintf('\ngenerating output png for docs.\n')
  print(fig1, 'output/verifPropeller.png','-dpng')
else
  fprintf('\n === NOT in docs workflow. ===\n')
end