Ingegneria Forum

@ Fazio
Si può utilizzare lo stesso programma di telaio tridimensionale assenando un solo impalcato in altezza (ad esempio alto 1 metro) e dichiarando il piano stesso estensionalmente rigido (master di piano baricentrico). Al master di piano si assegnano forze orizzontali tali da generare nelle sezioni di incastro alla base i momenti Mx,My desiderati (Fx=My*1=My:  Fy=Mx*1=Mx). L'eventuale sforzo normale va ripartito nelle 3 sottosezioni rettangolari o in base all'effettivo sforzo presente o in proporzione alle rispettive aree sezionali.
Risolto lo schema tridimensionale ci ritroviamo nelle sezioni di incastro delle tre mensole gli sforzi di verifica relativi alle 3 sottosezioni rettangolari.
Se nel master di piano applichiamo anche un momento torcente di piano il calcolo ci consegnerà anche lo sforzo flettente che fa equilibrio al torcente applicato (il famoso bimomento di Vlasov). In più se la sezione a parete sottile è chiusa e pluriconnessa otteniamo anche i tagli che riferiti alle singole sottosezioni realizzano l'equilibrio con il torcente di piano.
Si può volere di più da uno schema a telaio!

Parte 9

Larghezza delle fibre e loro baricentro
Bene, rispetto a quanto detto nella parte 8, dovrò rivedere gli aspetti deformativi della sezione, in modo da tenere conto anche di sforzi normali di trazione agenti sulla sezione, accogliendo così i suggerimenti di Alberto e Renato.
La definizione della lunghezza di ogni singola fibra della (eventuale) parte compressa di calcestruzzo può però essere affrontato, in quanto esso dipende da sole considerazioni geometriche.

Con un po' di geometria analitica bisognerà valutare le intersezioni della retta y=yf (con yf la ordinata della fibra di interesse) con il perimetro di uno dei poligono.
Prima bisognerà però svolgere alcune considerazioni “topologiche”.
La situazione standard infatti è quella della figura a sinistra della volta scorsa, che ripropongo:



Un po' più articolato è quanto accade, invece, nella parte destra della stessa figura. Infine, ecco il caso che mi pare il più “complesso” che possa capitare (ma accetto smentite):



In ogni caso si vede che la retta y=yf interseca un singolo poligono, sempre in un numero pari di punti (n), staccando un numero dispari di segmenti (n-1). Si vede anche come la lunghezza complessiva della fibra sia la somma delle lunghezze dei segmenti dispari.
Bene, vediamo adesso le coordinate delle intersezioni tra retta y=yf e poligono. Prendiamo un singolo segmento, come in figura:



Il segmento i / i+1 ha una equazione y=mx+q, dove:

m=(y[i+1]-y[i ])/(x[i+1]-x[i ])
q=y[i ]-m*x[i ]

Mettendo a sistema le due equazioni:

y=mx+q
y=yf

si ottiene per xf:

xf=(yf-q)/m

Adesso sostituendo ai vari simboli le espressioni, ricaviamo questa espressione per xf un po' articolata:

xf=[(yf-y[i ])*(x[i+1]-x[i ])+(y[i+1]-y[i ])*x[i ]]/(y[i+1]-y[i ])

Ma che mi consente di non dover gestire il caso particolare di segmento del poligono perfettamente verticale (che avrebbe m, coefficiente angolare, infinito). Certo, la formula cadrebbe in difetto nel caso di segmento orizzontale (y[i+1]-y[i ] diventa nullo), ma due rette orizzontali....si incontrano all'infinito comunque.
Per rendere l'algoritmo della funzione un po' più veloce basterà non processare, dunque, tutti i segmenti orizzontali, e tutti i segmenti che hanno punto iniziale e finale contemporaneamente più in alto o più in basso di yf.
Adesso quindi, facendo il giro, segmento per segmento del poligono, troverò tutte le coordinate xf di intersezione del poligono con la retta y=yf.
Faccio notare, ma magari lo avrete già notato guardando la figura con la “greca” e comunque riguardatela, che le coordinate xf non vengono ricavate in maniera crescente, ma a “casaccio”. Può accadere anche nel caso semplice di normale doppia intersezione. Dipende infatti da quale è il primo lato del poligono.
Pertanto dopo aver ricavato tutte le xf bisognerà procedere al loro ordinamento crescente, in modo da poter procedere quindi alla definizione della lunghezza della fibra, che ricordo consiste nel sommare le lunghezze dei segmenti dispari (quindi xf2-xf1 e xf4-xf3 e xf6-xf5, ecc.).

