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 <fabien.gend...@univ-lr.fr> 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
>> <thibault.bridellel...@gmail.com> 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 <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...
>>
>> <PastedGraphic-5.png>
>>
>> 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("<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.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("<Enter to finish>… ")
>>
>>
>> Sorry for all my messages...
>>
>> Fabien
>>
>>
>>> Le 17 sept. 2018 à 21:52, Thibault Bridel-Bertomeu
>>> <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 <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 à 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 <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> a écrit :
>>>>>
>>>>> 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)
>>>>>> <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>
>>>>>>> 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
>>>>>>> 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 ]
>>>>>
>>>>
>>>> _______________________________________________
>>>> fipy mailing list
>>>> fipy@nist.gov
>>>> http://www.ctcms.nist.gov/fipy
>>>> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>>>> <PastedGraphic-2.png><PastedGraphic-1.png><PastedGraphic-2.png>_______________________________________________
>>>>
>>>> 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 ]
>>> _______________________________________________
>>> 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 ]
>> _______________________________________________
>> 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 ]
_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
