0

I have tried to simplify my set of PDEs: Set of equations.

In this new set, p and m are my variables, where all the other values are constant. The m and p values in fc²m|m|/2DA²p are treated ass mean values on each cell. Furthermore I want to fix the inlet mass flow (m(0, t)) and the outlet pressure (p(L,t)), so I have defined it as constrains.

I have two questions: my code seems to diverge, for each time step it either gives me a very high value of m in the pipeline, while the next one gives me a very negative value and so on. Also, it seems that simply define inlet mass flow and outlet pressure gives an unreal solution.

Any thoughts and suggestions? Thank you in advance, dear community, once again.

code below:

import fipy as fi
import matplotlib.pyplot as plt
import numpy as np
import scipy as sci
from matplotlib import animation
from fipy.variables.cellVariable import CellVariable

#1. Domain
L = 101
dx = .1
mesh = fi.Grid1D(nx = L, dx=dx)

#2. Parameters values (Arbitrary) 
Lambda = 0.0001 # Friction factor
D = .4 # Pipe diameter
T = 350 # Gas Temperature
c = 380
A = 3.14 * (D **2) / 4

#3. Variables
## Pressure and mass flow
p = fi.CellVariable(name="pressure", 
                    hasOld=True, 
                    mesh=mesh, 
                    value=0.)
p.setValue(10000000.) # This is the initial condition

m = fi.CellVariable(name="mass flow", 
                    hasOld=True, 
                    mesh=mesh, 
                    value=0.)
m.setValue(300.) # This is the initial condition

#4. Boundary conditions
p.constrain (2000000., where = mesh.facesRight)
p.constrain (15., where = mesh.facesLeft)
m.constrain (500., where = mesh.facesLeft)
m.constrain (500., where = mesh.facesRight)

#5. PDE
eq1 = fi.TransientTerm(var=p) == fi.ConvectionTerm(coeff = [-(c**2) / A], var=m)
eq2 =  fi.TransientTerm(var=m) == fi.ConvectionTerm(coeff = [-A], var=p) - (Lambda*(c**2)*m.cellVolumeAverage*abs(m.cellVolumeAverage)/(2*D*A*p.cellVolumeAverage))

eqn = (eq1 & eq2)

timeStepDuration = .001
steps = 50

#This make a chart of the results
for step in range(steps):
    plt.plot(np.linspace(0, L, nx), m.value)
    #plt.xlim(0,1.1)
    #plt.ylim(0, 20)
    plt.show()
    
    eqn.sweep(dt=timeStepDuration)
    p.updateOld()
    m.updateOld()
petesing
  • 13
  • 3
  • I don't really know, as fluid flow really not my bailiwick, but `eq1` is purely hyperbolic and `eq2` is largely(?) so. FiPy is bad at hyperbolic equations and something like [CLAWPACK](https://www.clawpack.org/) might work better. On the other hand, the SIMPLE algorithm demonstrated in our [Stokes](https://www.ctcms.nist.gov/fipy/examples/flow/generated/examples.flow.stokesCavity.html) and [reactive wetting](https://www.ctcms.nist.gov/fipy/examples/reactiveWetting/generated/examples.reactiveWetting.liquidVapor1D.html) examples, is an approach that FiPy can do a reasonable job with. – jeguyer Apr 23 '21 at 19:24
  • Simplifying the equations this way may have actually made them harder for FiPy to solve. – jeguyer Apr 23 '21 at 19:25
  • Thank you. Do you think it would be possible to adapt the Stokes example for a 1D domain? I am having some troubles ("The coefficient must be rank 0 for a rank 0 solution variable.", for instace) – petesing Apr 25 '21 at 00:22
  • Do not confuse the rank with the dimension. A 1D domain can have rank-0 (scalar) fields, rank-1 (vector) fields, rank-2 (tensor) fields, ... Write your equations in general form (gradients, divergences, etc., instead of partial derivatives wrt x), and then solve them on whatever domain you wish. – jeguyer Apr 25 '21 at 01:56
  • I am sorry, I don't think I get what you are saying, sir. I know you have been helping me a lot, but I will post what I've done to adapt it to a 1D domain, since the example in FiPy's website shows a 2D case – petesing Apr 25 '21 at 19:15

0 Answers0