The size effect in the failure of a hybrid adhesive joint of a metal with a fiber-polymer composite, which has been experimentally demonstrated and analytically formulated in preceding two papers, is here investigated numerically. Cohesive finite elements with a mixed-mode fracture criterion are adopted to model the adhesive layer in the metal-composite interface. A linear traction-separation softening law is assumed to describe the damage evolution at debonding in the adhesive layer. The results of simulations agree with the previously measured load-displacement curves of geometrically similar hybrid joints of various sizes, with the size ratio of 1:4:12. The effective size of the fracture process zone is identified from the numerically simulated cohesive stress profile at the peak load. The fracture energy previously identified analytically by fitting the experimentally observed size effect curves agrees well with the fracture energy of the cohesive crack model obtained numerically by optimal fitting of the test data.