Non-commutative algebra

I recommend you try the NCAlgebra package if you are going to do these kinds of computations. It is a mature package that has been under development for many years. The function you are looking for is NCExpand.

If you want something lighter-weight than NCAlgebra you can do things like

NonCommutativeMultiply[H___, Plus[a_, addends_], T___] :=  NonCommutativeMultiply[H, a, T] + 
    NonCommutativeMultiply[H, Plus[addends], T]

NonCommutativeMultiply[H___, a1, b1, T___] := NonCommutativeMultiply[H, c1 + d1, T]

Then (a1 + a2 + a3)**(b1 + b2 + b3) evaluates to

NonCommutativeMultiply[c1] + NonCommutativeMultiply[d1] + a1 ** b2 + 
a1 ** b3 + a2 ** b1 + a2 ** b2 + a2 ** b3 + a3 ** b1 + a3 ** b2 + a3 ** b3

where the only "issue" is that the singletons are wrapped with a NonCommutativeMultiply.

If you add

NonCommutativeMultiply[H___, a2, a1, T___] := NonCommutativeMultiply[H, e1, T]

then, at least with $Version="12.3.0 for Mac OS X x86 (64-bit) (May 10, 2021)", your first paren example gives e1**b1 and the other gives a2 ** c1 + a2 ** d1.

As of Version 14.3, this basic functionality is now built into Mathematica with the noncommutative algebra features. There is first of all the default

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

which is the free algebra on any number of generators. In order to expand noncommutative expressions, we can use:

NonCommutativeExpand[(a1 + a2 + a3) ** (b1 + b2 + b3)]
(* a1 ** b1 + a1 ** b2 + a1 ** b3 + a2 ** b1 + a2 ** b2 + a2 ** b3 + a3 ** b1 + a3 ** b2 + a3 ** b3 *)

In order simplify an expression relative to some known relations, we can use NonCommutativePolynomialReduce. This one's a little trickier to use, because the canonical form of the expression relative to the relations$-$which is the remainder of the polynomial reduction$-$depends on the order of the variables, so this requires some by-hand work.

Nonetheless, here are some examples from the OP. Suppose we know the relation a2 ** a1 == e1; we define a set of relations as

rels = {a2 ** a1 - e1};

If the expression we are reducing is

expr = a2 ** a1 ** b1;

then we construct a list of the relevant variables via

vars = NonCommutativeVariables[{rels, expr}];

and compute the polynomial reduction

NonCommutativePolynomialReduce[expr, rels, vars]
(* {{#1 ** b1 &}, e1 ** b1} *)

whose remainder e1 ** b1 is exactly the simplification we're looking for. note that in reversing the order of the variables, we get

NonCommutativePolynomialReduce[expr, rels, Reverse@vars]
(* {{0 &}, a2 ** a1 ** b1} *)

which doesn't effect the same simplification, but that's because the canonical form of the expression depends on the order of the variables.

As another example from the OP,

rels = {a1 ** b1 - (c1 + d1)};
expr = (a1 + a2 + a3) ** (b1 + b2 + b3);
vars = NonCommutativeVariables[{rels, expr}];
NonCommutativePolynomialReduce[expr, rels, vars][[2]]
(* c1 + d1 + a1 ** b2 + a1 ** b3 + a2 ** b1 + a2 ** b2 + a2 ** b3 + a3 ** b1 + a3 ** b2 + a3 ** b3 *)

All of these reductions are done relative to the default NonCommutativeAlgebra described above. We can change things by adding symbolic scalars (explicit complex numbers are automatically scalars) if we need to and then explicitly including the algebra in the call, e.g.

rels = {a2 ** a1 - e1};
expr = a2 ** (s a1) ** b1;
vars = NonCommutativeVariables[{rels, expr}, alg]
NonCommutativePolynomialReduce[expr, rels, vars, alg][[2]]
(* {a1, a2, b1, e1} *)
(* s e1 ** b1 *)