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 > <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... > > <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 <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 <fabien.gend...@univ-lr.fr >> <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 à 21:29, Thibault Bridel-Bertomeu >>> <thibault.bridellel...@gmail.com <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 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 >>> <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("<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 >>>> <mailto: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 <mailto: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 >>>>>> <mailto: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 <mailto:fipy@nist.gov> >>>>>> http://www.ctcms.nist.gov/fipy <http://www.ctcms.nist.gov/fipy> >>>>>> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy >>>>>> <https://email.nist.gov/mailman/listinfo/fipy> ] >>>>> >>>>> >>>>> _______________________________________________ >>>>> fipy mailing list >>>>> fipy@nist.gov <mailto:fipy@nist.gov> >>>>> http://www.ctcms.nist.gov/fipy <http://www.ctcms.nist.gov/fipy> >>>>> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy >>>>> <https://email.nist.gov/mailman/listinfo/fipy> ] >>>> >>> >>> _______________________________________________ >>> fipy mailing list >>> fipy@nist.gov <mailto:fipy@nist.gov> >>> http://www.ctcms.nist.gov/fipy <http://www.ctcms.nist.gov/fipy> >>> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy >>> <https://email.nist.gov/mailman/listinfo/fipy> ] >> >>> <PastedGraphic-2.png><PastedGraphic-1.png><PastedGraphic-2.png>_______________________________________________ >> >>> >>> fipy mailing list >>> fipy@nist.gov <mailto:fipy@nist.gov> >>> http://www.ctcms.nist.gov/fipy <http://www.ctcms.nist.gov/fipy> >>> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy >>> <https://email.nist.gov/mailman/listinfo/fipy> ] >> >> _______________________________________________ >> fipy mailing list >> fipy@nist.gov <mailto:fipy@nist.gov> >> http://www.ctcms.nist.gov/fipy <http://www.ctcms.nist.gov/fipy> >> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy >> <https://email.nist.gov/mailman/listinfo/fipy> ] >> _______________________________________________ >> fipy mailing list >> fipy@nist.gov <mailto:fipy@nist.gov> >> http://www.ctcms.nist.gov/fipy <http://www.ctcms.nist.gov/fipy> >> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy >> <https://email.nist.gov/mailman/listinfo/fipy> ] > > _______________________________________________ > fipy mailing list > fipy@nist.gov <mailto:fipy@nist.gov> > http://www.ctcms.nist.gov/fipy <http://www.ctcms.nist.gov/fipy> > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy > <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 ]
