· · Meccanica Computazionale · Nuovo Metodo di Integrazione Numerica
dell'Ing. Lamberto Bertoli · · LA NUOVA INTEGRAZIONE NUMERICA PER GLI ELEMENTI FINITI · Lamberto Bertoli · · Gli elementi finiti e la loro integrazione · LA NUOVA QUADRATURA DIRETTA · L’integrazione di funzioni di 2 variabili reali · L’integrazione delle funzioni di 3 variabili reali
· Conclusioni · · ·
·
· I fondamenti scientifici del metodo proposto dall’Ing. Lamberto
Bertoli, i cui contenuti fondamentali sono discussi in questo articolo, sono
stati presentati dai docenti Carmelo Majorana, Stefano Odorizzi e Renato
Vitaliani dell’Università di Padova al “Fourth seminar about method and
variational methods” tenutosi a Plzen (Praga) tra il 12 e il 15 maggio
1981. Successivamente la trattazione fu pubblicata in
versione ridotta sulla rivista “Advanced in Engineering Software” (C.
Majorana, S. Odorizzi, R. Vitaliani, “Shortened quadrature rules for
finite elements”, Advanced in Engineering Software, 1982, Vol. 4, N.2) e,
come selected paper, nel volume “Software in Engineering Problems”. Le soluzioni del problema dell’integrazione numerica, a
partire dai casi più semplici, oggetto della tesi di laurea, fino ai più sofisticati,
che sono stati affrontati nelle ricerche degli ultimi anni, sono state
pubblicate dall’Ing. Bertoli nel testo: “Quadratura
diretta degli integrali multipli”, 2006, Ed. Libreria Internazionale
Cortina, Padova. Gli elementi finiti e la loro integrazione Il metodo degli elementi finiti pone il suo fondamento
nel principio di minima energia che Jacques Bernoulli (1654-1705) intuì alla
fine del ‘600 riflettendo sulla forma assunta da una fune inestensibile tesa
alle sue estremità. Egli comprese che il baricentro del sistema in equilibrio
doveva necessariamente situarsi alla minima altezza rispetto a qualunque
altra configurazione congruente coi vincoli ai suoi estremi. Da questo assunto, il grande studioso concluse pertanto
che l’equilibrio della catena impone la sua minima energia potenziale
rispetto a qualunque forma assumibile dalla medesima. Nel caso di una fune
estensibile, oltre all’energia potenziale è presente l’energia elastica di
deformazione, per cui la complessiva è ottenuta aggiungendo quest’ultima
all’energia potenziale gravitazionale. Joseph Louis Lagrange, ritenuto il più
grande matematico della seconda metà del XVIII secolo, comprese le
potenzialità di questo principio, e sviluppò il calcolo variazionale già
anticipato dalle ricerche pionieristiche di Bernoulli. L’energia di una fune
sospesa dipende infatti dalla sua configurazione e quando raggiunge un
minimo, ne consegue che per ogni cambiamento infinitesimo della sua
curvatura, la variazione di energia deve essere nulla. La forma assunta dalla
fune in equilibrio deve quindi necessariamente soddisfare la condizione che
per qualunque spostamento virtuale dei suoi punti, l’energia deve rimanere
costante. Il principio dei lavori virtuali che viene così frequentemente
utilizzato nella Scienza delle Costruzioni, si fonda proprio sul principio di
minima energia, e ne rappresenta una delle più importanti applicazioni. Molti fisici matematici dopo le fruttuose ricerche di Lagrange proseguirono gli studi sul
calcolo variazionale, la cui diffusione si è estesa dopo la nascita degli
elaboratori elettronici. Il metodo degli elementi finiti si fonda su questo
tipo di calcolo e deve il suo crescente interesse alla sua generalità, in
quanto è in grado di affrontare anche i problemi in cui si rende improponibile
il calcolo analitico, che al contrario deve limitarsi allo studio di un
numero limitato di situazioni idealizzate. Sostituendo ad un sistema continuo un insieme di
elementi di dimensioni finite, il principio di minima energia mantiene la sua
validità assicurando che la configurazione assunta da questo modello in
condizioni di equilibrio riprodurrà la condizione di energia costante per
qualunque spostamento virtuale. Attribuendo a ciascun elemento un numero
limitato di gradi di libertà determinato da punti collocati sui suoi bordi, è
possibile determinarne l’energia di deformazione sulla base degli spostamenti
assunti da ciascuno dei suoi punti campione collimanti con quelli degli
elementi contigui. La funzione che esprime lo spostamento di ogni singolo
punto dell’elemento una volta noto quello del punto di riferimento, sarà
costituita da polinomi in grado di annullarsi in tutti gli altri punti base e
dalla derivazione di queste funzioni è possibile calcolare la deformazione
assunta dal materiale in tutto l’elemento. Da quest’ultima si risale
all’energia attraverso il modulo di elasticità del materiale, ed integrando
su tutto l’elemento si può infine determinare la matrice di rigidezza
relativamente a ciascun grado di libertà stabilito per ogni punto campione. Il problema principale del metodo degli elementi finiti
è proprio il calcolo di questa matrice ottenuta con l’integrazione numerica
estesa a tutto l’intervallo delle tensioni e delle deformazioni indotte dagli
spostamenti unitari di ciascuno dei suoi gradi di libertà. Questa
computazione è infatti così impegnativa da assorbire normalmente oltre la
metà del tempo di calcolo di un programma strutturale. La procedura comunemente usata si basa sulla quadratura
di Gauss, che rappresenta la più rapida integrazione teoricamente possibile
per funzioni di una sola variabile reale. Tale metodo fu studiato dal suo
autore unicamente per funzioni di un’unica variabile, perché a quel tempo non
era ancora di attualità la loro quadratura in un insieme a più dimensioni. Gli elementi finiti invece richiedono normalmente
l’integrazione di funzioni di 2, 3 e 4 variabili, per cui la tecnica di Gauss
ha dovuto essere generalizzata ricorrendo a quegli stessi polinomi di
Legendre su cui si basa la quadratura del grande matematico. Questo metodo tuttavia non affronta in modo diretto
l’integrazione di funzioni dipendenti da più grandezze, e quindi richiede un
numero di punti crescente esponenzialmente col numero di variabili in gioco.
Per questo motivo l’integrazione di Gauss-Legendre utilizza un numero di
punti antieconomico rispetto a quello rigorosamente necessario, producendo
così un aumento del tempo di calcolo, soprattutto negli elementi a 3
dimensioni e in quelli in cui è presente la dimensione temporale. LA NUOVA
QUADRATURA DIRETTA I notevoli tempi di calcolo di programmi sempre più
sofisticati, hanno imposto la necessità di affrontare in modo diretto
l’integrazione di funzioni definite negli elementi normalizzati di forma
quadrata nel caso di 2 variabili e di forma cubica negli elementi a 3
dimensioni. Tale quadratura è stata affrontata dallo scrivente per
la prima volta nel 1977 nella sua tesi di laurea integrando polinomi fino al
7° grado per funzioni di 2, 3 e 4 variabili, e in seguito ampliata
estendendone la ricerca a polinomi di grado superiore. Il nuovo metodo
consente di ottenere un risparmio di tempo, a parità di precisione, del 25%
per funzioni definite nel piano, del 50% per polinomi in 3 dimensioni, e del
77% per la quadratura di polinomi in 4 variabili. L’integrazione
di funzioni di 2 variabili reali Essa è stata ottenuta a partire dall’espansione di
Taylor delle funzioni di due variabili contenenti, oltre ai termini del tipo
xn e ym, anche
le componenti miste del tipo xpyq. Tale sommatoria
comprende la presenza di monomi antimetrici rispetto alla variabile x, del
tipo x2n+1, e alla variabile y della forma y2m+1, i cui
effetti si compensano nell’elemento normalizzato di forma quadrata collocando
l’origine degli assi nel suo centro, gli assi coordinati paralleli ai lati
del dominio, e disponendo i punti campione in modo simmetrico rispetto ai due
assi. Anche i termini antimetrici x2n+1 y2m e x2p
y2q+1 vengono annullati perché assumono valori opposti nei
semipiani simmetrici rispetto agli assi coordinati. Di conseguenza tutte le componenti dispari
dell’espansione di Taylor vengono eliminate nei calcoli collocando i punti
campione in posizioni simmetriche rispetto agli assi coordinati, in quanto
ogni valore positivo calcolato nel punto base di un quadrante viene
compensato da un uguale valore negativo dovuto al suo simmetrico situato o
nel semipiano opposto, o nel quadrante antimetrico. Per esempio, il termine x3
y2 produce valori positivi nel semipiano positivo delle x, ma
altrettanti valori opposti nel semipiano negativo. Invece il termine x3
y3 produce valori uguali nel I e III quadrante che vengono
eliminati dai risultati opposti ottenuti nel II e IV quadrante. Solo i termini pari del tipo x2m y2n
producono valori identici in ogni quadrante, e perciò è sufficiente
considerarne gli effetti indotti limitandoci al I quadrante individuato dai
valori positivi delle ascisse e delle ordinate. I punti campione su cui va calcolato il valore della
funzione sono in grado di risolvere identicamente sia i termini x2m y2n
che le componenti x2n y2m, a condizione che essi siano
disposti simmetricamente rispetto alla bisettrice del I e III quadrante,
collocazione che riduce il numero di incognite e di equazioni risolutrici nelle
quadrature di grado più elevato. Si deve precisare tuttavia che la quadratura rapida dei
polinomi di grado uguale o inferiore al quarto è stata possibile rinunciando
alla simmetria rispetto alla bisettrice, perché solo in questo modo è stato
limitato al minimo teorico il numero di punti campione. Questa circostanza ha
tuttavia comportato un incremento delle difficoltà di calcolo, dovendo
considerare un maggior numero di termini dell’espansione di Taylor al fine di
una risoluzione corretta del problema. Dal grado 6° in poi, la quadratura diretta è stata
invece risolta unicamente con una disposizione simmetrica dei punti sia
rispetto agli assi coordinati, sia rispetto alle loro bisettrici, rendendo
così possibile l’eliminazione dei termini x2m y2n una
volta risolta l’integrazione dei termini x2n y2m. Con queste premesse, la quadratura di un qualunque
polinomio di grado uguale o inferiore ad un numero prestabilito, richiede che
i punti campione siano scelti in modo da integrare esattamente i termini x2n,
y2m e le componenti miste x2py2q
fino al raggiungimento di un grado uguale a quello fissato a priori. Infatti,
se i punti campione, di cui non si conosce ancora la posizione e i relativi
coefficienti di peso, sono in grado di integrare esattamente i singoli monomi
di grado pari dell’espansione di Taylor, allora saranno in grado di integrare
qualsiasi combinazione lineare di questi ultimi, risolvendo così la
quadratura di un qualunque polinomio di grado uguale, ma anche inferiore, a
quello prestabilito. Ogni polinomio è infatti per definizione una somma di
monomi, per cui è possibile generarne uno qualsiasi attribuendo ai
coefficienti dei singoli termini i valori da noi scelti. La quadratura di un
qualsiasi polinomio è quindi possibile se i punti campione vengono collocati
in posizioni tali da essere in grado di integrare ciascuno dei suoi
componenti di grado pari. Ognuno dei termini indipendenti del polinomio produce
un’equazione non lineare in cui al primo membro si calcola la sommatoria,
estesa a tutti i punti, dei prodotti del valori assunti dal monomio
integrando in ogni punto, moltiplicati per il proprio coefficiente di peso.
Per esempio, il termine x2 y6 dovrà essere calcolato in
ciascuno dei punti base incogniti, e tali valori xi2 yi6
saranno moltiplicati per i coefficienti di peso Hi relativi ai
medesimi punti e ancora incogniti. Al secondo membro si colloca il valore assunto
dall’integrale del monomio x2 y6 nell’intervallo delle
x e delle y compreso fra zero e uno. Tale valore nel caso del monomio x2
y6 è pari a 1/21 ed è sempre un valore frazionario per qualsiasi
termine considerato. Estendendo l’integrazione a tutte le componenti di
grado pari, con l’esclusione di quelle simmetriche rispetto alla bisettrice
del I quadrante, si ottiene con questa procedura un sistema di equazioni non
lineari in cui le incognite sono rappresentate dalle coordinate (xi,yi)
e dai coefficienti di peso Hi di ciascun punto. Si deve poi precisare che nelle incognite del problema
compaiono pure alcuni punti campione situati nelle bisettrici degli assi, le
cui incognite si limitano alla sola coordinata xi e al
coefficiente di peso Hi, in quanto yi coincide con xi. Il coefficiente di peso Hi ha un significato
non solo algebrico, ma anche geometrico. Esso rappresenta infatti il volume
sotteso dalla superficie generata da un polinomio in (x, y), avente la
proprietà di assumere un valore unitario nel punto Pi e di
annullarsi in tutti gli altri. Tale polinomio è di grado uguale o anche
inferiore a quello del polinomio di grado massimo integrabile esattamente, ed
è costituito da un insieme di termini di grado sia pari che dispari, che
riportano a zero il valore assunto dal polinomio in tutti i punti base,
tranne che in quello in cui deve essere calcolato il valore della funzione
integranda. Se ora vogliamo quadrare una generica funzione in due
variabili di cui sono noti i valori assunti dalla medesima nei punti
prestabiliti, il computo prevede che dobbiamo sommare i prodotti Hi f(xi,yi)
valutati in ogni punto Pi. Questa operazione equivale alla sostituzione della
funzione integranda col polinomio che assume nei punti base gli stessi valori
della medesima, attribuendo così all’integrale ignoto della funzione data, il
corrispondente integrale del polinomio che intercetta i medesimi valori della
funzione nei punti campione. E’ chiaro che quanto più elevato sarà il numero di
punti di riferimento, quanto più il polinomio interpolante si confonderà con
la funzione integranda anche nel contorno dei medesimi punti, per cui col
crescere del numero di punti base, si incrementa pure la precisione della
quadratura. La differenza fra il valore dell’integrale esatto e di quello
così calcolato viene chiamata resto, e dipende dalle derivate parziali
(n+1)-esime della funzione integranda calcolate in un punto interno al
dominio di definizione. Come già precisato, il polinomio interpolante è a sua
volta la somma dei polinomi ottenuti moltiplicando per f(xi,yi)
ciascuno dei polinomi che assumono un valore nullo in tutti i punti, tranne che
in Pi, dove valgono 1. Questi nuovi polinomi ottenuti
moltiplicando per f(xi,yi) quelli generanti con la loro
cubatura il coefficiente di peso Hi, descrivono una superficie che
ha la proprietà di intercettare nel punto Pi lo stesso valore della
funzione integranda, mentre interseca il piano (x, y) in tutti gli altri
punti base, dove assume un valore nullo. In questo modo, con la sommatoria
dei termini Hi f(xi,yi), non si fa altro che
calcolare il volume del polinomio che riproduce esattamente nei punti
prestabiliti tutti i valori assunti dalla funzione integranda. La quadratura diretta delle funzioni di due variabili
non è altro che una generalizzazione nel piano (x,y) del metodo di Gauss, in
quanto affronta la risoluzione di un sistema di equazioni non lineari in cui
sia i coefficienti di peso, sia le coordinate dei punti campione
costituiscono le incognite del problema. La equazioni non lineari descritte dall’integrazione
gaussiana vengono riprese nella quadratura diretta, ampliando tuttavia sia la
loro quantità, sia il numero di incognite, che ora non si limitano più alle
posizioni xi dei punti base di una retta, ma anche alle coordinate
yi dei punti del piano di
cui devono essere stabilite entrambe le coordinate (xi,yi
). Questo comporta che a parità di grado, il numero di incognite e quindi di
equazioni, è superiore rispetto a quanto si presenti nella trattazione di una
sola variabile, anche se l’incremento non è così rapido come con l’utilizzo
dei polinomi di Legendre. Si presenta inoltre una nuova complicazione che ha
ritardato di molto la soluzione del problema che altrimenti avrebbe potuto
essere già risolto fin dagli esordi del metodo degli elementi finiti. Essa
riguarda la condizione tassativa di accettare unicamente quelle soluzioni che
riproducono unicamente valori positivi dei coefficienti di peso Hi.
Valori negativi sono indice di una scelta infelice delle coordinate dei punti
base e rendono instabile la soluzione creando problemi nei programmi di
calcolo. Nel metodo di Gauss, i coefficienti di peso sono per
loro natura tutti positivi, e questo è reso implicito dalle proprietà
intrinseche dei polinomi di Legendre su cui si fonda la quadratura. Ne
consegue che una volta impostate le equazioni costitutive, qualunque sia il
loro numero, i sistemi di equazioni che ne derivano hanno sempre soluzioni
positive che possono essere calcolate mediante un’equazione algebrica
associata di grado dipendente sia dal numero di punti base situati nel
semiasse positivo delle ascisse che eventualmente da quello collocato nella
loro origine. I metodi dell’analisi numerica non hanno perciò
difficoltà a risolvere il sistema di equazioni ottenuto dalla quadratura di
Gauss, perché il problema viene trasformato nella risoluzione di un’equazione
algebrica di una sola variabile che restituisce tante soluzioni positive
quante il numero di punti campione prestabilito. Diverso è il caso della risoluzione delle equazioni
relative all’espansione di Taylor di una funzione in due variabili. Queste
equazioni non generano automaticamente valori positivi dei coefficienti di
peso e basta anche la presenza di un solo Hi negativo per
invalidare l’intera quadratura. Questo problema si aggrava con l’aumento del
grado del polinomio interpolante e diventa di ardua difficoltà a partire dai
polinomi di 12° grado per funzioni di due variabili e già di 8° grado per
polinomi in 3 variabili. La causa di questo fatto è dovuta alle condizioni
sempre più sottili che devono essere rispettate, ma soprattutto individuate,
richieste nella collocazione dei punti base del I quadrante, allo scopo di
originare in seguito sequenze di risultati sempre positivi. Questi punti interferiscono in tutte le equazioni
risolutrici e sono presenti non solo nell’integrazione dei termini x2my2n,
ma anche nella quadratura dei monomi del tipo x2p e y2q
che si incontrano nello sviluppo in serie di Taylor. Le funzioni del piano governate da una sola variabile
come x2p o y2q sono delle superfici rigate, ottenute
quindi da un insieme di infinite rette parallele e perpendicolari ai piani
coordinati, intersecanti la parabola di grado corrispondente definita nel
piano (x, z) o (y, z). Queste superfici generano una cubatura nel piano (x,
y), e il nostro interesse riguarda l’integrale definito nel dominio di
definizione compreso tra 0 e 1 degli assi x e y. E’ chiaro che esse assumono
in ogni punto del piano e degli assi dei valori ben precisi che sono
calcolabili e successivamente moltiplicabili per i coefficienti di peso dei
punti base. I punti campione collocati sugli assi cartesiani ignorano i
termini del tipo x2ny2p e quindi non sono presenti
nelle equazioni relative alla loro integrazione. I punti situati sugli assi
sono invece rilevanti unicamente nell’integrazione dei termini del tipo x2p
e y2q e sono necessari per la quadratura di un polinomio in
più variabili. Se i punti Pi(xi yi)
non sono collocati nel numero e nella posizione più adeguata, saranno proprio
i coefficienti H0i dei punti situati sugli assi coordinati a
decadere più facilmente in quei valori negativi che invalidano l’intera
soluzione, evidenziando così la criticità delle scelte di base. Poiché il
loro valore viene calcolato in modo sequenziale solo dopo la risoluzione
delle relazioni dipendenti unicamente dai punti Pi(xi,yi),
sono proprio quelle contenenti solamente i termini x2my2n
a determinare la riuscita delle equazioni in cui compaiono anche i punti
appartenenti agli assi cartesiani. L’integrazione
delle funzioni di 3 variabili reali In questo tipo di integrazione sono presenti nuovi
termini del tipo x2my2nz2p che non
esistevano nei polinomi in due variabili e che ora richiedono l’introduzione
di nuovi punti campione nello spazio (x,y,z). Ancora una volta il problema
può essere semplificato ricorrendo a simmetrie che consentono di isolare le
equazioni nel diedro formato dai semipiani positivi delle variabili x, y e z. Mentre i problemi più sofisticati nel piano vengono
risolti da coppie di punti con coordinate simmetriche rispetto alla
bisettrice del primo quadrante del tipo PiA(xi,yi)
e PiB(yi, xi), nell’integrazione dei
polinomi più complessi in 3 variabili, si deve soddisfare una condizione di
triplice simmetria con insiemi di punti del tipo PiA(xi,yi,yi),
PiB(yi,xi,yi) e PiC(yi,yi,xi).
Queste terne sono in grado di integrare correttamente tutti i termini del
tipo x2my2nz2p quando sono associate ad
altri punti base situati sulla diagonale che congiunge l’origine degli assi
col vertice di coordinate x, y e z pari a (1,1,1). In questo modo il problema
si riconduce all’integrazione dei polinomi in 2 variabili, anche se ora il
numero di punti da considerare è notevolmente superiore, e le equazioni
risolutrici dipenderanno da un più elevato numero di incognite. Disporre di
una terna di punti nello spazio equivale a collocare 2 punti simmetrici
rispetto alla bisettrice del piano (x,y), oltre ad un terzo punto sulla
medesima retta. La quadratura prevede la contemporanea risoluzione
delle equazioni generate dai termini del tipo x2my2nz2p
e delle relazioni generate dalla proiezione nel piano (x,y) dei medesimi
punti base nello spazio (x,y,z). Pure nel problema spaziale devono infatti
essere risolte le equazioni in 2 variabili generate dai termini x2my2n,
e questa soluzione può essere ottenuta ricorrendo alle proiezioni sul piano
(x,y) sia delle terne PiA, PiB e PiC, che
dei punti situati sulla diagonale del cubo che pure vengono proiettati sulla
bisettrice del quadrato. La soluzione è resa possibile dall’aggiunta di un
certo numero di punti sulla medesima bisettrice del piano (x,y) che
consentono così di integrare esattamente anche i termini del tipo x2my2n.
Infine devono ancora una volta essere disposti altri punti sugli assi
coordinati, al fine di poter integrare le funzioni del tipo x2n, y2m,
e z2p. Anche nello spazio la principale difficoltà consiste
nel posizionare nel modo più opportuno la terna di punti simmetrici rispetto
ai tre piani cartesiani, il cui scopo non è solo quello di consentire
l’integrazione dei termini x2my2nz2p, ma
anche dei monomi x2my2n definiti nel piano (x,y). Solo una
disposizione corretta sarà in grado di produrre unicamente valori positivi
dei coefficienti di peso e quindi di rendere accettabile la soluzione. La validità della configurazione scelta dipenderà
dall’individuazione delle sottili e non sempre evidenti condizioni che consentono
di indurre una cascata di valori tutti positivi dei coefficienti di peso Hi,
sia nei punti situati nel piano (x,y), che sugli assi coordinati. Per questo
motivo l’integrazione di un polinomio in 3 variabili è, a parità di grado,
più impegnativa della quadratura delle funzioni in 2 variabili, ma la maggior
ostilità dei calcoli viene compensata da un più generoso risparmio di punti
rispetto alle funzioni di 2 variabili. Conclusioni In questo articolo sono stati presentati alcuni dei
principi costitutivi del nuovo metodo di integrazione che consente, a parità
di precisione, un risparmio di punti del 25% per funzioni di 2 variabili, del
50% per polinomi di 3 variabili, e del 77% per funzioni di variabili. Come già anticipato, il procedimento qui presentato non
è altro che una generalizzazione del metodo che Gauss studiò per
l’integrazione numerica delle funzioni di variabile reale. Le equazioni
risolutrici sono quindi simili a quelle determinate dalla quadratura
gaussiana, cui si aggiungono nuove relazioni che comprendono termini misti
del tipo x2my2n per le funzioni di 2 variabili, mentre
nei polinomi in 3 dimensioni sono stati pure considerati i monomi della forma
x2my2nz2p.
|
||||||||||||||||
ingegneriastrutturale.net -
Tutti i Diritti Riservati |