Hello. I posted this previously, but no one was able to help me out so I am trying again. Here are the problem details...

I am modeling a channel that has three species in it. The first is nitrogen and is constant concentration (initially 80%)- so I set that as the primary fluid. The second is oxygen and has decreasing concentration (initially 19%) - so I am using sorsca.f user defined subroutine to simulate the removal of this scalar. The third is water vapor and has increasing concentration (initially 1%) - I am also using sorsca.f to simulate this mass source.

However, something must not be working correctly because the concentrations of the scalars are not changing. Here is the sorsca.f code:

C FOR O2...

IF(IS.EQ.1) THEN

S1P=(-0.032*1500)/(4*96487*0.0015) S2P=0

ENDIF

C FOR H2O...

IF(IS.EQ.2) THEN

S1P=(1.802*0.0000001*1500)/(2*96487*0.0015)

S2P=0

ENDIF

RETURN

END

These source terms should be specifying the addition or removal of the scalars 1 and 2 respectively. They are constant (kg/s*m^3) based on the governing equations which I am not varying anything within just yet.

Do I have something wrong in my code or is it in my model? Thanks in advance for you input and help.

- kevin

Let me Know ,What are the boundary conditions you gave?

There could be two reasons. One is the sorsca is not activated, the second is your IS is not the scalar number passed by the main code.

Here are my settings...

>Define Boundary Conditions

>Define Boundary Regions

>Inlet

>Default Conc. = 1, Velocity = 0.5 m/s

>Scalar Boundaries

>Inlet

>Scalar1(O2) Conc. = 1

>Scalar2(H2O) Conc. = 1

>Thermophysical Models and Props.

>Liquids and Gasses

>Molecular Props.

>Density = Ideal -f(T,P),

Mol Weight = 28.013 kg/kmol (N2)

>Turb. = off

>Temp. = off

> Bouancy = off

>Initialization

>P = 2 atm, T = 350 K, Values = constant,

code initialization exclude pressure

>Monitoring and Ref.

>P = 2 atm, T = 350 K

>Additional Scalars

>Initialization

>Scalar1(O2) Initial Conc. = .19,

values constant

>Scalar2(H2O) Initial Conc. = .01,

values constant

Alright here is the full code...I had merely shown the working portion of it before. What do you mean sorsca is not activated and that the IS is not the scalar number passed on by the main code? I do not follow. Thanks for the help.

- Kevin

C************************************************* ************************

SUBROUTINE SORSCA(S1P,S2P) C Source-term for scalar species C************************************************* ************************ C--------------------------------------------------------------------------* C STAR RELEASE 3.150 * C--------------------------------------------------------------------------*

INCLUDE 'comdb.inc'

COMMON/USR001/INTFLG(100)

INCLUDE 'usrdat.inc'

DIMENSION SCALAR(50)

EQUIVALENCE( UDAT12(001), ICTID )

EQUIVALENCE( UDAT03(001), CON )

EQUIVALENCE( UDAT03(002), TAU )

EQUIVALENCE( UDAT03(009), DUDX )

EQUIVALENCE( UDAT03(010), DVDX )

EQUIVALENCE( UDAT03(011), DWDX )

EQUIVALENCE( UDAT03(012), DUDY )

EQUIVALENCE( UDAT03(013), DVDY )

EQUIVALENCE( UDAT03(014), DWDY )

EQUIVALENCE( UDAT03(015), DUDZ )

EQUIVALENCE( UDAT03(016), DVDZ )

EQUIVALENCE( UDAT03(017), DWDZ )

EQUIVALENCE( UDAT03(019), VOLP )

EQUIVALENCE( UDAT04(001), CP )

EQUIVALENCE( UDAT04(002), DEN )

EQUIVALENCE( UDAT04(003), ED )

EQUIVALENCE( UDAT04(004), HP )

EQUIVALENCE( UDAT04(006), P )

EQUIVALENCE( UDAT04(008), TE )

EQUIVALENCE( UDAT04(009), SCALAR(01) )

EQUIVALENCE( UDAT04(059), U )

EQUIVALENCE( UDAT04(060), V )

EQUIVALENCE( UDAT04(061), W )

EQUIVALENCE( UDAT04(062), VISM )

EQUIVALENCE( UDAT04(063), VIST )

EQUIVALENCE( UDAT04(007), T )

EQUIVALENCE( UDAT04(067), X )

EQUIVALENCE( UDAT04(068), Y )

EQUIVALENCE( UDAT04(069), Z )

EQUIVALENCE( UDAT09(001), IS )

C FOR O2...

IF(IS.EQ.1) THEN

S1P=(-0.032*500)/(4*96487*0.0015)

S2P=0

ENDIF

C FOR H2O...

IF(IS.EQ.2) THEN

S1P=(1.802*0.0000001*500)/(2*96487*0.0015)

S2P=0

ENDIF

RETURN

END

