Rocket nozzles must be short and light but make the flow as parallel as possible at the exit. Further demands exist in the atmosphere. The eggcup shape answers that but needs numerical optimisation.
The thermodynamic state itself isn't algebraic, and a single point takes an effort to a PC. Dozens of chemical equilibria shift with T and P, notably 2CO + O2 <-> 2CO2. Condensation and solidification over the expansion, typically of Al2O3, worsen everything. The heat capacity depends heavily on the temperature. And limited heat exchanges between droplets and the gas would be useful to model. One common assumption is that the chemical composition shifts up to the throat but is frozen downstream as the expansion gets snappy.
Even in the axisymmetric shape, combining the flow computation with the thermodynamic state in finite elements would be a very heavy computation. An improvement, supposedly banal, is to pre-compute a table of the speed modulus versus one variable like T, P or any convenient one, and use some interpolation. The computation load remains heavy.
My proposal, and as usual I didn't check how common it is:
Don't solve finite elements for the speed vector, as this demands to find simultaneously the thermodynamic state.
Solve for the product of density and speed modulus. This conservation of mass is independent of the thermodynamic state. Its divergence is zero. It is even the gradient of a scalar harmonic potential if useful.
When this product is known everywhere in the nozzle, the thermodynamic state can be injected from the inverted precomputed table that gives P, T, composition an so on as a function of the modulus of the product. Solving the finite elements gets easy for a PC.
Solving for the product of speed and density applies also to other non-trivial flows, like the expansion of saturated steam. It's already useful for hypersonic flows of gases.
If any necessary, the 2D nozzle problem can become a 1D problem by adding a set of virtual flow sources (here tori) outside the nozzle walls and solving the intensity of the vitual sources that zeroes the flow through every wall section. One source behind each wall section, close enough for a decently conditioned matrix, far enough for even flow.
The huge but hollow matrix of finite elements becomes a small but full one. If solving by iterations (less useful nowadays), the small matrix can be approximated by a band-diagonal one.
A 3D problem also becomes a 2D one, with point sources behind area elements like triangles. The mass flow from each source is just R-2. The full matrix can be approximated by a multidiagonal one. This is "known" as I did it three decades ago, for a non-axisymmetric case, and other people may have done it before.
Virtual flow sources apply where turbulence is unimportant and the virtual flow behind a wall doesn't matter, for instance at a turbine blade or a fast aircraft.
Marc Schaefer, aka Enthalpy