Parte 7
Il codice SLE
In pratica per quanto riguarda la verifica SLE abbiamo già finito. Il problema è stato sviscerato nella sua parte teorica, rimane solamente da presentare un listato.
Ho provato a scrivere una funzione quanto più indipendente dal resto del codice che ci sta attorno. A parte le strutture dati che definiscono la geometria della sezione, le uniche variabili globali sono xg e yg, ovvero le coordinate x,y del baricentro della sezione omogenea.
E poiché le combinazioni cui una sezione viene sottoposta sono ormai numerose, la funzione accetta in ingresso una tripletta di sollecitazioni e restituisce la corrispondente tripletta deformativa. Sarà poi la funzione chiamante (che magari avrà un ciclo di chiamata per ognuna delle combinazioni) a decidere cosa fare di questo output (calcolare le tensioni, stampare una porzione di tabulato, disegnare la parte compressa, determinare il massimo valore di compressione calcestruzzo o acciaio tra più combinazioni, ecc.)
Nel solito header file bisognerà definire altre 2 strutture dati (lo so, ci sono affezionato e forse esagero)
struct deform
{
float u[3];
};
struct sollecitazioni_esterne
{
float n; /* Sforzo normale (positivo di compressione) */
float mx; /* Momento flettente asse x (positivo se comprime fibre +y) */
float my; /* Momento flettente asse y (positivo se comprime fibre +x) */
int tipo; /* 0=SleQp / 1=SleFr / 2=SleRa / 3=Slu */
};
Ed il prototipo della funzione:
struct deform verifica_SLE(struct sollecitazioni_esterne);
Infine, ecco la funzione. Intanto leggetela, alla fine il commento.
struct deform verifica_SLE(struct sollecitazioni_esterne soll)
{
int register k,np;
float norma0,norma1=1.0;
float s1,s2;
struct deform def;
struct poligono_sezione polic[4];
struct geometria_sezione geo;
/* Se tutte le sollecitazioni sono nulle ritorna la struttura def fatta da 0.0 */
if (soll.n==0.0 && soll.mx==0.0 && soll.my==0.0) { def.u[0]=0.0; def.u[1]=0.0; def.u[2]=0.0; return(def); }
def.u[0]=1.0; def.u[1]=0.0; def.u[2]=0.0;
do
{
float parz;
norma0=norma1;
/* Determina la geometria della sezione “reagente” */
for (np=0;np<=4;np++)
{
/* Se il poligono reagisce anche a compressione esso rimane sempre interamente reagente */
/* quindi "poli" viene copiato così come è in "polic" e la procedura passa al poligono successivo */
if (poli[np].traz==1) { polic[np]=poli[np]; continue; }
/* Prima determina i veri punti che definiscono il contorno "reagente" */
/* Calcola la tensione nel punto precedente quello di interesse */
s1=def.u[0]+(poli[np].y[poli[np].numv-1]-yg)*def.u[1]+(poli[np].x[poli[np].numv-1]-xg)*def.u[2];
polic[np].numv=0; polic[np].omog=poli[np].omog;
for (k=0;k<=poli[np].numv-1;k++)
{
int km1=k-1;
if (k==0) km1=poli[np].numv-1;
/* Calcola la tensione nel punto voluto per prendere decisioni */
s2=def.u[0]+(poli[np].y[k]-yg)*def.u[1]+(poli[np].x[k]-xg)*def.u[2];
if (s1*s2<=0.0)
{
polic[np].x[polic[np].numv]=poli[np].x[k]-s2*(poli[np].x[k]-poli[np].x[km1])/(s2-s1);
polic[np].y[polic[np].numv]=poli[np].y[k]-s2*(poli[np].y[k]-poli[np].y[km1])/(s2-s1);
polic[np].numv++;
}
if (s2>0.0)
{
polic[np].x[polic[np].numv]=poli[np].x[k];
polic[np].y[polic[np].numv]=poli[np].y[k];
polic[np].numv++;
}
/* Fa diventare s2 la tensione del punto precedente per il ciclo successivo */
s1=s2;
}
}
/* Con i nuovi contorni chiama la funzione che determina le caratteristiche geometriche della sezione reagente */
geo=calcola_caratt_geometriche (polic);
/* Risolve l'equazione a 3 incognite */
parz=geo.ix*geo.area-pow(geo.scx,2);
def.u[2]=((soll.my*geo.area-soll.n*geo.scy)*parz-(soll.mx*geo.area-soll.n*geo.scx)*(geo.ixy*geo.area-geo.scy*geo.scx))/((geo.iy*geo.area-geo.scy*geo.scy)*parz-pow(geo.ixy*geo.area-geo.scy*geo.scx,2));
def.u[1]=((soll.mx*geo.area-soll.n*geo.scx)-(geo.ixy*geo.area-geo.scx*geo.scy)*def.u[2])/parz;
def.u[0]=soll.n/geo.area-def.u[1]*geo.scx/geo.area-def.u[2]*geo.scy/geo.area;
norma1=sqrt(pow(def.u[0],2)+pow(def.u[1],2)+pow(def.u[2],2));
}
while (fabs(norma1-norma0)>1e-3);
return (def);
}
Rispetto a quanto visto nella parte precedente, come si può vedere, mi sono limitato a riportare la formula che individua il punto a tensione 0.0 all'interno di un lato del poligono per definire il nuovo poligono “reagente”, ed inoltre ho aggiunto il caso in cui uno o più dei poligoni invece che parzializzarsi possa reagire anche a trazione (è il valore di traz, tra i vari dati della struttura dati che definisce ogni singolo poligono). E' una cosa comoda quando si vuol calcolare una sezione metallica composta, o altro.
Il tutto è concentrato tra l'istruzione do e l'istruzione while.
In sintesi:
1) Con le deformazioni dell'iterazione precedente si definisce la reale geometria della sezione reagente (quando si entra per la prima volta nel ciclo le deformazioni valgono 1,0,0 come da riga precedente che le inizializza in tal senso).
2) Si chiama la funzione calcola_caratt_geometriche che determina i valori da inserire nella matrice di rigidezza della sezione
3) Si risolve il sistema lineare 3x3 (senza scomodare alcun metodo, Gauss-Siedel, Kholesky, ecc. Ho riportato materialmente le espressioni definitive avendo risolto prima “a mano” il sistema in maniera simbolica, dunque tre sole formuline)
4) Si calcola la norma della nuova soluzione e la si confronta con la norma della precedente soluzione (nell'istruzione while)
5) Se l'esito della condizione è positivo la funzione termina restituendo la struttura deform alla funzione chiamante, altrimenti il ciclo viene ripetuto fino a convergenza.
Faccio notare che nel determinare le “tensioni” per definire la porzione reagente di poligono, non ho affatto adoperato il fattore omog. Non è necessario infatti che le tensioni siano quelle vere; importante è che il rapporto s2/(s2-s1) sia quello corretto (e lo è visto che parliamo dello stesso poligono ed il fattore omog sarebbe comune), insomma una moltiplicazione in meno.
Altra cosa che faccio notare è che potrebbe capitare, in sezioni costituite da più poligoni, che nella parzializzazione della sezione un intero poligono potrebbe essere tutto in trazione e quindi scomparire. Che succede in quel caso?
Nulla. Perchè, se è vero che nessun vertice viene definito nelle sue coordinate x e y, il numero dei vertici del poligono reagente polic[np].numv impostato inizialmente a 0 dalla funzione rimane tale. E la funzione calcola_caratt_geometriche non si offende affatto di un simile evento, infatti con la riga:
for (k=0;k<=polic[np].numv-1;k++)
semplicemente l'intero ciclo non viene svolto in quanto la condizione k<=polic[np].numv-1 è subito vera. Ovvero come se il poligono non esistesse.
Terza cosa che faccio notare è che una volta usciti dalla funzione inesorabilmente si perdono anche le caratteristiche statiche della sezione reagente “faticosamente” trovate dall'algoritmo (in quanto variabili locali valide solamente all'interno della funzione). Se queste fossero di interesse poco male. Si estenderebbe la struttura dati struct deform andando a definire all'interno di essa, oltre alla tripletta deformativa, anche scx, scy, ix, iy, ixy. Con una semplice assegnazione di queste variabili poco prima del return (def), ovvero def.scx=geo.scx; ecc. ci si potrà “portare indietro” anche queste informazioni.
Infine una osservazione: Se considerate che questa funzione, comprese le righe di commento, ha circa 70 righe, e che l'altra funzione realmente necessaria, quella di calcolo delle caratteristiche geometriche, ne ha poco meno di 50, possiamo considerare che è realmente “poca cosa” la sola parte risolutiva di un programma che verifica le sezioni in c.a. agli SLE.
Un'ultima porzione di codice, tanto per gradire. La funzione chiamante sarebbe più o meno così:
for (nc=0;nc<=numcomb-1;nc++)
{
if (soll[nc].tipo<=2) def[nc]=verifica_SLE(soll[nc]);
/* Calcola le tensioni vertice per vertice della sezione */
for (np=0;np<=4;np++)
{
omog=fabs(poli[np].omog);
for (k=0;k<=poli[np].numv-1;k++)
{
tens=omog*(def[nc].u[0]+def[nc].u[1]*(poli[np].y[k]-yg)+def[nc].u[2]*(poli[np].x[k]-xg));
if (tens<0.0 && poli[np].traz==0) scarta il risultato
else gestisce questo valore per farne ciò che si vuole
}
}
/* Calcola le tensioni armatura per armatura */
for (k=0;k<=arm.numarm-1;k++)
{
tens=arm.omogarm*(def[nc].u[0]+def[nc].u[1]*(arm.y[k]-yg)+def[nc].u[2]*(arm.x[k]-xg));
ed anche in questo caso si può farne ciò che si vuole
}
}
Dove numcomb sarebbe il numero di combinazioni cui è sottoposta la sezione. Le scritte dopo il calcolo della tensione (es. scarta il risultato) non fanno parte ovviamente della funzione ma sono state messe lì da me per indicare la possibilità che si vuole.
Ho dimenticato di dire finora, che il codice fin qui presentato dovrebbe essere puro ANSI C e pertanto compilabile da qualsiasi compilatore C, C++.
Vista la presenza di funzioni tipo fabs (valore assoluto di numeri float) e sqrt (radice quadrata) sarà necessario, nel solito header file, inserire il comando:
#include "math.h"
Ovvero includere la libreria matematica standard del C.
E con questo penso di poter chiudere questa prima parte, passando tosto, con le altre prossime parti, alle verifiche SLU.