Hi,
in the addition is the code of the simulation. The file is a little bit larger.
I have there a lot of correlations inside the file.
And the Problem is a little bit larger, because I have there also a fluid
outside the tube and inside the tube a metal foam. But my Problem is the same,
that I don´t know how to write the flux boundary between the sperated parts
(fluid,tupe and foam)
Thanks for your reply
Alexander
-------- Original-Nachricht --------
> Datum: Wed, 5 May 2010 11:44:53 -0400
> Von: William Gathright <[email protected]>
> An: Multiple recipients of list <[email protected]>
> Betreff: Re: flux boundary
>
> Hi Alexander,
> I would suggest that we need a little more information about your
> problem to be as helpful as possible. Are you working in 1D or 2D?
> Are both the tube and the fluid included in your simulation domain, or
> just the tube?
>
> The reason I ask is that the boundary conditions are designed to be
> applied at the external faces of your domain. For example, if your
> domain features a tube surrounded by a fluid, the boundary conditions
> should only pertain to the fluid (since the fluid is the only thing
> touching the boundary of your domain).
>
> If you want a flux that is calculated in every cell of your system,
> you'll get that by solving your heat equations. For that, you may
> want to try declaring "Heat" or "Temperature" to be your CellVariable
> instead of "Fluid" or "Tube". Again, with a little more information
> about your problem and your code we may be able to provide more
> detailed help.
>
> - Will Gathright
> Rensselaer Polytechnic Institute
> Department of Materials Science and Engineering
>
> On Wed, May 5, 2010 at 11:14 AM, Alexander Holbach <[email protected]>
> wrote:
> >
> > Hi,
> > I am using the PDE solver in a system.
> > But I have problems with my boundarys. The system shall simulate a heat
> transfer from tube to Fluid. So I wrote two systems with two seperated PDE
> solvers.
> >
> > I'd like to know how I may implement the following boundary condition
> over the right and left faces:
> > the equations are:
> >
> > …
> > tube = CellVariable(name="temperature [C]", mesh=mesh,value = 0.)
> > tube.equation = (ImplicitDiffusionTerm(coeff=k_r))
> > tube.equation.solve(var=tube,
> > boundaryConditions=boundaryConditionsWand,
> > solver=DefaultAsymmetricSolver(tolerance=1.e-100,
> iterations=10000))
> >
> > …
> > and
> > …
> >
> > Fluid = CellVariable(name="temperature [C]", mesh=mesh,value = 0.)
> > Fluid.equation = ImplicitDiffusionTerm(coeff=k_fg)==
> rho*cp*ExplicitUpwindConvectionTerm(coeff=((v_f,)))
> > Fluid.equation.solve(var=Fluid,
> > boundaryConditions=boundaryConditionsFluid,
> > solver=DefaultAsymmetricSolver(tolerance=1.e-100,
> iterations=10000))
> >
> > and now I want to implement the boundary condition, that there is a Flux
> in each mesh cell from the tube to the Fluid.
> > What I can write is:
> > …
> >
> > FixedFlux(faces=mesh.getFacesLeft(),
> value=h_rf*a_rf*(tube(1,0,6000)-Fluid(1,0,6000))),
> > FixedFlux(faces=mesh.getFacesRight(),
> value=h_rf*a_rf*(tube(1,0,6000)-Fluid(1,0,6000))))
> > …
> >
> > but this is a fix Flux and I want that the Flux will be calculate in
> each cell of the system.
> >
> > Is this possible?
> >
> > Thank for your help
> > Alexander
> > --
> > GRATIS für alle GMX-Mitglieder: Die maxdome Movie-FLAT!
> > Jetzt freischalten unter http://portal.gmx.net/de/go/maxdome01
> >
> >
>
--
GRATIS für alle GMX-Mitglieder: Die maxdome Movie-FLAT!
Jetzt freischalten unter http://portal.gmx.net/de/go/maxdome01
"""
Created on Wed 5. 3. 10:56:36 2010
@author: holbach_a
"""
import fipy as fp
from fipy import *
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from fipy.variables.cellVariable import CellVariable
from fipy.boundaryConditions.fixedValue import FixedValue
from fipy.boundaryConditions.fixedFlux import FixedFlux
from fipy.viewers import make
from fipy.meshes.grid2D import Grid2D
from fipy.terms.implicitDiffusionTerm import ImplicitDiffusionTerm
from fipy.terms.exponentialConvectionTerm import ExponentialConvectionTerm
from scipy.special import gamma
from fipy.meshes.numMesh.mesh2D import Mesh2D
from fipy.viewers.tsvViewer import TSVViewer
##################################################################
diff=0.000000000000011 # Genauigigkeit der Lösung
steps = 20.
T_w_start=70.
T_s_start=25.
#initials:
T_ef = 15. # Temperatur des Fluids am Eingang [°C]
T_hwe = 80. # Temperatur des Heizw. am Eingang [°C]
# constants:
v_fw = 1 # Geschwindigkeit des Fluids im Mantel [m/s]
v = 100 # Volumenstrom des Fluids [l/h]
L = 0.10 # LÀnge [m]
D = 0.02 # Durchmesser [m]
eps = 0.95 # PorösitÀt [-]
D_m = 0.01 # Manteldickedicke Rohr [m]
D_r = 0.002 # Wanddicke Rohr [m]
k_f = 0.58 # WÀrmeleitfÀhigkeit Fluid [W/(m*K)]
[-]
k_s = 260. # WÀrmeleitfÀhigkeit Schaum(Metall) [W/(m*K)]
d_p = 0.003 # Durchmesser der Pore [m]
d_f = 0.0003 # Faser Durchmesser [m]
rho = 1000. # Dichte des Fluids [kg/m^3]
cp = 4190. # WÀrmekapazitÀt des Fluid [J/(kg*K)]
mu = 0.001 # dynamische ViskositÀt [Pa*s]
C_d = 0.1 # thermischer Dispersionskoeffizient [-]
r = 0.09 # Konstante (aus [Calmidi 2000]) [-]
lam = 0.42
C = 0.52
#a_rs = 0.0002 # spezifische OberflÀche [1/m]
h_rs = 5000 # WÀrmeÌbergang Rohr/Schaum [W/(m^2*K)]
k_r = 50 # WÀrmeleitfÀhigkeit Rohr [W/(m*K)]
######################################################################
v_0 = v/(1000*3600*pi*(D/2)**2)
print(v_0,'m/s')
v_f=v_0/eps
# Berechnung der FlÀche zwischen Rohr und Schaum und der FlÀche zwischen Wand
und Rohr
a_wr = pi*(D+D_r)*L
a_rs = pi*D*L*(1-eps)
a_rf = pi*D*L*eps
# Berechnung des WÀrmeÌbergangs Mantel Rohr nach Gnielinski
Re=v_f*d_p*rho/mu
Pr=cp*mu/k_f
Re_m=v_fw**pi/2*D_m*rho/mu
Nu_t=0.037*Re_m**(0.8)*Pr/(1+2.443*Re_m*(-1)*(Pr**(2/3)-1))
Nu_l= 0.664*np.sqrt(Re_m)*(Pr)**(1/3)
h_wr=(0.3+np.sqrt(Nu_l**2+Nu_t**2))*k_f/(D_m*pi/2)
# Berechnung des WÀrmeÌbergangs Rohr/Fluid mit konstantem heat Flux nach Shah
h_rf = 1.953*(Re*Pr*D/L)**(1/3)*k_f/(D)
# Iteration:
# mesh and lenght of the grid
nx = 100
ny = nx
dx_1 = D_m/2
dy_1 = L/3
#dx_2 = D_m/5
#dy_2 = L/2
class Cylinderizer(Mesh2D):
def __init__(self, mesh):
Mesh2D.__init__(self,
mesh.getVertexCoords(),
mesh.faceVertexIDs,
mesh.cellFaceIDs)
def _getFaceAreas(self):
return Mesh2D._getFaceAreas(self) * self.getFaceCenters()[0]
def getCellVolumes(self):
return Mesh2D.getCellVolumes(self) * self.getCellCenters()[0]
mesh = Cylinderizer(Grid2D(dx=D_m/100, dy=L/100, nx=nx, ny=ny))
#mesh = Cylinderizer(Grid2D(dx=D_m*3/200, dy=L/100, nx=200/3, ny=100) +
Grid2D(dx=D_m*2/500, dy=L/100, nx=500/2, ny=100))
# defintion of the variable
Mantel = CellVariable(name="temperature [C]", mesh=mesh,value = 0.)
T_m=Mantel
# Boudary Conditions
RHSBC = (((1,),) * mesh.getFacesLeft()).getDivergence()
boundaryConditionsMantel = (
FixedValue(faces=mesh.getFacesTop(), value=T_hwe),
FixedValue(faces=mesh.getFacesRight(), value=T_hwe),
FixedFlux(faces=mesh.getFacesRight(), value=0),
FixedFlux(faces=mesh.getFacesLeft(), value=-h_wr*a_wr*(T_hwe-T_w_start)))
# equations
Mantel.equation = ImplicitDiffusionTerm(coeff=k_f)==
rho*cp*ExplicitUpwindConvectionTerm(coeff=((0,v_fw)))
# initial Solver
Mantel.equation.solve(var=Mantel,
boundaryConditions=boundaryConditionsMantel,
solver=DefaultAsymmetricSolver(tolerance=1.e-100, iterations=10000))
#####################################################################
# mesh and lenght of the grid
nx = 100
ny = nx
dx = D_r/100
dy = L/100
'''
class Cylinderizer(Mesh2D):
def __init__(self, mesh):
Mesh2D.__init__(self,
mesh.getVertexCoords(),
mesh.faceVertexIDs,
mesh.cellFaceIDs)
def _getFaceAreas(self):
return Mesh2D._getFaceAreas(self) * self.getFaceCenters()[0]
def getCellVolumes(self):
return Mesh2D.getCellVolumes(self) * self.getCellCenters()[0]
'''
mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
#mesh = Cylinderizer(Grid2D(dx=dx, dy=dy, nx=nx, ny=ny) + (Grid2D(dx=dx, dy=dy,
nx=nx, ny=dy) + ((0,), ((D+D_r)/2,))))
# defintion of the variable
Wand = CellVariable(name="temperature [C]", mesh=mesh,value = 0.)
T_w=Wand
# Boudary Conditions
boundaryConditionsWand = (
FixedValue(faces=mesh.getFacesRight(), value=T_w_start),
FixedFlux(faces=mesh.getFacesRight(), value=h_wr*a_wr*(T_hwe-T_w_start)),
FixedFlux(faces=mesh.getFacesLeft(),
value=-h_rs*a_rs*(T_w_start-T_s_start)),
)
# equations
Wand.equation = (ImplicitDiffusionTerm(coeff=k_r))
# initial Solver
Wand.equation.solve(var=Wand,
boundaryConditions=boundaryConditionsWand,
solver=DefaultAsymmetricSolver(tolerance=1.e-100, iterations=10000))
##########################################################################
T_es = T_w (1,0,5000)
##########################################################################
# mesh and lenght of the grid
nz = 100
nr = nz
dz = L/100
dr = D/200
mesh = CylindricalGrid2D(dr=dr, dz=dz, nr=nr, nz=nz)
# defintion of the variable
Schaum = CellVariable(name="temperature [C]", mesh=mesh,value = 0.)
Fluid = CellVariable(name="temperature [C]", mesh=mesh,value = 0.)
T_s=Schaum
T_f=Fluid
# Berechnug der Prandtl-Zahl
Pr=cp*mu/k_f
# Berechnung der Reynolds-Zahl
v_f=v_0/eps
Re=v_f*d_p*rho/mu
# Berechnung des WÀrmeÌbergangkoeffizienten Schaum/Fluid h_sf
h_sf=k_f/d_p*0.62*lam*Pr**(0.36)*Re**(0.5) #VDI
#h_sf_2=k_f/d_p*C*Re**(0.5)*Pr**(0.37) #Zukauskas
#Paper
# Berechnung der spezifischen FlÀche a_sf (aus[Calmidi 200])
a_sf=3*pi*d_f/d_p**2
# Berechnung von (b/l)
x_1=-r+np.sqrt(r**2+4*(1-eps)*np.sqrt(3)/2*(2-r*1+4/np.sqrt(3))/3)
x_2=(4-2*r*(1+4/(np.sqrt(3))))/3
x_3=x_1/x_2
# Berechnung der WÀrmeleitfÀhigkeit des Schaums k_se (aus [Calmidi 2000])
m_1=2/np.sqrt(3)
n_1=r*(x_3)/((1+(x_3))*k_s/3)
o_1=(3*(1-r))/(2*k_s)
p_1=((np.sqrt(3)/2)-(x_3))/(4*r/(3*np.sqrt(3))*(x_3)*k_s)
k_se=(m_1*(n_1+o_1+p_1))**(-1)
# Berechnung der WÀrmeleitfÀhigkeit des Fluids k_fg (aus [Calmidi 2000])
m_2=2/np.sqrt(3)
n_2=r*(x_3)/(k_f+(1+(x_3))*(-k_f)/3)
o_2=((1-r)*x_3)/(k_f+2/3*x_3*(-k_f))
p_2=((np.sqrt(3)/2)-(x_3))/(k_f+4*r/(3*np.sqrt(3))*(x_3)*(-k_f))
m_3=2/np.sqrt(3)
n_3=r*(x_3)/(k_f+(1+(x_3))*(k_s-k_f)/3)
o_3=((1-r)*x_3)/(k_f+2/3*x_3*(k_s-k_f))
p_3=((np.sqrt(3)/2)-(x_3))/(k_f+4*r/(3*np.sqrt(3))*(x_3)*(k_s-k_f))
v_f=v_0/eps
v = np.array([v_f,0])
k_fe=(m_2*(n_2+o_2+p_2))**(-1)
k_e=(m_3*(n_3+o_3+p_3))**(-1)
k_d=k_e*C_d*v_f/v_0
k_fg=k_fe+k_d
# Boudary Conditions
boundaryConditionsFluid = (
FixedFlux(faces=mesh.getFacesTop(), value=0),
FixedValue(faces=mesh.getFacesBottom(), value=T_ef),
FixedFlux(faces=mesh.getFacesRight(), value=0),
FixedFlux(faces=mesh.getFacesLeft(), value=0),
)
boundaryConditionsSchaum = (
FixedFlux(faces=mesh.getFacesTop(), value=0),
FixedValue(faces=mesh.getFacesRight(), value=T_s_start),
FixedValue(faces=mesh.getFacesLeft(), value=T_s_start),
FixedFlux(faces=mesh.getFacesLeft(), value=h_rs*a_rs*(T_w_start-T_s_start)),
FixedFlux(faces=mesh.getFacesRight(), value=h_rs*a_rs*(T_w_start-T_s_start))
)
# equations
Schaum.equation = (ImplicitDiffusionTerm(coeff=k_se)
- ImplicitSourceTerm(h_sf*a_sf) +h_sf*a_sf*T_f)
Fluid.equation = ImplicitDiffusionTerm(coeff=k_fg)==
rho*cp*ExplicitUpwindConvectionTerm(coeff=((v_f,)))+
ImplicitSourceTerm(h_sf*a_sf) - h_sf*a_sf*T_s
# initial Solver
Schaum.equation.solve(var=Schaum,
boundaryConditions=boundaryConditionsSchaum,
solver=DefaultAsymmetricSolver(tolerance=1.e-100, iterations=10000))
Fluid.equation.solve(var=Fluid,
boundaryConditions=boundaryConditionsFluid,
solver=DefaultAsymmetricSolver(tolerance=1.e-100, iterations=10000))
#################################################################
#################################################################
residual = 10
for step in range(steps):
#while residual > diff:
# print Fluid_alt
boundaryConditionsMantel = (
FixedValue(faces=mesh.getFacesTop(), value=T_hwe),
FixedValue(faces=mesh.getFacesRight(), value=T_hwe),
FixedFlux(faces=mesh.getFacesRight(), value=0),
FixedFlux(faces=mesh.getFacesLeft(),
value=-h_wr*a_wr*(Mantel(1,0,6000)-Wand(1,0,6000))))
boundaryConditionsWand = (
FixedValue(faces=mesh.getFacesRight(), value=Mantel(1,0,5000)),
FixedFlux(faces=mesh.getFacesRight(),
value=h_wr*a_wr*(Mantel(1,0,6000)-Wand(1,0,6000))),
FixedFlux(faces=mesh.getFacesLeft(),
value=-h_rs*a_rs*(Wand(1,0,6000)-Schaum(1,0,6000))-h_rf*a_rf*(Wand(1,0,5000)-Schaum(1,0,5000))),)
boundaryConditionsSchaum = (
FixedFlux(faces=mesh.getFacesTop(), value=0),
# FixedValue(faces=mesh.getFacesRight(), value=T_w_neu),
# FixedValue(faces=mesh.getFacesLeft(), value=T_w_neu),
FixedFlux(faces=mesh.getFacesLeft(),
value=h_rs*a_rs*(Wand(1,0,6000)-Schaum(1,0,6000))),
FixedFlux(faces=mesh.getFacesRight(),
value=h_rs*a_rs*(Wand(1,0,6000)-Schaum(1,0,6000))))
boundaryConditionsFluid = (
FixedFlux(faces=mesh.getFacesTop(), value=0),
FixedValue(faces=mesh.getFacesBottom(), value=T_ef),
# FixedFlux(faces=mesh.getFacesRight(), value=0),
FixedFlux(faces=mesh.getFacesLeft(), value=0),
FixedFlux(faces=mesh.getFacesLeft(),
value=h_rf*a_rf*(Wand(1,0,6000)-Schaum(1,0,6000))),
FixedFlux(faces=mesh.getFacesRight(),
value=h_rf*a_rf*(Wand(1,0,6000)-Schaum(1,0,6000))))
nx = 100
ny = nx
dx = D_m/100
dy = L/100
mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
Mantel = CellVariable(name="temperature [C]", mesh=mesh,value = 0.)
Mantel.equation = ImplicitDiffusionTerm(coeff=k_f)==
rho*cp*ExplicitUpwindConvectionTerm(coeff=((0,v_fw)))
Mantel.equation.solve(var=Mantel,
boundaryConditions=boundaryConditionsMantel,solver=DefaultAsymmetricSolver(tolerance=1.e-100,
iterations=10000))
T_m=Mantel
nx = 100
ny = nx
dx = D_r/100
dy = L/100
mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
RHSBC = (((1,),) * mesh.getFacesRight()).getDivergence()
Wand.equation = (ImplicitDiffusionTerm(coeff=k_r))
Wand.equation.solve(var=Wand, boundaryConditions=boundaryConditionsWand,
solver=DefaultAsymmetricSolver(tolerance=1.e-100, iterations=10000))
T_w=Wand
nz = 100
nr = nz
dz = L/100
dr = D/200
mesh = CylindricalGrid2D(dr=dr, dz=dz, nr=nr, nz=nz)
Schaum.equation = (ImplicitDiffusionTerm(coeff=k_se)-
ImplicitSourceTerm(h_sf*a_sf) +h_sf*a_sf*T_f)
Fluid.equation = ImplicitDiffusionTerm(coeff=k_fg)==
rho*cp*ExplicitUpwindConvectionTerm(coeff=((v_f,)))+
ImplicitSourceTerm(h_sf*a_sf) - h_sf*a_sf*T_s
Schaum.equation.solve(var=Schaum,
boundaryConditions=boundaryConditionsSchaum,solver=DefaultAsymmetricSolver(tolerance=1.e-100,
iterations=10000))
Fluid.equation.solve(var=Fluid,
boundaryConditions=boundaryConditionsFluid,solver=DefaultAsymmetricSolver(tolerance=1.e-100,
iterations=10000))
T_s=Schaum
T_f=Fluid
Fluid_neu=Fluid(1,0,9999)
print (Fluid_neu,'Fluidtemperatur [C]')
residual = Fluid.equation.sweep(var=Fluid,
boundaryConditions=boundaryConditionsFluid,solver=DefaultAsymmetricSolver(tolerance=1.e-100,
iterations=10000))
##################################################################
# initial plot
if __name__ == '__main__':
viewer = make(vars=Fluid, datamin=15., datamax=25.)
title('fluid')
xlabel('radius [m]')
ylabel('lenght [m]')
xticks((0, D/2), color = 'k', size = 10)
viewer.plot()
TSVViewer(vars=Fluid).plot(filename="Fluidtemperatur.tsv")
'''
viewer = make(vars=Mantel, datamin=58, datamax=65.)
title('shell')
xlabel('size [m]')
ylabel('lenght [m]')
xticks((0,D_m), color = 'k', size = 10)
viewer.plot()
viewer = make(vars=Wand, datamin=58, datamax=80.)
title('tube')
xlabel('size [m]')
ylabel('lenght [m]')
xticks((0,D_r), color = 'k', size = 10)
viewer.plot()
viewer = make(vars=Schaum, datamin=15., datamax=40.)
title('foam')
xlabel('radius [m]')
ylabel('lenght [m]')
xticks((0, D/2), color = 'k', size = 10)
viewer.plot()
