-
.
Riporto i codici delle varie funzioni di BESSEL di vario tipo ed ordine, riprese dalla libreria di Leo Volpi. CODICEFunction BesselJx(x, Optional n)
'********************************************************************
'Funzione di Bessel di prima specie e di ordine n, Jn(x)
'********************************************************************
Dim n1#, BJ0#, DJ0#, BJ1#, DJ1#, BY0#, DY0#, BY1#, DY1#
Dim NM#, BJ#(), dj#(), by#(), dY#()
If IsMissing(n) Then n1 = 0 Else n1 = CDbl(n)
If n1 <= 1 Then
Call JY01A(x, BJ0, DJ0, BJ1, DJ1, BY0, DY0, BY1, DY1)
If n1 = 0 Then BesselJx = BJ0 Else BesselJx = BJ1
Else
Call JYNA(n1, x, NM, BJ, dj, by, dY)
BesselJx = BJ(n1)
End If
End Function
Function BesselYx(x, Optional n)
'********************************************************************
'Funzione di Bessel di seconda specie e di ordine n, Yn(x)
'********************************************************************
Dim n1#, BJ0#, DJ0#, BJ1#, DJ1#, BY0#, DY0#, BY1#, DY1#
Dim NM#, BJ#(), dj#(), by#(), dY#()
If IsMissing(n) Then n1 = 0 Else n1 = CDbl(n)
If n1 <= 1 Then
Call JY01A(x, BJ0, DJ0, BJ1, DJ1, BY0, DY0, BY1, DY1)
If n1 = 0 Then BesselYx = BY0 Else BesselYx = BY1
Else
Call JYNA(n1, x, NM, BJ, dj, by, dY)
BesselYx = by(n1)
End If
End Function
Function BesseldJ(x, Optional n)
'********************************************************************
'Derivata della Funzione di Bessel di prima specie e di ordine n, J'n(x)
'********************************************************************
Dim x1#, n1#, BJ0#, DJ0#, BJ1#, DJ1#, BY0#, DY0#, BY1#, DY1#
Dim NM#, BJ#(), dj#(), by#(), dY#()
If IsMissing(n) Then n1 = 0 Else n1 = CDbl(n)
x1 = CDbl(x)
If n1 <= 1 Then
Call JY01A(x1, BJ0, DJ0, BJ1, DJ1, BY0, DY0, BY1, DY1)
If n1 = 0 Then BesseldJ = DJ0 Else BesseldJ = DJ1
Else
Call JYNA(n1, x1, NM, BJ, dj, by, dY)
BesseldJ = dj(n1)
End If
End Function
Function BesseldY(x, Optional n)
'********************************************************************
'Derivata della Funzione di Bessel di seconda specie e di ordine n, Y'n(x)
'********************************************************************
Dim x1#, n1#, BJ0#, DJ0#, BJ1#, DJ1#, BY0#, DY0#, BY1#, DY1#
Dim NM#, BJ#(), dj#(), by#(), dY#()
If IsMissing(n) Then n1 = 0 Else n1 = CDbl(n)
x1 = CDbl(x)
If n1 <= 1 Then
Call JY01A(x1, BJ0, DJ0, BJ1, DJ1, BY0, DY0, BY1, DY1)
If n1 = 0 Then BesseldY = DY0 Else BesseldY = DY1
Else
Call JYNA(n1, x1, NM, BJ, dj, by, dY)
BesseldY = dY(n1)
End If
End Function
Function BesselIx(x, Optional n)
'********************************************************************
'Funzione di Bessel modificata di prima specie e di ordine n, In(x)
'********************************************************************
Dim n1#, BI0#, DI0#, BI1#, DI1#, BK0#, DK0#, BK1#, DK1#
Dim NM#, bi#(), Di#(), BK#(), DK#()
If IsMissing(n) Then n1 = 0 Else n1 = CDbl(n)
If n1 <= 1 Then
Call IK01A(x, BI0, DI0, BI1, DI1, BK0, DK0, BK1, DK1)
If n1 = 0 Then BesselIx = BI0 Else BesselIx = BI1
Else
Call IKNA(n1, x, NM, bi, Di, BK, DK)
BesselIx = bi(n1)
End If
End Function
Function BesselKx(x, Optional n)
'********************************************************************
'Funzione di Bessel modificata di seconda specie e di ordine n, Kn(x)
'********************************************************************
Dim n1#, BI0#, DI0#, BI1#, DI1#, BK0#, DK0#, BK1#, DK1#
Dim NM#, bi#(), Di#(), BK#(), DK#()
If IsMissing(n) Then n1 = 0 Else n1 = CDbl(n)
If n1 <= 1 Then
Call IK01A(x, BI0, DI0, BI1, DI1, BK0, DK0, BK1, DK1)
If n1 = 0 Then BesselKx = BK0 Else BesselKx = BK1
Else
Call IKNA(n1, x, NM, bi, Di, BK, DK)
BesselKx = BK(n1)
End If
End Function
Function BesseldI(x, Optional n)
'********************************************************************
'Derivata della Funzione di Bessel modificata di prima specie e di ordine n, I'n(x)
'********************************************************************
Dim x1#, n1#, BI0#, DI0#, BI1#, DI1#, BK0#, DK0#, BK1#, DK1#
Dim NM#, bi#(), Di#(), BK#(), DK#()
If IsMissing(n) Then n1 = 0 Else n1 = CDbl(n)
x1 = CDbl(x)
If n1 <= 1 Then
Call IK01A(x1, BI0, DI0, BI1, DI1, BK0, DK0, BK1, DK1)
If n1 = 0 Then BesseldI = DI0 Else BesseldI = DI1
Else
Call IKNA(n1, x1, NM, bi, Di, BK, DK)
BesseldI = Di(n1)
End If
End Function
Function BesseldK(x, Optional n)
'********************************************************************
'Derivata della Funzione di Bessel modificata di seconda specie e di ordine n, K'n(x)
'********************************************************************
Dim x1#, n1#, BI0#, DI0#, BI1#, DI1#, BK0#, DK0#, BK1#, DK1#
Dim NM#, bi#(), Di#(), BK#(), DK#()
If IsMissing(n) Then n1 = 0 Else n1 = CDbl(n)
x1 = CDbl(x)
If n1 <= 1 Then
Call IK01A(x1, BI0, DI0, BI1, DI1, BK0, DK0, BK1, DK1)
If n1 = 0 Then BesseldK = DK0 Else BesseldK = DK1
Else
Call IKNA(n1, x1, NM, bi, Di, BK, DK)
BesseldK = DK(n1)
End If
End Function
In questa prima parte di codice l'ordine viene fornito come input e quindi non c'è la distinzione delle funzioni con il loro ordine, ma solo la distinzione per tipo e specie.
Assumono quindi la funzione di "selettore". In realta i valori della funzioni vengono calcolate dalla funzioni che richiamano dopo aver fatto la selezione.
Per esempio considerando la prima funzione, icodice controlla se il parametro n passato è 0 1 o mancante oppure >1.
nel primo caso rinvia alla funzione Call JY01A nel secondo caso rinvia alla funzione JYNA
E' notevole notare che queste due funzioni calcolano sia la funzione J e sia la funzione Y quindi sia quella di prima specie e sia quella di seconda specie.
Vi ricordo la seguente distinzione:
J e Y bessel ordinarie (DJ e DY derivate bessel ordinarie)
I e K bessel modificate (DI e DK derivate bessel modificate)
poi:
J ed I funzioni di bessel di prima specie
Y e K funzioni di bessel di seconda specie
l'ordine invece viene stabilito dal parametro "n" passato alla funzione"
---
il parametro deve essere un intero positivo (o anche nullo nel caso di ordine 0)
Da notare che il parametro può essere qualsiasi reale. Riporto a proposito quanto scritto su Wikipedia (e sui documenti matematici)
" The most important cases are when alpha is an integer or half-integer. Bessel functions for integer alpha are also known as cylinder functions or the cylindrical harmonics because they appear in the solution to Laplace's equation in cylindrical coordinates. Spherical Bessel functions with half-integer alpha are obtained when solving the Helmholtz equation in spherical coordinates.
(Alpha= n)
----
Nel caso di parametro pari a mezzo-intero si tratta delle funzioni sferiche di bessel.
Edited by afazio - 11/11/2023, 10:33 -
.la formulazione analitica di queste curve dovrebbe essere la seguente (estratto dall'art. originale di goyal & chopra):
...
Hai il link da dove poter scaricare l'articolo originale? -
.la formulazione analitica di queste curve dovrebbe essere la seguente (estratto dall'art. originale di goyal & chopra):
c'è una somma infinita, ma è chiaro che sono sufficienti i primi termini.
Esatto.
La formula è proprio quella (che ho anche riportato).
La somma infinita ai fini pratici la si può arrestare al terzo termine. La mia funzione ne considera 50.
Sulle funzioni di Bessel esiste una discreta confusione specie sulle denominazioni delle stesse in Excel.
Le funzioni di Bessel sono le soluzioni (parte delle soluzioni) di due equazioni differenziali del secondo ordine (che non scrivo) aventi un parametro reale n.
Una prima equazione è l'equazione di Bessel ordinaria, mentre l'altra che considera variabili immaginarie, viene detta equazione di Bessel modificata.
Da qui nasce la distinzione tra funzione di Bessel (ordinarie) e funzioni di Bessel modificate.
Le soluzioni delle equazioni di Bessel sono del tipo:
equazione ordinaria:
y( x ) = a*Jn( x ) +b*Yn( x )
equazione modificata:
y( x ) = c*In( x ) +d*Kn( X )
Quindi sono delle combinazioni lineari delle funzioni Jn( x ) e Yn( x ) per le ordinarie e In( x ) e Kn( x ) per le modificate.
Jn( x ) e In( x ) vengono dette di prima specie di ordine n
Yn( x ) e Kn( x ) vengono dette di seconda specie in ordine n
Esistono poi le funzioni derivate che prendono lo stesso nome delle funzioni di partenza
Il nome delle funzioni di Bessel è stato mantenuto in excel, quindi
BESSELK( x,0 ) rappresenta la funzione di Bessel modificata di seconda specie e di ordine 0
e cosi via
Edited by afazio - 11/11/2023, 10:56 -
.
Ecco qui il risultato. Direi accettabile. I colori che vedete sono quelli settati per le serie di dati che occultano i colori originari delle curve.
D'altra parte se ci pensiamo bene questo diagramma potrebbe essere anche sostituito da una semplice tabella numerica.
A noi invece interessa proprio la codificazione del diagramma per poterlo usare in qualsiasi nostra procedure automatizzata. -
.
Il diagramma fornito dagli autori per la determinazione delle masse aggiunte è il seguente allegato (a colori).
IN ascissa il rapporto tra massa e massa_infinita variabile da 0 a 1, in ordinata il rapporto z/Ho variabile da 0 a 1
La funzione VBA deve essere tale da riprodurre le medesime curve su un diagramma a dispersioni di punti.
Ho pensato pertanto di ritagliare con la massima precisione possibile, il diagramma di Goyal-Chopra, di incollarlo come sfondo in un grafico dispersione di punti in excel, settare gli assi da 0 a 1 (entrambi) e plottarvi sopra le curve ricavate utilizzando la funzione VBA
Ho ottenuto la perfetta corrispondenza tra le curve. -
.
IL codice VBA che ho scritto, sfrutta le funzioni presenti in Excel per il calcolo delle funzioni di Bessel mediante le istruzioni:
WorksheetFunction.BesselK(x, n)
Finchè il codice resta in ambiente Excel tutto continua a funzionare, ma se si vuole portare il codice in altri ambienti mediante un porting, è necessario avere a disposizione il codice che calcola le varie funzioni di Bessel, per poterlo tradurre.
Ricordavo di avere scaricato tempo fa un file add-in per excel prodotto dal Fox-team di Leo Volpi, che riportava il codice di un centinaio di funzioni speciali tra cui proprio le funzioni di Bessel.
Sfruttando questi codici ho voluto calcolare le masse aggiunte i sia con le funzioni di excel e sia con le funzioni di Leo Volpi ottenendo i medesimi risultati.
Il codice di Leo Volpi lo troverete all'interno del file che pubblicherò alla fine della discussione.
Intanto da una ricerca odierna su codici relativi alle funzioni di Bessel ho trovato la seguente pagine che voglio linkare.
https://people.math.sc.edu/Burkardt/f_src/...l_functions.f90 -
.Visto che il "cuore" della faccenda è simile ai serbatoi, appena è possibile mi piacerebbe adattare un mio script a questa casistica e fare qualche modellino.
Dubito che mi chiederanno mai una simile opera, però è un argomento che stuzzica quindi seguo con interesse.
Già che ci sei, volevo invece chiederti sulla possibilità in Enexys delle "section cut".
Ho visto ma non ancora sperimentato che forse è possibile avendo modellato una parete (lineare) con elementi shell, di poter determinare a valle dell'analisi le sollecitazioni risultanti sull'insieme di elementi che costituiscono le pareti.
Non sono riuscito a fare la stessa cosa se le pareti hanno sezione curvilinea, nello specifico per le sezioni circolari cave schematizzate con elementi shell.
Esiste la possibiltià o un qualche script in Enexsys che consente di eseguire le "section cut" presenti in altri software?
Non mettere limiti al cielo.
Potrei chiedertelo io come forma di collaborazione. -
.
La Funzione Goyal()
La funzione che ci fornisce il fattore adimensionale m/(rho*PI*ro²) è la seguenteCODICEOption Explicit
Function Goyal(ro As Double, ho As Double, z As Double) As Double
Dim PI As Double
Dim i As Integer 'contatore per la sommatoria
Dim rh As Double ' rapporto ro/Ho
Dim alfaM As Double ' coefficiente per funzione di Bessel
Dim Em As Double ' fattore con le funzioni di bessel
Dim Chopra As Double
PI = 4 * Atn(1)
rh = ro / ho
Chopra = 0
For i = 1 To 50
alfaM = (2 * i - 1) * PI / 2
Em = WorksheetFunction.BesselK(alfaM * rh, 1) / (WorksheetFunction.BesselK(alfaM * rh, 0) + WorksheetFunction.BesselK(alfaM * rh, 2))
Chopra = Chopra + (-1) ^ (i - 1) / (2 * i - 1) ^ 2 * Em * Cos(alfaM * z / ho)
Next i
Goyal = 16 / PI ^ 2 * ho / ro * Chopra
End Function
E' sufficiente fornire per ogni nodo immerso della torre i valore fissi di ro e di Ho e il valore variabile di z (del nodo) per ottenere il diagramma pseudoparabolico della variazione delle masse aggiunte in altezza.
Ricordiamo che sono ancora masse per unità di lunghezza ossia kg/m. Occorre trasformarle in masse concentrate e per questo considereremo la lunghezza di pertinenza di ciascun nodo come la semisomma delle due lunghezza sottostate e sovrastante.
Prossimamente il procedimento che ho usato per verificare che i valori ottenuti con la formulazione coincidono con quelli eventualmente letti dal diagramma. -
.
Aggiungeremo le "Masse Aggiunte" come masse concentrate nei nodi (lumped mass) in cui abbiamo discretizzato la torre anzichè come masse distribuite in altezza.
--
Molti programmi non hanno la possibilità di dare in input delle masse, ma consentono solo input di forze. Altri come il SAP2000 o anche Strauss consentono l'input di masse.
Occorre tenere però presente che le masse idrodinamiche aggiunte non hanno effetti gravitazionali ma sono semplicemente masse aggiunte alle masse sismiche che producono quindi azioni sismiche orizzontali.
Per intenderci queste masse vanno considerate come masse associate con coefficiente di partecipazione unitario per la valutazione degli effetti dell'azione sismica (vedi 2.5.3 del DM2018), mentre non intervengono nelle varia combinazioni sia di tipo statico e sia di tipo sismico.
Personalmente come strumento di analisi avevo a disposizione Enexys e per risolvere la questione dell'input delle masse non consentito, dopo consulto con il caro Alex Drake, ho dovuto trasformare le "masse idrodinamiche aggiunte" in carichi gravitazionali (moltiplicandoli per g) e quindi inserirli come carichi gravitazionali partecipanti ai fini degli effetti sismici e assumere per essi coefficienti parziali di combinazione nulli. -
.
Consideriamo lo schema della torre della immagine allegata al presente messaggio.
siano a titolo di esempio:
H = altezza complessiva della torre = 60 m
Ho = quota rispetto al basamento del massimo livello di regolazione = 50.00 m
ro = raggio della sezione (vuoto per pieno) della torre = 4.0 m
Per semplicità in questo esercizio si assume una torre cilindrica. Nel caso di sezione diversa occorre prima ricavare la sezione equivalente.
z= la generica quota rispetto al basamento della torre in cui vogliamo applicare la massa aggiunta
m= valore della massa aggiunta alla quota z
Supponiamo di voler determinare la massa aggiunta da collocare alla quota z=30,00m
I parametri che entrano in gioco nella formula di Goyal e Chupra sono:
ro/Ho=4.00/50.00 =0.08
Notiamo che, purtroppo, la curva relativa a questo valore non è stata diagrammata (vi sono quella relativa a 0.05 e quella relativa a 0.10)
il rapporto z/Ho vale 30.0/50.0 =0.6
Entreremmo nel diagramma proprio in corrispondenza di questo valore, orizzontalmente fino ad intercettare la curva relativa a ro/Ho e quindi verticalmente fino a leggere il valore della massa adimensionalizzata da applicare.
La massa è adimensionalizzata in rapporto alla massa denominata "massa per torre infinita" il cui valore è:
Minf= rho*pi.greco*ro²
quindi basta moltiplicare il valore letto per quest'ultimo valore per ottenere la massa in kg/m da aggiungere alla quota prefissata.
Edited by afazio - 5/11/2023, 09:38 -
.
Vi sottopongo la seguente considerazione relativamente alle frequenze di strutture immerse in un "fluido". -
.
Naturalmente in una procedura automatizzata o semi-automatica, non possiamo pensare di leggere i valori adimensionalizzati dal grafico proposto in funzione della quota e di altri parametri.
E poi se ci capita un rapporto ro/Ho = 0.55 dovremmo anche interpolare tra le curva relativa a 0.50 e quella relativa a 0.60. Una roba non proponibile se non esclusivamente per esercizio didattico.
E' necessaria una qualche automazione che ci fornisca direttamente il valore richiesto per fissati valori degli altri parametri.
Nelle prossime puntate dopo aver analizzato i parametri che entrano nelle formula vedremo come fare. -
.
Eccole!. Rispuntano le somme infinite di prodotti che coinvolgono le funzioni di Bessel di ordini e specie varie. Come accadeva nei serbatoi cilindrici e rettangolari.
In realtà le sommatorie che teoricamente vanno da 1 a infinito, possono arrestarsi ai fini pratici anche al terzo termine.
Il diagramma di cui si fa cenno nel precedente messaggio è il seguente:
Edited by afazio - 4/11/2023, 13:12 -
.
Per tenere conto delle azioni idrodinamiche su una torre di presa in un lago, la teoria più accreditata ad oggi è quella elaborata da Goyal e Chupra del 1989. e basata sul concetto di "masse idrodinamiche aggiunte".
L’espressione che fornisce le “masse aggiunte” è stata ottenuta sulla base della soluzione analitica dell’equazione di Laplace che in coordinate polari assume la forma:
(∂^2 p)/(∂r^2 )+1/r ∂p/∂r+1/r^2 (∂^2 p)/(∂ϑ^2 )+(∂^2 p)/(∂z^2 )=0
In cui
r,ϑ e z sono le coordinate cilindriche con z misurato a partire dalla base della torre
p(r,ϑ,z) è la pressione del liquido nel punto.
Per una torre cilindrica, la soluzione dell’equazione di Laplace che fornisce il valore della “massa aggiunta” alla quota generica z, è data dalla seguente espressione riportata nella seguente immagine: -
.
Qualcuno di voi ricorderà questa formula o qualcosa di assai simile nel calcolo delle azioni sismiche dell'acqua a tergo dei muri di sostegno.
Ricordo vagamente che questa formula era utilizzata sui programmi per muri della sts.
Altri ravvederanno in questa formula una qualche similitudine con la formula prevista dalla attuale normativa sulla progettazione e costruzione delle dighe.
Il Bar dell'Ingegneria
Posts written by afazio |
