deformable object squeezing

It's straightforward to model very slight deformation of the sphere (i.e., slightly flattening its surface). It's also pretty straightforward to model the case when the sides are already fully flattened and meeting at the corners. It's trickier to bridge the gap.

The small-deformation case has an exact solution for a single side (Case 2 here): the deformation $alpha$ (retaining the nomenclature used in the link) is $$alpha=frac{(3pi)^{2/3}}{2}P^{2/3}left(frac{1-sigma^2}{pi E}right)^{2/3}D^{-1/3}=frac{1}{2D^{1/3}}left(frac{3P}{ E^star}right)^{2/3},$$

where $P$ is the applied force, $sigma$ is sphere Poisson's ratio, $E$ is the sphere Young's modulus, $D$ is the sphere diameter, and $E^star=frac{E}{1-sigma^2}$ is the reduced Young's modulus.

To add a second side (the other side), we apply symmetry to obtain Case 3 in the link.

To add two more sides connected by an orthogonal axis (i.e., a sphere very slightly flattened on four sides), we can refer to generalized Hooke's Law and guess that we need to apply the same force as calculated from the two-side case, plus an additional force to suppress Poisson expansion (Example 5 in the link). In other words, $E$ has to be replaced by a modified biaxial modulus $E^prime=frac{E}{1-nsigma}$. (For a cuboid, $n=1$; I'm not sure what value to use for a sphere and request guidance. For small deformations of a sphere, $n$ may in fact be essentially zero, with no correction required.) In other words, squeezing an object along axis 1 causes lateral strain in axes 2 and 3, so if you want to obtain a certain strain along axis 2, you have to consider these Poisson effects.

(Two minor points here: (1) If Poisson's ratio were negative, then this form of the biaxial modulus would predict easier contraction, which accords with intuition (because the axes would be precontracted from the orthogonal stress; this behavior is unusual in engineering materials, but it's a good check). (2) If we proceed to apply compression along all three axes, then the modified elastic modulus—the triaxial modulus—for a cube is $E^prime=frac{E}{1-2sigma}$, which is one-third the bulk modulus, so that's consistent as well.)

When the four faces are fully compressed to meet at the new corners, then we might idealize the object as a rectangular cuboid and again apply generalized Hooke's Law to give

$$varepsilon_{11}=varepsilon_{22}=frac{s_{11}(1-2sigma)}{E},$$

where $varepsilon_{11}$ is further strain in the 1 direction and $s_{11}=s_{22}$ is additional stress applied by the constriction channel on each of the four faces; i.e., we apply the biaxial modulus because the strains in the 1 and 2 directions are being controlled and the stress in the 3 direction is idealized as zero.

For the tricky intermediate case of large deformation of the sphere before it assumes its square cross section, you could look into extensions to the small-deformation relations (e.g., here) or use large-displacement finite element analysis. You may also find this report on packing deformable spheres useful.

Edit: See also this paper on the packing of elastic spheres.