# Re: Boussinesq Equations

```Thanks Thibault,
You’re right, the top ans the bottom walls are adiabatic and I tried to put
velocity conditions.
What do you think about these conditions :```
```
# Boundary Conditions
T.constrain(Tmax, mesh.facesLeft)
T.constrain(Tmin, mesh.facesRight)
xVelocity.constrain(0, mesh.facesLeft)
xVelocity.constrain(0, mesh.facesRight)
zVelocity.constrain(0, mesh.facesTop)
zVelocity.constrain(0, mesh.facesBottom)

Thanks,

Fabien

> Le 17 sept. 2018 à 21:29, Thibault Bridel-Bertomeu
> <thibault.bridellel...@gmail.com> a écrit :
>
> Good evening Fabien,
>
> I have not been following the development or the answers given to you by
> other members, so, apologies if I stumble on something that’s already been
> told.
>
> I think the figure you send with your question is the answer. From what I can
> see, at the top and the bottom of the domain, there is some extrapolation
> going on, or Neumann sort of boundary condition : the temperature gradient is
> zero. In physical terms, those two walls look adiabatic to me.
>
> I don’t have my laptop anywhere close but I would encourage you to browse the
>
> As for the mathematical aspect of your question: when it comes to the
> Navier-Stokes equations, you have to impose boundary conditions on all
> boundary patches of the domain. The boundary conditions can change depending
> on what you want to happen, but you have to have one. Also it is surprising
> to me that only imposing the temperature on the left and right patches is
> enough: I’m not sure about fipy core behavior, but with Navier-Stokes, you
> gotta impose temperature and velocity to have a well-posed problem !
>
> Cheers and bon courage !
>
> T
> Le lun. 17 sept. 2018 à 21:09, fgendr01 <fabien.gend...@univ-lr.fr
> <mailto:fabien.gend...@univ-lr.fr>> a écrit :
> Hi Jon and the community,
>
> Since a few weeks, I develop a small code to solve a 2D problem of water
> transport in a simple mesh.
> My problem is the resolution of the Navier-Stokes equations with the
> Boussinesq approximation in 2D with a thermal gradient.
> I try to solve these equations :
>
> A member (Jonathan) helped me and wrote a part of the code but the code seems
> not to be stable and I don’t know why !
> Can you help me ?
> There must be a problem with the boundary conditions.
> I want to impose a fixed temperature on the left wall (303K) and a fixed
> temperature on the right wall (283K). Should I put a condition on the upper
> and lower walls?
> Here is my code :
>
> # -*- coding: utf-8 -*-
> from fipy import *
> import matplotlib.pyplot as plt
> # Parameter
> L = 1.
> N = 100.
> dL = L/N
> T0 = 293.
> Tmax = 303.                                   # Temperature max
> Tmin = 283.                                   # Temperature min
> alpha = 0.0002                                # Thermal dilatation
> lamda = 0.6                                   # Thermal conductivity (W/m/K)
> c = 4186.                                             # Specific Heat of sea
> water (J/kg/K)
> ro0 = 1023.                                   # average density of sea water
> g = 10.                                               # Gravity (m/s2)
> DT = lamda/(c*ro0)                            # Thermal diffusivity (m2/s)
> # Mesh
> mesh = Grid2D(nx=N, ny=N, dx=dL, dy=dL)
> # Variables
> T = CellVariable(mesh=mesh, name = 'T', value = T0)
> 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
> T.setValue(T0)
> # Boundary Conditions
> T.constrain(Tmax, mesh.facesLeft)
> T.constrain(Tmin, mesh.facesRight)
> # Viewer init conditions
> viewer = None
> if __name__ == '__main__':
>       viewer = Viewer(vars=T, datamin=Tmin, datamax=Tmax)
>       viewer.plot()
>       raw_input("<Enter to continue>...")
> # Boussinesq equations
> eqX = (TransientTerm(var=xVelocity) + ConvectionTerm (var=xVelocity,
> coeff=velocity) == 0)
> eqZ = (TransientTerm(var=zVelocity) + ConvectionTerm (var=zVelocity,
> coeff=velocity) == alpha*g*(T-T0))
> eqT = (TransientTerm(var=T) + ConvectionTerm (var=T, coeff=velocity) ==
> DiffusionTerm(var=T, coeff=DT))
> eq = eqX & eqZ & eqT
> # Solving Boussinesq equations
> timeStepDuration = 0.1*dL**2/(2*DT)
> 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 finish>… ")
>
> I’d like to obtain something like this :
>
>
> Thanks again for your help, maybe Jon again ;-)
>
> Fabien
>
>
>> Le 16 août 2018 à 07:20, fgendr01 <fabien.gend...@univ-lr.fr
>> <mailto:fabien.gend...@univ-lr.fr>> a écrit :
>>
>> Thanks a lot Jon.
>> I will try to apply your remarks but I will certainly still need your help
>> again ;-) to develop the code.
>>
>> Thanks again.
>>
>> Fabien
>>
>>> Le 15 août 2018 à 19:46, Guyer, Jonathan E. Dr. (Fed)
>>> <jonathan.gu...@nist.gov <mailto:jonathan.gu...@nist.gov>> a écrit :
>>>
>>> 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
>>>> <mailto: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 <mailto:fipy@nist.gov>
>>>> http://www.ctcms.nist.gov/fipy <http://www.ctcms.nist.gov/fipy>
>>>> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy
>>>> <https://email.nist.gov/mailman/listinfo/fipy> ]
>>>
>>>
>>> _______________________________________________
>>> fipy mailing list
>>> fipy@nist.gov <mailto:fipy@nist.gov>
>>> http://www.ctcms.nist.gov/fipy <http://www.ctcms.nist.gov/fipy>
>>> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy
>>> <https://email.nist.gov/mailman/listinfo/fipy> ]
>>
>
> _______________________________________________
> fipy mailing list
> fipy@nist.gov <mailto:fipy@nist.gov>
> http://www.ctcms.nist.gov/fipy <http://www.ctcms.nist.gov/fipy>
>   [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy
> <https://email.nist.gov/mailman/listinfo/fipy> ]
> <PastedGraphic-2.png><PastedGraphic-1.png><PastedGraphic-2.png>_______________________________________________
> fipy mailing list
> fipy@nist.gov <mailto:fipy@nist.gov>
> http://www.ctcms.nist.gov/fipy <http://www.ctcms.nist.gov/fipy>
>  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy
> <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 ]
```