I don't understand how you're getting the strain energy from the eigenvalue solver. It is my understanding that the eigenvalue solver will only give you the strain energy at a particular eigenfrequency. We're interested in the strain energy from the beam deformation at frequencies below the first eigenmode.
We have been encountering an error in COMSOL for a while of the form "Failed to find a solution.; The relative residual (###) is greater than the relative tolerance.; Returned solution is not converged.; -Feature: Stationary Solver 1 (sol1/s1); -Error: Failed to find a solution." The error has occurred when attempting to find a stationary solution in models with boundary loads and no fixed constraints (preventing an edge from moving). We wanted to determine what the error was caused by to allow us to use the stationary solver, or at least to confirm that the error was not indicative of a problem in our model. To this end, I designed a few very simple COMSOL models in which I was able to reproduce the behaviour and attempted to determine the root of the issue.
I first constructed a somewhat similar model to ours using a cylinder of fused silica with all of the default values and a normal meshing. I applied a boundary load of 1N on one of the faces and ran a stationary solver. As expected, the solver failed to converge since it had no boundary condition preventing it from accelerating continuously. Applying a force of 1N on the opposing face resulted in the same error as above, which replicates the previous error since the system is failing to converge in a case where it should. I decided also to make an even simpler 2-D model of a square. Applying 1N forces to opposing sides on the square again returned the error above.
Both of these models were able to be evaluated using at least an eigenfrequency solver as noted on the primary model in the previous eLog. I looked on the COMSOL forums and read through some more of their documentation and saw the suggestion in response to this error to use a time-dependent solver and simply view times after the system will have settled to a stationary state (#2 https://community.cmc.ca/docs/DOC-1453). I attempted this on the test models and both of the time dependent solutions converged without error to their expected solutions (compression between the faces on which forces were being applied). This may be a sub-optimal computational method though as even in the simple cylinder case with 6133 elements and a simple force profile, the solution took several minutes to run. For the cylindrical model, I evaluated the strain energy using both an eigenfrequency and time dependent solver and obtained the same result using both of the solvers. The eigenfrequency solver evaluates much more quickly than the time dependent solver, and in the primary model as I noted in my previous eLog, the strain energy obtained using the eigenfrequency and frequency domain solvers agreed, so it seems that the best manner in which to proceed is to use the eigenfrequency solver to compute the strain energy.
The source of the error is still unknown, but given these tests, it seems very unlikely to be indicative of a problem in our model. We still have a very significant disagreement between the simulated results and the calculated values. I am going to spend the next day or so looking through both the COMSOL model and the analytic calculation and checking them for errors which could cause this discrepancy. I will then start reading the documentation on Livelink for Matlab and try to implement it.