To provide a more accurate and efficient theoretical tool for optimizing the design of biomimetic muscle-driven soft robots with large deformations, we propose an open-source multiphysics finite element framework in this paper. Within this framework, the directional arrangement and anisotropic properties of biomimetic muscles are modeled using the Holzapfel-Ogden model, while the electromechanical coupling characteristics of electrically driven muscles, such as dielectric elastomers and biological tissues, are described by the Aliev-Panfilov model. The soft robots, composed of various materials, are divided into distinct domains, each governed by its own constitutive equation. Furthermore, a unified continuum and variational multiscale (VMS) formulation based on the Gibbs free energy function is employed to effectively address the volumetric locking issue in modeling the near-incompressible behavior of these robots. The framework’s effectiveness and practicality are demonstrated through the application case study of a bionic jellyfish robot.