Modeling correlation (and covariance) matrices is a challenging problem due to the large dimensionality and positive-definiteness constraint. Here, we propose a novel Bayesian framework based on modeling the correlation matrix as a product of vectors on unit spheres. The covariance matrix is then modeled by using its decomposition into correlation and variance matrices. This approach allows us to induce flexible prior distributions for covariance matrices by proposing a wide range of distributions on spheres (e.g. the squared-Dirichlet distribution). Additionally, our method can be easily extended to dynamic settings in order to model real-life spatio-temporal processes with complex dependence structures (e.g., brain signals during cognitive tasks). To handle the intractability of the resulting posterior, we introduce a novel sampling algorithm called Adaptive Spherical Hamiltonian Monte Carlo. Using an example of Normal-Inverse-Wishart problem, a simulated periodic process, and an analysis of local field potential data, we demonstrate the validity and effectiveness of our proposed framework for (dynamic) modeling of covariance and correlation matrices.