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
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 : >>> >>> Hi again Fabien, >>> >>> Well, I would say that
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 à 21:29, Thibault Bridel-Bertomeu >>> mailt
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 >> timeStepDuration = 0.1*dL**2/(2*DT) >> ste
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.arithmeticFaceValue >velocity[..., mesh.exteriorFac
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=velocity) == >> -DiffusionTerm(var=dT, coeff=D)) >> eq = eqX & eqZ & eqT >&g
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 ]