why are the results different from moist air domain(MA) using only dry air and gas domain(G)? - simscape

Hello, I used Simscape to model with the constant volume chamber, local restriction, and controlled heat flow rate source blocks. I used the semi-ideal gas and now I'm trying to switch to the moist air domain to include gases with different properties. I modeled it the same way as before and set the humidity ratio to 0 to ensure there is no water vapor. However, the results are significantly different. Are there any conditions I might have overlooked?

[Yifeng Tang]

Could you share the simple model to demonstrate your observations? My experience is that it's usually some parameters, likely default ones, that are not consistent between the two models.

[정수 이]

I appreciate the quick response. I'll simply send over the two modeling files. I used default properties for water vapor. While attempting to create a simple model, I found that using the dry air properties from the MA domain in the Gas domain yielded the same values. However, when I input air properties ranging from 200K to 35000K into both models simultaneously, they still produce different results.

Answer 1

[Yifeng Tang]

Hello,

Given the range of temperature you'd like to implement, I believe you are seeing the difference in the way the specific heat is used & derived. Let me explain in a bit of details.

In the Gas domain when using the semiperfect specification for gas properties, you'll be asked to provide, as functions of temperature, the specific enthalpy, specific heat (cp), dynamic viscosity, and thermal conductivity. In the MA domain, you'll be asked for the same list, except for the specific heat. hmm ... that's because in the MA domain, the cp is "derived" from the specific enthalpy array. The MA domain doesn't store the actual cp array, and instead it uses a 3rd order polynomial function of temperature to fit to the specific enthalpy array and take the analytical derivative. The resulting coefficient are carried within the domain properties. The code for this fit may be found in the source code of "Moist Air Properties (MA)", line 663-672 and 697-701.

Now, with a very large temperature range, a 3rd order polynomial may be insufficient for the h array. So the resulting cp array can be different from what you have for the Gad domain. I believe this is most likely the source of error. The following script shows the difference in the polynomial fit vs the original data.

T_TLU = [200 300 400 500 600 700 800 900 1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 7500 8000 8500 9000 9500 10000 11000 12000 13000 14000 15000 16000 17000 18000 19000 20000 21000 22000 23000 24000 25000 26000 27000 27800 29000 30000 31000 32000 33000 34000 35000];
h_TLU = [-103000 -2630 98200 200000 304000 411000 519000 630000 743000 1340000 1980000 2710000 3770000 5490000 7370000 8710000 9900000 11600000 14400000 18900000 25100000 31900000 37300000 40800000 43400000 45600000 47900000 53900000 63100000 76700000 94900000 116000000 136000000 152000000 164000000 173000000 179000000 185000000 190000000 196000000 204000000 215000000 229000000 248000000 272000000 300000000 329000000 358000000 386000000 410000000 432000000 450000000];
cp_TLU = [996 1000 1010 1030 1050 1070 1100 1120 1140 1230 1340 1670 2720 3990 3210 2360 2710 4260 7100 10900 13700 12600 8760 5790 4520 4390 4920 7300 11100 15800 20000 21200 18500 14000 9940 7280 5830 5290 5570 6730 8960 12400 16800 21600 25900 28700 29500 28400 26000 22800 19700 17100];

p = polyfit(T_TLU,h_TLU,3);

h_Fit = p(1)*T_TLU.^3 + p(2)*T_TLU.^2 + p(3)*T_TLU.^1 + p(4);
cp_Fit = p(1)*T_TLU.^2*3 + p(2)*T_TLU.^1*2 + p(3);

figure(101); clf; hold on; box on;
plot(T_TLU,h_TLU)
plot(T_TLU,h_Fit)
xlabel("Temperature (K)")
ylabel("Specific Enthalpy (J/kg)")
legend("Reference Array","Polynomial Fit",location="northwest")

figure(102); clf; hold on; box on;
plot(T_TLU,cp_TLU)
plot(T_TLU,cp_Fit)
xlabel("Temperature (K)")
ylabel("Specific Heat (J/K/kg)")
legend("Reference Array","Polynomial Fit",location="northwest")

If using MA for your future work is critical, please consider either of the following:

(1) wait for R2024b release. I've heard about a different way to generate the cp array that will be more robust. This does requires that your license is up to date by the time R2024b pre-release or formal release become available (in a few month?).

(2) if your license is up to date now, you may reach out to your MathWorks sales rep and mention this post. We can share a prototype of the MA domain that uses the provided enthalpy arrary to calculate cp using a differencing scheme.