5.3 Damping#
Models used to represent dissipation at the local material level and at the global system level should typically be different. Simple viscous behavior is very often not appropriate to describe material damping while a viscous model is appropriate in the normal mode model format (see details in Ref. [34]). This section discusses typical damping models and discusses how piece-wise Rayleigh damping is implemented in SDT.
5.3.1 Viscous damping in the normal mode model form#
In the normal mode form, viscous damping is represented by the modal damping matrix Γ which is typically used to represent all the dissipation effects at the system level.
Models with modal damping assume that a diagonal Γ is sufficient to represent dissipation at a system level. The non-zero terms of Γ are then usually expressed in terms of damping ratios Γjj = 2 ζj ωj. The damping ratio ζj are accepted by most SDT functions instead of a full Γ. The variable name damp is then used instead of ga in the documentation.
For a model with modal damping, the matrices in (6.81) are diagonal so that the contributions of the different normal modes are uncoupled and correspond exactly to the spectral decomposition of the model (see cpx for the definition of complex modes). The rational fraction expression of the dynamic compliance matrix (transfer from the inputs {u} to displacement outputs {y}) takes the form
| (5.5) |
where the contribution of each mode is characterized by the pole frequency ωj, damping ratio ζj, and the residue matrix Tj (which is equal to the product of the normal mode output shape matrix {c φj} by the normal mode input shape matrix {φjTb}).
Modal damping is used when lacking better information. One will thus often set a uniform damping ratio (ζj=1% or damp = 0.01) or experimentally determined damping ratios that are different for each pole (po=ii_pof(po,3); damp=po(:,2);).
Historically, modal damping was associated to the proportional damping model introduced by Lord Rayleigh which assumes the usefulness of a global viscously damped model with a dynamic stiffness of the form
| (5.6) |
While this model indeed leads to a modally damped normal mode model, the α and β coefficients can only be adjusted to represent physical damping mechanisms over very narrow frequency bands. The modal damping matrix thus obtained writes
| (5.7) |
which leads to damping ratios
| (5.8) |
Mass coefficient α leads to high damping ratios in the low frequency range. Stiffness coefficient β leads to a damping ratio linearly increasing with the frequency.
Using a diagonal [Γ] can introduce significant errors when normal mode coupling through the spatial distribution of damping mechanisms is possible. The condition
| (5.9) |
proposed by Hasselman [35], gives a good indication of when modal coupling will not occur. One will note that a structure with a group of modes separated by a few percent in frequency and levels of damping close to 1% does not verify this condition. The un-coupling assumption can however still be applied to blocks of modes [19].
A normal mode model with a full Γ matrix is said to be non-proportionally damped and is clearly more general/accurate than the simple modal damping model. The SDT leaves the choice between the non-proportional model using a matrix ga and the proportional model using damping ratio for each of the pole frequencies (in this case one has ga=2*diag(damp.*freq) or ga=2*damp*diag(freq) if a scalar uniform damping ratio is defined).
For identification phases, standard approximations linked to the assumption of modal damping are provided by (id_rc, id_rm and res2nor), while id_nor provides an original algorithm of the determination of a full Γ matrix. Theoretical aspects of this algorithm and details on the approximation of modal damping are discussed in [19]).
5.3.2 Viscous damping in finite element models#
Standard damped finite element models allow the incorporation of viscous and structural damping in the form of real C and complex K matrices respectively.
fe_mk could assemble a viscous damping matrix with user defined elements that would support matrix type 3 (viscous damping) using a call of the form
fe_mk(MODEL,'options',3) (see section 7.16 for new element creation). Viscous damping models are rarely appropriate at the finite element level [34], so that it is only supported by celas and cbush elements. Piece-wise Rayleigh damping where the viscous damping is a combination of element mass and stiffness on element subsets
| (5.10) |
is supported as follows. For each material or group that is to be considered in the linear combination one defines a row entry giving GroupId MatId AlphaS BetaS (note that some elements may be counted twice if they are related to a group and a material entry). One can alternatively define ProId as a 5th column (useful for celas element that have no matid). Note that each line is separately accounted for, so that duplicated entries or multiple references to same GroupId, MatId or ProId will also be combined. For example
model=demosdt('demogartfe'); model=stack_set(model,'info','Rayleigh', ... [10 0 1e-5 0.0; ... % Elements of group 10 (masses) 9 0 0.0 1e-3; ... % Elements of group 9 (springs) 0 1 0.0 1e-4; ... % Elements with MatId 1 0 2 0.0 1e-4]); % Elements with MatId 2 % Note that DOF numbering may be a problem when calling 'Rayleigh' % See sdtweb simul#feass for preferrred assembly in SDT c=feutilb('Rayleigh',model); figure(1);spy(c);dc=fe_ceig(model,[1 5 20 1e3]);cf=feplot(model,dc);
Such damping models are typically used in time integration applications. Info,Rayleigh entries are properly handled by Assemble commands.
You can also provide model=stack_set(model,'info','Rayleigh',[alpha beta]).
Note that in case of Rayleigh damping, celas element viscous damping will also be taken into account.
5.3.3 Hysteretic damping in finite element models#
Structural or hysteretic damping represents dissipation by giving a loss factor at the element level leading to a dynamic stiffness of the form
| (5.11) |
The name loss factor derives from the fact that η is equal to the ratio of energy dissipated for one cycle Ed=∫0Tσ є′ dt by 2π the maximum potential energy Ep=1/2E.
If dissipative materials used have a loss factor property, these are used by Assemble commands with a desired matrix type 4. If no material damping is defined, you can also use DefaultZeta to set a global loss factor to eta=2*DefaultZeta.
Using complex valued constitutive parameters will not work for most element functions. Hysteretic damping models can thus be assembled using the Rayleigh command shown above (to assemble the imaginary part of K rather than C or using upcom (see section 6.5). The following example defines two loss factors for group 6 and other elements of the Garteur FEM model. Approximate damped poles are then estimated on the basis of real modes (better approximations are discussed in [36])
Up=upcom('load GartUp'); cf=feplot(Up); Up=fe_case(Up,'parReset', ... 'Par k','Constrained Layer','group 6', ... 'Par k','Main Structure','group~=6');% type cur min max vtype par = [ 1 1.0 0.1 3.0 1 ; ... 1 1.0 0.1 3.0 1 ]; Up=upcom(Up,'ParCoef',par);
% assemble using different loss factors for each parameter B=upcom(Up,'assemble k coef .05 .01'); [m,k]=upcom(Up,'assemble coef 1.0 1.0'); Case=fe_case(Up,'gett');
% Estimate damped poles on real mode basis def=fe_eig({m,k,Case.DOF},[5 20 1e3]); mr=def.def'mdef.def; % this is the identity cr=zeros(size(mr)); kr=def.def'kdef.def+i*(def.def'Bdef.def); dr=fe_ceig({mr,cr,kr,[]});dr.def=def.def*dr.def;dr.DOF=def.DOF; cf.def=dr
Note that in this model, the poles λj are not complex conjugate since the hysteretic damping model is only valid for positive frequencies (for negative frequencies one should change the sign of the imaginary part of K).
Given a set of complex modes you can compute frequency responses with res2xf, or simply use the modal damping ratio found with fe_ceig. Continuing the example, above one uses
Up=fe_case(Up,'Dofload','Point loads',[4.03;55.03], ... 'SensDof','Sensors',[4 55 30]'+.03); Sens=feutilb('placeindof',def.DOF,fe_case(Up,'sens')); Load=fe_load(Up); ind=find(dr.data(:,1)>5); % flexible modes% Standard elastic response with modal damping f=linspace(5,60,2048); d1=def; d1.data(7:20,2)=dr.data(ind,2); nor2xf(d1,Up,f,'hz iiplot "Normal" -reset -po');
% Now complex modes RES=struct('res',[],'po',dr.data(ind,:),'idopt',idopt('new')); RES.idopt.residual=2;RES.idopt.fitting='complex'; for j1=1:length(ind); % deal with flexible modes Rj=(Sens.cta*dr.def(:,ind(j1))) * ... % c psi (dr.def(:,ind(j1)).'*Load.def); % psi^T b RES.res(j1,:)=Rj(:).'; end
% Rigid body mode residual RES.res(end+1,:)=0; for j1=1:6; Rj=(Sens.ctadef.def(:,j1))(def.def(:,j1)'*Load.def); RES.res(end,:)=RES.res(end,:)+Rj(:).'; end res2xf(RES,f,'hz iiplot "Res2xf"');
damp=dr.data(ind,2); d2=def;d2.data(7:20)=sqrt(real(d2.data(7:20).^2)).sqrt(1+idamp*2); nor2xf(d2,Up,f,'hz iiplot "Hysteretic"'); iicom('submagpha');
Note that the presence of rigid body modes, which can only be represented as residual terms in the pole/residue format (see section 5.6), makes the example more complex. The plot illustrates differences in responses obtained with true complex modes, viscous modal damping or hysteretic modal damping (case where one uses the pole of the true complex mode with a normal mode shape). Viscous and hysteretic modal damping are nearly identical. With true complex modes, only channels 2 and 4 show a visible difference, and then only near anti-resonances.
To incorporate static corrections, you may want to compute complex modes on bases generated by fe2ss, rather than simple modal bases obtained with fe_eig.
The use of a constant loss factor can be a crude approximation for materials exhibiting significant damping. Methods used to treat frequency dependent materials are described in Ref. [37].