Però in effetti quello che ci serve per i calcoli successivi non è solamente la larghezza della fibra (necessaria per definire un N), ma anche il baricentro di questa serie di segmenti (necessari, trovato l'N voluto, a definire i momenti di rottura della sezione).
Pertanto la funzione che calcola la lunghezza della fibra dovrà restituirmi anche la coordinata xb di tale punto.
Lo sappiamo che in C le funzioni restituiscono un solo valore alla volta. Per cui si potrebbero scrivere 2 funzioni, ognuna con un ben preciso valore restituito. Sarebbe un peccato perchè le 2 funzioni farebbero al 95% lo stesso lavoro. Si potrebbe inserire un flag tra i parametri di chiamata di una singola funzione, in modo da far restituire o la lunghezza della fibra, o il suo baricentro. Ma questo significa chiamare 2 volte la funzione e far calcolare due volte all'elaboratore le stesse cose. La soluzione, quella che come al solito mi piace, ormai lo avrete capito, è quella di una struttura dati. Così:

struct lungh_fibra
{
 float lungh;
 float xb;
};
La funzione completa invece sarà:

struct lungh_fibra calcola_lungh_fibra(struct poligono_sezione,float); /* Prototipo funzione */

struct lungh_fibra calcola_lungh_fibra (struct poligono_sezione polic,float yf)
{
 int register j,k;
 int ni=0;    /* Numero intersezioni con il poligono */
 float xf[20]; /* Coordinate intersezioni con il poligono */
 float provv,scy=0.0;
 struct lungh_fibra fibra;

 for (k=0;k<=polic.numv-1;k++)
     {
      int kp1=k+1;

      if (k==polic.numv-1) kp1=0;
     
      if (polic.y[k]>yf && polic.y[kp1]>yf) continue; /* Segmento tutto più in alto della fibra */
      if (polic.y[k]<yf && polic.y[kp1]<yf) continue; /* segmento tutto più in basso della fibra */
      if (polic.y[k]==polic.y[kp1]) continue; /* Segmento orizzontale */
     
      xf[ni]=((yf-polic.y[k])*(polic.x[kp1]-polic.x[k])+(polic.y[kp1]-polic.y[k])*polic.x[k])/(polic.y[kp1]-polic.y[k]);
      ni++;
     }
 
 if (ni!=0) /* Controlla se effettivamente ci sono intersezioni */
    {
     /* Ottenute tutte le intersezioni provvede a metterle in ordine crescente */
     for (k=0;k<=ni-1;k++)
         {
          for (j=k+1;j<=ni-1;j++) if (xf[j]<xf[k]) { provv=xf[j]; xf[j]=xf[k]; xf[k]=provv;  }
         }

     /* Infine calcola le lunghezze dei vari tratti */
     fibra.lungh=0.0;
     for (k=0;k<=ni-1;k+=2)
         {
          fibra.lungh+=xf[k+1]-xf[k];
          scy+=(xf[k+1]-xf[k])*(xf[k+1]+xf[k])/2;
         }
     /* e definisce il baricentro della fibra nella sua totalità */
     fibra.xb=scy/fibra.lungh;
     
     return (fibra);
    }
 else { fibra.lungh=0.0; fibra.xb=0.0; return(fibra); }
}
La funzione chiamante provvederà a moltiplicare la lunghezza della fibra per la larghezza della fibra stessa (abbiamo detto un parametro gestibile dall'utente, oppure una variabile globale del programma – 1 mm? 2 mm?), in modo da avere un'area e dunque a moltiplicarla per la tensione ricavata dal diagramma costitutivo del materiale in funzione della effettiva deformazione della fibra.

Per il momento, mi fermo qui, rimandando a dopo pasqua altri sviluppi. Nel frattempo spero non vi siate accapigliati “alla Vlasov” per decidere se una sezione si mantiene piana o meno a seguito della deformazione.

Per rispondere a chi ha tirato fuori la teoria di Vlasov, essa si applica alla torsione di sezioni (forse) aperte a pareti sottili.
In pratica in questi casi oltre a tensioni tangenziali nello stesso piano della sezione si destano nella sezione anche delle tensioni a farfalla agenti in maniera normale alla sezione. Tensioni dunque che si sommano a quelle della flessione (che si potranno ricavare con l'ipotesi di sezione piana?). La tensione a farfalla ha comunque un diagramma strano, nel senso che si inverte sia in senso verticale che orizzontale. Vado a memoria, perchè mi pare, sempre da letteratura tecnica, che sia possibile ricavare un “indicatore geometrico” che consente di stabilire se il “fenomeno Vlasov” ha una influenza importante o meno sulla sezione che si sta studiando.
Per le sezioni in cemento armato le realizzazioni “sottili” sono in genere elementi prefabbricati del tipo tegoli “alari” e simili proposti da vari prefabbricatori. Orbene, per gli spessori usuali (minimo 6-7 cm) e le altre dimensioni della sezione, alla fine anche queste sezioni non sono definibili “sottili” alla Vlasov.

@Gilean

Non ho capito in cosa consisterebbe la "messa in rilievo" del topic visto che, bontà vostra, mi sembra che sia stato messo al numero uno fin quasi dall'inizio (sfondo marroncino pallido nella mia visualizzazione).

Parte 10

Ritornando sui propri passi – Le deformazioni a rottura
Visti gli interventi di Renato ed Alberto Ferrari rettifico alcune cose della parte 8 riguardanti la deformazione a rottura della sezione.
Approfitto della figura che Renato ha gentilmente postato per rubargliela e riaggiustare il discorso.



La figura può essere letta in due modi:

1) Ho una certa situazione deformativa, da cui ricavo la posizione dell'asse neutro;
2) Ho una certa posizione dell'asse neutro, da cui ricavo la situazione deformativa.

E' questa seconda via quella che perseguirò.

Definiamo un po' di cose:

yn: ordinata dell'asse neutro;
ymax: massima ordinata tra tutti i vertici della sezione;
yamin: minima ordinata tra tutte le barre di armatura;
?a: deformazione in corrispondenza di ymax;
?b: deformazione in corrispondenza di yamin;
?yd: deformazione massima dell'acciaio (valore negativo se di trazione);
?cu: deformazione massima nel calcestruzzo (valore positivo se di compressione)
?c0: deformazione in cui termina il tratto parabolico del “classico” diagramma parabola-rettangolo

Conoscendo le deformazioni in ymax ed yamin, per l'ipotesi di conservazione della planarità, si ricava la deformazione in qualsiasi altra fibra della sezione. Infatti per ogni fibra posta alla generica quota y, si avrà:

?=?a+(?b-?a)*(ymax-y)/(ymax-yamin)

Dove ?b e ?a vanno inseriti con il loro effettivo segno nella formula.

Le situazioni che possono riscontrarsi guardando la figura sono solamente tre, ovvero:

a) Posizione dell'asse neutro: +?>yn>(??yd?*ymax+?cu*yamin)/(?cu+??yd?)
In questo caso:
?a=?yd*(yn-ymax)/(yn-yamin)
?b=?yd (valore negativo)

