Fabien -
I think the code below should get you going.
The changes I made were:
- `xVelocity` and `zVelocity` changed to rank-0 CellVariables. FiPy *always*
solves at cell centers.
- Created a rank-1 FaceVariable to hold the velocity. The cell components must
be interpolated and manually inserted into this field at each sweep.
- Changed the sign of the DiffusionTerm in the heat equation; it's
unconditionally unstable otherwise
- Slowed down the time steps to better see the evolution
The result is not identical to the image you showed, but I think that's down to
boundary conditions, sign of gravity, etc.
- Jon
# -*- coding: utf-8 -*-
from fipy import *
import matplotlib.pyplot as plt
# Parameter
L = 1.
N = 50.
dL = L/N
alpha = 0.0002 # Thermical dilatation
landa = 0.6 # Thermical conductivity
ro0 = 1023. # Average volumic Mass
g = 10. # Gravity
# Mesh
mesh = Grid2D(nx=N, ny=N, dx=dL, dy=dL)
# Variables
dT = CellVariable(name = 'dT', mesh = mesh, value = 0.)
xVelocity = CellVariable(mesh=mesh, name='Xvelocity', value = 0.)
zVelocity = CellVariable(mesh=mesh, name='Zvelocity', value = 0.)
velocity = FaceVariable(mesh=mesh, name='velocity', rank=1)
# Init Condition
x = mesh.cellCenters[0]
dT.setValue(1.)
dT.setValue(-1., where = x > L/2)
# Viewer
viewer = None
if __name__ == '__main__':
viewer = Viewer(vars=dT, datamin=-1., datamax=1.)
viewer.plotMesh()
raw_input("<Enter to continue>...")
# Boussinesq equations
D = landa/ro0
eqX = (TransientTerm(var=xVelocity) + ConvectionTerm (var=xVelocity,
coeff=velocity) == 0)
eqZ = (TransientTerm(var=zVelocity) + ConvectionTerm (var=zVelocity,
coeff=velocity) + alpha*g*dT == 0)
eqT = (TransientTerm(var=dT) + ConvectionTerm (var=dT, coeff=velocity) ==
DiffusionTerm(var=dT, coeff=D))
eq = eqX & eqZ & eqT
# Solving Boussinesq equations
timeStepDuration = 1 * dL**2 / (2 * D)
steps = 50
sweeps = 5
for step in range(steps):
for sweep in range(sweeps):
eq.sweep(dt=timeStepDuration)
velocity[0] = xVelocity.arithmeticFaceValue
velocity[1] = zVelocity.arithmeticFaceValue
velocity[..., mesh.exteriorFaces.value] = 0.
if viewer is not None:
viewer.plot()
plt.pause(0.1)
raw_input("<Enter to finsh>… ")
> On Aug 14, 2018, at 11:22 AM, fgendr01 <fabien.gend...@univ-lr.fr> wrote:
>
> Hello,
>
> I’m a PhD Student in La Rochelle University (France) and I’m working on
> Boussinesq Approximation.
> I’m starting with fipy and I’d like to make a little algorithm to solve these
> equations in 2D (x,z) :
> <PastedGraphic-1.png>
> In my problem T = T0 + dT.
> (ux, uz) is the velocity.
>
> So, here is my algorithm :
> # -*- coding: utf-8 -*-
> from fipy import *
> import matplotlib.pyplot as plt
> # Parameter
> L = 1.
> N = 50.
> dL = L/N
> alpha = 0.0002 # Thermical dilatation
> landa = 0.6 # Thermical conductivity
> ro0 = 1023. # Average volumic Mass
> g = 10. # Gravity
> # Mesh
> mesh = Grid2D(nx=N, ny=N, dx=dL, dy=dL)
> # Variables
> dT = CellVariable(name = 'dT', mesh = mesh, value = 0.)
> xVelocity = FaceVariable(mesh=mesh, name='Xvelocity', value = 0., rank=1)
> zVelocity = FaceVariable(mesh=mesh, name='Zvelocity', value = 0., rank=1)
> # Init Condition
> x = mesh.cellCenters[0]
> dT.setValue(1.)
> dT.setValue(-1., where = x > L/2)
> # Viewer
> viewer = None
> if __name__ == '__main__':
> viewer = Viewer(vars=dT, datamin=-1., datamax=1.)
> viewer.plotMesh()
> raw_input("<Enter to continue>...")
> # Boussinesq equations
> D = landa/ro0
> velocity = ((xVelocity), (zVelocity))
> eqX = (TransientTerm(var=xVelocity) + ConvectionTerm (var=xVelocity,
> coeff=velocity) == 0)
> eqZ = (TransientTerm(var=zVelocity) + ConvectionTerm (var=zVelocity,
> coeff=velocity) + alpha*g*dT == 0)
> eqT = (TransientTerm(var=dT) + ConvectionTerm (var=dT, coeff=velocity) ==
> -DiffusionTerm(var=dT, coeff=D))
> eq = eqX & eqZ & eqT
> # Solving Boussinesq equations
> timeStepDuration = 10 * dL**2 / (2 * D)
> steps = 50
> for step in range(steps):
> eq.solve(dt=timeStepDuration)
> if viewer is not None:
> viewer.plot()
> plt.pause(0.1)
> raw_input("<Enter to finsh>… »)
>
> My error : all the input arrays must have same number of dimensions
>
> I’d like to obtain something like this :
> <PastedGraphic-3.png>
>
> If someone can help me, I will be very happy.
>
> Thanks a lot.
>
> Fabien G.
> _______________________________________________
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
