Computational Fluid Dynamics (CFD) is a powerful tool which can help with the geometry optimization of continuous milli-scale reactors, which often are highly complex devices. Attempting to perform this optimization by manually modifying and testing geometry configurations can however be tedious and computationally inefficient. Addressing this problem, we present a framework in which the CFD software COMSOL Multiphysics is coupled with the multi-objective Bayesian Optimization algorithm TSEMO (Thompson sampling efficient multiobjective optimization), implemented in MATLAB. The mixing element geometry of a Miprowa Lab millireactor is parameterized, and the framework automatically executes CFD simulations to minimize areas of stagnating flow and maximize the mixing performance. The framework is able to find Pareto-optimal reactor variations, and can easily be adapted for other devices and objectives.