Come suggerito da Renato, ad ?yd viene assegnato il valore -0.0675 (pari al 90% di 7.5% deformazione a rottura di un acciaio B450C – vedi §4.1.2.1.2.3 e Tab. 11.3.Ib)

b) Posizione dell'asse neutro (??yd?*ymax+?cu*yamin)/(?cu+??yd?)?yn>yamin
In questo caso:
?a=?cu
?b=?cu*(yn-yamin)/(yn-ymax)

c) Posizione dell'asse neutro yamin?yn>-?
In questo caso il diagramma ruota attorno al punto di deformazione ?c0 posto a quota yc0=yamin+?c0*(ymax-yamin)/?cu, e quindi:
?a=?c0*(ymax-yn)/(yc0-yn)
?b=?c0*(yamin-yn)/(yc0-yn)

Osservate come io abbia dato un valore numerico alla sola deformazione massima ?yd, lasciando le altre grandezze “indefinite”. Perchè a questo punto scatta quanto detto da Alberto Ferrari circa i differenti valori di ?c0 ed ?cu per i calcestruzzi con classe di resistenza superiore a C50/60 (vedi §4.1.2.1.2.2).
Questo fatto, come accennato in un post precedente, comporta delle ulteriori elaborazioni da fare prima di definire la corretta deformazione della sezione.

Intanto occorrerà modificare la struttura dati dei poligoni, inserendo anche tre nuovi valori:

struct poligono_sezione
{
 float x[100];   /* coordinata x del vertice */
 float y[100];   /* coordinata y del vertice */
 int numv;       /* numero di vertici del poligono */
 float omog;     /* coefficiente di omogeneizzazione */
 int traz;       /* 0=non reagente a trazione; 1=reagente a trazione */
 char classe[8]; /* NUOVO Classe di resistenza del cls */
 float fd;       /* Resistenza massima di calcolo */
 int dominio;    /* 0=parabola-rettangolo (cls) 1=bilatera (acciaio) */
 int fase;       /* A futura memoria per gestire fasi di realizzazione */
 float epsc0;    /* NUOVO deformazione di fine parabola del dominio */
 float epscu;    /* NUOVO deformazione ultima del dominio */
};
Io immagino che ci sia una porzione di programma che in funzione di un input dell'Utente (la classe di resistenza) calcoli con le formule di normativa i valori corretti di epsc0 ed epscu.

Poi bisognerà rettificare la funzione che avendo in entrata la deformazione restituisce il valore alfa da moltiplicare per fd per ottenere la tensione di rottura. La funzione presentata nella parte 8, se ricordate, presentava dei valori fissi, che ovviamente adesso non vanno più bene.
Per la parte parabolica del diagramma, facendo passare la parabola per i punti 0;0 e ?c0;1 ed imponendo la tangenza orizzontale in questo ultimo punto, si ottengono i parametri della parabola: a=-1/?co^2 e b=2/?c0

Pertanto la funzione, aggiungendo il parametro epsc0 in ingresso, diventa:

float parabola_rett(float,float); /* Prototipo della funzione */

float parabola_rett(float def,float epsc0)
{
 float alfa;

 if (def>=epsc0) alfa=1.0;
 else alfa=-pow(def,2)/pow(epsc0,2)+2*def/epsc0;

 return (alfa);
}
A questo punto sarei pronto, rettificando tutto quanto ho rettificato, a presentare la funzione che, dato un certo valore yn di ordinata dell'asse neutro determina le deformazioni nei punti di maggiore e minore “altezza” della sezione.

Non senza prima però aggiungere un ulteriore tassello alla questione.
Infatti, se la sezione è costituita da più poligoni, e se alcuni di questi poligoni sono costituiti da calcestruzzo con classe di resistenza superiore a C50/60, e quindi con valori di ?cu ed ?c0 differenti dai soliti 3.5 per mille (minore di tale valore) e 2 per mille (maggiore di tale valore), chi ci garantisce che la deformazione in ymax del 3.5 per mille (poiché per esempio ymax è di un vertice appartenente ad un poligono “normale”) sia compatibile con la massima deformazione possibile di un vertice del poligono “special”, magari posto pochi centimetri più in basso di ymax?
Nella figura successiva viene presentato proprio il caso di cui parlo.



