Re: Boussinesq Equations
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
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
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
Hi Fabien, I looked over your problem and I haven't been able to debug it. I tried various combinations of constraints and I'm not sure why it's becoming nonphysical (dT >1 or dT < -1). It would require getting into the code and looking at the numbers in the offending cells. One thing that worries me is that the equations don't have pressure or the continuity equation, which seems strange. I would be rather suspicious of this set of equations because there is no strong coupling between the x and z momentum and also because there they are pure convection without diffusion, which requires special schemes to solve correctly and which FiPy doesn't handle very well. Cheers, Daniel On Thu, Sep 20, 2018 at 8:09 AM fgendr01 wrote: > > Hi Jon, Thibaut and the community ! > Thanks for your help. > I change my parameters back and my problem seems to be stable now. > But I think, there’s still a problem. > On my simulation, there is a continuity problem at the top left and bottom > right. I do not understand why a bubble of average temperature appears! It’s > a bit strange ! > I tested with several boundary conditions but it does not change anything. > Yet I impose a temperature of 1 to the left and -1 to the right and adiabatic > walls at the top and bottom. Do you have an idea ? > Here is some screenshots : -- 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
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 :
>>>
>>&g
Re: Boussinesq Equations
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 <mailto:fabien.gend...@univ-lr.fr>> 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 > <mailto:fabien.gend...@univ-lr.fr>> 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 à 2
Re: Boussinesq Equations
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
>>
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)
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 <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("...")
> # 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.arithmeticFa
Re: Boussinesq Equations
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=
Re: Boussinesq Equations
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 ]
