Re: Boussinesq Equations

[Thibault Bridel-Bertomeu]

Hi again Fabien,

Well, I would say that physically it makes sense. With such boundary
condition, I think the problem is well posed. Now it is all about the
numerics.
Have you tried ? What does the solver yield ?

Cheers
Thibault
Le lun. 17 sept. 2018 à 21:49, fgendr01  a
écrit :

> 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)
> T.faceGrad.constrain(0, mesh.facesTop)
> T.faceGrad.constrain(0, mesh.facesBottom)
> 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 examples to find an adiabatic or zero-gradient boundary condition !
>
> 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  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("...")
>> # 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("… ")
>>
>> 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  a écrit :
>>
>> T

Re: Boussinesq Equations

[fgendr01]

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)
T.faceGrad.constrain(0, mesh.facesTop)
T.faceGrad.constrain(0, mesh.facesBottom)
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 
>  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 
> examples to find an adiabatic or zero-gradient boundary condition ! 
> 
> 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  > 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("...")
> # 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("… ")
> 
> 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 > > a