Chiaramente quindi, la funzione che definisce le deformazioni limite della sezione, dovrà tenere conto del problema.
Correggere le cose sembra semplice a prima vista (la linea rossa della figura rimette a posto le cose), però, riflettendo, la deformata che “comanda” a questo punto è quella del poligono “forte” (quello più scuro della figura), e questo dovrà rispettare le condizioni a) b) o c) comunque (massima tensione di rottura in un lembo piuttosto che nell'altro).
Però gli effettivi casi che possono realmente verificarsi sono il b) (in quanto al lembo superiore si ha già una deformazione limite, il materiale cls) oppure il c).

Allora, per il caso b) le nuove deformazioni, indicate con ?'a ed ?'b in figura, valgono:

?'a=?'cu*(ymax-yn)/(y'max-yn)
?'b=?'cu*(yn-yamin)/(yn-y'max)

Invece nel caso c) si modifica anche il punto di rotazione della sezione, ovvero:

y'c0=yamin+?'c0*(y'max-yamin)/?'cu, e quindi:
?'a=?'c0*(ymax-yn)/(y'c0-yn)
?'b=?'c0*(yamin-yn)/(y'c0-yn)

Questi controlli di “congruenza” vanno estesi a tutti i poligoni che costituiscono la sezione, ritoccando di volta in volta i numeri. Certamente è uno ed uno solo il poligono che “comanda”.

Infine ecco una nuova struttura dati:

struct deform_ultime_sezione
{
 float epsa;
 float epsb;
};
Ed ecco la funzione (preceduta dal solito prototipo), che restituisce proprio una struttura dati come quella di cui sopra:

struct deform_ultime_sezione calcola_deform_ultime(struct poligono_sezione *,float,int,float,float);

/* Dati di ingresso della funzione ------------------------ */
/* ymax=massima ordinata della sezione in cls               */
/* npm=numero del poligono in cui si ha la massima ordinata */
/* yamin=minima ordinata delle armature                     */
/* yn=ordinata asse neutro                                  */
struct deform_ultime calcola_deform_ultime(struct poligono_sezione *polic,float ymax,int npm,float yamin,float yn)
{
 int register k,np;
 float y1[5],yc0,dpa;
 struct deform_ultime_sezione def;

 /* Calcola le massime ordinate di ogni singolo poligono */
 for (np=0;np<4;np++)
     {
     y1[np]=-1000000.0;
      for (k=0;k<=polic[np].numv-1;k++)
          {
           if (polic[np].y[k]>y1[np]) y1[np]=polic[np].y[k];
          }
     }

 /* Definisce le deformazioni ultime della sezione */
 if (yn>(fabs(epsyd)*ymax+polic[npm].epscu*yamin)/(polic[npm].epscu+fabs(epsyd)))
    {
     def.epsa=epsyd*(yn-ymax)/(yn-yamin);
     def.epsb=epsyd;
    }
 if (yn<=(fabs(epsyd)*ymax+polic[npm].epscu*yamin)/(polic[npm].epscu+fabs(epsyd)) && yn>yamin)
    {
     def.epsa=polic[npm].epscu;
     def.epsb=polic[npm].epscu*(yn-yamin)/(yn-ymax);
    }
 if (yn<=yamin)
    {
     yc0=yamin+polic[npm].epsc0*(ymax-yamin)/polic[npm].epscu;
     def.epsa=polic[npm].epsc0*(ymax-yn)/(yc0-yn);
     def.epsb=polic[npm].epsc0*(yamin-yn)/(yc0-yn);
    }

 /* Controlla che per tutti i contorni valga quanto ricavato */
 for (np=0;np<=4;np++)
     {
      if (np==npm || polic[np].numv==0) continue;

      dpa=def.epsa+(def.epsb-def.epsa)*(ymax-y1[np])/(ymax-yamin);

      if (dpa>polic[np].epscu)
         {
          if (yn<=(fabs(epsyd)*y1[np]+polic[np].epscu*yamin)/(polic[np].epscu+fabs(epsyd)) && yn>yamin)
             {
              def.epsa=polic[np].epscu*(ymax-yn)/(y1[np]-yn);
              def.epsb=polic[np].epscu*(yn-yamin)/(yn-y1[np]);
             }
          if (yn<=yamin)
             {
              yc0=yamin+polic[np].epsc0*(y1[np]-yamin)/polic[np].epscu;
              def.epsa=polic[np].epsc0*(ymax-yn)/(yc0-yn);
              def.epsb=polic[np].epsc0*(yamin-yn)/(yc0-yn);
             }
         }
     }

 return(def);
}
Come vedete non ho fatto altro che implementare tutte le formule prima ricavate. La variabile epsyd è una variabile di tipo globale. Ovviamente si presuppone che prima di chiamare per la prima volta questa funzione il programma abbia già provveduto a calcolare preventivamente ymax ed yamin.

E' stato complicato, non difficile. Speriamo bene per i prossimi “mattoni”.

Parte 12

La determinazione dell'asse neutro
Adesso bisogna scrivere una funzione che, utilizzando tutte le altre già riportate, sia in grado di identificare il corretto posizionamento dell'asse neutro, in modo che lo sforzo normale totale generato nella sezione sia pari ad un Nd sforzo normale dato.
Una volta ricavata la posizione dell'asse neutro sarà un gioco da ragazzi determinare gli mxr ed myr, momenti flettenti di rottura, avendo a disposizione anche la posizione delle risultanti dello sforzo normale della parte compressa e tesa della sezione (e queste le abbiamo perchè le funzioni finora riportate svolgono questo doppio calcolo: risultante e sua posizione nel piano x,y).
E' chiaro che la ricerca della posizione dell'asse neutro deve avvenire per via iterativa. Si stabilisce una posizione dell'asse neutro, si ricavano le risultanti di compressione e trazione, e si controlla che la loro somma algebrica sia pari ad Nd. Se così non è si sposta l'asse neutro e si prova nella nuova configurazione geometrica.
Ovviamente bisognerà rendere intelligente una simile ricerca, senza sprecare inutilmente la potenza dell'elaboratore, considerando poi che la posizione effettiva dell'asse neutro andrebbe da +∞ a -∞.

Come al solito non invento nulla, ed ho implementato un meccanismo banale, una ricerca della posizione dell'asse neutro “a passo variabile”. Sono arci sicuro che esistono dei metodi migliori di quello che propongo, e su questo attendo fiducioso.

Con un semplice esempio vi spiego in cosa consisterebbe il metodo “a passo variabile”.
Immaginiamo che Nd sia 10000 e che da un primo tentativo di posizionamento dell'asse neutro, la somma algebrica di N compressione ed N trazione nella sezione sia pari a 7000.
Risulta evidente che bisogna aumentare la parte compressa di cls per arrivare ad Nd, modifico pertanto il posizionamento dell'asse neutro, posizionandolo più in basso di una quantità “delta”. Da questa nuova posizione, ricalcolando il tutto ottengo una risultante nella sezione di 11000. Ecco quindi che avendo superato Nd mi accorgo che la posizione effettiva dell'asse neutro si trova in mezzo al tratto “delta”. A questo punto riduco il passo “delta” (ad esempio della metà, ad un terzo, ecc.) e “torno indietro” sui miei passi (lo inverto di segno), ed ogni volta che accade che la differenza tra Nd ed N calcolato della sezione vari di segno riduco ulteriormente il passo “delta” di variazione dell'asse neutro (ed il suo segno). Il tutto fino ad ottenere una differenza minima tra Nd ed N calcolato, oppure un passo delta realmente microscopico.
La funzione in questione deve restituire i due valori di momento resistente ricavati. Ecco quindi nascere l'esigenza della ennesima struttura dati:

struct risultante_n_finale
{
 struct risultante_n ncf; /* Risultante degli sforzi di compressione nella sezione */
 struct risultante_n ntf; /* Risultante degli sforzi di trazione nella sezione */
};
Ed ecco infine la funzione che implementa questo meccanismo, con l'immancabile prototipo:

struct risultante_n_finale determina_asse_neutroSLU(struct poligono_sezione *,struct armature_sezione,float);
struct risultante_n_finale determina_asse_neutroSLU(struct poligono_sezione *polic,struct armature_sezione armc,float Nd)
{
 int register k,np;
 int npm;
 float ymax=-1000000.0,yamin=+1000000.0,yn,ntot,delta=10;
 struct risultante_n_finale nfin;
 struct deform_ultime_sezione deform;
 struct risultante_n nc,nt;
 
 /* Per prima cosa determina ymax della sezione */
 for (np=0;np<=4;np++)
     {
      for (k=0;k<=polic[np].numv-1;k++)
          {
           if (polic[np].y[k]>ymax) { ymax=polic[np].y[k]; npm=np; }
          }
     }
 /* Poi determina yamin tra tutte le armature */
 for (k=0;k<=armc.numarm-1;k++)
     {
      if (armc.y[k]<yamin) yamin=armc.y[k];
     }
 /* Quindi definisce un primo posizionamento dell'asse neutro */
 yn=yamin+0.8*(ymax-yamin);
 
 /* Inizia il ciclo iterativo alla ricerca della effettiva posizione dell'asse neutro */
 do
  {
   deform=calcola_deform_ultime(polic,ymax,npm,yamin,yn);
   /* Calcola lo sforzo normale nella parte compressa della sezione */
   nc=risult_compr(polic,armc,ymax,yamin,yn,deform);
   nt=risult_traz(armc,ymax,yamin,yn,deform);
   ntot=nc.c+nt.n;
   
   if (delta*(Nd-ntot)>0.0 && fabs(Nd-ntot)>1.0) { delta=-delta/5; yn+=delta; continue; }
   if (fabs(Nd-ntot)>1.0) yn+=delta;
  }
 while (fabs(Nd-ntot)>1.0 && fabs(delta)>1e-6)
 
 nfin.ncf.n=nc.n; nfin.ncf.x=nc.x; nfin.ncf.y=nc.y;
 nfin.ntf.n=nt.n; nfin.ntf.x=nt.x; nfin.ntf.y=nt.y;

 return (nfin);
}
Come si vede la parte iterativa (il do.....while) è contenuta in pochissime righe (visto che il grosso del lavoro viene svolto comunque da risult_compr e risult_traz e tutte le altre funzioni ad esse collegate, precedentemente esposte). E che la parte “a passo variabile” è contenuta in un solo rigo che merita qualche spiegazione.
La variabile delta è la quantità di cui si alza o si abbassa l'asse neutro. Essa oltre a diminuire via via nelle varie iterazioni deve anche poter cambiare di segno.
Pertanto ci si potrà trovare in queste condizioni:
1) delta è positiva (significa che l'asse neutro si alza via via facendo diminuire lo sforzo normale di compressione); se l'espressione (Nd-ntot) è anch'essa positiva vuol dire che ntot deve aumentare, e quindi delta, oltre a diminuire in valore deve anche modificare il proprio segno (deve aumentare la parte compressa della sezione ovvero abbassarsi l'asse neutro).
2) delta è negativa (significa che l'asse neutro si abbassa ad ogni iterazione); se l'espressione (Nd-ntot) è anch'essa negativa vuol dire che ntot deve diminuire, e quindi delta diminuirà e si modificherà di segno (deve diminuire la parte compressa della sezione ovvero alzarsi l'asse neutro).

Badate bene che questo è vero anche se Nd è uno sforzo di trazione e quindi di segno negativo.

Ed alla fine sono partito da un delta da 10 cm in 10 cm che via via diminuisce di 5 volte (quindi 2 cm, poi 0.4 cm, ecc.).
E questo è tutto per adesso. Alla prossima ed ultima parte.

A dire il vero avevo pensato inizialmente di posizionare l'asse neutro di 1° tentativo in corrispondenza dell'assetto 2.
Poi però riflettendo, ho pensato, iterazione per iterazione, cosa cambia?

Io non saprei dirti il numero di iterazioni per arrivare al mio 1 kg di tolleranza che avrei imposto. Immagino che molto dipenda dalla fortuna (ovvero da quanto parta prossimo al valore "vero" dell'asse neutro).
In ogni caso a livello concettuale non ci sono differenze marcate tra quanto proponi (ogni volta esamino cosa succede dividendo a metà la variabile che ho chiamato"delta"), e quanto ho proposto io, ovvero di dividere per 5 tale valore iniziale e "spazzolare" in questo range.

Ciò che non mi è chiaro dalla tua spiegazione è perchè di "auto-limiti" all'assetto 2-3. Ciò che proponi non sarebbe valido anche per gli altri possibili assetti?
O viceversa, se non ricado in quell'ambito, che dovrei fare?

Caro Renato, mi hai fatto tornare indietro nel tempo.
Quando ragazzotto di periferia mi recavo allo sportello della banca con un assegno in mano nella speranza che me lo cambiassero senza tante storie.
Possibile mai che io non sia mai riuscito nell'impresa senza inconvenienti?
Ogni volta il cassiere, mentre lo guardavo attonito, mi costruiva una montagna virtuale di impedimenti e di eccezioni per cui io, con quell'assegno, con quella richiesta, mi trovavo nel posto sbagliato al momento sbagliato, ecc.
Salvo poi, improvvisamente, e con un guizzo degno di fantastisca scuola attoriale, ecco illuminarsi in volto, e TROVARE la soluzione che mi avrebbe consentito di fare finalmente l'operazione, e che effettivamente faceva. Grazie, grazie, grazie, rispondevo.
Ma un pò di amaro in bocca queste scene me lo lasciavano. Perchè se la soluzione c'è, perchè fare tante storie?
Si è trattato di un rodaggio comunque istruttivo. Nella vita ho poi avuto modo di incontrare parecchi "costruttori" di montagne, per fortuna sempre con traforo incorporato.

Per tornar seri. Condivido ciò che dici. Se sto tra la posizione limite 2 e la posizione limite 3, ecco che il problema si semplifica assai. E condivido anche che nel 90% dei casi "reali" è proprio in quell'ambito che mi troverò.
Ma tu devi scusarmi se mi permetto. Io posso darti ragione una volta, non tutte le volte.
Specie se in una di queste volte, poi, ti contraddici.
Proprio tu, quando nella parte 8 ho parlato della configurazione deformata a rottura della sezione, mi hai bacchettato dicendomi che mi ero limitato ad una sola delle possibile configurazione di rottura.
Ed allora ti ho dato assolutamente ragione. Tanto è vero che nella parte 10 ti ho risposto, partendo proprio dalla tua osservazione, sviluppando (a meno che non ci siano errori) tutte le possibilità che possono riscontrarsi.
Pertanto, proprio perchè il programma che propongo è generico e generale, posso solamente rammaricarmi che se certe ipotesi si verificassero le cose si semplificano. Ma io devo (e così ho fatto) considerare anche le cose che si possono verificare solamente nel 10% del casi.

E quindi stavolta non ti do ragione.

Parte 13

Il dominio di rottura - finalmente
Siamo praticamente alla fine di questa esposizione. Manca infatti solamente la funzione che materialmente determini il dominio di rottura della sezione per dato Nd.
Come si era detto proprio ad inizio trattazione delle verifiche SLU, tale dominio verrà ricavato da coppie di mrx,mry calcolate non già facendo ruotare l'asse neutro, ma facendo ruotare fittiziamente la sezione stessa.
In questo modo sarà possibile sfruttare tutte le routine che ho via via proposto e che appunto avevano questa impostazione, ovvero di considerare orizzontale l'asse neutro (e con porzione compressa della sezione al di sopra dell'asse neutro).
Chi ha avuto la voglia di seguirmi fin qui, ed è stato particolarmente attento, avrà notato che tutte le funzioni che dovevano “manipolare” i dati geometrici della sezione oppure delle armature non sono state scritte per operare direttamente sulle strutture dati di input della sezione stessa (poli ed arm), ma esse accettano in ingresso sempre le strutture dati polic e armc dove la c finale sta per “copia”. Ciò consente di dare in pasto alle varie funzioni le coordinate della sezione ruotata rispetto all'input dell'utente. E' chiaro ovviamente che quando le funzioni restituiscono le coordinate della risultante di sforzi normale di compressione o trazione, queste ultime devono essere riportate al primigenio sistema di riferimento per determinare i momenti di rottura (che ad unico sistema devono essere riferiti).

Ci servono quindi le formule per la trasformazione delle coordinate dei vertici e delle armature costituenti i vari poligoni della sezione.
Nell'immagine inferiore, lo schema geometrico dell'operazione più generale possibile, ovvero la roto-traslazione:



Il problema è risolvibile per via geometrica-trigonometrica e permette di ricavare le “nuove” coordinate X,Y dei vertici (in maiuscolo), semplicemente conoscendo le coordinate della origine del nuovo sistema di riferimento, rispetto al vecchio, e la rotazione degli assi.
Le formule generali per la roto-traslazione sono le seguenti:

X= x*cos(alfa)+y*sen(alfa)-x0
Y=-x*sen(alfa)+y*cos(alfa)-y0

Se l'operazione è di semplice rotazione del sistema di riferimento ecco che nelle formule di cui sopra scompaiono i fattori x0 ed y0 solamente.

Quindi propongo l'ennesima struttura dati:

struct dominio_rottura
{
 float nd;
 float mrx[24];
 float mry[24];
};
Infine la funzione, preceduta dal solito prototipo:

struct dominio_rottura calcola_dominio_SLU(float);
struct dominio_rottura calcola_dominio_SLU(float Nd)
{
 int register k,np,step;
 struct poligono_sezione polic[5];
 struct armature_sezione armc;
 struct risultante_n_finale nfin;
 struct dominio_rottura dom;
 float xc,yc,xt,yt;

 dom.nd=Nd;

 for (step=0;step<=23;step++)
     {
      float alfa=step*0.261799; /* alfa in radianti – 0.261799=15° */
      float c=cos(alfa),s=sin(alfa);

      /* Prima copia la sezione in altra struttura dati e poi la ruota di alfa */
      for (np=0;np<=4;np++)
          {
           polic[np]=poli[np];
           for (k=0;k<=poli[np].numv-1;k++)
               {
                polic[np].x[k]= poli[np].x[k]*c+poli[np].y[k]*s;
                polic[np].y[k]=-poli[np].x[k]*s+poli[np].y[k]*c;
               }
          }

      /* Copia le armature in altra struttura dati e poi le ruota di alfa */
      armc=arm;
      for (k=0;k<=arm.numarm-1;k++)
          {
           armc.x[k]= arm.x[k]*c+arm.y[k]*s;
           armc.y[k]=-arm.x[k]*s+arm.y[k]*c;
          }

      /* Chiama la funzione per definire l'asse neutro della sezione ruotata di alfa */
      nfin=determina_asse_neutroSLU(polic,armc,Nd);

      /* Ripristina il sistema di riferimento delle risultanti */
      c=cos(-alfa); s=sin(-alfa);
      xc= nfin.ncf.x*c+nfin.ncf.y*s;
      yc=-nfin.ncf.x*s+nfin.ncf.y*c;

      xt= nfin.ntf.x*c+nfin.ntf.y*s;
      yt=-nfin.ntf.x*s+nfin.ntf.y*c;

      dom.mrx[step]=nfin.ncf.n*(yc-yg)+nfin.ntf.n*(yt-yg);
      dom.mry[step]=nfin.ncf.n*(xc-xg)+nfin.ntf.n*(xt-xg);
     }

 return (dom);
}
A parte l'ingresso con unico parametro, Nd, da segnalare l'utilizzo delle variabili globali xg ed yg, baricentro della sezione omogenea (che ho calcolato precedentemente senza considerare la presenza di armature nella sezione). Polo da cui ricavare i momenti di rottura.

E questo sembrerebbe tutto. Nel senso che effettivamente la funzione restituisce 24 punti del dominio di rottura della sezione per dato Nd (mi è piaciuto così, ma nulla toglie di rendere modificabile da un ipotetico utente questa assunzione).
Al solito la funzione chiamante (ancora più “grande” delle scatole una dentro l'altra che finora abbiamo visto) potrà eseguire questa funzione tutte le volte che serve per tutte le combinazioni di carico cui è sottoposta la sezione, o magari procedere ad Nd via via crescenti per ottenere più “fette di salame” sovrapposte in modo da poter diagrammare in 3D la “safe ball” come la chiamano gli inglesi, ecc.

Ma proprio a questo proposito dobbiamo chiarire un ultimissimo aspetto, accennato in qualche commento, ma ancora da affrontare. Ovvero che non tutti gli Nd sono “buoni” e che non da tutti è possibile ricavare la “fetta di salame”. Ci saranno in effetti un valore di sforzo normale di compressione ed uno di trazione “limite” oltre i quali il dominio di rottura non esiste. Pertanto la funzione che chiama calcola_dominio_SLU dovrà prima capire se Nd è “buono” o meno.
L'Nd massimo di trazione, oltre il quale la sezione non è in grado di spingersi, è dato dal prodotto dell'area complessiva delle barre di armature per fyd. Viceversa l'Nd massimo di compressione è dato dalla sommatoria di tutte le aree dei singoli poligoni per il rispettivo fd, e dall'area totale della barre per fyd, e questo poiché al minimo la sezione è uniformemente deformata del 2 per mille – alcuni contorni di cls particolarmente “buono” come abbiamo visto possono giungere a valori di deformazione limite più alti (vedi formule di normativa) – e con quella deformazione l'acciaio con modulo elastico 2100000 daN/cm² è già snervato (infatti 3913/2100000=0.00186).
Ma in ogni caso e per ogni evenienza è bene fare qualche controllo visto che i numeri possono anche essere differenti da quelli previsti (fyd maggiori di acciai migliori).

In pratica la parte di programma che “chiama” calcola_domino_SLU dovrebbe prima determinare Nd-trazione limite ed Nd-compressione limite, in modo da non “disturbare” la funzione se l'Nd della combinazione è al di fuori di questi limiti. Ovvero:

/* Definisce il massimo valore nd di compressione per la sezione */
/* Determina dapprima il minimo valore di epsc0 tra tutti i poligoni della sezione */
epsc0min=1000.0;
for (np=0;np<=4;np++) if (poli[np].epsc0<epsc0min) epsc0min=poli[np].epsc0;

/* Quindi determina l'area e la forza di ogni singolo poligono */
nd_comprlim=0.0;
for (np=0;np<=4;np++)
    {
     area=0.0;
     for (k=0;k<=poli[np].numv-1;k++)
         {
          if (k==poli[np].numv-1) kp1=0; else kp1=k+1;
          area+=(poli[np].x[kp1]-poli[np].x[k])*(poli[np].y[kp1]+poli[np].y[k])/2;
         }
      alfa=parabola_rett(epsc0min,poli[np].epsc0);
      nd_comprlim+=area*alfa*poli[np].fd;
    }
/* Inserisce anche le barre di armatura */
for (k=0;k<=arm.numarm-1;k++)
    {
     alfa=elast_plast_indef(epsc0min);
     nd_comprlim+=arm.af[k]*alfa*fyd;
    }

/* Definisce il massimo valore di nd di trazione per la sezione */
nd_trazlim=0.0;
for (k=0;k<=arm.numarm-1;k++) nd_trazlim-=arm.af[k]*fyd;

/* Ed infine chiama la funzione di costruzione del dominio per ogni combinazione */
for (k=0;k<=numcomb-1;k++)
    {
     if (soll[k].tipo!=3) continue;
     if (soll [k].n>nd_trazlim && soll[k].n<nd_comprlim) dominio=calcola_dominio_SLU(soll[k].n);
     else /* Messaggio di avviso, semplice salto della combinazione, quel che si vuole */
    }
Come si vede ho sempre utilizzato funzioni già definite precedentemente che mi hanno consentito di tenere conto di tutte le possibilità che possano verificarsi. Per le armature compresse ho fatto calcolare il valore alfa dalla funzione elast_plast_indef (ma dovrebbe essere sempre 1.0) inserendo come parametro in ingresso il valore di epsc0min nel caso in cui esso fosse più piccolo del valore di plasticizzazione dell'acciaio (cosa che non ho fatto per le armature tese, per ovvie ragioni).
Per la costruzione del dominio di rottura faccio in modo che vengano processate le sole combinazioni SLU (vedi il controllo di soll[k].tipo che deve valere 3 per proseguire nella elaborazione del ciclo, soll che è una struttura dati definita precedentissimamente) ed inoltre che il valore di n della combinazione stia all'interno dei valori limite.

Più o meno è tutto.
E' arrivato il momento di dire davvero fine a tutta questa trattazione. Non so voi, ma io mi sento pronto adesso, dopo i dieci anni che vi dicevo all'inizio, di rimetterci le mani e completare il programma.

Sono consapevole di aver lasciato qualcosa per strada (ad esempio la possibilità, mai sfruttata, di dichiarare il tipo di diagramma a rottura di un singolo poligono, oppure la possibilità di reagire anche a trazione per le verifiche SLU, le fasi costruttive, ecc.). D'altra parte spero non ci siano errori nelle varie routine proposte. In ogni caso mi impegno fin da adesso a correggerle nel caso in cui mi accorgessi di qualcosa che non funziona.