You implied that both H2O and O2 both have inlet conc. 1.0, which can't be right. Unless you are also using fluinj, your sorsca will not act as a mass source or sink, just a source term in the scalar (mass fraction) equations. The sources over all the scalars should sum to zero, which they don't seem to. In any case, there is nothing with H in it which could react with the O2 to create H2O. Are you modelling a reacting flow, or the injection of water vapour into air? If it's the latter, you might do better to discard sorsca and just use fluinj to inject H2O.

Richard,

I had originally set the concentrations of O2 and H2O to .19 and .01, respectively, but the end result from that was a constant concentration equal to those initial concentrations throughout the channels. Something goofy seemed to be goin on in that the initial concentration of the fluid (in my case N2) seemed to have to be equal to that of one of the initial concentrations of the scalars. I would attempt to change the value and "apply" but it would change it on its own back to whatever value the scalar was at. See below:

>Define Boundary Conditions

>Define Boundary Regions

>Inlet

>Default Conc. = [***ends up here as well?], Velocity = 0.5 m/s

>Scalar Boundaries

>Inlet

>Scalar1(O2) Conc. = [*whatever I put here or below]

>Scalar2(H2O) Conc. = [**which ever is greater?]

So I must not understand how to initialize this mixture gas flow. There seems to me to be two places in the model settings to specify concentrations. Under scalar boundaries, as well as under additional scalars initialization. Do you know how to set this up properly?

Regarding the chemical reaction, I am modeling only the mass transfer that is associated with a reaction occuring outside my channels that is effectively drawing out O2 and introducing H2O. Hence the consumption of oxygen and the addition of water vapor.

Regarding the subroutine usage, I was instructed that I should only have to use sorsca.f to accomplish this model. You mentioned:

"Unless you are also using fluinj, your sorsca will not act as a mass source or sink, just a source term in the scalar (mass fraction) equations. The sources over all the scalars should sum to zero, which they don't seem to."

I do not follow. When you have more time, perhaps you could expound a little.

Much thanks for the help, sincerely.

- Kevin

I don't know why the initial and boundary conditions are mixed up. You should just be able to specify the concentrations of O2 and H2O at the inlet and in the initial field, and the N2 (the background fluid) is inferred from that.

My point about sorsca is that it is only for redistributing mass between the scalars, not introducing it. So if you had a chemical reaction where X kg/s of fuel was reacting with Y kg/s of oxidant to produce Z kg/s of product, you would have three scalars (fuel, oxidant and products) with three sources (X, Y and Z) which are related by X + Y = Z (conservation of mass).

On the other hand, you are extracting O2 and replacing it with H2O, so there is not conservation of mass. Does the extraction take place within the domain or at the boundaries? If it's within the domain I think you will need a sorsca and a fluinj; if it's on the boundaries you should be able to manage with just boundary conditions.

ok, please just add a write statement in this code to see if it is called or not.

I know the code is being called and utilized because I have run trials with different calues for S1P and S2P and the results have shown variations in pressure that seem reasonable.

Richard,

This reaction effectively takes place throughout the domain - I need to simulate removal and injection in kg/m^3*s - so in every cell. I believe fluinj.f can only be used for the primary fluid - how would I use it for the additional scalars that I wish to increase/decrease?

Regarding sorsca.f the user guides says "it specifies additional source terms for the scalar variable equations" So I thought that I could use this to model my removal/injection. You are saying that this is not the case? What then can I do to simulate this model? Is there some other subroutine that will accomplish this? Perhaps scalfn.f? That is the only other that I find that is even close to relevant. Thanks again.

- Kevin

You can inject/remove scalars with fluinj: that is what its scinj argument specifies. The sample coding shows how to inject scalar(1) at 0.1 kg/m3/s.

In your case you need to add one scalar and remove another, so a combination of fluinj and sorsca is needed. In the fluinj you need to specify fluxi to be the net injection rate S(H2O)+S(O2) and scinj(H2O) to be 1.0. This accounts for the "extra" mass, and then you use sorsca to redistribute the remaining mass, so for H2O you need a source term of -S(O2) and for O2 just S(O2).

So your net source for H2O will be S(H2O)+S(O2)-S(O2)=S(H2O) and for O2 just S(O2), exactly what you wanted. It would be much easier if sorsca could do all this in one but, as I say, it can only redistribute mass, so you need the fluinj to inject/remove mass.

Or at least I think that's how you do. Give it a try. You can use the conservation checks in the info file to make sure it's behaving.

Richard,

I reformatted my subroutines as you suggested using both fluinj and sorsca, but I am still getting the same results...0.19 concentration for O2 across the whole model and 0.01 concentration for H2O across the whole model. Any other ideas? Model settings perhaps? However, in the info file I am getting this:

*** WARNING #052 *** INITIAL RESIDUAL BELOW ROUND-OFF ERROR LIMIT SOLUTION ,EQ: SC1

*** WARNING #052 *** INITIAL RESIDUAL BELOW ROUND-OFF ERROR LIMIT SOLUTION ,EQ: SC2

Could this have anything to do with it? You also mentioned conservation checks...what part of the info file relates to that? Gracias!

- Kevin

