Here's how to use the `NonCommutativeAlgebra` functionality to implement an algebra with a finite number of generators and relations, illustrated by implementing the oscillator algebra.  This will address Questions 1 and 2 in the OP (I'll wrap up at the end explaining why this won't work for a countably infinite number of generators or relations).

First define a set of generators, which are the raising and lowering operators $\hat{a}$ and $\hat{a}^{\dagger}$ for the algebra:

    gens = {a, aD};

The defining relation for the algebra is that the commutator is the identity, i.e., $\hat{a}\hat{a}^{\dagger} - \hat{a}^{\dagger}\hat{a} -1 = 0,$ so we define a set of relations as

    reps = {Commutator[a,aD] - 1};

This is a noncommutative Groebner basis:

    NonCommutativeGroebnerBasis[rels, gens]
    (* {-1 + a**aD - aD**a} *)

so it can be used to reduce noncommutative polynomials in the generators to a canonical form.  Since I have specified the generators in the order `{a, aD}`, this canonical form will always be a noncommutative polynomial in which each monomial is of the form $(\hat{a}^{\dagger})^m\hat{a}^n$.  Reversing the order of the generators reverses the order of these two factors in the canonical form.

Before I show these details, I need to define the algebra.  The default algebra is 

    NonCommutativeAlgebra[]
    (* NonCommutativeAlgebra[<|"Multiplication" -> NonCommutativeMultiply, "Addition" -> Plus, "Unity" -> 1, "Zero" -> 0, "CommutativeVariables" -> {}, "ScalarVariables" -> {}|>] *)

most of which is good.  However, the only assumed scalars are explicit complex numbers, so I'll add some scalar variables that I know I'll need for later calculations for the harmonic oscillator algebra:

    alg = NonCommutativeAlgebra[<|"ScalarVariables" -> {ℏ, m, ω}|>];

This is the minimum needed to do calculations with the oscillator algebra.  As an example, let's look at the operator $\hat{x}^2$, which is

    xx = (Sqrt[ℏ/(2 m ω)] (a + aD)) ** (Sqrt[ℏ/(2 m ω)] (a + aD));

To write this in canonical form, we reduce the polynomial relative to the generators, relations, and algebra using `NonCommutativePolynomialReduce`.  Note that this function uses `NonCommutativeExpand` behind the scenes, so we don't need to call that function explicitly.

    NonCommutativePolynomialReduce[xx, rels, gens, alg]
    (* {{(ℏ #1)/(2 m ω) &}, ℏ/(2 m ω) + (ℏ GeneralizedPower[NonCommutativeMultiply, a, 2])/(2 m ω) + (ℏ GeneralizedPower[NonCommutativeMultiply, aD, 2])/(2 m ω) + (ℏ aD ** a)/(m ω)} *)

The second term is the remainder, and this is exactly the canonical form of $\hat{x}^2$: notice that the lowering operators are all on the right in each term.  Note that if we reverse the order of the generators, we get the reverse with all the lowering operators on the *left*:

    NonCommutativePolynomialReduce[xx, rels, gens, alg]
    (* {{-((ℏ #1)/(2 m ω)) &}, -(ℏ/(2 m ω)) + (ℏ GeneralizedPower[NonCommutativeMultiply, a, 2])/(2 m ω) + (ℏ GeneralizedPower[NonCommutativeMultiply, aD, 2])/(2 m ω) + (ℏ a ** aD)/(m ω)} *)


This does everything we need it to!  But, we can add functionality with our rules.  For instance, let's suppose we wanted to implement a vacuum expectation value.  We can "abuse" `NonCommutativePolynomialReduce` to make this happen.

    Clear@vacuumExpectation
    vacuumExpectation[{{__Function}, expr_}] := NonCommutativePolynomialReduce[expr, gens, gens, alg]
    vacuumExpectation[expr_] := vacuumExpectation[NonCommutativePolynomialReduce[expr, gbrel, gens, alg]][[2]]


The second part of this definition takes an expression which can be expanded into a noncommutative polynomial in our generators and reduces it relative to the generators and relations, leaving a sum of terms that are in normal order (lowering operators all to the right).  Since the vacuum expectation value will render all non-constant terms zero in this canonical form, we can call `NonCommutativePolynomialReduce` *again*, but this time replace the relations with just the generators, which renders any non-constant term in the canonical form equal to zero, leaving just the constant term. This is what the first part of the definition does.

As an example, we can do

    vacuumExpectation[xx]
    (* ℏ/(2 m ω) *)

which should be the right answer!  Note that to ensure things are more robust, it probably makes more sense to define `vacuumExpectation` in such a way that you feed in the generators and such explicitly.  However, the trick of calling `NonCommutativePolynomialReduce` again with the generators in place of the relations already assumes a set of generators, relations, and canonical form.  So perhaps we just need to be careful.

---

I'm working on implementing a function that will compute general matrix elements in the oscillator eigenstates, but that can be fiddly, so perhaps I'll add it later.

---

Final note to answer Question 3:  This functionality cannot be used for an infinite number of generators and/or relations, because it relies on the explicit computation of a Groebner basis, and then the explicit computation of the reduction of the polynomial relative to this basis.  This can only be done, of course, if the numbers of generators and relations are both finite.