Consulting

Results 1 to 9 of 9

Thread: Problem with excel macro

  1. #1
    VBAX Newbie
    Joined
    Jul 2009
    Posts
    2
    Location

    Problem with excel macro

    Please help me

    [VBA]
    Private Sub CommandButton1_Click()
    ' Calcolo delle condizioni all'uscita del canale in rame
    Dim M1, M2, M3, M4, M5, M6 As Double
    Dim M3new, M5new As Double
    Dim L1, L2, L3, L4 As Double
    Dim L23, L56 As Double
    Dim RE As Double
    Dim eps As Double

    ' Dati geometrici
    ' Aree
    A1 = 0.0257 'Ingresso canale
    A2 = 0.00623 'Parte a sezione rettangolare 361
    A3 = 0.00951 'Area immediatamente a monte delle pale, calcolata su un piano parallelo alle stesse
    A4 = 0.00586 'Area di gola delle pale
    A5 = 0.00622 'Area immediatamente a valle delle pale, calcolata su un piano parallelo alle stesse
    A6 = 0.00361 'Area minima nel canale (piu' o meno all'altezza del gomito) Dai miei calcoli e' 0.00427!
    A7 = A1 'Area di uscita del canale
    ' Altezza canale (sezione di prova)
    H = 0.045 '45 mm
    ' Lunghezza canale (lungo l'ascissa curvilinea presa sull'asse) per calcolo perdite distribuite
    L23 = 0.172 * 1# 'Lunghezza tratto isotermo a monte delle pale
    L56 = 0.29 * 1# 'Lunghezza tratto isotermo a valle delle pale
    ' Diametri idraulici
    Per1 = 4# * A1 / 3.14159
    D1 = 4# * A1 / Per1
    Per23 = 2# * (H + A2 / H)
    D23 = 4# * A2 / Per23
    Per56 = 2# * (H + A6 / H)
    D56 = 4# * A6 / Per56
    D7 = D1
    ' Rugositą assoluta [m]
    eps = 0.00002 'Rugositą tipica di un tubo in rame
    '''' Fine dati geometrici ''''''
    '========================================================================== ==============
    ' Proprietą del fluido
    gam = 1.33
    r = 288#
    nmax = Range("E2")
    'print*,"Numero dati da elaborare?"
    'read*,nmax
    For niter = 1 To nmax
    ' Dati all'ingresso del canale
    'print*,"G, p1,T1"
    'read*, G,p1,T1
    G = Cells(1 + niter, 1)
    p1 = Cells(1 + niter, 2)
    T1 = Cells(1 + niter, 3)

    ' Calcolo tratto isentropico 1-2
    ro1 = p1 / (r * T1)
    u1 = G / (ro1 * A1)
    M1 = u1 / (gam * r * T1) ^ 0.5
    'print*,"M1",M1,ro1,p1,u1
    ' Rapporto Area/Area critica nella sez. 1
    A1sAc = area(M1, gam)
    ' Calcolo area critica tratto isentropico 1-2
    Ac12 = A1 / A1sAc
    If (A2 < Ac12) Then
    ' Se A2<Ac il flusso e' strozzato nella sez. 2
    MsgBox "Condizioni soniche raggiunte nella sez. 2: A2,Ac12,A2,Ac12"
    End If
    ' Grandezze totali in ingresso
    p01 = p1 * (1# + 0.5 * (gam - 1#) * M1 ^ 2) ^ (gam / (gam - 1#))
    T01 = T1 * (1# + 0.5 * (gam - 1#) * M1 ^ 2)
    ' Condizione di isentropia
    p02 = p01
    T02 = T01
    ' Area/area critica sez. 2
    A2sAc = A2 / Ac12
    ' Calcolo delle condizioni nella sezione 2
    StartMach = 0.1
    ' Calcolo del Mach 2 dalla relazione Area/Mach
    M2 = calcM(A2sAc, gam, StartMach)
    ' Calcolo delle condizioni 2
    p2 = p02 * (1# + 0.5 * (gam - 1#) * M2 ^ 2) ^ (-gam / (gam - 1#))
    T2 = T02 * (1# + 0.5 * (gam - 1#) * M2 ^ 2) ^ (-1#)
    u2 = M2 * (r * gam * T2) ^ 0.5
    'print*,"M2,p2,T2,u2"
    'print*,M2,p2,T2,u2
    ' Verifica conservazione della portata
    'print*,"TEST G2=",p2*u2*A2/(r*T2)
    '============ Fine tratto 1-2 ============
    '======= Perdita di carico per restringimento di sezione =======
    ro2 = p2 / (r * T2)
    coef = 0.1 'E' un po' piu' grande per tenere conto delle perdite distribuite
    ' densitą e velocita' costanti, solo perdita di pressione
    p2cor = p2 - 0.5 * coef * ro2 * u2 ^ 2
    ' u2cor = Sqr(2# * (1 - coef)) * u2 'Perdita di energia cinetica
    u2cor = u2
    T2cor = T2 * p2cor * u2cor / p2 / u2
    ' T2cor = T2
    ' u2cor = u2 * T2cor / T2 * p2 / p2cor
    corMach = u2cor / (gam * r * T2cor) ^ 0.5
    'print*,"corMach,p2cor,T2cor,u2"
    'print*,corMach,p2cor,T2cor,u2
    p2 = p2cor
    T2 = T2cor
    u2 = u2cor
    M2 = corMach
    ' Verifica conservazione della portata
    ' print*,"TEST G2cor=",p2*u2*A2/(r*T2)
    '========= Perdite distribuite tratto 2-3 : ipotesi di flusso isotermo con attrito ==========
    ' Calcolo del numero di Reynolds
    vis2 = visc(T2)
    RE23 = (ro2 * u2 * D23) / vis2
    RE = RE23
    ' print*,"Reynolds tratto 2-3",Re23
    ' Calcolo del coefficiente d'attrito

    f = amoody(RE23, eps / D23)

    ' Calcolo iterativo del rapporto di pressione

    'print*,M2,gam,f,L23,D23
    beta = calcbeta(M2, gam, f, L23, D23)
    'print*,beta
    p3 = beta * p2
    u3 = u2 / beta
    T3 = T2
    M3 = u3 / (r * gam * T3) ^ 0.5
    'print*,"M3,p3,T3,u3"
    'print*,M3,p3,T3,u3
    'print*,"TEST G3=",p3*u3*A2/(r*T3)
    '======== Fine perdite distribuite =========
    ' Calcolo tratto isentropico 3-5 a cavallo delle pale. L'area di gola e' A4
    ro3 = p3 / (r * T3)
    u3new = G / ro3 / A3 'A3 č l'area su un piano parallelo alle pale, quindi e' maggiore di A2
    M3new = u3new / (gam * r * T3) ^ 0.5
    ' Area/Area critica
    A3sAc = area(M3new, gam)
    ' Area critica tratto isentropico 3-5
    Ac35 = A3 / A3sAc
    If (Ac35 > A4) Then
    ' print*,"Area critica",Ac35,
    '& "maggiore area di gola A4",A4
    ' print*,"Flusso supersonico nel divergente!!"
    StarMach = 3#
    Else
    ' print*,"Flusso subsonico nel divergente"
    StartMach = 0.001
    End If

    ' Grandezze totali sez. 3
    p03 = p3 * (1# + 0.5 * (gam - 1#) * M3new ^ 2) ^ (gam / (gam - 1#))
    T03 = T3 * (1# + 0.5 * (gam - 1#) * M3new ^ 2)
    ' Condizione di isentropia
    p05 = p03
    T05 = T03
    ' Area/Area critica sezione 5
    A5sAc = A5 / Ac35
    ' Calcolo delle condizioni nella sezione 5 subito a valle delle pale
    ' Calcolo Mach dalla relazione Area/MAch
    M5 = calcM(A5sAc, gam, StartMach)
    ' Condizioni sez. 5
    p5 = p05 * (1# + 0.5 * (gam - 1#) * M5 ^ 2) ^ (-gam / (gam - 1#))
    T5 = T05 * (1# + 0.5 * (gam - 1#) * M5 ^ 2) ^ (-1#)
    u5 = M5 * (r * gam * T5) ^ 0.5
    'print*,"M5,p5,T5,u5"
    'print*,M5,p5,T5,u5
    ' Verifica conservazione della portata
    'print*,"TEST G5=",p5*u5*A5/(r*T5)
    '======== Fine tratto 3-5 =========
    '======== Perdita di carico per gomito all'altezza delle palette =========
    ro5 = p5 / (r * T5)
    coef = 0.7 'tiene conto anche delle perdite distribuite del tratto 5-6
    ' densitą e velocita' costanti, solo perdita di pressione
    p5cor = p5 - 0.5 * coef * ro5 * u5 ^ 2 'perdita di pressione
    ' u5cor = Sqr(2# * (1# - coef)) * u5 'perdita di energia cinetica
    u5cor = u5
    T5cor = T5 * p5cor * u5cor / p5 / u5
    ' T5cor = T5
    ' u5cor = u5 * T5cor / T5 * p5 / p5cor
    corMach = u5cor / (gam * r * T5cor) ^ 0.5
    'print*,"corMach,p5cor,T5cor,u5"
    'print*,corMach,p5cor,T5cor,u5
    p5 = p5cor
    T5 = T5cor
    u5 = u5cor
    M5 = corMach

    ' Verifica conservazione della portata
    ' print*,"TEST G5cor=",p5*u5*A5/(r*T5)
    '======== Perdita di carico distribuita tratto 5-6 a valle delle palette : ipotesi di flusso isotermo =========
    ' Calcolo del numero di Reynolds
    ro5 = p5 / (r * T5)
    u5new = G / (ro5 * A6) 'prendo l'area trasversale, minore dell'area // alle pale
    M5new = u5new / (gam * r * T5) ^ 0.5
    vis5 = visc(T5)
    Re56 = (ro5 * u5 * D6) / vis5
    RE = Re56
    ' print*,"Reynolds tratto 5-6",Re56
    ' Calcolo del coefficiente d'attrito

    f = amoody(Re56, eps / D56)
    ' Calcolo iterativo del rapporto di pressione
    beta = calcbeta(M5new, gam, f, L56, D56)
    p6 = beta * p5
    u6 = u5new / beta
    T6 = T5
    M6 = u6 / (r * gam * T6) ^ 0.5
    'print*,"M6,p6,T6,u6"
    'print*,M6,p6,T6,u6
    ' Verifica conservazione della portata ERRORE!!!!
    'print*,"TEST G6=",p6*u6*A6/(r*T6)
    '======== Fine perdite distribuite =========

    '========= Perdita di carico dovuta al distacco di vena a valle del primo gomito ========
    ro6 = p6 / (r * T6)
    coef = 1#
    p6cor = p6 - 0.5 * coef * ro6 * u6 ^ 2 'perdita di pressione
    ' u6cor = Sqr(2# * (1# - coef)) * u6 'perdita di energia cinetica
    u6cor = u6
    T6cor = T6 * p6cor * u6cor / p6 / u6
    ' T6cor = T6
    ' u6cor = u6 * T6cor / T6 * p6 / p6cor
    corMach = u6cor / (gam * r * T6cor) ^ 0.5
    'print*,"corMach,p6cor,T6cor,u6"
    'print*,corMach,p6cor,T6cor,u6
    p6 = p6cor
    T6 = T6cor
    u6 = u6cor
    M6 = corMach
    '========================================================================== ======
    '== ==
    '== Attenzione: a causa del distacco di vena a valle del gomito, ==
    '== il successivo allargamento di sezione non ha praticamente alcun effetto ==
    '== sulla pressione, che resta praticamente costante (zona di ricircolo) ==
    '== per questa ragione, il calcolo del divergente a partire dall'ipotesi di ==
    '== isentropia del flusso e' totalmente inappropriata. E' piu' accurato ==
    '== ipotizzare che le condizioni del flusso restino invariate fino all'uscita ==
    '== ==
    '========================================================================== ======
    ' Verifica conservazione della portata
    'print*,"TEST G6cor=",p6*u6*A6/(r*T6)
    Worksheets("Risultati").Range("A25") = G
    Worksheets("Risultati").Cells(2 + niter, 1) = G
    Worksheets("Risultati").Range("B25") = p1
    Worksheets("Risultati").Cells(2 + niter, 2) = p1
    Worksheets("Risultati").Range("C25") = T1
    Worksheets("Risultati").Cells(2 + niter, 3) = T1
    Worksheets("Risultati").Range("D25") = p6
    Worksheets("Risultati").Cells(2 + niter, 4) = p6
    Worksheets("Risultati").Range("E25") = M1
    Worksheets("Risultati").Cells(2 + niter, 5) = M1
    Worksheets("Risultati").Range("F25") = M2
    Worksheets("Risultati").Cells(2 + niter, 6) = M2
    Worksheets("Risultati").Range("G25") = M3
    Worksheets("Risultati").Cells(2 + niter, 7) = M3
    Worksheets("Risultati").Range("H25") = M5
    Worksheets("Risultati").Cells(2 + niter, 8) = M5
    Worksheets("Risultati").Range("I25") = M6
    Worksheets("Risultati").Cells(2 + niter, 9) = M6
    portrid = G * Sqr(T01) / p01 * 1000
    Worksheets("Risultati").Range("L25") = portrid
    Worksheets("Risultati").Cells(2 + niter, 10) = portrid
    rapp = p01 / p6
    Worksheets("Risultati").Range("M25") = rapp
    Worksheets("Risultati").Cells(2 + niter, 11) = rapp

    'write(10,*) G,p1,T1,p6,T6,M1,M2,M3,M5,M6
    'write(20,*) G*sqrt(T01)/p01*1000,p01/p6
    Next niter
    End Sub
    '======================
    '===== FUNZIONI =======
    '======================
    Function area(M, gam)
    ' Relazione Area/Mach
    gamp1 = gam + 1#
    gamm1 = gam - 1#
    area = (2# / gamp1 * (1# + gamm1 / 2# * M ^ 2)) ^ (gamp1 / (2# * gamm1)) / M
    End Function

    '===============================
    Function calcM(AsAc, gam, StartMach)
    ' Calcolo inverso del Mach dalla legge delle aree

    gamp1 = gam + 1#
    gamm1 = gam - 1#

    ' print*,"A/A*",AsAc
    M = StartMach
    dM = 10000000000#

    ' Calcolo iterativo con Newton-Raphson
    Do

    f = (2# / gamp1 * (1# + gamm1 / 2# * M ^ 2)) ^ (gamp1 / (2# * gamm1)) / M - AsAc
    fp = -(2# / gamp1 * (1# + gamm1 / 2# * M ^ 2)) ^ (gamp1 / (2# * gamm1)) / M ^ 2 + (2# / gamp1 * (1# + gamm1 / 2# * M ^ 2)) ^ (gamp1 / (2# * gamm1) - 1)
    dM = -f / (fp + 0.00000001)
    M = M + dM
    ' print*,M,dM
    Loop Until (Abs(dM) > 0.00001)
    calcM = M
    End Function
    '====================================================================

    Function visc(T)
    ' Calcolo della viscositą dalla formula di Sutherland

    T0 = 288.15
    S = 110#
    visc = 0.0000178 * (T / T0) ^ 1.5 * (T0 + S) / (T + S)
    End Function
    '============================
    Function amoody(RE, eos)
    ' Calcolo del coefficiente d'attrito dalla formula di Colebrook
    ' Valore stimato dalla formula di Swamee-Jain


    f = 0.25 / ((Log(eos / 3.72 + 5.74 / RE ^ 0.9)) / 2.302585) ^ 2

    ' Calcolo iterativo
    usrf = 1# / Sqr(f)
    Do
    usrf1 = usrf
    usrf = -2# * (Log(2.51 * usfr1 / RE + eps / 3.72)) / 2.302585
    ' print*,usrf,usrf1
    Loop Until Abs(usrf1 - usrf) > 0.000001

    amoody = usrf ^ (-2#)

    End Function
    '==========================
    Function calcbeta(M, gam, f, L, D)
    ' Calcolo del rapporto di pressione per flusso isotermo


    '10 continue
    beta = 0#
    beta1 = 1#
    nit = 0
    Do
    nit = nit + 1
    beta = beta1
    B = Log(beta ^ (-1#))
    A = 1# - 2# * gam * (M ^ 2) * (0.5 * B + 2# * f * L / D)
    If A > 0.8 Then beta1 = A
    If A < 0.8 Then beta1 = 0.8
    ' Limitatore inserito per evitare problemi in caso di divergenza

    Loop Until (Abs(beta - beta1) > 0.0001)
    'print*,sqrt(beta)
    calcbeta = Sqr(beta)

    End Function
    '========================================================================== ===================
    [/VBA]

  2. #2
    Administrator
    VP-Knowledge Base
    VBAX Grand Master mdmackillop's Avatar
    Joined
    May 2004
    Location
    Scotland
    Posts
    14,478
    Location
    Welcome to VBAX.
    Please explain the problem.
    MVP (Excel 2008-2010)

    Post a workbook with sample data and layout if you want a quicker solution.


    To help indent your macros try Smart Indent

    Please remember to mark threads 'Solved'

  3. #3
    VBAX Newbie
    Joined
    Jul 2009
    Posts
    2
    Location
    thiis the part of code that gives problem:
    cccccc Calcolo del coefficiente d'attrito dalla formula di Colebrook ccc
    function amoody(Re,eps)
    cc Valore stimato dalla formula di Swamee-Jain
    f=0.25/(log10(eos/3.72+5.74/Re**0.9)**2)
    cc Calcolo iterativo
    usrf=1./sqrt(f)
    do while(abs(usrf1-usrf).gt.1e-6)

    usrf1=usrf
    usrf=-2.*log10(2.51*usfr1/Re+eps/3.72)
    c print*,usrf,usrf1
    enddo

    amoody=usrf**(-2.)
    return
    end

  4. #4
    VBAX Newbie
    Joined
    Oct 2011
    Posts
    3
    Location
    I am newbie of this game.I feel good.We can often talk about games

  5. #5
    VBAX Newbie
    Joined
    Jun 2012
    Posts
    4
    Location
    so complicated

  6. #6
    Ciao Andrea,
    come possiamo aiutarti??

  7. #7
    long code ... difficult to understand !

  8. #8
    complicated code

  9. #9
    VBAX Regular
    Joined
    Feb 2013
    Posts
    10
    Location
    Uff!

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •