Hello I am trying to simulate crawling motion of a bilayer beam using Abaqus. Kind of similar to this: https://youtu.be/os-Z3ucoFpA

I am using multiphysics enabled (diffusion, deformation, thermal) user-defined element subroutine to simulate this. But the main problem so far is coming from contact formulation and boundary conditions. The model looks like as the attached image.

I defined hard frictional contact between the bottom surface and an analytical rigid surface. Gravitational body force is working as y-axis constrain. Using only friction in x-axis gives me numerical singularity issues because of rigid body motion. So I decided to constrain some nodes in the model. In first diffusion loading step, N2 was constrained in x-direction. In subsequent cooling and heating step, N1 and N3 were constrained, respectively, to obtain crawling motion. N1 and N3 are nodes slightly above the contact surface as Abaqus has conflict with applying BC on contact surface. First loading step works fine and it maintains N2 boundary condition, but N1 and N3 boundary conditions are being ignored in subsequent steps. With these BCs imposed, it still has rigid body motion (see attached animation please). I used keywords for changing boundary condition. It works with temperature and diffusion but not with displacement. I am also attaching important section from the input file. I am not sure what exactly I am missing here. What would be an appropriate way of modeling contact and BC together? Any help will be appreciated.