Re: Boussinesq Equations

[Guyer, Jonathan E. Dr. (Fed)]

Fabien -

I don't know where the problem lies, but my recommendation would be to 
systematically take the problem back toward the code I previously gave you. 
While that code didn't produce the image you were looking for, the behavior 
seemed consistent and it wasn't unstable.

As far as I can tell, in addition to changing the initial condition and the 
boundary conditions, you've increased the Peclét number by three orders of 
magnitude (by including c in DT) and increased the thermal dilatation by factor 
of 10 (by changing to a temperature delta from a dimensionless value +/- 1 to 
an absolute temperature difference). By changing these parameters back, you may 
regain the earlier stability and get some insight as to what the problem is.

- Jon

> On Sep 18, 2018, at 3:10 AM, fgendr01  wrote:
> 
> Yes I tried without kinematic boundary conditions and there is no difference !
> Maybe the scheme is not stable.
> Other ideas ?
> 
> Fabien
> 
>> Le 17 sept. 2018 à 22:27, Thibault Bridel-Bertomeu 
>>  a écrit :
>> 
>> Have you tried without the kinematic boundary conditions ? Just constrain T 
>> and dT/dn, what happens then ? 
>> 
>> Le lun. 17 sept. 2018 à 22:17, fgendr01  a écrit :
>> Alas, it doesn’t work :-(
>> I tried to change the CFL condition but there is not real differences…
>> 2 screenshots at t = 0 and after 20s...
>> 
>> 
>> 
>> If someone have the solution ?
>> 
>> here is the code with last modifications :
>> 
>> # -*- 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)#Tmax to the 
>> Left
>> T.constrain(Tmin, mesh.facesRight)   #Tmin to the right
>> T.faceGrad.constrain(0, mesh.facesTop)   #Adiabatic wall to the 
>> Top
>> T.faceGrad.constrain(0, mesh.facesBottom)#Adiabatic wall to the 
>> Bottom
>> xVelocity.constrain(0, mesh.facesLeft)
>> xVelocity.constrain(0, mesh.facesRight)
>> zVelocity.constrain(0, mesh.facesTop)
>> zVelocity.constrain(0, mesh.facesBottom)
>> # Viewer
>> 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.001*dL**2/(2*DT)#CFL Condition
>> 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("… ")
>> 
>> 
>> Sorry for all my messages...
>> 
>> Fabien
>> 
>> 
>>> Le 17 sept. 2018 à 21:52, Thibault Bridel-Bertomeu 
>>>  a écrit :
>>> 
>>> 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)
>>> 

Re: Boussinesq Equations

[fgendr01]

Yes I tried without kinematic boundary conditions and there is no difference !
Maybe the scheme is not stable.
Other ideas ?

Fabien

> Le 17 sept. 2018 à 22:27, Thibault Bridel-Bertomeu 
>  a écrit :
> 
> Have you tried without the kinematic boundary conditions ? Just constrain T 
> and dT/dn, what happens then ? 
> 
> Le lun. 17 sept. 2018 à 22:17, fgendr01  > a écrit :
> Alas, it doesn’t work :-(
> I tried to change the CFL condition but there is not real differences…
> 2 screenshots at t = 0 and after 20s...
> 
> 
> 
> If someone have the solution ?
> 
> here is the code with last modifications :
> 
> # -*- 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) #Tmax to the 
> Left
> T.constrain(Tmin, mesh.facesRight)#Tmin to the right
> T.faceGrad.constrain(0, mesh.facesTop)#Adiabatic wall to the 
> Top
> T.faceGrad.constrain(0, mesh.facesBottom) #Adiabatic wall to the 
> Bottom
> xVelocity.constrain(0, mesh.facesLeft)
> xVelocity.constrain(0, mesh.facesRight)
> zVelocity.constrain(0, mesh.facesTop)
> zVelocity.constrain(0, mesh.facesBottom)
> # Viewer
> 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.001*dL**2/(2*DT) #CFL Condition
> 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("… ")
> 
> 
> Sorry for all my messages...
> 
> Fabien
> 
> 
>> Le 17 sept. 2018 à 21:52, Thibault Bridel-Bertomeu 
>> mailto:thibault.bridellel...@gmail.com>> a 
>> écrit :
>> 
>> 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 
>>> mailto: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