Fit of structural damping

Fit of structural damping#

Please note that this feature is not easy to use, and some iterations must be foreseen in order to end at satisfactory result.

The aim of this feature is to develop a method to fit the damping parameters in a HAWC2 model in such a way that desired damping ratios are obtained for specified eigenmodes. Further, it is the aim that the formulation can be used in both HAWC2 and HAWCStab2.

The method is described below in Section Formulation. It fits element stiffness matrices and saves them to file so that they can be used for bodies using the damping method “damping_file <damping file> ;” (where the <damping file> is generated by the method). This requires a special block inside the htc file which must be placed after the “new_htc_structure” block:

Obl.

Command name

Explanation

*

begin damping_fit ;

First line in damping fit block.

*

damping_file

Name of damping file. This file name MUST match the <damping file> used for the “damping_file” method in the “main_body” block.

twin_bodies

The two body names given as arguments will share damping properties. This is used e.g. to specify the same damping for all the blades on the rotor.

1. Body name for 1st twin.

2. Body name for 2nd twin.

*

cmd_solver

1. Command executed by HAWC2 which does the actual fitting, (e.g. python.exe damping_fit.py ;). If you rely on a virtual Python environment, make sure to activate this first before running HAWC2 within this environment. Within this Python environment the numpy and scipy packages are required to be installed.

*

mode

Repeated line for each mode to be fitted:

1. Mode number.

2. Damping ratio.

-

*

end damping_fit ;

Last line in damping fit block.

The example below shows the setup for fitting the damping of a single blade with requested damping ratios specified for the first six modes. 0.5% damping ratio (i.e. approx. 3% log.decr.) is requested for modes 1 and 2, 1% for modes 3 and 4, and 2% for modes 5 and 6. Note the use of the damping file (blade.dmp) in two locations which links the damping fit only to include blade damping in the fit. If other bodies were present in the example, the damping specified for those bodies would enter the total damping fit, but only the damping parameters for the blade will change the total damping.


;---------------------------------------------------------------------------
begin new_htc_structure;
;---------------------------------------------------------------------------
  begin main_body;
 name    blade1 ;    
 type    timoschenko ;
 nbodies   10 ;
 node_distribution  c2_def;
  damping_file blade.dmp ;
  begin timoschenko_input ;
  filename ./data/DTU_10MW_RWT_Blade_st.dat;
  set 1 1 ;        set subset
  end timoschenko_input;
 begin c2_def;       Definition of centerline (main_body coordinates)
  nsec 27 ;
  sec	 1	0.00000E+00	7.00600E-05	4.44089E-16	-1.45000E+01	;
    ..
    ..
  sec	27	-8.98940E-02	-3.33685E+00	8.63655E+01	3.42796E+00	;
  end c2_def ;                 
  end main_body;
;---------------------------------------------------------------------------
  begin orientation;
  begin base;
   body  blade1;
   inipos    0.0 0.0 0.0 ;     initial position of node 1
   body_eulerang 0.0 0.0 0.0;
  end base;
  end orientation;
;---------------------------------------------------------------------------
 begin constraint;  
  begin fix0; fixed to ground in translation and rotation of node 1
   body blade1;
  end fix0;
  end constraint;
end new_htc_structure;
;---------------------------------------------------------------------------
begin damping_fit ;
  damp_file blade.dmp ;
  cmd_solver C:\Users\anmh\Anaconda3\Scripts\conda.exe run python damping_fit.py ;
  mode 1 0.005 ;	Damping ratio of mode 1
  mode 2 0.005 ;			etc.
  mode 3 0.01 ;
  mode 4 0.01 ;
  mode 5 0.02 ;
  mode 6 0.02 ;
end damping_fit ;
;---------------------------------------------------------------------------

A fully functional example is available for download at https://​gitlab​.windenergy​.dtu​.dk​/HAWC2Public​/examples​/​-​/tree​/master​/hawc2​/structure​/damping​_fit​/IEA​_15MW​_RWT.

Formulation#

The linearised HAWC2 EOMs are of the usual 2nd order form:

(1)#\[\mat{M}\,\ddot{\vec{x}} + \mat{C}\,\dot{\vec{x}} + \mat{K}\,{\vec{x}} = \vec{0}\]

