The time dependent behaviours like creep and stress relaxation of materials like polymer, bitumen, rubber, and soft tissues are idealised using viscoelastic constitutive relations. There are many frameworks like integral, rate type, and fractional order representations which relates the stress and strain explicitly or implicitly to model the viscoelastic response. However, a general numerical technique capable of solving boundary value problems from any of these representations is lacking in literature. In this talk a methodology to formulate a 2D rectangular element to solve plane problems involving viscoelastic constitutive relation is developed. Stress fields satisfying equilibrium equations using Airy’s potentials which are expressed as a linear combination of C2 basis functions. Strain field is derivedfrom a continuous displacement field obtained from a linear combination of C0 displacementbasis functions. The appropriate linear combination of these stress and displacement basisfunctions are determined such that the resulting stress and strain fields satisfy the constitutive relation subjected to the satisfaction of the constraints arising from the traction anddisplacement boundary conditions. Since a viscoelastic constitutive relation involves, stress,strain and their rates, the basic unknowns of stress potentials and displacement fields or therate of basic unknowns can be considered as optimization variables for minimizing the errorin satisfying the constitutive relation. Two Algorithms are proposed based on this choiceof optimization variable. The accuracy and efficiency of the proposed algorithms are studied through different boundary value problems involving various forms of the viscoelasticconstitutive relations – Linear integral, Rate type, an explicit non linear, a fractional order integral form – for two loading histories. Using the developed rectangular element viscoelasticbeam bending problem is solved for the same set of different constitutive relations. It is found that the algorithm involving rate of the basic unknowns as the optimization variable could solve boundary value problem involving any viscoelastic constitutive relation and is accurate and computationally efficient.