|
| Tuple[float, np.ndarray] | _element_geometry (np.ndarray v) |
| float | _triple_integral_factor (int i, int j, int k) |
| np.ndarray | _mass_r_element (float area, np.ndarray r_vertices) |
| np.ndarray | _stiffness_r_element (float area, np.ndarray grad, np.ndarray r_vertices) |
| np.ndarray | _azimuthal_element (float area, np.ndarray r_vertices) |
| Tuple[sparse.csr_matrix, sparse.csr_matrix, sparse.csr_matrix] | assemble_passive (MeridionalMesh mesh, np.ndarray c_sound, int m) |
| np.ndarray | _boundary_mass_r_edge (float r_a, float r_b) |
| sparse.csr_matrix | assemble_boundary_admittance (MeridionalMesh mesh, np.ndarray boundary_nodes, float coefficient) |
| np.ndarray | point_sampling_vector (MeridionalMesh mesh, float x0, float r0) |
| np.ndarray | disk_load_vector (MeridionalMesh mesh, np.ndarray column_nodes) |
FEM assembly of the passive (no flame, no boundary impedance) axisymmetric
thermoacoustic operator: the ``K``, ``Km``, ``M2`` matrices of Eq. (10)/(11) in
``docs/stability/thermoacoustic_global_stability_paper.md``.
K = int_Omega c^2 (grad p . grad phi) r dOmega stiffness
Km = int_Omega (c^2 / r) p phi dOmega azimuthal term (x m^2)
M2 = int_Omega p phi r dOmega mass
With these, the passive discrete eigenproblem (flame off, rigid walls: C=0, F=0) is
(K + m^2 Km) p = -lambda^2 M2 p
a *linear* generalized eigenvalue problem in mu = -lambda^2 (real, >= 0 for a lossless
rigid-wall cavity, since K, Km, M2 are all real symmetric positive-(semi)definite) —
see ``eigen/passive.py``.
Why K and M2 have exact closed forms but Km does not
-----------------------------------------------------
On a P1 (linear) triangle, the shape functions N_i are the barycentric coordinates
L_i, and grad(N_i) is *constant* over the element. The radial coordinate r(x) is
itself a linear (degree-1) function of position, since it is a coordinate. So:
- K's integrand, c^2 * (grad N_i . grad N_j) * r, is (const) * (linear in r) ->
integrable exactly by the 1-point centroid rule: int_T r dA = Area * r_centroid.
- M2's integrand, N_i * N_j * r, is a *cubic* polynomial (degree 2 from N_i N_j,
degree 1 from r). Triangle integrals of monomials in barycentric coordinates have
a standard closed form (see ``_mass_r_element`` below), so M2 is also exact.
- Km's integrand, N_i * N_j / r, is *not* polynomial (r appears in the denominator)
so no closed form exists; it is evaluated with a 6-point Gauss quadrature rule
(exact for polynomials up to degree 4), which is far more resolution than the
(non-polynomial) 1/r factor needs for engineering accuracy, and cheap since it
only touches the m>=1 (transverse) matrix.
Mean-flow fields (``c``, i.e. sound speed) are supplied per-node and averaged over the
three vertices of an element to get a single c^2 per element — piecewise-constant per
element, consistent with the P1 discretization.
| Tuple[sparse.csr_matrix, sparse.csr_matrix, sparse.csr_matrix] engine.stability_hifi.acoustics.assembly.assemble_passive |
( |
MeridionalMesh | mesh, |
|
|
np.ndarray | c_sound, |
|
|
int
| m ) |
Assemble the passive (flame-off, rigid-wall) K, Km, M2 matrices (Eq. 10, 11).
``c_sound`` is the sound-speed field [m/s], given EITHER per-node (shape
``(n_nodes,)`` — the mean flow varies smoothly, e.g. V1's uniform field or a real
CFD/parametric temperature profile; each element's c^2 is the mean of its 3 vertex
values, i.e. the P1 linear interpolant evaluated at the centroid) OR per-element
(shape ``(n_tri,)`` — the material is piecewise-constant across a genuine interface
that a shared node cannot represent, e.g. V2's temperature-jump duct; each element's
c^2 is used directly, with no averaging across the interface). See
``_element_c_squared`` above for why this distinction matters physically.
``m`` is the azimuthal wavenumber; ``Km`` is returned *without* the m^2 factor
(caller multiplies when assembling ``K + m^2 * Km``), so the m=0 matrix (which has
no 1/r term at all) can be skipped cheaply by callers that only need
longitudinal/radial modes.
Returns ``(K, Km, M2)`` as real symmetric sparse CSR matrices of size n_nodes^2.
Axis-regularity Dirichlet elimination for m>=1 is the caller's job (``eigen/passive.py``),
since it depends on how the reduced system is solved, not on the assembly.
Definition at line 145 of file assembly.py.
References _azimuthal_element(), _element_geometry(), _mass_r_element(), and _stiffness_r_element().
| np.ndarray engine.stability_hifi.acoustics.assembly.disk_load_vector |
( |
MeridionalMesh | mesh, |
|
|
np.ndarray | column_nodes ) |
The g_k vector for a compact-in-x, UNIFORM-across-the-cross-section flame at the
axial column given by ``column_nodes`` (node indices there, SORTED BY RADIUS —
exactly what ``MeridionalMesh.nodes_at_x`` returns; the flame location must
therefore fall on an exact mesh column, e.g. by building the mesh with
``two_zone_duct_mesh`` even when there is no real material jump, purely to get an
exact shared column at the flame's axial location).
``g_i = int_r N_i(x_f, r) * r dr`` along that column: same r dOmega measure as K,
M2, C (see the module-level note above for why that consistency matters), computed
per edge via the same barycentric-monomial exact-integral trick as
``_boundary_mass_r_edge``, but for a LOAD VECTOR (int N_i * 1 * r dr) rather than a
mass MATRIX (int N_i N_j r dr): for an edge from r_a to r_b, the exact two-node
contribution is ``L*[r_a/3+r_b/6, r_a/6+r_b/3]`` (L = r_b - r_a).
Sanity check used while developing this: summed over the whole column, ``g``'s
total is exactly ``int_0^R r dr = R^2/2`` — reproduced below to machine precision,
confirming the exact-integral bookkeeping is right.
Definition at line 309 of file assembly.py.
1.15.0