The solution to the undamped eigenvalue problem (\(\mat{C=0}\)) are defined by the eigenvectors \(\mat{\Gamma}\) and diagonal eigenfrequency matrix \(\mat{\Omega}\). The eigensolution fulfills the identity \(\mat{M}\mat{\Gamma}\mat{\Omega}^2=\mat{K}\mat{\Gamma}\). By using the eigenvectors as basis, \(\vec{x}\) can be transformed as \(\vec{x}(t) = \mat{\Gamma}\,\vec{\alpha}(t)\). By using the above relations, (1) can be manipulated as:

(2)#\[\ddot{\vec{\alpha}} + \mat{\Gamma}^{-1}\mat{M}^{-1}\mat{C}\mat{\Gamma}\,\dot{\vec{\alpha}} + \mat{\Omega^2}\,{\vec{\alpha}} = \vec{0}\]

Note that the undamped part of (2) is a diagonal system, and that the total set of equations can de uncoupled if the damping matrix part is also a diagonal matrix. If we choose this matrix as \(2\,\mat{\zeta}\mat{\Omega}\) (\(\mat{\zeta}\) is a diagonal matrix), then the system damping matrix \(\mat{C}\) can be calculated as

(3)#\[\mat{C}=2\,\mat{\Gamma}\mat{M}\mat{\zeta}\mat{\Omega}\mat{\Gamma}^{^-1}\]

Unfortunately, such a damping matrix cannot directly be used in neither HAWC2 nor in HAWCStab2, so something else must be done. Instead, the damping of one mode at a time is formulated as function of element damping matrices:

(4)#\[\ddot{\vec{\alpha_i}} + (\mat{\gamma_i}^{T}\mat{M}\mat{\gamma_i})^{-1}(\mat{\gamma_i}^{T}\mat{C}\mat{\gamma_i})\,\dot{\vec{\alpha_i}} + \Omega_i^2\,{\vec{\alpha_i}} = \vec{0}\]

where all variables with sub-script i relate to the i’th eigenmode. The damping coefficient in ((4)) must then fulfill the equation

(5)#\[(\mat{\gamma_i}^{T}\mat{M}\mat{\gamma_i})^{-1}(\mat{\gamma_i}^{T}\mat{C}\mat{\gamma_i}) = 2\zeta_i\,\Omega_i\]

The system damping matrix, \(\mat{C}\), is assembled based on element damping matrices \(\mat{c}_j\) (for the j’th element), where the element damping matrices are defined as having the same eigenvectors as the element stiffness matrices. By using this formulation, the structure only dissipates energy when it is deformed and not during rigid body motion.

\[\mat{c}_j = \mat{v}_j\mat{x}_j \mat{v}_j^T\]

where \(\mat{v}_j\) is the eigenvectors of the j’th stiffness matrix (or rather the six eigen vectors that have non-zero eigenvalues) and \(\mat{x}_j\) is a \(6\times6\) diagonal matrix containing the unknown damping parameters.

For each eigenmode that needs to be fitted, (5) provides one equation that needs to be fulfilled, and the unknowns are the six element damping parameters, \(diag(\mat{x}_j)\), for all elements. Further, in order for the element damping matrices to be positive semi-definite, \(\mat{x}_j \geq 0\) for all diagonal components \(\mat{x}_j\) and for all elements (all j). The equations in ((5)) are linear in \(\vec{x}\) and the (one of many!) solution is found by solving the optimization problem

(6)#\[\min_{wrt. \vec{x}}\left ( \left |\mat{W}\left ( \mat{A}\,\vec{x} - \vec{\zeta}\right ) \right |^2 \right ),\,\,s.t.\,\, \vec{x} \geq 0\]

where \(\vec{x}\) are the element damping parameters for all elements collected in a vector, \(\vec{\zeta}\) are the (user-specified) damping ratios prescribed for the individual modes, \(\mat{A}\) are the coefficients to \(\vec{x}\) in accordance with ((5)), and \(\mat{W}\) is a diagonal weighting matrix which is included in order to weigh the individual eigenmodes in the optimization.

HAWC2 implementation#

Currently, the solution of ((6)) is handled by an external call to a python script outside of HAWC2. This means that HAWC2 calculates the matrices and vectors in ((6)) and exports those to a binary file (currently named ’dfit_a.bin’). This binary file is then handled in a Python script that reads the system, solves ((6)) for \(\vec{x}\) and writes back the solution to file (currently named ’dfit_x.bin’). This solution is then read back into HAWC2 and the resulting element damping matrices are calculated and written to file for subsequent use in HAWC2 simulations.

Note that even though only a few of the total number of eigenmodes have prescribed damping ratios (specified in the htc-file), all eigenmodes are included in the outputted binary file, however, the components in \(\mat{W}\) associated with non-prescribed modes are all set to zero.

