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 > <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 http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
