Liu and Thorne demonstrated in their paper that the optimal set of boundary conditions for a problem like the one I am working on is the balance the force applied to one of the mirror faces by an opposing force of equal magnitude distributed throughout the body of the mirror.  When we first attempted to implement this boundary condition, we obtained an error from COMSOL indicating that a solution could not be generated as "The relative error is greater than the relative tolerance."  I attempted to find a solution to this problem so that we can run our simulations using the correct boundary conditions.

I tried a number of settings changes and tweaks suggested on the COMSOL forums and in official documentation on the error (https://community.cmc.ca/docs/DOC-1453), but without success.  Most of the suggestions were not relevant to our model, and confusingly when the number of elements in the mesh was increased from 13628 to 24570 the relative residual increased contrary to our initial attempts in http://nodus.ligo.caltech.edu:8080/COMSOL/3.  It seems like the next action to take is to attempt to generate a simpler model and attempt reproduce and fix the same error.

I also decided to run both an eigenfrequency and frequency domain solver on our model using improved meshing with 24570 elements.  I first ran these solvers using the rigid back constraint and both of them returned Umax=1.49931*10^-10 which gives S_x(100 Hz)=1.97673*10^-40 m^2/Hz which agrees quite well with the solution obtained from our previous simulations which gave S_x(100 Hz)=2.0157*10^-40 m^2 / Hz.  This still significantly differs from the analytically obtained value of 7.80081*10^-40 m^2 / Hz using Liu and Thorne's calculation.

I then changed the boundary condition to use Liu and Thorne's gravitational body load boundary condition and both of these methods returned Umax=1.17164*10^-9 J which corresponds to S_x(100 Hz)=1.54472*10^-39 m^2/HZ.  This only differs from the analytic result by a factor of 1.98, which is a significant improvement over the previous simulations, but still a much more significant difference than I would hope for given that they now use the same boundary conditions and it seems unlikely that using a stationary state solver will give a different result when both the eigenfrequency and frequency domain solvers agree to 6 digits.

edit: Taking the advice of the error documentation, I also ran a time dependent solver which returned the same result as the eigenfrequency and frequency domain solvers.