kill(all);

dati recuperati da modale FEM "lamella_fatto.mfd", 
modi propri normalizzati a massa modale unitaria.
 
Si riprendendo le soluzioni del problema della risposta armonica in un oscillatore ad un grado di
libertà smorzato e forzato, associando al numero di modi propri scelti, quindi 5, la variabile m.

m:5;

Si definiscono gli estremi della risposta in frequenza 0Hz e 1000Hz (sul fem ho definito tra 1Hz e
1000Hz),
questi saranno i valori estremi dei grafici sottoriportati.

fmin:0$
fmax:1000$

Lo smorzamento frazione critico è definito come sul fem, quindi 1%

zeta:0.01$

Dopo aver estratto dal fem la modale, si campionano da essa le frequenze proprie, naturalmente 5.

freq:[
83.37,
524.7,
828.9,
1488,
1608
];

Nel caso in esame si ha un unico g.d.l che viene caricato dall'azione esterna, di esso si prende lo 
spostamento normalizzato in massa modale unitaria e  si ottengono i seguenti spostamenti 
lungo l'asse z al nodo caricato, misurati in millimetri.

disp:[
863.197,
-865.160,
0,
-870.81,
0
];

Tali modi probabilmente non potranno essere eccitati da quella forzante, perchè essa batte su un nodo
del modo, detto nodo l'opposto di un ventre. Nelle forme modali il grado di libertà che ha spostamento
nullo viene detto nodo, mentre il gdl che ha spostamento massimo viene detto ventre.

Si definisce ora l'azione esterna sull'unico g.d.l. caricato che e' una forza quota parte coseno (1.2N),
quota parte seno (0 N);

fr:1.2$
fi:0$

Per estrarre come risultato qualcosa di terzo rispetto allo spostamento in z
del nodo in esame,si campiona una quantità di interesse in output, che sarà la tensione
al top layer del punto centrale della lamella zona incastro, quindi sarà una tensione in direzione
assiale; 

Secondo i 5 modi propri avrò cinque valori di tale tensione, misurati in MPa.

s11:[
-11673.2,
-65554.5,
0,
166612,
0
];

PER SEMPLICITA' DI TRATTAZIONE SI TORNA SUL PDF PER CONTINUARE L'ANALISI DELLA RISPOSTA
ARMONICA TRAMITE SOVRAPPOSIZIONE MODALE

Avendo 5 modi si imposta un ciclo for sui modi, calcolando per ogni l-esimo modo il fattore scalare 
complesso. 
Come prima cosa si creano due strutture per raccogliere i risultati:
si costruisce una lista vuota, es. A = lista vuota, poi si definisce A come la concatenazione 
di se stessa e un termine, procedo il ciclo dopo si definisce A come concatenazione di se stessa 
e di un altro termine.
Come risultati si raccolgono lo spostamento dell'oscillatore ad un grado di libertà ur e ui 
(che di seguito verranno chiamati xii e xir).
Questa lista nasce vuota e andra' ad accumulare modo proprio per modo proprio il fattore
di amplificazione dello stesso.

A questo punto si cicla sui singoli oscillatori secondo i singoli modi propri.
(xir[l] + %i xii[l]) che va a modularlo. 
Si utilizza nella sua definizione la risposta dell'oscillatore monodimensionale equivalente a quel modo.
Si definiscono xir e xii tramite le formule trovate nel foglio di calcolo precedente
(sdof_finale_versione_studenti.wmxm %o13)
La frequenza naturale è definita come 2*pigreco*frequenza i-esima.
Il carico esterno è definito come una quota parte coseno e una quota parte meno seno: la parte coseno qr
è il rapporto della forzante fr=1.2N
e dello spostamento monitorato al gdl caricato secondo l'i-esimo modo proprio, qi è dato dal rapporto di
fi=0N per lo stesso spostamento 
(quindi qi=0). 
f è la frequenza che scalerà da 0Hz e 1000Hz, r è il rapporto tra la frequenza della forzante e quella
dell'oscillatore.
Noti questi parametri si calcola xir e xii e li accodo alle liste xir_list (che contiene la parte reale
o modulata coseno) e xii_list 
(che contiene la parte immaginaria o modulata in meno seno).

xir_list:[]$
xii_list:[]$
for l : 1 thru m do (
    qr:fr*disp[l],
 qi:fi*disp[l],
    omega: 2*%pi*freq[l],
    r:f/freq[l],
    xir: (2*qi*r*zeta-qr*r^2+qr)/(4*omega^2*r^2*zeta^2+omega^2*r^4-2*omega^2*r^2+omega^2),
    xii:-(2*qr*r*zeta+qi*r^2-qi)/(4*omega^2*r^2*zeta^2+omega^2*r^4-2*omega^2*r^2+omega^2),
    xir_list:append(xir_list,[xir]),
    xii_list:append(xii_list,[xii])
);

In alternativa è possibile implementare nel ciclo for le formule in fondo al 
paragrafo sulla risposta armonica per sovrapposizione modale
Di seguito si estraggono direttamente modulo e fase.