The Python script is listed below from which the individual binary file formats can be deduced, if needed. This script is part of the distributed HAWC2 files.


#-------------------------------------------------------------------------
# -*- coding: utf-8 -*-
"""
Damping fit for HAWC2

This script finds the parameters for structural HAWC2 damping type
"damp_file", based on the mode damping matrix, A, calculated by HAWC2
and written to file "dfit_a.bin".

Each row of the A matrix corresponds to a mode shape in the HAWC2 model, 
ordered in increasing order of eigenfrequency, i.e. first row corresponds to 
the mode with the lowest eigenfrequency. By multiplication with the damping 
parameter vector, x, gives the damping rati0 vector, d = (A*x).

The purpose of this script is then to find the best fit of x which gives the
specified damping ratio for the individual modes using the constraint that 
x>0 for all x.

The A matrix contains all modes, and not all modes can be fitted for any 
damping level. Normally the first (say 10) modes are of interest. This is 
handled by the weighting vector, w, below. See code below for further details.

"""

import numpy as np
import struct
from scipy.optimize import nnls

def damping_fit():
  
  wmin = 1.e-6
  
  # Read A matrix from file
  f=open('dfit_a.bin','rb')
  (nr,nc) = struct.unpack('ii',f.read(8))
  ntot = nr*nc
  data = np.zeros(ntot,dtype=np.dtype('f8'))
  for i in range(ntot):
    (data[i],) = struct.unpack('d',f.read(8))
  A = np.reshape(data,[nr,nc], order='F')

  # Read target damping for optimization
  d = np.zeros(nr,dtype=np.dtype('f8'))
  for i in range(nr):
    (d[i],) = struct.unpack('d',f.read(8))
  w = np.zeros(nr,dtype=np.dtype('f8'))
  for i in range(nr):
    (w[i],) = struct.unpack('d',f.read(8))
    if w[i] == 0.0:
      w[i] = wmin
  f.close()  
    
  # Solve
  res = nnls(np.matmul(np.diag(w),A), np.matmul(np.diag(w),d))

  # Write results to file  
  f = open('dfit_x.bin','wb')
  f.write(np.array([nc,1],dtype='i4'))
  f.write(res[0])
  f.close()

  # Check solution
  dfit = np.matmul(A,res[0])
  print(('*'+'{0:1s}'*36+'*').format('*'))
  print('*{:^36s}*'.format('Damping fit result'))
  print(('*'+'{0:1s}'*36+'*').format('*'))
  print('*{:^8s}{:^14s}{:^14s}*'.format('Mode','Target','Fit'))
  for i in range(nr):
    if w[i] > wmin:
      print('*{:^8d}{:^14.3e}{:^14.3e}*'.format(i+1,d[i],dfit[i]))
  print(('*'+'{0:1s}'*36+'*').format('*'))
      
  return res

#-----------
# DO IT....
#-----------
res = damping_fit()
#-------------------------------------------------------------------------

Usage considerations#

Consistently use matched input and damping file#

The result of the structural damping fitting procedure is a main body element damping matrix file that will match the user defined damping for the relevant modes. This file is specific for a given combination of nodes, number of bodies and structural input (st-file). If any changes are made in either of these inputs the element damping matrix file will have to be redefined based on the procedure outlined here. Users are especially cautioned to carefully track that the number of bodies used for generating the damping fit is also the same number of bodies used in subsequent simulations.

Tune on full or main body only models#

It is more likely to obtain a good damping fit for many frequencies when tuning the damping for a HAWC2 model containing only the main body of interest. When a model with several main bodies is used (tower, blades, etc) the optimisation problem becomes inherently more difficult to solve. When using multiple main bodies, make sure to verify that the targeted damping ratios in the damping_fit section relate to the total systems modes (1st and 2nd modes likely to be the tower, etc), as opposed to when using a model that only contains the main body of interest. The user is responsible for tracking which mode number relates to which body. For example, fitting the damping for tower modes while only adjusting the damping coefficients related to the blades is not likely to give meaningful results. It is therefore recommended to only list/target mode numbers of the body at interest, and leave out the others (especially rigid body modes) in the damping_fit section.

Number of modes to target#

When fitting to a low number of modes a very good result can be expected. The more modes a user attempts to fit a damping value to, the more difficult the trade-off becomes. In those cases an advanced user could consider changing the weights \(\mat{W}\) in the example script damping_fit.py (defined as w, see above) to obtain a specific trade-off in which some modes are allowed to differ more compared to others with respect to the requested target values.