The warning about initial residual being below the round-off error would occur if the boundary conditions and initial conditions (or current field) were in perfect agreement, and there were no source terms. It's the solver's way of saying it can't converge the field any more than it is already. So yjl's question about whether you have actually activated and are using your subroutines is valid. It's not enough just to compile them in, you have to activate them in prostar too.

You can activate the conservation checks in the Output Controls -> Monitor Numeric Behavior panel in prostar. The numbers come out every iteration in the info file. You need to look at MSOR for each of your scalars, and make sure that they are what you want them to be (kg/s).

Hey...

Yea I goofed that up! I fixed the initial residual being below the round-off error problem. I had to activate the subroutines like you said, I had previously saved a copy of my model and started on a new one and then went back and forgot to set up the subroutines settings properly. Brain fart. Thanks for pointing that out (again) So now...

This is from my run file:

-> USER FORTRAN CODING USED

DEFINED SOURCE TERM FOR SPECIES O2 IN FLUID 1

DEFINED SOURCE TERM FOR SPECIES STEAM_H2O IN FLUID 1

FLUID INJECTION/SUCTION

Also, I have been able to produce varying output results (relative pressure for instance) by changing the values used in my subroutines - so I do believe that they are now being activated, compiled, and utilized in the solution.

However, the model is still not functioning as I would expect it to. The concentrations of the scalars still seems to be remaining constant.

Also, these errors are now appearing in the info file: (This from iteration 1, but they are present throughout)

*** WARNING #042 *** PROCESS NOT FULLY CONVERGED -- EQ.,NIT,RESI: SC1 101 0.36594E+01

*** WARNING #045 *** SOLUTION DIVERGES; RESCUE NOT ATTEMPTED; RETURNING OLD SOLUTION

*** WARNING #042 *** PROCESS NOT FULLY CONVERGED -- EQ.,NIT,RESI: SC2 101 0.33123E+03

*** WARNING #045 *** SOLUTION DIVERGES; RESCUE NOT ATTEMPTED; RETURNING OLD SOLUTION

Something else must be wrong I am assuming. Any ideas anyone?

- Kevin

Good progress ... The scalar fields are not changing because of the warnings you quote, especially the "returning old solution" bit. Are the MSOR values in the info file what you expect? They should be exactly right, as I think star just sums up what you specify, no massaging before printing it out. Does reducing the under-relaxation factor for the scalars help? (I assume it's a steady-state sim).

The MSOR (species mass sources) values in the info file are:

-5.0300E-08 for O2 and

5.6608E-08 for H20 (in kg/s)

From my subroutines I have:

FLUXI=0.00348745

SCINJ(2)=1.0

from fluinj.f and:

C FOR O2...

IF(IS.EQ.1) THEN

S1P=-0.02763755

S2P=0

ENDIF

C FOR H2O...

IF(IS.EQ.2) THEN

S1P=0.02763755

S2P=0

ENDIF

from sorsca.f.

I am not sure if those MSOR values are correct or not because I thought the mass source terms were in kg/m^3*s. Has star converted them to a per cell basis using the volume of my cell/s (aka kg/s)?

The volume of my cells are all 2.7e-11 m^3 for each.

O2 Consumption = 0.02763755 kg/m^3*s

H20 Production = 0.031125 kg/m^3*s

If so, I don't think that they are mathcing up right. Why would this be? What should the relation between the MSOR values and my numbers above (consumption and production) be?

I tried reducing the under-relaxation factor for the scalars from 1 to 0.7 for both scalars and that got rid of the error messages I spoke of last.

I was actually able to see some sort of scalar concentration gradient for both of the scalars in question. However, the gradients only last for 1/40 of the length of the model channel and then they are constant for the remainder (majority!) of the model.

Something is working a little better now...but something must still be ary. Let me know what you make of this...

Thanks again.

- Kevin

The under-relaxation factors sound to have done the trick. Since you're applying the fluinj and sorsca in all the cells, the MSORs (kg/s) will be (specified source kg/m3/s)*(Volume of domain m3). Hopefully that matches up. The last thing to make sure of is that you converge the equations completely - the terms in the conservation check should balance, i.e. (Inlet flux of scalar) + MSOR = (Outlet Flux of scalar).

Richard,

I checked the numbers for MSOR and they are close.

Info file has -5.0300E-08 for O2 and I calculated -5.03755E-08.

Info file has 5.6608E-08 for H2O and I calculated 5.67322E-08.

Are these close enough or is this cause for concern?

Regarding your check for converging the equations completely, I am not getting good results. For O2 I have:

4.2468E-07 inlet and 4.2491E-07 outlet from the info file. These have a difference of -2.30E-10.

For H2O I have:

2.2344E-08 inlet and 2.3487E-08 outlet from the info file. These have a difference of 1.143E-08.

These numbers are nowhere near where you say they should be (equal to MSOR) Any ideas? Let me know...

- Kevin

Sorry, I made one error...For O2 I have:

4.2468E-07 inlet and 4.2491E-07 outlet from the info file. These have a difference of 2.30E-10.

Which does not make sense that it is not negative, being that oxygen is being CONSUMED.

and the plot thickens...