| Author | Message | |||
     Roman Fursenko (roman) New member Username: roman Post Number: 1 Registered: 02-2009 |
Dear Mr. Nelson, I try to solve the problem of 1D flow through the heated channel. The model consists of mass conservation and momentum equations. I set mass flow rate at the inlet (x=-Lx), and dx(u)=0 at the outlet (x=Lx). I have also added the “artificial viscosity” in mass conservation equation in order to smooth solution. The following script was failed at time about 0.099. And solution oscillations at earlier time were observed. Could you, please, give any suggestions about this problem. Thanks in advance, Roman. !SCRIPT: TITLE 'New Problem' COORDINATES cartesian1 VARIABLES p, u DEFINITIONS Tg0=300 Tgb=2000 v0=10.0 Lx=3 Tg=Tg0 +(Tgb-Tg0)*exp(-5*(x)**2) wM=29 R=8.31 p0=0.1013 ro=p*wM/(R*Tg) ro0=p0*wM/(R*Tg0) mu=1E-1 muro=1E-5 Q=ro0*v0 !Mass flow rate INITIAL VALUES u=v0 p=p0 EQUATIONS p: dt(ro)+dx(u*ro)=muro*dxx(ro) u: dt(u)+u*dx(u)=mu*dxx(u)-dx(p)/ro boundaries REGION 1 START(-Lx) point value(u)=Q/ro line TO(Lx) point natural(u)=0 TIME 0 TO 20 MONITORS PLOTS for t = 0 by 0.01 to 1.0 by 0.1 to endtime GRID(x) ELEVATION(ro) FROM (-Lx) to (Lx) ELEVATION(u) FROM (-Lx) to (Lx) ELEVATION(p) FROM (-Lx) to (Lx) ELEVATION(ro*u) FROM (-Lx) to (Lx) ELEVATION(Tg) FROM (-Lx) to (Lx) END | |||
| Robert G. Nelson (rgnelson) Moderator Username: rgnelson Post Number: 1217 Registered: 06-2003 |
You start out with uniform flow velocity and density, then at time zero the temperature of the gas at the center of the tube is instantaneously raised to 2000 degrees. This is a very fast detonation, and sends strong shock waves through the gas. The computation mesh is not dense enough to model these sharp fronts, and oscillation results. So, the first step is to increase the mesh density. SELECT NGRID=200 seems to do the job. The second step is to slow the detonation speed. If this were a tube with a heater, then you could not heat the gas instantly. It would take some time. So ramp the central temperature over a time long compared to the sound transit time through the pipe. This will ameliorate some of the shocks. Also, I am curious how you plan to cool the gas again after it pases through the central section. Cooling coils wrapped next to the heater coil? If the Gaussian temperature shape was intended as a heat input, then you will need a temperature equation to track the heating of the gas. This would allow the gas to exit the pipe at the heated temperature. The attached script does a few of the things I mentioned.
| |||
| Roman Fursenko (roman) New member Username: roman Post Number: 2 Registered: 02-2009 |
Thank you very much for the prompt and helpful reply! I agree with you about the necessity of the equation on temperature. The script I have attached in the first message was just a test, and so on I have simplified the problem statment. Thanks again. | |||
| Roman Fursenko (roman) Junior Member Username: roman Post Number: 3 Registered: 02-2009 |
Dear Mr. Nelson, Need your help again. The problem statement is the same but I have changed the constants on the more reasonable. Due to great value of the initial pressure (~10E6) I wrote pressure=p0+p. Where p0 is the initial pressure variation p should be determinated. As the first step (see script below) I have set the temperature Tg equal to the constant value Tg0 and started from the constant initial conditions representing exact stationary solution. Even in this trivial case I have got the oscillations of the solution. Would you be so kind to give any suggestions about this. In the case of non-constant temperature the oscillations amplitude are much greater and program is failed. Even significant increasing of the ngrid parameter does not help. Thanks in advance. !SCRIP TITLE 'New Problem' COORDINATES cartesian1 select ngrid=200 VARIABLES ro, u DEFINITIONS Tg0=300 Tgb=2000 v0=10.0 Lx=3 Tg=Tg0 !Tg0 +uramp(t,t-0.1)*(Tgb-Tg0)*exp(-5*(x)^2) wM=29 R=8.31E7 p0=1.013E6 p=ro*Tg*R/wM-p0 ro0=p0*wM/(R*Tg0) mu=1E-1 muro=1E-5 Q=ro0*v0 !Mass flow rate INITIAL VALUES u=v0+0*(x+Lx) ro=ro0 EQUATIONS ro: dt(ro)+dx(u*ro)=muro*dxx(ro) u: dt(u)+u*dx(u)=mu*dxx(u)-dx(p)/ro boundaries REGION 1 START(-Lx) point value(u)=Q/ro point natural(u)=0 line TO(Lx) point natural(u)=0 TIME 0 TO 20 MONITORS for cycle=10 GRID(x) ELEVATION(ro) FROM (-Lx) to (Lx) ELEVATION(u) FROM (-Lx) to (Lx) ELEVATION(p) FROM (-Lx) to (Lx) ELEVATION(ro*u) FROM (-Lx) to (Lx) ELEVATION(Tg) FROM (-Lx) to (Lx) END |
