I don't know how you did your simulation, but you should at least include some kind of damping term in the motion of the birds, otherwise every configuration starting out-of-equilibrium will keep oscillating forever. Maybe add a term like $\overrightarrow{F_{\mathrm{viscous}}} = - \alpha \overrightarrow{v}$. If you are not only interested in the equilibrium position but also the dynamics of the flock for small perturbations, you can also add a random Langevin force $\overrightarrow{F_{\mathrm{Langevin}}}(t) = \eta(t)$, with $\left\langle \eta(t) \eta(t') \right\rangle = A \delta(t-t')$.
For the exact form of the interaction force between the birds, I don't know if it is possible to capture all of the complexity of the phenomenon simply by adding pairwise, conservative force to the problem, but we can play around and see if the results are good enough in the end.
If you want a "homogeneous" density of birds, you can try $$\overrightarrow{F}_{\overrightarrow{x}\to\overrightarrow{x'}}= \frac{k}{N} \left(\overrightarrow{x} - \overrightarrow{x'}\right) - \beta \frac{\left(\overrightarrow{x} - \overrightarrow{x'}\right)}{\left|\overrightarrow{x} - \overrightarrow{x'} \right|^2},$$
where $\overrightarrow{F}_{\overrightarrow{x}\to\overrightarrow{x'}}$ is the force exerted by a bird at position $\overrightarrow{x}$ on a bird at position $\overrightarrow{x'}$, $N$ is the number of birds in the flock (a normalization factor), and $k$ and $\beta$ are positive constants whose relative value will change the average spacing between birds at equilibrium. You also need to fix the relative value of the Langevin force and viscous force to somewhat realistic values so that you reach an equilibrium fast enough but you can still see the small perturbations in position. I will come with a more detailed explanation for the form of the force and a simulation later.
Reason for the choice of the force:
The first term in $\overrightarrow{F}_{\overrightarrow{x}\to\overrightarrow{x'}}$, $\frac{k}{N} \left(\overrightarrow{x} - \overrightarrow{x'}\right)$, simply correspond to the force of an harmonic oscillator of strengh $k$ centered on the center of mass of the flock. So, assuming the center of mass is centered around $0$, the force is simply $\overrightarrow{F}(\overrightarrow{x}) = -k \overrightarrow{x}$.
The reason for the second term is actually a little bit subtle. By analogy with the 3D Coulomb force, you can show that the field $\overrightarrow{f}(\overrightarrow{x})$ created by a density of birds $\rho(\overrightarrow{x}')$ follows a 2D analog to Gauss law in electromagnetism, which is $\overrightarrow{\nabla} \cdot \overrightarrow{f} = 2 \pi \beta \rho$. You can then show that a constant density of birds $\rho(r) = N/(\pi R^2)$ for $0 < r < R$, and $\rho(r>R) = 0$ is creating a repulsive force of the form $\overrightarrow{f}(\overrightarrow{x}) = \beta N \overrightarrow{x}/R^2$ (for $\left| \overrightarrow{x} \right| < R$), which exactly compensates the harmonic oscillator force when $R = \sqrt{\beta N/k}$. This is of course only true for a continuous "density" of birds, which is not the case in reality but for $N \gg 1$, you can assimilate $\rho$ to the average density of birds. In short, for $N \gg 1$, the birds at equilibrium should arrange themselves in a somewhat homogeneous configuration within a radius $R = \sqrt{\beta N/k}$.
Simulation
I implemented the problem in Python using all the forces that I described earlier, with $k = \beta = 1$, $A = 3$ and $\alpha = 2$. I integrate in time steps of $dt = 0.01$ from $t = 0$ to $t = T = 10$. The mass $m$ of each bird is also set to $1$. Initially, the position of each bird is drawn randomly and uniformly in a square smaller than the expected radius, so that you can see the out-of-equilibrium dynamics.
Here is what the trajectories look like on a typical run
And here is a nice gif of the time evolution of the bird flock:
As you can see, despite some random fluctuations in position, the birds keep moving around an "equilibrium" configuration where the number of birds per surface area is approximately constant within a given radius.
If you're worried that the birds get too close to each other and it is unrealistic, try decreasing the $A$ parameter (controlling the amplitude of the random motion of birds). See below a gif for $A = 0.5$.
Code used
import numpy as np
import matplotlib.pyplot as plt
plt.rc('font', size=20)
N = 160 #number of birds
sizex = 20 #size of the image for plotting
sizey = 20
k = 1.0 #strength of the (long-range) attractive force
beta = 1.0 #strength of the (short-range) repulsive force
alpha = 2.0 #"viscosity"
A = 3.0 #intensity of the random erratic movement of each bird
dt = 0.01 #timestep
T = 10 #final time
expected_radius = (beta*N/k)**0.5 #expect radius of the bird flock by analogy with the continuous case
X_pos_0 = (np.random.random(N)-0.5)*0.5*sizex #draw the initial X position of birds randomly
X_pos_0 -= np.mean(X_pos_0) #center it on 0
Y_pos_0 = (np.random.random(N)-0.5)*0.5*sizey #same in the Y direction
Y_pos_0 -= np.mean(Y_pos_0)
X_pos = np.zeros((N, int(T/dt)+1)) #array of X position of each bird
Y_pos = np.zeros((N, int(T/dt)+1)) #array of Y position of each bird
X_speed = np.zeros((N, int(T/dt)+1)) #array of X velocity of each bird
Y_speed = np.zeros((N, int(T/dt)+1)) #array of Y velocity of each bird
X_pos[:,0], Y_pos[:,0] = X_pos_0, Y_pos_0
@np.vectorize
def repulsive_force_x(X,Y): #function used to compute the pairwise forces at each timestep of the algorithm
if X==0 and Y==0:
return 0.0
else:
return -beta*X/(X**2+Y**2)
@np.vectorize
def repulsive_force_y(X,Y):
if X==0 and Y==0:
return 0.0
else:
return -beta*Y/(X**2+Y**2)
time = np.linspace(0, T, int(T/dt)+1)
for i in range(int(T/dt)):
langevin_force_x = np.random.normal(0, A, N)/dt**0.5 #the Langevin force has to be normalized like this in the limit dt -> 0
langevin_force_y = np.random.normal(0, A, N)/dt**0.5
viscous_force_x = -alpha*X_speed[:, i]
viscous_force_y = -alpha*Y_speed[:, i]
attractive_force_x = -k*(X_pos[:,i]-np.mean(X_pos[:,i]))
attractive_force_y = -k*(Y_pos[:,i]-np.mean(Y_pos[:,i]))
X1, X2 = np.meshgrid(X_pos[:,i], X_pos[:,i])
distance_matrix_X = X1-X2 #array containing all pairwise distances x-x' along X
Y1, Y2 = np.meshgrid(Y_pos[:,i], Y_pos[:,i])
distance_matrix_Y = Y1-Y2 #array containing all pairwise distances y-y' along Y
force_X = np.sum(repulsive_force_x(distance_matrix_X, distance_matrix_Y), axis=1)
force_Y = np.sum(repulsive_force_y(distance_matrix_X, distance_matrix_Y), axis=1)
X_speed[:, i+1] = X_speed[:, i] + dt*(langevin_force_x + viscous_force_x + attractive_force_x + force_X) #v(t+dt) = v(t) + a(t)*dt (m=1)
Y_speed[:, i+1] = Y_speed[:, i] + dt*(langevin_force_y + viscous_force_y + attractive_force_y + force_Y)
X_pos[:, i+1] = X_pos[:, i] + dt*X_speed[:, i+1] #x(t+dt) = x(t) + v(t+dt)*dt (symplectic integration)
Y_pos[:, i+1] = Y_pos[:, i] + dt*Y_speed[:, i+1]
if i%(int(T/dt)//10) ==0: #progression of the calculation
print(10*i//(int(T/dt)//10), '%')
cmap_cool = plt.get_cmap('cool')
#plot initial and final positions only
fig, ax = plt.subplots(1, figsize=(15,12))
theta = np.linspace(0, 2*np.pi, 200)
ax.plot(expected_radius*np.cos(theta), expected_radius*np.sin(theta), 'black', label="Expected radius")
ax.plot(X_pos[:, 0], Y_pos[:, 0], 'o', label="Initial position", color=cmap_cool(0)) #initial position in blue
ax.plot(X_pos[:, -1], Y_pos[:, -1], 'o', label="Final position", color=cmap_cool(0.99)) #final position in purple
ax.axis("scaled")
ax.set_xlim(left=-24, right=24)
plt.legend(loc='best')
#plot trajectories too
fig, ax = plt.subplots(1, figsize=(15,12))
theta = np.linspace(0, 2*np.pi, 200)
ax.plot(expected_radius*np.cos(theta), expected_radius*np.sin(theta), 'black', label="Expected radius")
ax.plot(X_pos[:, 0], Y_pos[:, 0], 'o', label="Initial position", color=cmap_cool(0)) #initial position in blue
ax.plot(X_pos[:, -1], Y_pos[:, -1], 'o', label="Final position", color=cmap_cool(0.99)) #final position in purple
ax.axis("scaled")
ax.set_xlim(left=-24, right=24)
for i in range(20):
for bird in range(N):
ax.plot([X_pos[bird, i*(int(T/dt)//20)], X_pos[bird, (i+1)*(int(T/dt)//20)]], [Y_pos[bird, i*(int(T/dt)//20)], Y_pos[bird, (i+1)*(int(T/dt)//20)]], '--', color=cmap_cool(i/20), lw=1)
plt.legend(loc='best')
Hope this helps.
