The goal of the current study was to develop a computational framework for modelling the coupled fluid-structure interaction problem of squeeze films often encountered in MEMS devices. Vibratory MEMS devices such as gyroscopes, RF switches, and 2D resonators often have a thin plate like structure vibrating transversely to a Fixed substrate, and are generally not perfectly vacuum packed. This results in a thin air film being trapped between the vibrating plate and the fixed substrate which behaves like a squeeze film offering both stiffness and damping to the vibrating plate. For accurate modelling of the squeeze film effect, one must account for the coupled fluid-structure interaction. The majority of prior works attempting to address the coupled problem either approximate the mode shape of the vibrating plate or resort to cumbersome iterative solution strategies to address the problem in an indirect way. In the current work, we discuss the development of a fully coupled finite element based numerical scheme to solve the 2D Reynolds equation coupled with the 3D plate elasticity equation in a single step. The squeeze film solver so developed has been implemented into a commercial FEA package NISA as part of its Micro-Systems module. Further, extending on a prior analytical work, the effect of variable ow boundaries for an all sides clamped plate on squeeze film parameters has been thoroughly investigated. The developed FEM based numerical scheme has been used to validate the results of the prior analytical study. The developed numerical scheme models the 2D Reynolds equation thus limiting the model to account for the effects of the fluid volume strictly confined between the structure and the substrate. To study the effect of surrounding fluid volume ANSYS FLOTRAN simulations have been performed by numerically solving the full 3D Navier Stokes equation in the extended fluid domain for the different flow boundary scenarios. Cut-off frequencies are established beyond which one can consider a 2D fluid domain without considerable loss of accuracy. First, a displacement based finite element formulation is presented for the 2D Reynolds equation coupled with the 3D elasticity equation. Both lower order 8 node and higher order 27 node 3D elements are developed. Only a single type of 3D element is used for modelling along with a 2D fluid layer represented by the \wet" face of the 3D structural domain. The results from our numerical model are compared with experimental data from literature for a MEMS cantilever. The results from the 27 node displacement based elements show good agreement with published experimental data. The results from the lower order 8 node displacement based elements however show huge errors even for relatively fine meshes due to locking issues in modelling high aspect ratio structures. This limits the implementation of the displacement based solver in commercial FE packages where the available mesh generators are generally restricted to lower order 3D elements. In order to overcome the limitations faced by lower order elements (primarily locking issues) in modelling high aspect ratio MEMS geometries, a coupled hybrid formulation is developed next. A thorough performance study is presented considering both the hybrid and displacement based elements for lower order 8 node and higher order 27 node ele- ments. The optimal element choice for modelling squeeze film geometries is determined based on the comparative studies. The effect of element aspect ratio for hybrid and displacement based elements are studied and the superiority of hybrid formulation over displacement based formulations is established for lower order 8 node elements. The coupled hybrid nite element formulation developed for lower order elements is implemented in the commercial FEA package NISA. The implementation scheme to integrate the developed coupled hybrid 8 node squeeze film solver into the commercial FEA package is discussed. The pre-integration analysis and subsequent requirement gaps are first investigated. Based on the gap analysis, certain GUI modifications are undertaken and parser programs are developed to re-format data according to NISA input requirements. Certain special features are included in the package to aid in post processing data analysis by MEMS designers such as \frequency sweep" and \node of interest" selection. As a case study for validation, we also present the modelling of a MEMS cantilever and show that the simulation results from our software are in good agreement with experimental data reported in the literature. Finally as a case study, an extension of a prior analytical work, which studies the effect of varying flow boundaries on squeeze film parameters, is discussed. Explanations are provided for the findings reported in the prior analytical work. The concept of using variation in flow boundaries as a frequency tuning tool is introduced. The analytical results are validated with the coupled numerical scheme discussed before, by considering imposed mode shape for an all sides clamped plate as prescribed displacement to the fluid domain. The simulated results are used to study the intricacies in squeeze film damping and stiffness variations with respect to spatial changes in the fluid flow boundary conditions. In particular, it has been shown that the boundary venting conditions can be used effectively to tune the dynamic response of a micromechanical structure over a fairly large range of frequencies and somewhat smaller range of squeeze film damping. Next, the effect of the surrounding fluid volume for various venting conditions is studied. ANSYS FLOTRAN is used to solve for the full 3D Navier Stokes equation over the extended fluid domain. Results from the extended domain study are used to determine cut-off frequencies beyond which one need not resort to an extended mesh study, and yet be within 5% accuracy of the full extended mesh model.