This paper presents a novel trust-region method for the optimization of multiple expensive functions. We apply this method to a biobjective optimization problem in fluid mechanics, the optimal mixing of particles in a flow in a closed container. The three-dimensional time-dependent flows are driven by Lorentz forces that are generated by an oscillating permanent magnet located underneath the rectangular vessel. The rectangular magnet provides a spatially non-uniform magnetic field that is known analytically. The magnet oscillation creates a steady mean flow (steady streaming) similar to those observed from oscillating rigid bodies. In the optimization problem, randomly distributed mass-less particles are advected by the flow to achieve a homogeneous distribution (objective function 1) while keeping the work done to move the permanent magnet minimal (objective function 2). A single evaluation of these two objective functions may take more than two hours. For that reason, to save computational time, the proposed method uses interpolation models on trust-regions for finding descent directions. We show that, even for our significantly simplified model problem, the mixing patterns vary significantly with the control parameters, which justifies the use of improved optimization techniques and their further development.