  Twisted bilayer graphene (TBG) exhibits a wide range of intriguing physical
properties, such as superconductivity, ferromagnetism, and superlubricity.
Depending on the twist angle, periodic moir\'e superlattices form in twisted
bilayer graphene, with inhomogeneous interlayer coupling and lattice
deformation. For a small twist angle (typically <2{\deg}), each moir\'e
supercell contains a large number of atoms (>10,000), making it computationally
expensive for first-principles and atomistic modeling. In this work, a finite
element method based on a continuum model is used to simulate the inhomogeneous
interlayer and intralayer deformations of twisted graphene flakes on a rigid
graphene substrate. The van der Waals interactions between the graphene layers
are described by a periodic potential energy function, whereas the graphene
flake is treated as a continuum membrane with effective elastic properties. Our
simulations show that structural relaxation and the induced strain localization
are most significant in a relatively large graphene flake at small twist
angles, where the strain distribution is highly localized as shear strain
solitons along the boundaries between neighboring domains of commensurate AB
stacking. Moreover, it is found that there exist many metastable equilibrium
configurations at particular twist angles, depending on the flake size. The
nonlinear mechanics of twisted bilayer graphene is thus expected to be
essential for understanding the strain distributions in the moir\'e
superlattices and the strain effects on other physical properties.
