Edit on GitHub

Aerodynamic non-linear cantilever beam example

Octave script

In this tutorial, the aerodynamic non-linear cantilever beam example is solved using ONSAS. The aim of this problem is to validate aerodynamic steady wind loads applied to a cantilever beam undergoing small strains. The aerodynamic force modification due to the beam deformation is considered (drag reconfiguration). Given the aforementioned characteristics and under the hypothesis of small displacements regime a semi-analytic solution is available.

The beam is submitted to a uniform air wind velocity field $V_a$, at 20 degrees and atmospheric pressure, along axis $y$. Due to revolution symmetry of the problem lift $c_l$ and torsional moment $c_m$ coefficients are null. A drag coefficient $c_d=1.2$ is extracted from this reference. The beam has a length $L$ and a circular solid cross section with diameter $d$ as it is shown in the following figure:

general dimensions sketch

Small displacements 2D case

Static semi-analytic solution

The wind load forces of a generic cross section can be derived within the quasi-steady-theory. Considering a cross section located at $x$, then the projected wind velocity into the transverse deformed plane is then $V_p=V_acos(\theta_z)$ (the axial drag is neglected). Subsequently a drag force per unit of length $F_d= \frac{1}{2} \rho d c_d ||V_p||^2$ with $\frac{V_p}{||V_p||}$ direction is applied. In order to link the force $F_d$ with the beam deflection, the uniform distributed force along $y$ is computed as $F_y=F_d.cos(\theta_z)$. This leads to the following third order differential equation:

  • \[EI_{zz} \frac{\partial ^3 \theta_z}{\partial x ^3} = -q_0 c_d\cos^3(\theta_z)\]

in which $q_0 = \frac{1}{2} \rho_f d ||V_a||^2$ and the air density is $\rho_f = 1.225$ kg/m$^3$.

Numerical solution

Before defining the structs, the workspace is cleaned and ONSAS directory is added:

close all, clear all ; addpath( genpath( [ pwd '/../../src'] ) );

material and geometrical parameters are defined:

E = 1e9 ;  nu = 0.3 ; rho = 1800 ; G = E / (2 * (1+nu)) ;
l = 10 ; d = l/100 ; J = pi * d ^ 4 / 32 ; Iyy = J / 2 ; Izz = Iyy ;

the fluid properties are:

rhoA = 1.225 ; nuA = 1.6e-5;

next the number of frame elements for all cases is set:

numElements = 10 ;

MEB parameters

materials

Since the example contains only one linear Euler Bernoulli element the fields of the materials struct will have only one entry. Although, the constitutive behavior law selected is Saint-Venant-Kirchhoff:

materials.hyperElasModel  = 'linearElastic' ;
materials.hyperElasParams = [ E nu ]        ;

note that the use of this linear elastic element guarantees the left hand side of the differential equation stated above.

elements

Two different types of elements are considered, node and frames. The nodes will be assigned in the first entry (index $1$) and the beam at the index $2$. The elemType field is then:

elements(1).elemType = 'node'  ;
elements(2).elemType = 'frame' ;

The node has not cross section geometry to assign (an empty array is automatically set). The solid circular cross section is preset in ONSAS, and to load it just use:

elements(2).elemCrossSecParams{1,1} = 'circle' ;
elements(2).elemCrossSecParams{2,1} = d        ;

Now the element aerodynamic properties are defined. First the drag coefficient function located at the folder's example is declared into userDragCoef field as:

boundaryConds

Only one welded (6 degrees of freedom are set to zero) boundary condition (BC) is considered:

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

mesh parameters

The coordinates of the mesh nodes are given by the matrix:

