Attribute VB_Name = "CellFLow" ' *******CellFlow, an n-dimensional CFD code******** 'G. Westendorp, 2003 'The dimensions x, y, z, ... are numbered from 1 to nD 'The first coordinate, or x-coordinate , 'corresponds to rows in a spreadsheet. 'A geometry of 1 row will result in a 1D problem 'The second coordinate, or y-coordinate is the row-direction of a spreadsheet. 'All higher dimensions are then specified 'Each layermust be terminated by a letter "L", optionally followed by a number 'the number is the hyperslice that ends after this layer. Default is 3. 'Below is an example of 2X2X2X2 hypercube ' x x ' x x 'L End of first 2-slice (Could have put L3 instead of L) ' x x ' x x 'L ' x x ' x x 'L4 End of first 3-slice ' x x ' x x 'L ' x x ' x x 'L ' x x ' x x 'L4 '***conventions used with indices*** 'i,j direction indices (x,y,z...) 'n Ordinal of point 'pt pseudo time, in case of state variables '***Index Order:*** 'ARRAY(direction,ordinal,pseudotimestep) 'coordinate to Ordinal 'n = dm(1) + ( dm(2)-1) * dnghbr(2)) + ( dm(3) - 1) * dnghbr(3)) '*********Change list since july 2008*************************** ' july 2008: 'Error corrected in '.out' file: Columns and rows incorrect 'heat flux as display variable 7 (was cell pitch) 'aug 2008: 'Export facility kinetic energy spectrum '*************************************************************** '***Global variables*** Option Explicit Public padnaam As String Private infilenaam As String Public mijnfilenaam As String Private outfilenaam As String Private ffilenaam As String Public GraphTitle As String Private dt As Single 'timestep Private ut As Long 'screen update interval Private st As Long 'autosave update interval Private ft As Long 'autobmp file create interval Private logt As Long 'logfile update interval Private tRamp As Long ' Ramp rate for (time) ramped boundaries Public nD As Integer 'number of dimensions Public nspecies As Integer 'number of chemical species Private nPoints As Long ' total number of points Private Cylinder As Integer ' if 1, then use cylinder symmetry Private cellen() As String Private ntimesteps As Long Private tt As Long 'actual timestep Private nHeaderRows As Long Private Rgas As Single Private Boltzman_4 ' fourth root of Boltzman Private rho_min As Single 'minimum density to prevent cavitation Public DispDim As Integer, WhatDisp As Integer, PleaseStop As Integer, DoInit As Integer Public DispZoom() As Single Public DispOffset() As Single Private nx As Long, layersize As Long Private specie As Integer '***Dynamic Arrays are commented on in DoRedim Sub '***Geometry Private nndm() As Integer 'number of points in dimension i. Private dd() As Single Private dh() As Single 'Extra dimension for cylider symmetry etc Private xx() As Single 'coordinates of cell Private Vol() As Single Private dNghbr() As Long '***State variables (Extensity) Private mv() As Single Private m() As Single Private Q() As Single Private mtotal() As Single ' sum of masses of species 'Intensity variables Private v() As Single Private v0() As Single Private p() As Single Private T() As Single Private Conc() As Single 'Parameters Input code Private R() As Single Private F() As Single Private g() As Single Private WP() As Integer Private WT() As Integer Private WV() As Integer Private WC() As Integer Private Lambda() As Single Private rho() As Single Private Cp() As Single Private eta() As Single 'Internal variables 'Private am() As Single Private ac() As Single Private aT() As Single Private av() As Single Private ap() As Single Public GlobalPar(1000) As Single Private Viewfactor() As Single Private Viewconnect() As Long Private Nviews As Long Private DoRadiation As Integer Public DemoMode As Integer Public MaxScale As Single Public Minscale As Single Public SliceNrs() As Integer Public WhichSlice(1) As Integer Dim foutje As Long 'to trap errors Dim act As Integer Sub DoReDim() 'Redim dynamic arrays 'Geometry 'ReDim nndm(nD) As Integer 'number of points in dimension i ReDim dd(1 To nD, nPoints) As Single ' distance between 2 points in dimension i ReDim xx(1 To nD, nPoints) As Single ' coordinates of cells (lower corner) ReDim dNghbr(1 To nD + 1) As Long 'indexdistance between 2 neighbours ReDim Vol(nPoints) As Single 'Volume of a cell ReDim dh(nPoints) '***State variables (Extensity) ReDim mv(1 To nD, nPoints, 5) As Single ReDim m(nPoints, nspecies, 5) As Single ReDim Q(nPoints, 5) As Single ReDim mtotal(nPoints) 'Intensity variables ReDim v(1 To nD, nPoints) ReDim v0(1 To nD, nPoints) ReDim p(nPoints) As Single ReDim T(nPoints) As Single ReDim Conc(nspecies, nPoints) As Single 'Parameters Input code ReDim R(nPoints) As Single 'R 'ReDim F(1 To nD, nPoints) As Single 'Force ReDim g(nD) As Single 'gravity ReDim WV(nD, nPoints) As Integer 'inhibit flow and mass transport in direciton i ' =0: No inhibition (default) '>=1: inhibited transport (=impregnable wall) '<=-1: Convection, but constant velocity (intensity variable held constant) '=-2: Constant velocity, ramping up over TR timesteps ReDim WP(nPoints) As Integer 'Constant pressure ' =0: No inhibition (default) '>=1: inhibited transport '<=-1: Constant pressure (intensity variable held constant) ReDim WT(nPoints) As Integer 'inhibit change in energy ' =0: No inhibition (default) '>=1: inhibited transport (=perfect insulation) '<=-1: Constant tempature (intensity variable held constant) ReDim WC(nPoints) As Integer 'chemical boundary conditions ' =0: No inhibition (default) '>=1: inhibited transport (=perfect insulation) '<=-1: Constant concentration (intensity variable held constant) ReDim Lambda(nPoints) As Single 'la ReDim rho(nPoints) As Single 'rho ReDim Cp(nPoints) As Single 'Cp ReDim eta(nPoints) As Single 'et 'Internal variables 'ReDim am(nD, nPoints) As Single ReDim ac(1 To nD, nPoints) As Single ReDim ap(1 To nD, nPoints) As Single ReDim aT(1 To nD, nPoints) As Single ReDim av(1 To nD, nPoints) As Single ReDim SliceNrs(nD) As Integer ReDim DispZoom(nD) As Single ReDim DispOffset(nD) As Single End Sub Public Sub GetinputRange() Dim masage As String Dim nr As Long Dim a As String Dim aa As String Dim lasttab As Long Dim ntab As Long Dim irow As Long Dim icol As Long Dim nlayers As Long Dim HyperSliceSize As Long Dim cell As String Dim i As Integer, Dimension As Integer Dim SliceNumber As Integer Dim endofSlicePack As Integer Dim waarpunt As Integer Dim filenaam As String 'padnaam = App.Path padnaam = CF_in.Dir1.Path infilenaam = CF_in.File1.FileName waarpunt = InStr(1, infilenaam, ".") If waarpunt > 0 Then mijnfilenaam = Mid$(infilenaam, 1, waarpunt - 1) Else mijnfilenaam = infilenaam End If outfilenaam = mijnfilenaam & "_out" GraphTitle = mijnfilenaam Open padnaam & "\" & infilenaam For Input As #1 nr = 0 nx = 0 Do Line Input #1, a If nr = 0 Then aa = a nr = nr + 1 CF_in.Text1.Text = nr Loop While Not EOF(1) Close #1 lasttab = 0 Do ntab = InStr(lasttab + 1, aa, Chr$(9)) lasttab = ntab nx = nx + 1 DoEvents Loop Until ntab = 0 CF_in.Text1.Text = " rows: " & nr & Chr$(13) & Chr$(10) & " columns: " & nx ReDim cellen(nr, nx) Open padnaam & "\" & infilenaam For Input As #1 For irow = 1 To nr Line Input #1, aa lasttab = 0 For icol = 1 To nx ntab = InStr(lasttab + 1, aa, Chr$(9)) If icol = nx Then cellen(irow, icol) = Mid$(aa, lasttab + 1, Len(aa) - lasttab) Else cellen(irow, icol) = Mid$(aa, lasttab + 1, ntab - lasttab - 1) End If lasttab = ntab Next icol Next irow Close #1 nspecies = 0 nHeaderRows = 0 nlayers = -1 For irow = 1 To nr cell = cellen(irow, 1) cell = UCase$(cell) Call findval(cell, "NS", nspecies) ' we need nspecies for redim If Mid$(cell, 1, 1) = "L" Then nlayers = nlayers + 1 If nHeaderRows = 0 Then nHeaderRows = irow End If Next irow If nlayers = 0 Then MsgBox ("You must define at least one layer by specifying a letter 'L' above and below the layer geometry") Else If Mid$(cell, 1, 1) <> "L" Then MsgBox ("You must terminate a set of rows with an 'L' to specify end of layer.") Else 'Sort out the maximum number of dimensions Call findval(cell, "L", nD) If nD < 3 Then nD = 3 If nD = 3 Then If nlayers = 1 Then nD = 2 'Only one slice -> 2D If nD = 2 Then 'Only one row -> 1D If nr - 2 = 1 Then nD = 1 End If End If End If ReDim nndm(nD) As Integer '*** check layersizes layersize = Int((nr - nHeaderRows) / nlayers) If layersize <= 1 Then MsgBox ("Layersize <= 1!") Else HyperSliceSize = layersize For i = 3 To nD nndm(i) = 0 SliceNumber = 1 endofSlicePack = 0 Do irow = (SliceNumber) * HyperSliceSize + nHeaderRows If irow > nr Then endofSlicePack = 1 Else If irow = nr Then endofSlicePack = 1 cell = cellen(irow, 1) ' cell = cellen(irow + 1, 1) cell = UCase$(cell) If Len(cell) = 0 Then cell = "Error" If Mid$(cell, 1, 1) <> "L" Then masage = "row :" & irow & Chr$(13) & Chr$(10) masage = masage & "The layers are not of equal size." & Chr$(13) & Chr$(10) masage = masage & "The 'L' symbol of layer " & Dimension & " was not found at the correct spot." masage = masage & Chr$(13) & Chr$(10) masage = masage & "-Remember the row with defaults at the top should also be selected" masage = masage & Chr$(13) & Chr$(10) masage = masage & "-Remember the last row with 'L' should also be selected" MsgBox (masage) ' End End If Call findval(cell, "L", Dimension) Dimension = Int(Dimension) If Dimension <= 0 Then Dimension = 3 If Dimension = i Then nndm(i) = nndm(i) + 1 SliceNumber = SliceNumber + 1 Else If Dimension = i + 1 Then nndm(i) = nndm(i) + 1 HyperSliceSize = HyperSliceSize * nndm(i) endofSlicePack = 1 Else masage = "Row :" & irow & Chr$(13) & Chr$(10) masage = "A slice nr of " & i & " or" & i + 1 & Chr$(13) & Chr$(10) masage = masage & " was not found at the correct spot." MsgBox (masage) End If End If End If Loop While endofSlicePack = 0 Next i nndm(1) = nx If nD >= 2 Then nndm(2) = layersize - 1 'Check for dimension last dimension of 1 layer (make nD-1) If nD > 1 Then If nndm(nD) = 1 Then nD = nD - 1 End If End If 'For i = 3 To nD ' nndm(i) = nndm(i) + 1 'Next i nPoints = 1 For i = 1 To nD nPoints = nPoints * nndm(i) Next i Call GiveInfo DoInit = 2 DoReDim End If End If End If End Sub Sub init() Dim i As Integer, ii As Integer, j As Integer Dim irow As Long Dim icol As Long Dim cell As String Dim icell As String Dim theEnd As Integer Dim lastdelim As Long Dim nextdelim As Long Dim nntimesteps As Long Dim xxx As Single Dim n As Long, nn As Long, iii As Long, x As Long, layer As Long, icolumn As Long Dim Nghbr As Long Dim pt As Integer Dim lLambda As Single Dim conctotal As Single Dim filenaam As String If nD = 0 Then MsgBox ("select an input geometry first") Else If WhatDisp < 1 Then WhatDisp = 1 'If DispDim < 1 Then DispDim = 1 'If DispDim > nD Then DispDim = nD Call resetstop 'neighbour structure dNghbr(1) = 1 For i = 1 To nD dNghbr(i + 1) = nndm(i) * dNghbr(i) Next i For i = 1 To nD g(i) = 0 'default no gravity Next i '******fill parameter arrays******************* '****pre-defaults (before reading header) rho(0) = 1# rho_min = 0.01 'to prevent cavitation Cp(0) = 1000# T(0) = 273.15 + 25# p(0) = 101300# Conc(0, 0) = 1 For specie = 1 To nspecies Conc(specie, 0) = 0 Next specie Lambda(0) = 0.02 eta(0) = 0.00001 R(0) = 0# dh(0) = 1# For i = 1 To nD WV(i, 0) = 0 dd(i, 0) = 0.01 Next i WT(0) = 0 ntimesteps = -1 dt = 0.01 ut = 25 logt = 0 st = 10000 ft = 10000000 'st = 10 Cylinder = 0 tRamp = 1 If DoInit <> -1 Then '*******************************ONLY RECALC******************** '***initial state is determined by p,T and Volume** MaxScale = 1000 Minscale = 600 For i = 1 To nD DispZoom(i) = 1 DispOffset(i) = 0 SliceNrs(i) = 0 Next i WhichSlice(0) = 1 WhichSlice(1) = 2 If nD >= 1 Then SliceNrs(1) = -2 If nD >= 2 Then SliceNrs(2) = -1 End If '***read header rows**** For irow = 0 To nHeaderRows - 1 For icol = 0 To nx - 1 ' cell = Inputrange(1).Offset(iy, ix) cell = cellen(irow + 1, icol + 1) cell = UCase$(cell) If Len(cell) > 0 Then theEnd = 0 lastdelim = 1 Do nextdelim = InStr(lastdelim, cell, ":") If nextdelim = 0 Then theEnd = 1 icell = Mid$(cell, lastdelim, Len(cell) - lastdelim + 1) Else icell = Mid$(cell, lastdelim, nextdelim - lastdelim) lastdelim = nextdelim + 1 End If If Len(icell) > 0 Then Call findval(icell, "DT", dt) Call findval(icell, "UT", ut) Call findval(icell, "ST", st) Call findval(icell, "FT", ft) Call findval(icell, "LT", logt) Call findval(icell, "CY", Cylinder) Call findval(icell, "TR", tRamp) ' Call findval(icell, "NS", nspecies) Call findval(icell, "NT", nntimesteps) If DoInit <> -1 Then 'display options Call findval(icell, "SA", MaxScale) Call findval(icell, "SI", Minscale) Call findval(icell, "SW", WhatDisp) Call findval(icell, "SD", DispDim) End If Call readvalues(icell, 0) i = -1: Call find2val(icell, "GG", i, xxx) If i > 0 Then g(i) = xxx ' gravity i = -1: Call find2val(icell, "GP", i, xxx) If i > 0 Then GlobalPar(i) = xxx End If xxx = -1: Call find2val(icell, "DD", i, xxx) If xxx > -1 Then If i = 0 Then For ii = 1 To nD dd(ii, 0) = xxx Next ii Else If Abs(i) <= nD Then dd(Abs(i), 0) = xxx End If End If End If Loop While theEnd = 0 End If Next icol Next irow Rgas = p(0) / (rho(0) * T(0)) '****Overwrite pre-defaults with defaults from header For n = 1 To nPoints rho(n) = rho(0) For specie = 0 To nspecies Conc(specie, n) = Conc(specie, 0) Next specie Cp(n) = Cp(0) Lambda(n) = Lambda(0) eta(n) = eta(0) T(n) = T(0) p(n) = p(0) For i = 1 To nD v(i, n) = 0 Next i R(n) = R(0) dh(n) = dh(0) For i = 1 To nD WV(i, n) = WV(i, 0) If getcoori(n, i) = nndm(i) - 1 Then WV(i, n) = 1 End If ' extra wall at boundary dd(i, n) = dd(i, 0) Next i WT(n) = WT(0) Next n '*** Now look for parameters in non-header cells. For n = 1 To nPoints nn = n - 1 layer = Int(nn / nx / (layersize - 1)) irow = Int(nn / nx) - layer * (layersize - 1) icolumn = (nn - layer * nx * (layersize - 1) - irow * nx) irow = layer * (layersize) + irow + nHeaderRows ' cell = Inputrange(1).Offset(irow, icolumn) cell = cellen(irow + 1, icolumn + 1) cell = UCase$(cell) If Len(cell) > 0 Then theEnd = 0 lastdelim = 1 Do nextdelim = InStr(lastdelim, cell, ":") If nextdelim = 0 Then theEnd = 1 icell = Mid$(cell, lastdelim, Len(cell) - lastdelim + 1) Else icell = Mid$(cell, lastdelim, nextdelim - lastdelim) lastdelim = nextdelim + 1 End If If Len(icell) > 0 Then Call readvalues(icell, n) xxx = -1: Call find2val(icell, "DD", i, xxx) If xxx > -1 Then If i = 0 Then For ii = 1 To nD dd(ii, n) = xxx Next ii Else If Abs(i) <= nD Then If i > 0 Then dd(Abs(i), n) = xxx Else 'if i = negative, then the entire hyperplane of i=const gets the same space interval ii = getcoori(n, Abs(i)) For nn = n To nPoints '***does NOT overrule previous points iii = getcoori(nn, Abs(i)) If iii = ii Then dd(Abs(i), nn) = xxx End If Next nn End If End If End If End If xxx = -1: Call find2val(icell, "HH", i, xxx) If xxx > -1 Then If Abs(i) <= nD Then If i < 1 Then dh(n) = xxx Else 'i the entire hyperplane of i=const gets the same height ii = getcoori(n, i) For nn = 1 To nPoints '***overrules previous points also iii = getcoori(nn, i) If iii = ii Then dh(nn) = xxx End If Next nn End If End If End If End If Loop While theEnd = 0 End If Next n 'If DoInit <> -1 Then ReDim m(nPoints, nspecies, 5) As Single For n = 1 To nPoints Vol(n) = dh(n) For i = 1 To nD Vol(n) = Vol(n) * dd(i, n) Next i conctotal = 0 For specie = 1 To nspecies conctotal = conctotal + Conc(specie, n) Next specie Conc(0, n) = 1 - conctotal ' renormalize concentrations '****************************************************ONLY RECALC******************** If DoInit <> -1 Then '*******************************ONLY RECALC******************** '***initial state is determined by p,T and Volume** tt = 0 rho(n) = p(n) / (Rgas * T(n)) mtotal(n) = Vol(n) * rho(n) For specie = 0 To nspecies m(n, specie, 0) = mtotal(n) * Conc(specie, n) Next specie Q(n, 0) = mtotal(n) * T(n) * Cp(n) For i = 1 To nD v(i, n) = v0(i, n) If WV(i, n) = -1 Then v(i, n) = v0(i, n) * rho(0) / rho(n) 'standardized velocity is specified End If If WV(i, n) = -2 Then v(i, n) = (1 / tRamp) * v0(i, n) * rho(0) / rho(n) 'standardized velocity is specified End If mv(i, n, 0) = mtotal(n) * v(i, n) Next i For specie = 0 To nspecies m(n, specie, 5) = m(n, specie, 0) Next specie Q(n, 5) = Q(n, 0) For i = 1 To nD mv(i, n, 5) = mv(i, n, 0) Next i For pt = 1 To 4 For specie = 0 To nspecies m(n, specie, pt) = 0 Next specie Q(n, pt) = 0 For i = 1 To nD mv(i, n, pt) = 0 Next i Next pt End If Next n '****************************************************(ONLY RECALC)******************** 'extra walls at far edges of the grid, and reinterpret wall=2 'wall=2 means the neighbour on the decementing side is also a wall For n = 1 To nPoints If WV(0, n) <> 0 Then For i = 1 To nD WV(i, n) = WV(0, n) Next i End If For i = 1 To nD If WV(i, n) = 2 Then If Int((n - 1) / dNghbr(i + 1)) = Int((n - dNghbr(i) - 1) / dNghbr(i + 1)) Then WV(i, n - dNghbr(i)) = 1 End If WV(i, n) = 1 End If ' If getcoori(n, i) = nndm(i) - 1 Then ' WV(i, n) = 1 ' End If Next i Next n 'transfercoefficients, coordinates For n = 1 To nPoints R(n) = R(n) * Vol(n) For i = 1 To nD If getcoori(n, i) < 1 Then xx(i, n) = 0 ' origin coordinates are zero Else xx(i, n) = xx(i, n - dNghbr(i)) + dd(i, n - dNghbr(i)) End If Nghbr = n + dNghbr(i) ' F(i, n) = F(i, n) * Vol(n) / dd(i, n) ' am(i, n) = Vol(n) / (dd(i, n) ^ 2) 'mass transfer coeff ac(i, n) = 1 / (dd(i, n)) 'convection transfer coeff (convert momentum to mass flux) ap(i, n) = Vol(n) / (dd(i, n)) 'pressure difference to dmv/dt transfer coeff lLambda = Lambda(n) ' If Int((n - 1) / dNghbr(i + 1)) = Int((Nghbr - 1) / dNghbr(i + 1)) Then ' lLambda = (1 / Lambda(n) + 1 / Lambda(Nghbr)) / 2 ' lLambda = 1 / lLambda ' End If aT(i, n) = lLambda * Vol(n) / (dd(i, n) ^ 2) 'thermal conduction transfer coeff For j = 1 To nD av(i, n) = eta(n) * Vol(n) / (dd(i, n) ^ 2) 'viscosity transfer coeff Next j Next i Next n ntimesteps = nntimesteps + tt If DoInit = 2 Then Close #3 If logt > 0 Then filenaam = mijnfilenaam & "_log" Open padnaam & "\" & filenaam For Output As #3 Close #3 End If End If DoInit = 1 ''ReDim dd(0) 'ReDim Lambda(0) 'ReDim eta(0) 'ReDim cellen(0) End If End Sub 'Sub findval(where As String, what As String, ByRef waarde As single) Sub findval(where As String, What As String, ByRef theValue) 'Find a numeric value in a string after a codestring 'first, replace commas with points Dim theEnd As Integer Dim pp As Long Dim codepos As Long theEnd = 0 pp = 1 Do pp = InStr(pp, where, ",") If pp = 0 Then theEnd = 1 Else Mid$(where, pp, 1) = "." pp = pp + 1 End If Loop While theEnd = 0 codepos = InStr(1, where, What) If codepos > 0 Then theValue = Val(Mid$(where, codepos + Len(What), Len(where) - codepos - Len(What) + 1)) End If End Sub Sub find2val(where As String, What As String, ByRef theValue1, ByRef theValue2) 'Find 2 numeric values in a string after a codestring. The 2 values are separated by semicolon 'first, replace commas with points Dim theEnd As Integer Dim pp As Long, codepos1 As Long, codepos2 As Long theEnd = 0 pp = 1 Do pp = InStr(pp, where, ",") If pp = 0 Then theEnd = 1 Else Mid$(where, pp, 1) = "." pp = pp + 1 End If Loop While theEnd = 0 codepos1 = InStr(1, where, What) If codepos1 > 0 Then codepos2 = InStr(codepos1 + Len(What), where, ";") ' waarde1 = Val(Mid$(where, codepos1 + Len(what), Len(where) - codepos - Len(what) + 1)) If codepos2 > 0 Then theValue1 = Val(Mid$(where, codepos1 + Len(What), codepos2 - codepos1 - Len(What) + 1)) theValue2 = Val(Mid$(where, codepos2 + 1, Len(where) - codepos2 - 1 + 1)) End If End If End Sub Sub findIndexed(where As String, What As String, MijnArray, Index) 'Find 2 numeric values (value1 and value2) in a string after a codestring. 'The 2 values are separated by semicolon 'the value1 is the first index of the array 'the value2 is the value that the array element gets 'the second index of the array is the one passed as 'Index' ' In other words: MijnArray(value1, Index) = value2 'first, replace commas with points Dim theEnd As Integer Dim foutje As Integer Dim pp As Long Dim codepos1 As Long Dim codepos2 As Long Dim value1 As Single Dim value2 As Single Dim i As Integer theEnd = 0 pp = 1 Do pp = InStr(pp, where, ",") If pp = 0 Then theEnd = 1 Else Mid$(where, pp, 1) = "." pp = pp + 1 End If Loop While theEnd = 0 pp = 1 For i = 0 To nD 'There could be more than one dimension specified in the string! codepos1 = InStr(pp, where, What) If codepos1 > 0 Then pp = codepos1 + 2 codepos2 = InStr(codepos1 + Len(What), where, ";") ' waarde1 = Val(Mid$(where, codepos1 + Len(what), Len(where) - codepos - Len(what) + 1)) If codepos2 > 0 Then value1 = Val(Mid$(where, codepos1 + Len(What), codepos2 - codepos1 - Len(What) + 1)) value2 = Val(Mid$(where, codepos2 + 1, Len(where) - codepos2 - 1 + 1)) foutje = 0 If value1 < LBound(MijnArray, 1) Then foutje = 1 If value1 > UBound(MijnArray, 1) Then foutje = 1 If Index < LBound(MijnArray, 2) Then foutje = 1 If Index > UBound(MijnArray, 2) Then foutje = 1 If foutje = 0 Then MijnArray(value1, Index) = value2 Else ' MsgBox ("Error reading :" & where) ' End End If Else ' MsgBox ("Error reading :" & where) ' End End If End If Next i End Sub Public Sub DoRestart() If nD = 0 Then MsgBox ("load an input geometry first") Else DoInit = 1 Call DoCellFLow End If End Sub Public Sub DoContinue() If nD = 0 Then MsgBox ("load an input geometry first") Else DoInit = -1 Call DoCellFLow End If End Sub Sub DoCellFLow() Dim pt As Integer, old As Integer, HasNghbr As Integer Dim n As Long Dim Nghbr As Long Dim bNghbr As Long Dim CNghbr As Long Dim iv As Long Dim n1 As Long, n2 As Long Dim b As Integer ' neigbour direction Dim i As Integer Dim j As Integer Dim ptf As Single Dim drctn As Single Dim delta_m As Single Dim T1 As Single, T2 As Single, Q12 As Single Dim Keq As Single Dim c_eq As Single Dim c_convert As Single Dim c_tau As Single Dim filenaam As String '*** '**** iterations init If nD = 0 Then MsgBox ("No geometry") Else old = 0 act = 5 PleaseStop = 0 Do tt = tt + 1 For pt = 1 To 4 If pt = 4 Then ptf = 1 * dt Else ptf = 0.5 * dt End If 'factors in Runge Kutta method DoEvents If pt = 1 Then 'create new state vector from previous data For n = 1 To nPoints If WP(n) = 0 Then For specie = 0 To nspecies m(n, specie, act) = m(n, specie, old) + dt * (m(n, specie, 1) + 2 * m(n, specie, 2) + 2 * m(n, specie, 3) + m(n, specie, 4)) / 6 Next specie End If ' If WT(n) = 0 And WP(n) = 0 Then If WT(n) = 0 Then Q(n, act) = Q(n, old) + dt * (Q(n, 1) + 2 * Q(n, 2) + 2 * Q(n, 3) + Q(n, 4)) / 6 End If For j = 1 To nD If WV(j, n) = 0 Then mv(j, n, act) = mv(j, n, old) + dt * (mv(j, n, 1) + 2 * mv(j, n, 2) + 2 * mv(j, n, 3) + mv(j, n, 4)) / 6 End If Next j For specie = 0 To nspecies m(n, specie, old) = m(n, specie, act) Next specie Q(n, old) = Q(n, act) For j = 1 To nD mv(j, n, old) = mv(j, n, act) Next j Next n Else 'construct the temporary state (act) for this pseudo timestep For n = 1 To nPoints If WP(n) = 0 Then For specie = 0 To nspecies m(n, specie, act) = m(n, specie, old) + ptf * m(n, specie, pt - 1) Next specie End If ' If WT(n) = 0 And WP(n) = 0 Then If WT(n) = 0 Then Q(n, act) = Q(n, old) + ptf * Q(n, pt - 1) End If For j = 1 To nD If WV(j, n) = 0 Then mv(j, n, act) = mv(j, n, old) + ptf * mv(j, n, pt - 1) End If Next j Next n End If Call CalcIntensity For n = 1 To nPoints For specie = 0 To nspecies m(n, specie, pt) = 0 Next specie Q(n, pt) = 0 For j = 1 To nD mv(j, n, pt) = 0 Next j For b = 1 To nD * 2 '***neighbour related calcs ' ********* c is transported quantity (C) per unit mass, ' multiply this by mass flux accross surface. ' The mass flux is momentum * transfercoeff(ac) ' So the coeff ac converts momentum to mass flux. ' .. -- ' : --->| left neigbour mflux>0: transfercoeff&mflux from left, c from left ' .. -- '***bnghbr = x-1, cnghbr = x-1, drctn > 0 ' .. -- ' : <-- | left neigbour mflux<0: transfercoeff&mflux from left, c from middle ' .. -- '***bnghbr = x-1, cnghbr = x, drctn < 0 ' ' -- .. ' | --> : right neigbour mflux>0: transfercoeff&mflux from middle, c from middle ' -- .. '***bnghbr = x, cnghbr = x, drctn < 0 ' ' -- .. ' | <-- : right neigbour mflux<0: transfercoeff&mflux from middle, c from right ' -- .. '***bNghbr = x, cNghbr = x+1, drctn > 0 i = 1 + Int((b - 1) / 2) 'dimension index If (b Mod 2) = 1 Then Nghbr = n - dNghbr(i) ' index of neighbour bNghbr = Nghbr 'index of edge (mass flux and transfer parameter) CNghbr = n 'for later use If Int((n - 1) / dNghbr(i + 1)) = Int((Nghbr - 1) / dNghbr(i + 1)) Then '**check if we are on the boundary (-> no neighbor on that side HasNghbr = 1 drctn = 1 * ac(i, bNghbr) Else ' drctn = 1 * ac(i, n) 'at edge, transport a copy of facing inner edge HasNghbr = 0 ' drctn = 0 'at edge, no transport End If Else Nghbr = n + dNghbr(i) ' index of neighbour bNghbr = n 'index of edge (mass flux and transfer parameter) CNghbr = Nghbr 'for later use If Int((n - 1) / dNghbr(i + 1)) = Int((Nghbr - 1) / dNghbr(i + 1)) Then '**check if we are on the boundary (-> no neighbor on that side HasNghbr = 1 drctn = -1 * ac(i, bNghbr) Else ' drctn = -1 * ac(i, n) 'at edge, transport a copy of facing inner edge HasNghbr = 0 ' drctn = 0 'at edge, no transport End If ' drctn = -1 * ac(i, bNghbr) End If '***index of convective source If HasNghbr = 1 Then If v(i, bNghbr) > 0 Then CNghbr = bNghbr Else CNghbr = CNghbr If (WV(i, bNghbr) <= 0) And (rho(CNghbr) > rho_min) Then 'convection For specie = 0 To nspecies m(n, specie, pt) = m(n, specie, pt) + drctn * (mv(i, bNghbr, act)) * Conc(specie, CNghbr) Next specie If WT(n) <= 0 Then Q(n, pt) = Q(n, pt) + drctn * (mv(i, bNghbr, act)) * Cp(CNghbr) * T(CNghbr) For j = 1 To nD If WV(j, n) <= 0 Then mv(j, n, pt) = mv(j, n, pt) + drctn * (mv(i, bNghbr, act)) * v(j, CNghbr) Next j End If If (b Mod 2) = 0 Then If WV(i, n) <= 0 Then mv(i, n, pt) = mv(i, n, pt) + (p(n) - p(Nghbr)) * ap(i, n) ' If WT(BNghbr) <= 0 Then Q(n, pt) = Q(n, pt) + p(Nghbr) * v(i, BNghbr) * ap(i, BNghbr) 'labour done by pressure on volume element End If 'thermal conduction If WT(n) <= 0 Then Q(n, pt) = Q(n, pt) + (T(Nghbr) - T(n)) * aT(i, bNghbr) 'viscosity For j = 1 To nD If WV(j, n) <= 0 Then mv(j, n, pt) = mv(j, n, pt) + (v(j, Nghbr) - v(j, n)) * av(i, bNghbr) Next j End If Next b If WV(0, n) = 0 Then For j = 1 To nD If WV(j, n) <= 0 Then mv(j, n, pt) = mv(j, n, pt) - R(n) * v(j, n) ' mv(j, n, pt) = mv(j, n, pt) + F(j, n) mv(j, n, pt) = mv(j, n, pt) + g(j) * mtotal(n) Next j End If If WP(n) < 0 Then 'corrections for mass conservation violation at constant pressure ' delta_m = (m(n, act) - m(n, old)) / ptf delta_m = 0 For specie = 0 To nspecies delta_m = delta_m - m(n, specie, pt) Next specie Q(n, pt) = Q(n, pt) + T(n) * delta_m * Cp(n) 'recalibrate Q to match change in mass For j = 1 To nD mv(j, n, pt) = mv(j, n, pt) + v(j, n) * delta_m 'recalibrate mv to match change in mass Next j For specie = 0 To nspecies m(n, specie, pt) = m(n, specie, pt) + Conc(specie, n) * delta_m Next specie End If If WC(n) = 1 Then ' chemical reaction 1 ' '***8 ata ' Keq = 0.0000000027 * Exp(18000# / T(n)) ' c_eq = 0.22 * Keq / (1 + Keq) '' '***1 ata '' Keq = 0.000000000016011 * Exp(20000# / T(n)) '' c_eq = 0.22 * Keq / (1 + Keq) ' ' c_tau = (1 / 0.05) * Exp(-8000# / T(n)) / Exp(-8000# / (273 + 700)) '***8 ata Keq = GlobalPar(11) * Exp(GlobalPar(12) / T(n)) c_eq = GlobalPar(13) * Keq / (1 + Keq) ' '***1 ata ' Keq = 0.000000000016011 * Exp(20000# / T(n)) ' c_eq = 0.22 * Keq / (1 + Keq) c_tau = (1 / GlobalPar(15)) * Exp(-GlobalPar(14) / T(n)) / Exp(-GlobalPar(14) / (273 + 700)) If c_tau > 0.5 / dt Then c_tau = 0.5 / dt ' ghsv = 5*3600-> tau = 0.2 sec c_convert = mtotal(n) * (Conc(1, n) - c_eq) * c_tau m(n, 1, pt) = m(n, 1, pt) - c_convert m(n, 0, pt) = m(n, 0, pt) + c_convert Q(n, pt) = Q(n, pt) - c_convert * 208 * 1000 * 1000 / (12 + 4) ' 208 kJ/mol is approximate value at 900 degrees, 50% CO conversion End If Next n 'radiation If DoRadiation = 1 Then For iv = 1 To Nviews n1 = Viewconnect(iv, 0) n2 = Abs(Viewconnect(iv, 1)) T1 = T(n1) * Boltzman_4 T2 = T(n2) * Boltzman_4 Q12 = Viewfactor(iv) * (T1 ^ 4 - T2 ^ 4) If WT(n1) <= 0 Then Q(n1, pt) = Q(n1, pt) - Q12 If WT(n2) <= 0 And n2 > 0 Then Q(n2, pt) = Q(n2, pt) + Q12 Next iv End If Next pt If (tt >= ntimesteps And ntimesteps >= 0) Then PleaseStop = 1 If (tt - 1) Mod st = st - 1 Then Call SaveState If (tt - 1) Mod ft = ft - 1 Then Call maakplaatje(0) If (tt - 1) Mod ut = 0 Or PleaseStop = 1 Then CF_in.Text1.Text = tt Call Dodisplay(WhatDisp, DispDim) DoEvents End If If logt > 0 Then If ((tt - 1) Mod logt = 0 Or PleaseStop = 1) Then filenaam = mijnfilenaam & "_log" Open padnaam & "\" & filenaam For Append As #3 'take point (3,16,16) n = 1 n = n + 4 * dNghbr(1) n = n + 16 * dNghbr(2) n = n + 16 * dNghbr(3) Print #3, v(1, n); Chr$(9); p(n) Close #3 End If End If DoEvents Call checkforstop Loop Until PleaseStop = 1 End If End Sub Sub checkforstop() Dim x As Integer Open "checkstop" For Input As #2 Input #2, x Close #2 If x = 1 Then Call exitorpause PleaseStop = 1 Call resetstop End If End Sub Sub resetstop() Open "checkstop" For Output As #2 Print #2, 0 Close #2 End Sub Sub CalcIntensity() Dim n As Long, j As Long Dim ramp As Single For n = 1 To nPoints ' ***construct intensity variables, '***and handle exensity variables with constant intensity conditions mtotal(n) = 0 For specie = 0 To nspecies mtotal(n) = mtotal(n) + m(n, specie, act) Next specie rho(n) = mtotal(n) / Vol(n) If WC(n) < 0 Then 'constant concentration For specie = 0 To nspecies m(n, specie, act) = mtotal(n) * Conc(specie, n) Next specie Else For specie = 0 To nspecies Conc(specie, n) = m(n, specie, act) / mtotal(n) Next specie End If If rho(n) < 0 Then MsgBox ("Stopped because of negative density in point " & n) PleaseStop = 1 End If If WT(n) < 0 Then 'constant temperature Q(n, act) = T(n) * (mtotal(n) * Cp(n)) Else T(n) = Q(n, act) / (mtotal(n) * Cp(n)) End If ' If WV(0, n) <= 0 Then For j = 1 To nD If WV(j, n) = 0 Then v(j, n) = mv(j, n, act) / mtotal(n) 'normal: calc intensity End If If WV(j, n) < 0 Then If WV(j, n) = -2 Then ramp = (tt / tRamp) If ramp > 1 Then ramp = 1 v(j, n) = ramp * v0(j, n) * rho(0) / rho(n) 'fixed standardized velocity Else v(j, n) = v0(j, n) * rho(0) / rho(n) 'fixed standardized velocity End If mv(j, n, act) = mtotal(n) * v(j, n) End If Next j ' End If If WP(n) < 0 Then 'constant pressure ' handled separately Else p(n) = Rgas * mtotal(n) * T(n) / Vol(n) End If Next n End Sub Sub Dodisplay(ByVal What As Integer, ByVal iDim As Integer) Dim n As Long, i As Integer Dim row As Long, col As Long, irows As Long, icols As Long Dim Nghbr As Long Dim waarde As Single, wwaarde As Single Dim maxwaarde As Single, minwaarde As Single Dim xx1 As Single, xx2 As Single, yy1 As Single, yy2 As Single Dim red As Integer, green As Integer, blue As Integer Dim Lgndwidth As Integer, lkey As Integer Dim xgeom As Single, ygeom As Single, xscreen As Single, yscreen As Single Dim gwscale As Single, Margin As Single, TitleHeight As Single Dim xxx1 As Single, xxx2 As Single, yyy1 As Single, yyy2 As Single Dim vvx As Single, vvy As Single, vvv As Single ReDim coordinate(nD) As Long If rho(0) > 0 Then ' check for empty arrays For i = 1 To nD coordinate(i) = SliceNrs(i) Next i icols = WhichSlice(0) irows = WhichSlice(1) Lgndwidth = 120 Margin = 5 TitleHeight = 25 'CF_out.Picture1.ForeColor = &HFFFFC0 CF_out.Picture1.ForeColor = &H800000 CF_out.Picture1.Left = 0 CF_out.Picture1.Top = 0 CF_out.Picture1.Width = CF_out.ScaleWidth CF_out.Picture1.Height = CF_out.ScaleHeight If DispZoom(icols) <= 0 Then DispZoom(icols) = 1 If DispZoom(irows) <= 0 Then DispZoom(irows) = 1 xgeom = (xx((icols), nPoints)) * DispZoom(icols) ygeom = (xx((irows), nPoints)) * DispZoom(irows) xscreen = (CF_out.Picture1.ScaleWidth - Lgndwidth - Margin) yscreen = CF_out.Picture1.ScaleHeight - 3 * Margin - TitleHeight gwscale = xscreen / xgeom If gwscale > yscreen / ygeom Then gwscale = yscreen / ygeom xx1 = (CF_out.Picture1.ScaleWidth - Lgndwidth) / 2 yy1 = Margin CF_out.Picture1.CurrentX = xx1 CF_out.Picture1.CurrentY = yy1 CF_out.Picture1.FontSize = 12 CF_out.Picture1.FontBold = True CF_out.Picture1.Print GraphTitle xx1 = (CF_out.Picture1.ScaleWidth - Lgndwidth) + Lgndwidth * 0.1 yy1 = TitleHeight + (0.2) / (10 + 3) * CF_out.ScaleHeight CF_out.Picture1.CurrentX = xx1 CF_out.Picture1.CurrentY = yy1 CF_out.Picture1.FontSize = 10 CF_out.Picture1.FontBold = True CF_out.Picture1.Print CF_dsp.Combo1.List(WhatDisp - 1) 'If WhatDisp = 2 Or WhatDisp = 3 Or WhatDisp = 6 Or WhatDisp = 7 Or WhatDisp = 10 Then If WhatDisp = 2 Or WhatDisp = 3 Or WhatDisp = 7 Or WhatDisp = 10 Then xx1 = (CF_out.Picture1.ScaleWidth - Lgndwidth) + Lgndwidth * 0.1 yy1 = yy1 + (0.5) / (10 + 3) * CF_out.ScaleHeight CF_out.Picture1.CurrentX = xx1 CF_out.Picture1.CurrentY = yy1 If DispDim > 0 Then CF_out.Picture1.Print "(" & Mid$(" X Y Zd4", DispDim * 2 - 1, 2) & ")" End If End If CF_out.Picture1.FontSize = 8 CF_out.Picture1.FontBold = False If WhatDisp <> 11 Then For lkey = 0 To 10 waarde = lkey / 10 xx1 = (CF_out.Picture1.ScaleWidth - Lgndwidth) + Lgndwidth * 0.2 yy1 = TitleHeight + (lkey + 1.5) / (10 + 3) * CF_out.Picture1.ScaleHeight xx2 = (CF_out.Picture1.ScaleWidth - Lgndwidth) + Lgndwidth * 0.8 yy2 = yy1 + 0.5 / (10 + 3) * CF_out.Picture1.ScaleHeight wwaarde = (3 * waarde) * 255 If wwaarde < 0 Then wwaarde = 0 If wwaarde > 255 Then wwaarde = 255 red = Int(wwaarde) wwaarde = (3 * waarde - 1) * 255 If wwaarde < 0 Then wwaarde = 0 If wwaarde > 255 Then wwaarde = 255 green = Int(wwaarde) wwaarde = (3 * waarde - 2) * 255 If wwaarde < 0 Then wwaarde = 0 If wwaarde > 255 Then wwaarde = 255 blue = Int(wwaarde) CF_out.Picture1.ForeColor = RGB(red, green, blue) CF_out.Picture1.Line (xx1, yy1)-(xx2, yy2), , BF waarde = Minscale + waarde * (MaxScale - Minscale) CF_out.Picture1.CurrentX = xx1 CF_out.Picture1.CurrentY = yy2 + 0.1 / (10 + 3) * CF_out.Picture1.ScaleHeight CF_out.Picture1.Print Format(waarde) Next lkey End If Open padnaam & "\" & outfilenaam For Output As #2 If WhatDisp = 11 Then CF_out.Picture1.Cls 'For n = 1 To nPoints For row = 0 To nndm(irows) - 1 For col = 0 To nndm(icols) - 1 coordinate(icols) = col coordinate(irows) = row n = 1 For i = 1 To nD n = n + coordinate(i) * dNghbr(i) Next i If What = 1 Then waarde = T(n) - 273.15 If What = 2 Then If DispDim > nD Then DispDim = nD End If If DispDim > 0 Then waarde = v(DispDim, n) Else DispDim = 0 waarde = 0 For i = 1 To nD waarde = waarde + v(i, n) ^ 2 Next i waarde = Sqr(waarde) End If End If If What = 3 And rho(0) > 0 Then If DispDim > nD Then DispDim = nD End If If DispDim > 0 Then waarde = v(DispDim, n) * rho(n) / rho(0) Else DispDim = 0 waarde = 0 For i = 1 To nD waarde = waarde + v(i, n) ^ 2 Next i waarde = Sqr(waarde) * rho(n) / rho(0) End If End If If What = 4 Then waarde = p(n) ' - p(0) If What = 5 Then waarde = rho(n) If What = 6 Then If DispDim > nspecies Then DispDim = nspecies End If If DispDim < 0 Then DispDim = 0 End If waarde = Conc(DispDim, n) ' dimension used here for specie End If If What = 7 Then ' heat flux to INCREMENTING neighbour waarde = 0 'default If DispDim > nD Then DispDim = nD End If If DispDim > 0 Then Nghbr = n + dNghbr(DispDim) ' index of neighbour If Int((n - 1) / dNghbr(DispDim + 1)) = Int((Nghbr - 1) / dNghbr(DispDim + 1)) Then If WT(n) <= 0 Then waarde = (T(n) - T(Nghbr)) * Lambda(n) / 1000# End If Else waarde = 0 For i = 1 To nD Nghbr = n + dNghbr(i) ' index of neighbour If Int((n - 1) / dNghbr(i + 1)) = Int((Nghbr - 1) / dNghbr(i + 1)) Then waarde = waarde + ((T(n) - T(Nghbr)) * Lambda(n) / 1000#) ^ 2 End If Next i waarde = Sqr(waarde) End If End If If What = 8 Then If DispDim > nD Then DispDim = nD End If If DispDim < 1 Then DispDim = 1 End If waarde = dd(DispDim, n) End If If What = 9 Then waarde = dh(n) If What = 10 Then If DispDim > nD Then DispDim = nD End If If DispDim < 1 Then DispDim = 1 End If waarde = getcoori(n, DispDim) End If If getcoori(n, icols) = 0 Then Print #2, Else Print #2, Chr$(9); End If Print #2, waarde; xx1 = (xx(icols, n) - DispOffset(icols)) * DispZoom(icols) yy1 = (xx(irows, n) - DispOffset(irows)) * DispZoom(irows) xx2 = xx1 + dd(icols, n) * DispZoom(icols) yy2 = yy1 + dd(irows, n) * DispZoom(irows) xx1 = Int((xx1 * gwscale) + 0.5) + 1 + Margin yy1 = Int((yy1 * gwscale) + 0.5) + 1 + TitleHeight xx2 = Int((xx2 * gwscale) + 0.5) - 1 + Margin yy2 = Int((yy2 * gwscale) + 0.5) - 1 + TitleHeight If WhatDisp <> 11 Then wwaarde = (3 * waarde - 3 * Minscale) / (MaxScale - Minscale) * 255 If wwaarde < 0 Then wwaarde = 0 If wwaarde > 255 Then wwaarde = 255 red = Int(wwaarde) wwaarde = (3 * waarde - 2 * Minscale - MaxScale) / (MaxScale - Minscale) * 255 If wwaarde < 0 Then wwaarde = 0 If wwaarde > 255 Then wwaarde = 255 green = Int(wwaarde) wwaarde = (3 * waarde - Minscale - 2 * MaxScale) / (MaxScale - Minscale) * 255 If wwaarde < 0 Then wwaarde = 0 If wwaarde > 255 Then wwaarde = 255 blue = Int(wwaarde) If xx1 >= Margin Then If xx2 <= (CF_out.Picture1.ScaleWidth - Lgndwidth) Then If yy1 >= TitleHeight Then If yy2 <= (CF_out.Picture1.ScaleHeight - Margin) Then CF_out.Picture1.DrawWidth = 1 CF_out.Picture1.ForeColor = RGB(red, green, blue) CF_out.Picture1.Line (xx1, yy1)-(xx2, yy2), , BF End If End If End If End If Else If Minscale <= 0 Then Minscale = 1 If MaxScale <= 0 Then MaxScale = 1 vvx = v(icols, n) * gwscale * DispZoom(icols) / Minscale vvy = v(irows, n) * gwscale * DispZoom(irows) / Minscale vvv = Sqr(vvx ^ 2 + vvy ^ 2) If vvv > MaxScale Then vvx = vvx * MaxScale / vvv vvy = vvy * MaxScale / vvv CF_out.Picture1.DrawWidth = 2 Else CF_out.Picture1.DrawWidth = 1 End If xxx1 = (xx1 + xx2) / 2 yyy1 = (yy1 + yy2) / 2 xxx2 = xxx1 + vvx yyy2 = yyy1 + vvy If xxx1 >= Margin Then If xxx2 <= (CF_out.Picture1.ScaleWidth - Lgndwidth) Then If yyy1 >= TitleHeight Then If yyy2 <= (CF_out.Picture1.ScaleHeight - Margin) Then CF_out.Picture1.ForeColor = RGB(0, 255, 255) CF_out.Picture1.Line (xxx1, yyy1)-(xxx2, yyy2) End If End If End If End If End If If WV(icols, n) > 0 Then If xx1 >= Margin Then If xx2 <= (CF_out.Picture1.ScaleWidth - Lgndwidth) Then If yy1 >= TitleHeight Then If yy2 <= (CF_out.Picture1.ScaleHeight - Margin) Then CF_out.Picture1.DrawWidth = 2 CF_out.Picture1.ForeColor = RGB(0, 0, 255) CF_out.Picture1.Line (xx2 + 1, yy1 - 1)-(xx2 + 1, yy2 + 1) End If End If End If End If End If If WV(irows, n) > 0 Then If xx1 >= Margin Then If xx2 <= (CF_out.Picture1.ScaleWidth - Lgndwidth) Then If yy1 >= TitleHeight Then If yy2 <= (CF_out.Picture1.ScaleHeight - Margin) Then CF_out.Picture1.DrawWidth = 2 CF_out.Picture1.ForeColor = RGB(0, 0, 255) CF_out.Picture1.Line (xx1 - 1, yy2 + 1)-(xx2 + 1, yy2 + 1) End If End If End If End If End If Next col Next row CF_out.Show Close #2 End If End Sub Sub DoeSpecrum() Dim nfreq As Long Dim i As Long Dim n As Long Dim freq As Long Dim ccos() As Double Dim ssin() As Double Dim ifreq() As Long Dim xfreq() As Double Dim ffreq() As Double Dim term As Double Dim Energy As Double nfreq = InputBox("How many frequency points?", , 100) ReDim ccos(nfreq, nD) ReDim ssin(nfreq, nD) ReDim ffreq(nfreq) ReDim ifreq(nD) ReDim xfreq(nD) For freq = 0 To nfreq For i = 1 To nD ccos(freq, i) = 0 ssin(freq, i) = 0 Next i Next freq For freq = 0 To nfreq ffreq(freq) = 0 For i = 1 To nD ifreq(i) = Int(Rnd(1) * nndm(i) / 2) xfreq(i) = ifreq(i) * 2 * 4 * Atn(1) / xx(i, nPoints) ' xfreq(i) = ifreq(i) * 2 * 4 * Atn(1) / (nndm(i) - 1) ffreq(freq) = ffreq(freq) + xfreq(i) ^ 2 Next i ffreq(freq) = Sqr(ffreq(freq)) For n = 1 To nPoints term = 0 For i = 1 To nD ' term = term + xfreq(i) * getcoori(n, i) term = term + xfreq(i) * xx(i, n) Next i For i = 1 To nD ccos(freq, i) = ccos(freq, i) + v(i, n) * Cos(term) ssin(freq, i) = ssin(freq, i) + v(i, n) * Sin(term) Next i Next n Next freq On Error Resume Next Open padnaam & "\" & "CF_spectrum" For Output As #1 foutje = Err.Number On Error GoTo 0 If foutje <> 0 Then MsgBox ("An error occured opening the file") Else For freq = 0 To nfreq Energy = 0 For i = 1 To nD Energy = Energy + ccos(freq, i) ^ 2 + ssin(freq, i) ^ 2 Next i Print #1, ffreq(freq); Chr$(9); Energy Next freq Close #1 Beep MsgBox ("Spectrum Done!") End If End Sub Sub maakplaatje(gwoptie As Integer) Dim mmijnfilenaam As String mmijnfilenaam = padnaam & "\" & mijnfilenaam & Format(Now(), "_yyyymmdd_hhmmss") & ".bmp" CF_out.Picture1.Cls CF_out.Picture1.AutoRedraw = True Call Dodisplay(WhatDisp, DispDim) Call SavePicture(CF_out.Picture1.Image, mmijnfilenaam) CF_out.Picture1.Cls DoEvents CF_out.Picture1.AutoRedraw = False If gwoptie = 1 Then MsgBox (mmijnfilenaam & " created") End Sub Sub readvalues(icell As String, point As Long) ' Call findval(icell, "HH", dh(point)) Call findval(icell, "CP", Cp(point)) Call findval(icell, "ET", eta(point)) Call findval(icell, "RH", rho(point)) Call findval(icell, "LA", Lambda(point)) Call findval(icell, "TT", T(point)) Call findval(icell, "PP", p(point)) Call findval(icell, "RR", R(point)) Call findIndexed(icell, "WV", WV, point) Call findIndexed(icell, "CC", Conc, point) Call findIndexed(icell, "VV", v0, point) Call findIndexed(icell, "FF", F, point) Call findval(icell, "WT", WT(point)) Call findval(icell, "WP", WP(point)) Call findval(icell, "WC", WC(point)) End Sub Sub SvpStop() PleaseStop = 1 End Sub Sub GiveInfo() Dim masage As String Dim i As Integer masage = "You have defined the cell structure: " & Chr$(13) & Chr$(10) For i = 1 To nD masage = masage & "D=" & i & " : " & nndm(i) & " Cells" & Chr$(13) & Chr$(10) Next i If DemoMode = 0 Then MsgBox (masage) End Sub Public Sub SaveState() Dim filenaam As String Dim n As Long Dim i As Integer If nD = 0 Then MsgBox ("select an input geometery first") Else filenaam = mijnfilenaam & "_state" On Error Resume Next Open padnaam & "\" & filenaam For Output As #1 foutje = Err.Number On Error GoTo 0 If foutje <> 0 Then MsgBox ("An error occured opening the file") Else Print #1, nPoints Print #1, tt Print #1, tt 'reserved Print #1, tt 'reserved For i = 1 To nD For n = 1 To nPoints Print #1, mv(i, n, 0); Chr$(9); If getcoori(n, 1) = nndm(1) - 1 Then Print #1, Next n Print #1, Next i For n = 1 To nPoints Print #1, Q(n, 0); Chr$(9); If getcoori(n, 1) = nndm(1) - 1 Then Print #1, Next n Print #1, For specie = 0 To nspecies For n = 1 To nPoints Print #1, m(n, specie, 0); Chr$(9); If getcoori(n, 1) = nndm(1) - 1 Then Print #1, Next n Print #1, Next specie For i = 1 To nD For n = 1 To nPoints Print #1, v(i, n) * rho(n) / rho(0); Chr$(9); If getcoori(n, 1) = nndm(1) - 1 Then Print #1, Next n Print #1, Next i For n = 1 To nPoints Print #1, T(n) - 273.15; Chr$(9); If getcoori(n, 1) = nndm(1) - 1 Then Print #1, Next n Print #1, For specie = 0 To nspecies For n = 1 To nPoints Print #1, Conc(specie, n); Chr$(9); If getcoori(n, 1) = nndm(1) - 1 Then Print #1, Next n Print #1, Next specie End If End If Close #1 End Sub Sub watdisp() Dim x As String Dim i As Integer, j As Integer CF_dsp.mMinscale = Minscale CF_dsp.mMaxScale = MaxScale 'CF_dsp.Dimension = DispDim For i = 1 To nD CF_dsp.zzoom(i - 1) = DispZoom(i) CF_dsp.Offset(i - 1) = DispOffset(i) CF_dsp.ComboSlice(i - 1).Clear CF_dsp.ComboSlice(i - 1).List(0) = "Columns" CF_dsp.ComboSlice(i - 1).List(1) = "Rows" For j = 0 To nndm(i) - 1 CF_dsp.ComboSlice(i - 1).List(2 + j) = "Slice " & j Next j CF_dsp.ComboSlice(i - 1).ListIndex = SliceNrs(i) + 2 Next i For i = nD + 1 To 4 CF_dsp.zzoom(i - 1) = "---" CF_dsp.Offset(i - 1) = "---" CF_dsp.ComboSlice(i - 1).Clear CF_dsp.ComboSlice(i - 1).List(0) = "---" CF_dsp.ComboSlice(i - 1).ListIndex = 0 Next i If WhatDisp >= 1 Then CF_dsp.Combo1.ListIndex = WhatDisp - 1 Call doedimencombo CF_dsp.Show End Sub Sub doedimencombo() Dim i As Integer, j As Integer, jj As Integer Const watstr = " X Y Zd4" i = WhatDisp CF_dsp.Dimen.Clear 'For j = 0 To CF_dsp.Dimen.ListItems.Count ' jj = CF_dsp.Dimen.ListItems.Count ' CF_dsp.Dimen.List(jj).Delete 'Next j If i = 2 Or i = 3 Or i = 7 Or i = 8 Or i = 10 Then CF_dsp.Dimen.List(0) = "Total" For j = 1 To nD CF_dsp.Dimen.List(j) = Mid$(watstr, (j - 1) * 2 + 1, 2) Next j If DispDim <= CF_dsp.Dimen.ListCount - 1 Then CF_dsp.Dimen.ListIndex = DispDim End If CF_dsp.Dimen.Visible = True ElseIf i = 6 Then For j = 0 To nspecies CF_dsp.Dimen.List(j) = "Species " & j + 1 Next j If DispDim <= CF_dsp.Dimen.ListCount Then CF_dsp.Dimen.ListIndex = DispDim End If CF_dsp.Dimen.Visible = True Else CF_dsp.Dimen.Visible = False End If End Sub Sub RecallState() Dim filenaam As String Dim n As Long, nnpoints As Long Dim i As Integer Dim pt As Integer Dim dummy As String If nD = 0 Then MsgBox ("select an input geometry first") Else filenaam = mijnfilenaam & "_state" On Error Resume Next Open padnaam & "\" & filenaam For Input As #1 foutje = Err.Number On Error GoTo 0 If foutje <> 0 Then MsgBox ("An error occured opening the file") Else DoInit = 1 Call init Input #1, nnpoints If nnpoints <> nPoints Then MsgBox ("incorrect number of cells in saved state") Else Input #1, tt Input #1, tt 'reserved Input #1, tt 'reserved For i = 1 To nD For n = 1 To nPoints Input #1, mv(i, n, 0) mv(i, n, 5) = mv(i, n, 0) Next n Input #1, dummy Next i For n = 1 To nPoints Input #1, Q(n, 0) Q(n, 5) = Q(n, 0) Next n Input #1, dummy For specie = 0 To nspecies For n = 1 To nPoints Input #1, m(n, specie, 0) m(n, specie, 5) = m(n, specie, 0) Next n Input #1, dummy Next specie For n = 1 To nPoints For pt = 1 To 4 For specie = 0 To nspecies m(n, specie, pt) = 0 Next specie Q(n, pt) = 0 For i = 1 To nD mv(i, n, pt) = 0 Next i Next pt Next n DoInit = 1 Call CalcIntensity End If End If End If Close #1 End Sub Function getcoori(point, direction) Dim ii As Long 'index coordinate. Counts from zero, (but points count from 1!) ii = (point - 1) Mod dNghbr(direction + 1) getcoori = 0 + Int(ii / dNghbr(direction)) End Function Public Sub StartRadiation() Dim iv As Long Dim A1 As Single Dim F12 As Single If nD = 0 Then MsgBox ("select an inputrange first") Else On Error Resume Next Open padnaam & "\" & mijnfilenaam & "_rad" For Input As #1 foutje = Err.Number On Error GoTo 0 If foutje <> 0 Then MsgBox ("An error occured opening the file") Else Input #1, Nviews Input #1, Nviews 'reserved Input #1, Nviews 'reserved DoRadiation = 1 Boltzman_4 = 0.0000000567 Boltzman_4 = Boltzman_4 ^ (1 / 4) ' use to rescale temerpatures ReDim Viewconnect(Nviews, 1) As Long ReDim Viewfactor(Nviews) As Single For iv = 1 To Nviews Input #1, Viewconnect(iv, 0) Input #1, Viewconnect(iv, 1) Input #1, F12 Input #1, A1 Viewfactor(iv) = F12 * A1 Next iv MsgBox (Nviews & " Viewfactors imported") End If Close #1 End If End Sub Public Sub exitorpause() Dim ja As Integer CF_einde.Show End Sub