Re: Boussinesq Equations

[Daniel Wheeler]

It certainly won't be that simple. You'll need to look through the
literature to figure it out. I suspect it will require re-introducing
a pressure into the equations and then transforming the continuity
equation into an equation for pressure.
On Thu, Oct 11, 2018 at 4:06 AM fgendr01  wrote:
>
> Thanks Daniel for your answer,
> I used the continuity equation (nabla vector velocity) to obtain the equation 
> in temperature.
> But actually, maybe we should use it as a constraint to my problem.
> But now, how can I create this constraint ? with a new equation ?
> I tried : velocity.divergence == 0 but it doesn’t work.
>
> Thank you,
>
> Fabien
>
>
>
>
> Le 10 oct. 2018 à 17:43, Daniel Wheeler  a écrit :
>
> Don't you still have a $\nabla . \vec{u} = 0$ equation though? It
> doesn't go away. That equation becomes like a constraint.
>
> https://www.comsol.com/multiphysics/boussinesq-approximation
>
> On Wed, Oct 10, 2018 at 5:58 AM fgendr01  wrote:
>
>
> Hi Daniel,
> Thank you for your answer.
> I thank you for trying to solve my problem.
> About my set of the equation here is my reasoning.
>
>
>
> --
> Daniel Wheeler
> ___
> 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 ]



-- 
Daniel Wheeler

___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: Boussinesq Equations

[fgendr01]

Thanks Daniel for your answer,
I used the continuity equation (nabla vector velocity) to obtain the equation 
in temperature.
But actually, maybe we should use it as a constraint to my problem.
But now, how can I create this constraint ? with a new equation ?
I tried : velocity.divergence == 0 but it doesn’t work.

Thank you,

Fabien




> Le 10 oct. 2018 à 17:43, Daniel Wheeler  a écrit :
> 
> Don't you still have a $\nabla . \vec{u} = 0$ equation though? It
> doesn't go away. That equation becomes like a constraint.
> 
> https://www.comsol.com/multiphysics/boussinesq-approximation
> 
> On Wed, Oct 10, 2018 at 5:58 AM fgendr01  wrote:
>> 
>> Hi Daniel,
>> Thank you for your answer.
>> I thank you for trying to solve my problem.
>> About my set of the equation here is my reasoning.
>> 
> 
> 
> -- 
> Daniel Wheeler
> ___
> 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 ]

Re: Boussinesq Equations

[Daniel Wheeler]

Don't you still have a $\nabla . \vec{u} = 0$ equation though? It
doesn't go away. That equation becomes like a constraint.

https://www.comsol.com/multiphysics/boussinesq-approximation

On Wed, Oct 10, 2018 at 5:58 AM fgendr01  wrote:
>
> Hi Daniel,
> Thank you for your answer.
> I thank you for trying to solve my problem.
> About my set of the equation here is my reasoning.
>


-- 
Daniel Wheeler
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

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 

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 :
>>
>> 

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 > > 

Re: Boussinesq Equations

[fgendr01]

Thanks a lot Jon.
I appreciate your help.
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) 
>  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("...")
> # 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("… ")
> 
> 
>> On Aug 14, 2018, at 11:22 AM, fgendr01  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) :
>> 
>> 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("...")
>> # 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("… »)
>> 
>> My error : all the input arrays must have same number of dimensions
>> 
>> I’d like to obtain something like this :
>> 
>> 
>> If someone can help me, I will be very happy.
>> 
>> Thanks a lot.
>> 
>> Fabien G.
>> ___
>> fipy mailing list
>> fipy@nist.gov
>> 

Re: Boussinesq Equations

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

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("...")
# 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("… ")


> On Aug 14, 2018, at 11:22 AM, fgendr01  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) :
> 
> 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("...")
> # 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("… »)
> 
> My error : all the input arrays must have same number of dimensions
> 
> I’d like to obtain something like this :
> 
> 
> 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 ]