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 ]
