Random Walk 'Snake Problem'

Question

Introduction

A snake is traveling around the infinite Cartesian grid. It starts at (0,0) and takes the first step to the left to (-1,0). It continues spiralling counterclockwise, so it visits (-1,-1) next. Then, as the snake arrives at (0,-1), it cannot continue to (0,0) since this place has already been visitted. In such situation the snake chooses a random direction and continues its journey.

Question

How many steps will the snake take on average before biting its own body?

Simulation

I wrote some python code to check the answer numerically and is seems to be about 56.261616. Is there a way to estimate the answer mathematically?

The histogram after 10mln runs:

An example journey:

import numpy as np
import matplotlib.pyplot as plt
import random as rd
start=[0,0]
right=[1,0]
up=[0,1]
left=[-1,0]
down=[0,-1]

directions={0:right, 1:up, 2:left, 3:down}
facing=0
hist_data=[]
for jj in range(10000):
    pos=start
    visits=[start]
    step_num=0
    while True:
        new_pos=[sum(x) for x in zip(pos,directions[facing])]
        if new_pos not in visits:
            pos=new_pos
            visits.append(pos)
            step_num+=1
            facing=(facing+1)%4
        else: 
            facing=(facing+rd.randint(1,3))%4
            new_pos_right=[sum(x) for x in zip(pos,directions[0])]
            new_pos_up=[sum(x) for x in zip(pos,directions[1])]
            new_pos_left=[sum(x) for x in zip(pos,directions[2])]
            new_pos_down=[sum(x) for x in zip(pos,directions[3])]
            if (new_pos_right in visits) and (new_pos_up in visits) and (new_pos_left in visits) and (new_pos_down in visits):
                hist_data.append(step_num)
                break

plt.figure()
for ii in range(len(visits)-1):
    plt.plot([visits[ii+1][0],visits[ii][0]],[visits[ii+1][1],visits[ii][1]],'b-')

plt.plot(0,0,'g.',markersize=25)
plt.plot(visits[-1][0],visits[-1][1],'r.',markersize=25)
plt.title('Snake got stuck after '+str(step_num)+' steps')

plt.figure()
w=1
plt.hist(hist_data,bins=np.arange(min(hist_data), max(hist_data) + w, w))
print(np.average(hist_data))
plt.show()