Cette thèse porte sur le développement et l’étude de méthodes numériques pour le problème de poroélasticité, modélisé par les équations de Biot. La discrétisation historique de ce système consiste à utiliser une méthode de type éléments finis pour traiter l’équation d’équilibre mécanique, et une méthode de type volumes finis pour traiter l’équation de conservation de la masse du fluide. Or, en géosciences, les procédures de maillage utilisées pour reproduire les propriétés géométriques des milieux telles que les hétérogénéités, discontinuités ou failles produisent des cellules incompatibles avec la discrétisation par éléments finis. En conséquence, on s’intéresse à la récente méthode des éléments virtuels qui paraît capable de traiter des maillages plus généraux, et qui permettrait donc l’utilisation d’un unique maillage pour les deux modélisations. Après avoir appliqué la méthode des éléments virtuels au problème d’élasticité, on effectue le couplage avec l'écoulement discrétisé par volumes finis. Selon l’approximation choisie pour cette méthode, deux schémas sont créés : le schéma VEM-TPFA, qui utilise l'approximation volumes finis la plus simple et permet d'en proposer une analyse mathématique, et le schéma VEM-MPFA qui se base sur une approximation plus complète de manière à traiter des maillages plus généraux. En plus de la discrétisation, on s’intéresse à la stratégie de résolution du problème couplé dans le but d’en améliorer l’efficacité et la robustesse, en particulier en parallèle. On considère ainsi la stratégie fixed-stress split qui permet de se ramener à la résolution de deux sous-problèmes, pour lesquels on exhibe des solveurs performants.