In the current study, the major aim was to develop computationally-effective numerical procedures with the capability of modelling the thermal and mechanical response of fibre-reinforced composites. The attempt was made with a focus on the applicability of the procedures to a variety of synthetic and natural fibres with various geometrical characteristics so that a materialindependent framework is developed. The work is divided into three main stages: developing methodologies, investigating the performance of the available analytical tools, and extending the material models by numerical techniques. The finite element computational framework was used and a library of general-purpose subroutines was developed that encompassed the methodologies. Then, finite element prototypes were challenged through a variety of thermal and mechanical analyses to set up the interaction between micro- and macro-scales. Representative volume elements (RVEs) were created for a variety of short/long and aligned/randomly-oriented fibres in thermal analyses. Fibre orientation tensors were used to carry out clustering and spectral analyses in order to characterise the morphology of the heterogeneities. According to the results of sensitivity analyses, although spectral analysis seems to be less sensitive to local variation of fibre direction, it correlates better with the effective properties. Furthermore, the element elimination technique was used to indirectly model the progression of damage. This technique was used in both discrete fibre modelling and homogenised elements. It was shown that by following the latter case, mesh-independent results could be obtained. Moreover, natural fibre-reinforced composites were modelled using discrete fibre elements and a new strength-updating scheme was proposed and implemented. The numerical results showed the detrimental effect of considering length-dependent strength properties in computational modelling. Namely, the strength of fibres must be updated for the remaining portions to obtain results that are closer to the experimental one. At the final stage of the study, the localised mesoscopic data was collected through the introduced concept of auxiliary maps. Auxiliary maps of volume fraction and fibre orientation data were obtained and it was shown that their resolution should be linked to the mesh density of the model. A semi-analytical model was created to demonstrate the performance of the purposed method. The proposed models were used in single-scale elastic analyses but are able to be extended to coupled multi-scale analyses.