mesh.nodesCoords = [ (0:(numElements))'*l/numElements  zeros(numElements+1,2) ] ;

The connectivity is introduced using the conecCell. Each entry of the cell contains a vector with the four indexes of the MEBI parameters, followed by the indexes of nodes that compose the element (node connectivity). For didactical purposes each element entry is commented. First the cell is initialized:

mesh.conecCell = { } ;

then the first welded node is defined with material (M) zero since nodes don't have material, the first element (E) type (the first entry of the elements struct), and (B) is the first entry of the the boundaryConds struct. For (I) no non-homogeneous initial condition is considered (then zero is used) and finally the node is assigned:

mesh.conecCell{ 1, 1 } = [ 0 1 1  1 ] ;

Next the frame elements MEB parameters are set. The frame material is the first material of materials struct, then $1$ is assigned. The second entry of the elements struct correspond to the frame element employed, so $2$ is set. Finally no BC is required for this element, then $0$ is used. Consecutive nodes build the element so then the mesh.conecCell is:

for i=1:numElements,
  mesh.conecCell{ i+1,1 } = [ 1 2 0  i i+1 ] ;
end

initial Conditions

Any non-homogeneous initial conditions are considered, then an empty struct is set:

initialConds = struct() ;

analysisSettings

First the wind velocity and the fluid properties are set into __ field of analysisSettings struct. This will apply a external wind loads for each element with elemTypeAero field into the elements struct. The name of the wind velocity function must located on the same example path is: 

analysisSettings.fluidProps = {rhoA; nuA; 'windVelNonLinearStatic'} ;

Inside that function a linear velocity $v_a = 30*t$ is declared. The final time will be set to $1$ in order to achieve 30 m/s. The geometrical non-linear effects are considered in this case to compute the aerodynamic force. As consequence the wind load forces are computed on the deformed configuration. The field geometricNonLinearAero into analysisSettings struct is then set to:

analysisSettings.geometricNonLinearAero = true;

note that if this boolean is not declared ONSAS will automatically assign it as true. since this problem is static, then a N-R method is employed. The convergence of the method is accomplish with ten equal load steps. The time variable for static cases is a load factor parameter that must be configured into the windVel.m function. A linear profile is considered for ten equal velocity load steps as:

analysisSettings.deltaT        =   0.1           ;
analysisSettings.finalTime     =   1             ;
analysisSettings.methodName    = 'newtonRaphson' ;

Next the maximum number of iterations per load(time) step, the residual force and the displacements tolerances are set to: 

analysisSettings.stopTolDeltau =   1e-6          ;
analysisSettings.stopTolForces =   1e-6          ;
analysisSettings.stopTolIts    =   40            ;

otherParams

The name of the problem and vtk format output are selected: 

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

ONSAS software is executed for the parameters above defined and the displacement solution of each load(time) step is saved into matUsSD matrix:

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

Verification

DifferentialEquations.jl (reconfiguration) solution.

DiffEq.jl solves the third order ordinary differential equation for this case by executing DiffEq.jl script. Then assembleJuliaSol.m script function is executed to build the julia solution with mesh and elements struct as:

% [dSolJulia] = assembleJuliaSol(elements,mesh) ;

Then the the relevance linear and angular displacements are extracted using:

% ydefJulia = dSolJulia(3:6:end)              ;
% thetaZdefJulia = dSolJulia(6:6:end)         ;
% xdefJulia = linspace(0,l,length(ydefJulia)) ;

Numeric solution

The numerical solution is computed:

xref    = mesh.nodesCoords(:,1)       ;
yref    = mesh.nodesCoords(:,2)       ;
zref    = mesh.nodesCoords(:,3)       ;
xdefNum = xref + matUsSD(1:6:end,end) ;
ydefNum = yref + matUsSD(3:6:end,end) ;
thetaZdefNum = matUsSD(6:6:end,end)   ;

Plot verification

The plot parameters are:

lw = 2 ; ms = 12 ;
labelTitle = [' Validating solution with ' num2str(numElements) ' elements' ] ;
axislw = 2 ; axisFontSize = 20 ; legendFontSize = 15 ; curveFontSize = 15 ;    
folderPathFigs = './output/figs/' ;
mkdir(folderPathFigs) ;

The linear $u_y$ displacements verification is plotted executing: 

fig1 = figure(1) ;
hold on, grid on
plot(xref      , ydefNum  , 'bo' , 'linewidth', lw, 'markersize', ms+5   );
% plot(xdefJulia , ydefJulia, 'b-' , 'linewidth', lw, 'markersize', ms     );
legend('y numeric SD', 'y semi-analytic SD' )
labx=xlabel(' x[m] ');    laby=ylabel('y[m]');
set(legend, 'linewidth', axislw, 'fontsize', legendFontSize, 'location','northWest' ) ;
set(gca, 'linewidth', axislw, 'fontsize', curveFontSize ) ;
set(labx, 'FontSize', axisFontSize); set(laby, 'FontSize', axisFontSize) ;
namefig1 = strcat(folderPathFigs, 'linDispSD.png') ;
print(fig1, namefig1,'-dpng') ;
plot check angular displacements

The angular $\theta_z$ displacements verification is plotted executing: 

fig2 = figure(2) ;
hold on, grid on
plot(xref,      thetaZdefNum,      'bo' , 'linewidth', lw,'markersize', ms+5   );
% plot(xdefJulia, thetaZdefJulia,    'b-' , 'linewidth', lw, 'markersize', ms    );
legend('\theta_z numeric SD', '\theta_z semi-analyitc SD')
labx=xlabel(' x[m] ');    laby=ylabel('Angle[rad]');
% title (labelTitle)
set(legend, 'linewidth', axislw, 'fontsize', legendFontSize, 'location','northWest' ) ;
set(gca, 'linewidth', axislw, 'fontsize', curveFontSize ) ;
set(labx, 'FontSize', axisFontSize); set(laby, 'FontSize', axisFontSize) ;
namefig2 = strcat(folderPathFigs, 'angDispSD.png') ;
print(fig2, namefig2,'-dpng')
plot check angular displacements

Large displacements 2D case

Now a large displacements 2D case is solved. The solution is computed using the co-rotational beam element formulation proposed in this reference

Numerical solution static case

MEBI parameters

materials

In order to reproduce large displacements results the materials struct is then changed to: 

materials.hyperElasModel  = '1DrotEngStrain' ;
materials.hyperElasParams = [ 1e6 nu ]       ;
materials.density         = rho              ;

elements

The element tangent matrices of the consistent inertial force vector are taking into account by the following boolean:

elements(2).massMatType = 'consistent' ;

otherParams

The name of this case problem is:

otherParams.problemName = 'nonLinearCantileverLDStatic' ;

ONSAS software is executed for the parameters above defined and the displacement solution of each load(time) step is saved into matUsLDStatic matrix:

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

Numerical solution dynamic case

Next a dynamic example considering large displacements motion is addressed to test the convergence of the dynamic solution disregarding any artificial damping.

analysisSettings

For such propose the wind velocity function name is now: 

analysisSettings.fluidProps = {rhoA; nuA; 'windVelNonLinearDynamic2D'} ;

Inside that function a ramp velocity profile $v_a(t) = 5*t*(t<6.6) + 5*t*(t>=6.6)$ is declared. This is an abrupt wind velocity load from 0 to $7$ m/s in $10$ s .

Regarding the integration time method scheme, a classic Newmark trapezoidal is set as: 

analysisSettings.deltaT     =  1        ;
analysisSettings.finalTime  =  200      ;
analysisSettings.methodName = 'newmark' ;

otherParams

The name of this case problem is:

otherParams.problemName = 'nonLinearCantileverLDDynamic2D' ;

ONSAS software is executed for the parameters above defined and the displacement solution for each time step is saved into matUsLDDynamic matrix:

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

Verification

The numerical dynamic convergence to the static solution is then verified . The degree of freedom selected for such task is $u_y(t)$ of node A.

Static solution.

Extract static numerical time history displacements $u_y$ of node A. First the selected degree of freedom is:

nodeA = numElements + 1                ;
relativeDofUyA = 3                     ;
dofUyA = (nodeA -1)*6 + relativeDofUyA ;

then node A $u_y$ time history accessed by:

UyAStaticSol = matUsLDStatic(dofUyA,:) ;

Dynamic solution.

Extract dynamic numerical solution as follows:

UyADynamicSol = matUsLDDynamic(dofUyA,:) ;

next, the time vector is given by:

timVecLD = linspace(0, analysisSettings.finalTime, size(matUsLDDynamic,2) ) ;

Verification Plot

Create folder to save figures

folderFigs = strcat('./output/', 'figs/') ;
mkdir(folderFigs) ;

The linear $u_y$ displacements verification of node A is finally plotted executing: 

fig3 = figure(3) ;
hold on,  grid on
% legend first point plot
plot(timVecLD(1), UyADynamicSol(1),...
     'color', 'b', 'linewidth', lw, 'linestyle', '-','markersize', ms, 'marker', 'o')
% static solution plot
plot(timVecLD   , UyAStaticSol(end)*ones(length(timVecLD)), 'k:' , 'linewidth', lw, 'markersize', ms     );
% markers plot
plot(timVecLD(1:8:end), UyADynamicSol(1:8:end),...
     'color', 'b', 'linewidth', lw, 'linestyle', 'none','markersize', ms, 'marker', 'o')
% continium line plot
plot(timVecLD, UyADynamicSol,...
     'color', 'b', 'linewidth', lw, 'linestyle', '-', 'marker', 'none')
legend('dynamic LD', 'static LD' )
labx=xlabel('x[m]');    laby=ylabel('y[m]');
set(legend, 'linewidth', axislw, 'fontsize', legendFontSize, 'location','northEast' ) ;
set(gca, 'linewidth', axislw, 'fontsize', curveFontSize ) ;
set(labx, 'FontSize', axisFontSize); set(laby, 'FontSize', axisFontSize) ;
namefig3 = strcat(folderPathFigs, 'uyA.png') ;
print(fig3, namefig3,'-dpng') ;
plot check angular displacements

Large displacements 3D case

A large displacements dynamic 3D case is presented as follows. This example is inspired on Vortex shedding exposed at Youtbue Video

MEBI parameters

materials

In order to reproduce large displacements results the materials struct is then changed to: 

materials.hyperElasParams = [ 1e8 nu ]       ;

analysisSettings

Regarding the integration time method scheme, a classic $\alpha-HHT$ method is employed. This method is more stable numerically than Newmark, the keen reader is refereed to this reference:

analysisSettings.methodName = 'alphaHHT';

the simulation time is defined such that: 

analysisSettings.deltaT     =  .2  ;
analysisSettings.finalTime  =  120 ;

The emulation of the vortex shedding vibration is generated by a synthetic wind velocity composed by two sinusoidal velocities. A low frequency $Vy_a$ along the mean flow direction $y$ and then a high frequency component $Vz_a$ along $z$. The high frequency component is selected to produce resonance effects between the flow and the beam, thus the high frequency velocity is selected equal to the first mode bending:

freqBendingFirstMode = (1.875)^2 * sqrt( materials.hyperElasParams(1) * Iyy / (materials.density * (pi * d^4 / 4) * l^4) ) ;
analysisSettings.fluidProps = {rhoA; nuA; 'windVelNonLinearDynamic3D'} ;

The velocity function componentes are assembled:

timeVecLD3d = linspace(0,analysisSettings.finalTime, ceil(analysisSettings.finalTime / analysisSettings.deltaT + 1) ) ;
windVelY = [] ; windVelZ = [] ;
for timeIndex = timeVecLD3d
    windVelVecTimeIndex = feval(analysisSettings.fluidProps{3}, 0, timeIndex) ;
    windVelY = [windVelY windVelVecTimeIndex(2) ] ;
    windVelZ = [windVelZ windVelVecTimeIndex(3) ] ;
end

otherParams

The name of this case problem is:

otherParams.problemName = 'nonLinearCantileverLD3D' ;

ONSAS software is executed for the parameters above defined and the displacement solution for each time step is saved into matUsLD3D matrix:

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

Plots

The wind velocity profile implemented is plotted executing: 

fig4 = figure(4) ;
hold on,  grid on
% legend first point plot
plot(timeVecLD3d(1), windVelY(1),...
     'color', 'b', 'linewidth', lw, 'linestyle', '-','markersize', ms, 'marker', 'o')
plot(timeVecLD3d(1), windVelZ(1),...
     'color', 'r', 'linewidth', lw, 'linestyle', '-','markersize', ms, 'marker', '^')
% markers plot
plot(timeVecLD3d(1:10:end), windVelY(1:10:end),...
     'color', 'b', 'linewidth', lw, 'linestyle', 'none','markersize', ms, 'marker', 'o')
plot(timeVecLD3d(1:17:end), windVelZ(1:17:end),...
     'color', 'r', 'linewidth', lw, 'linestyle', 'none','markersize', ms, 'marker', '^')
% continium line plot
plot(timeVecLD3d, windVelY,...
     'color', 'b', 'linewidth', lw, 'linestyle', '-', 'marker', 'none')
plot(timeVecLD3d, windVelZ,...
     'color', 'r', 'linewidth', lw, 'linestyle', '-', 'marker', 'none')
legend('Va_y', 'Va_z' )
labx=xlabel('t[s]');    laby=ylabel('V_a[m/s]');
set(legend, 'linewidth', axislw, 'fontsize', legendFontSize, 'location','northEast' ) ;
set(gca, 'linewidth', axislw, 'fontsize', curveFontSize ) ;
set(labx, 'FontSize', axisFontSize); set(laby, 'FontSize', axisFontSize) ;
namefig4 = strcat(folderPathFigs, 'windVel3D.png') ;
axis([0,50])
print(fig4, namefig4,'-dpng') ;
plot check angular displacements

Then $u_y$ of node A is computed using: 

UyADynamicSol3D = matUsLD3D(dofUyA,:) ;

analogosuly $u_z$ of A node is: 

UzADynamicSol3D = matUsLD3D(dofUyA + 2,:) ;

Open figure and plot

fig5 = figure(5) ;
hold on,  grid on
% legend first point plot uy 
plot(timeVecLD3d(1), UyADynamicSol3D(1),...
     'color', 'b', 'linewidth', lw, 'linestyle', '-','markersize', ms, 'marker', 'o')
% legend first point plot uz
plot(timeVecLD3d(1), UzADynamicSol3D(1),...
     'color', 'r', 'linewidth', lw, 'linestyle', '-','markersize', ms, 'marker', '^')
% markers plot uy
plot(timeVecLD3d(1:13:end), UyADynamicSol3D(1:13:end),...
     'color', 'b', 'linewidth', lw, 'linestyle', 'none','markersize', ms, 'marker', 'o')
% continium line plot uy
plot(timeVecLD3d, UyADynamicSol3D,...
     'color', 'b', 'linewidth', lw, 'linestyle', '-', 'marker', 'none')
% markers plot uz
plot(timeVecLD3d(1:23:end), UzADynamicSol3D(1:23:end),...
     'color', 'r', 'linewidth', lw, 'linestyle', 'none','markersize', ms, 'marker', '^')
% continium line plot uz
plot(timeVecLD3d, UzADynamicSol3D,...
     'color', 'r', 'linewidth', lw, 'linestyle', '-', 'marker', 'none')
legend('U_y node A', 'U_z node A' )
labx=xlabel('x[m]');    laby=ylabel('Dispalcements[m]');
set(legend, 'linewidth', axislw, 'fontsize', legendFontSize, 'location','northEast' ) ;
set(gca, 'linewidth', axislw, 'fontsize', curveFontSize ) ;
set(labx, 'FontSize', axisFontSize); set(laby, 'FontSize', axisFontSize) ;
namefig5 = strcat(folderPathFigs, 'uA3D.png') ;
print(fig5, namefig5,'-dpng') ;
plot check angular displacements

Finally a GIF to illustrate the motion amplitude is subsequently presented:

plot check angular displacements