  The Yarkovsky effect plays an important role in the motions of small
celestial bodies. Increasingly improving observations bring the need of
high-accuracy modelling of the effect. Using a multiphysics software COMSOL, we
model the diurnal Yarkovsky effect in three dimensions and compare the results
with that derived from the widely adopted theoretical linear model. We find
that the linear model presents a high accuracy for spherical asteroids in most
cases. The ranges of parameters in which the relative error of the linear model
is over 10\% are explored. For biaxial ellipsoidal asteroids (particularly
oblate ones), the linear model systematically overestimates the transverse
Yarkovsky force by $\sim$10\%. The diurnal effect on triaxial ellipsoids is
periodic for which no linear model is available. Our numerical calculations
show that the average effects on triaxial ellipsoids are stronger than that on
biaxial ellipsoids. We also investigate the diurnal effect on asteroids of real
shapes and find it be overestimated by the linear model averagely by 16\%, with
a maximum up to 35\%. To estimate the strength of Yarkovsky effect directly
from the shape, we introduce a quantity "effective area" for asteroids of any
shapes, and find a significant linear relationship between the Yarkovsky
migration rate and the effective area. This brings great convenience to the
estimation in practice.
