The AAA algorithm was proposed in [Nakatsukasa/Sete/Treftehen '18] as a rational approximation tool that computes approximants as rational functions represented in barycentric form. The latter represents a numerically stable format that also imposes interpolation implicitly. The recent review paper [Nakatsukasa/Trefethen '26] surveys various improvements, generalizations, and applications of the AAA algorithm in the last 8 years. By harmoniously blending interpolation (as for non-iterative Loewner matrix methods) and least-squares fitting (as in the vector fitting approach), AAA aims to find a rational approximant by iteratively adjusting the barycentric form of the fitted model based on greedily selected interpolation points and a least-squares fit. In summary, AAA is applied iteratively, obtaining rational approximations of increased order, until a desired tolerance on the l2 approximation error is achieved. The model enrichment is performed by adding support points at the locations where the approximation error is largest. In particular, one neither needs to fix the surrogate model order nor (even more importantly) the support points in advance, since both are automatically chosen by the algorithm. Originally, the AAA algorithm was developed for approximation of scalar-valued functions, although recent extensions deal with vector- or matrix-valued functions as well (such as set-valued AAA, tangential AAA, or Block-AAA). We aim to extend the Block-AAA algorithm in [Gosea/Guettel '21] to more than two dimensions (going beyond 2D matrices to, e.g., multi-dimensional tensors). The reason for choosing this approach over other vector-valued extensions is that the latter typically fix the poles of the approximant for various parameter values, while Block-AAA offers more flexibility by allowing the location of poles to vary for both space and parameter values. This is particularly relevant for the class of problems we target here; the application of interest stems from structural dynamics, in which simulations and data analysis are challenged by large multi-dimensional data. Apart from two or three spatial directions, the problem coordinates include the dimension in frequency actuation and, possibly, dimensions to model uncertainties. If stored as a multidimensional array, these data quickly exceed storage capacities. Hence, efficient approximative representations are needed for data handling and, respectively, for simulations as a surrogate model. The proposed method aims precisely at attenuating this burden by harnessing the efficiency of canonical/hierarchical tensor decompositions, encoded explicitly in the generalized barycentric form of the approximant. The application studied here concerns modeling vibrating plates or vibroacoustic-coupled phenomena. Additionally, we consider uncertainty in the material parameters by modeling Young’s modulus as a log-normal random field. Various numerical results are presented to support the theoretical claims.