Computation and estimation of transfers

Contents Functions

3 Computation and estimation of transfers

The objective of this session is to illustrate equations given in the modal basis (??) and (3). For this, one first analyzes the frequency response then simulation the test process using time transients and some signal processing steps.

To help through various steps, one will use different MATLAB toolboxes

Structural Dynamics Toolbox (sdt) for FEM computations
Control Toolbox (control) for time response simulation (lsim function to integrate linear time invariant systems, alternatives would have been calling ode45 under MATLAB or fe_time with SDT)
Signal Processing Toolbox (signal) for the generation of anti-aliasing filters (butter function)

To get the data

On your PC, download the folder VibAM.zip from https://www.sdtools.com/contrib/VibAM21.zip (zipped file with keyword sent by email). Once unzipped go to this folder in MATLAB (pwd must give the folder where you unzipped) and do startup.
To obtain a SDT license key for your PC use sdtcheck('sitereq','base-fem')

3.1 Computation of transfer functions in the frequency domain

The objective of this part is to understand the influence of terms in equation (2) for the computation of transfers as a sum of modal contributions (see course notes section 2.2.4)

H(ω) = [c] [−Mω²+K]⁻¹ [b] ≈

∑

j=1

[c]{φ_j}{φ_j}^T [b]

−ω ²+iωζ_jω_j+ω_j²

(2)

and thus to analyzed equations of motion in the modal basis

[Is²+[^\ 2ζ_jω_j _\ ]s+[^\ ω_j² _\ ]]{q_r}=[φ^T b]{u} et {y} = [cφ] {q_r} (3)

A FEM mesh of the considered engine cover was done earlier and can be loaded using

 load cc_model
 cf=feplot; cf.model=model; cf.def=def;fecom('ColorDataEvalY');
 model
 def

The data contains two data structures

>> model % model contains the model
model = 
    Node: [2031x7 double]  % sdtweb('node') 
     Elt: [2021x18 double] % sdtweb('elt')
>> def  % contains the modes
def = 
      def: [10195x30 double]
      DOF: [10195x1 double]
     data: [30x1 double]
      fun: [1 2 0 8 0 0]
    label: 'Normal modes'
      lab: {30x1 cell}

model gives nodes, elements, properties ... def contains the first 30 modes.

Identify the fields of def in formula (2) and in the spatial view of the mesh shown on your screen

The following commands allow the generation of transfer functions. Fill in the ??? to obtain an input at 426y (DOF 426.02) and an output at 1976y (DOF 1976.02) (one uses the low level calling format sys=nor2ss(def,damp,InDOF,OutDOF,'Hz acc')).

 load cc_model
 cf=feplot; cf.model=model; cf.def=def;fecom('ColorDataEvalY');
 damp=.01;
 f=linspace(0,5000,2048)'; % frequencies
 sys=nor2ss(def,damp, ??? , ??? ,'Hz acc');
 qbode(sys,f*2*pi,'iiplot "simul1" -po');
 ci=iiplot;iicom('sub1 1'); % pointer to responses
 ci.Stack{'simul1'} % Data structure containing response

Comment the effect of changing damping damp (ζ_j). Note : to overlay multiple simulations use

 sys1=nor2ss(def, ??? , ??? , ??? ,'Hz acc');
 qbode(sys1,f*2*pi,'iiplot "simul1" -po');
 sys2=nor2ss(def, ??? , ??? , ??? ,'Hz acc');
 qbode(sys2,f*2*pi,'iiplot "simul2" -po');

The Rayleigh damping model assumes C=α M + β K, for α=0 choose β giving 1% damping for the first flexible mode. Comment the (poor) transfer obtained (indication damp must be a vector proportional to def.data).
Comment the effect of choosing different input locations. Compare responses for input at 911.02 and 426.02). Indication to strengthen your argumentation: compare modal controlabilities φ_j^Tb with the code below and see amplitudes of various modes at these DOF
```
 % commandability computation 
 [(1:size(def.def,2))' def.def'*fe_c(def.DOF,[426;911]+.02)']
 % Show node numbers
 fecom('ColorDataEvalY -alpha.3');%fecom('showline');
 fecom('textnode 426 911','FontSize',16,'Color','r')
```
Propose a method to select an input location allowing a good controlability of all modes (implementation is not necessary)
Use equations of motion in a truncated modal basis (3) to estimate the effect of adding a 10 g mass at nodes 426 then 911. The model uses MM units (mass in kg, displacement mm, force milli-N). Rationale : a mass is an external force proportional to acceleration at the same location (collocated transfer). What is the error made by only applying this mass in the y direction ?

3.2 Time simulation and FRF estimation

From the FEM model, one builds a state-space model using SDT commands, then generate the transient and draws the response using cc_simul

 load cc_model
 cf=feplot; cf.model=model; cf.def=def;
 uf.t=linspace(0,.5,2048)';  % time vector
 uf.u0=uf.t*0;uf.u0(find(uf.t<=1e-4))=1;  % ch1 (input =impact 0.1 ms);
 uf.window='None';  uf.noise=0; uf.filt=[];
 damp=.001; 
 uf.sys=nor2ss(def,damp,426.02,1976.02,'Hz acc'); % State space model
 uf=cc_simul('simul',uf);cc_simul('plot',uf);

Run the script above. Questions:

The initial simulation has an aliasing problem. Describe the effects.
Verify that aliasing disappears for a higher sampling frequency. What modeling approximation allows this effect to occur ?
Describe the implementation of the simulation done in cc_simul with the simul command at lines 45 to 77, as a block diagram for example. (You can open the code at the right line using sdtweb cc_simul('Without_filter')). You can describe fields of the uf structure and their use.

One can use an anti-aliasing filter using the following call

 load cc_model; cf=feplot; cf.model=model; cf.def=def;
 uf.t=linspace(0,.5,2048)';  % time vector
 fmax=1500;[num,den]=butter(16,fmax*2*pi,'s');uf.filt=tf(num,den);
 uf.u0=uf.t*0;uf.u0(find(uf.t<=1e-4))=1;  % ch1 (input impact 0.1 ms);
 uf.window='None';  uf.noise=0;damp=.001; 
 uf.sys=nor2ss(def,damp,426.02,1976.02,'Hz acc'); % state-space model
 uf=cc_simul('simul',uf);cc_simul('plot',uf)

Describe the usefulness of the anti-aliasing filter in this simulation.

If you have a little time

analyze the effect of windowing using an exponential window
```
 uf.window='exponential 0 10 3'; % 0 pts at 0, 10 pts at 1, end at exp(-3)
 uf=cc_simul('simul',uf);cc_simul('plot',uf)
```
Indication : look at the low frequency response and at the resonances. If you really want to look, it is possible to estimate the damping induced by the exponential window.

analyze the effect of a sine sweep (often used for product qualification tests). How are modes seen in the response ?

 uf.t=linspace(0,.5,8000)';  % time vector
 fmin=200;fmax=5000;
 u=@(t)cos(2*pi*(fmin+(fmax-fmin)/2*t/uf.t(end)).*t);  uf.u0=u(uf.t);

End of lab work 2