xir_list2:[]$
xii_list2:[]$
for l : 1 thru m do (
    qr:fr*disp[l],
    qi:fi*disp[l],
    omega: 2*%pi*freq[l],
    r:f/freq[l],
    a:1-r^2,
    b:2*zeta*r,
    xir:(a*qr+b*qi)/omega^2/(a^2+b^2),
    xii:(a*qi-b*qr)/omega^2/(a^2+b^2),
    xir_list2:append(xir_list2,[xir]),
    xii_list2:append(xii_list2,[xii])
);

In questo modo si ottiene una lista che contiene i termini che andranno a 
modulare i singoli modi propri;

scostamento numerico tra le due soluzioni

wxplot2d(
append(
    (xir_list-xir_list2),
    (xii_list-xii_list2)
),
[f,fmin,fmax]
);

nullo...

mostro i fattori modulanti associati ai modi propri

for l : 1 thru m do(
    disp(
        [
            xi[l],
            ratsimp(
                xir_list[l]+%i*xii_list[l]
            )
        ]
    )
),ratprint=false;

Si nota che il 3 modo proprio è nullo, ciò a causa del fatto che in questo modo il punto che prende il
carico in direzione z
non si sposta, quindi è energeticamente disaccoppiato, pertanto il 3 modo proprio viene modulato per
zero, cioè non viene 
eccitato. Analogamente accade per il  5 modo proprio.

Si procede all'accumulo dei contributi tensionali associati ai vari modi andando a 
ricostruire lo stato tensionale al nodo e al layer scelto come rappresentativo. 
Questo stato tensionale viene rappresentato secondo due quantità:
s11r, quota modulata in coseno di omega*t, e s11i, quota modulata in meno seno di omega*t;
entrambe nascono 0 poi via via si accumulano i contributi modo per modo.
Si costruisce un secondo ciclo con i che va da 1 ad n:
- disp(I) scrive iterativamente il numero del ciclo in output;
- wxplot2d grafica, per ogni modo proprio e con frequenza che va da 0Hz a 1000Hz, i fattori di
amplificazione
  parte reale e parte immaginaria moltiplicati per s11, entità di stress di un punto specifico di un
modo proprio
  normalizzato a massa modale unitaria.

s11r:0$
s11i:0$
for l : 1 thru m do(
    disp(l),
    wxplot2d(
        s11[l]*[xir_list[l],xii_list[l]],
        [f,fmin,fmax],
        [legend,"Re","Im"]
    ),
    [s11r,s11i]:[s11r,s11i]+s11[l]*[xir_list[l],xii_list[l]]
);

ANALISI GRAFICI:
    1 MODO : si ha un picco di risonanza a 83Hz, tale risonanza presenta ampiezza limitata a causa dello
smorzamento zeta = 0.1, 
                      si raggiungono circa 2200MPa. La parte dominante risulta essere quella immaginaria.
    2 MODO : si ha un picco di risonanza a circa 500Hz, la tensione raggiunta rappresenta quella al nodo
di centro incastro, circa -320MPa.
    3 MODO: Piatto a zero perchè non eccitato.
    4 MODO: la risonanza è molto in là e la tensione è di circa 2MPa, quindi trascurabile.
    5 MODO: Piatto a zero perchè non eccitato.

valore statico, da FEM statico 39.2043MPa

Andando a valutare la quota parte reale della tensione vengono -39.79...MPa, il meno è indicativo del
fatto che
la tensione è in controfase rispetto alla forza. 

[s11r,s11i],f=0,numer;

plotto i risultati

Si plotta quindi la risposta cumulativa, somma della risposta secondo il primo modo proprio e il secondo
modo proprio, 
il quarto non è visibile.

wxplot2d(
[s11r,s11i,sqrt(s11r^2+s11i^2)],
[f,fmin,fmax],
[legend,"Re","Im","Mag"]
);

zoom nell'intorno della prima risonanza
Tre curve: parte reale, parte immaginaria e magnitudo (sqrt(Re^2 + Im^2))

wxplot2d(
[s11r,s11i,sqrt(s11r^2+s11i^2)],
[f,freq[1]-5,freq[1]+5],
[legend,"Re","Im","Mag"]
);

andamento al di fuori delle risonanze 
(zoom su valori limitati di amplificazione)

wxplot2d(
[s11r,s11i,sqrt(s11r^2+s11i^2)],
[f,fmin,fmax],[y,-100,100],
[legend,"Re","Im","Mag"]
);

magnitudo su scala logaritmica 
(il secondo elemento nullo nella legenda è introdotto come workaround al messaggio di errore 'The value
"Mag" is not of type LIST', vabbé..)

wxplot2d(
[sqrt(s11r^2+s11i^2)],
[f,fmin,fmax],
[legend,"Mag",""],
[logy]
);

Per concludere si considera la risposta in frequenza per l'esempio della lamella fatta su Marc. Si torni
sul pdf per le conclusioni.