Does hyperbolic-solver really use a nonlinear solver?!
6 views (last 30 days)
My PDE is time-dependent and has a coefficient c , which is depending on the gradient of the solution u (due to the usage of Green-Lagrange-strain measure).
Following the documentation „a nonlinear solver is available for the nonlinear elliptic PDE, the parabolic and hyperbolic equation solvers also solve nonlinear and time-dependent problems.“ For elliptic problems it is necessary to explicitly use pdenonlin instead of assempde . But according to the documentation „the parabolic and hyperbolic functions call the nonlinear solver automatically “.
I don't use the PDE-Tool-GUI, but implemented my own method of line with the pdetool-command-line functions. I looked up the codes of the hyperbolic / parabolic -solver and it's subroutines, and I can nowhere find any call of a nonlinear-solver in case of solution-dependent coefficients!
The results of my calculations also strongly imply, that the nonlinearity induced by c is not treated correctly, since the results are the same as in the linear case. So if anyone can tell me, if it is possible to solve nonlinear transient PDEs with the pdetool-command-line functions, I would be very happy, since I am at my wit's end for another error source producing the unexpected results.
Bill Greene on 29 Jan 2014
Unfortunately, I believe you have uncovered a bug in the hyperbolic function that occurs when you have a nonlinear equation but no Dirichlet boundary conditions. I believe there is a simple fix. If you replace the line
uFull = u;
in toolbox/pde/pdehypdf.m with
uFull = u(1:nu);
your example will run to completion. The line to replace is line 16 in MATLAB R2013b.
However, the results may not be what you intended. I'm guessing that your problem is just a rigid body rotation of the block. If so, the rotated block is larger than the original block which is obviously not correct for the rigid body case.
Before saying more, I should first say that I have wanted to derive a c-matrix for a Lagrangian formulation of 2D elasticity but have been somewhat discouraged by all the algebra required. So I don't have such formulation but believe it is possible to derive one and suspect it is more complicated than what you have in cCoeff_rotation.
I did a simple test of your cCoeff_rotation function (file attached) where I created a simple block, defined the displacement vector for a 45-degree rotation, and calculated the internal forces using cCoeff_rotation. This vector is not zero as expected.
You probably have already looked at this page:
but there is a key equation there that is easy to overlook. Under the line "The residual vector ρ(U) can be easily computed as" you will see that the residual vector is calculated from a product of K and U. This is different from the typical nonlinear mechanics notation where K is a tangent stiffness matrix.
Hope this documentation page and the fix will point you in the right direction.
More Answers (4)
The first thing to check is that your version of MATLAB is R2012b or newer. The nonlinear hyperbolic (and parabolic) solvers have been available only since then.
Since you didn't post your code for computation of the c-coefficient, I'm speculating about another possibility. The hyperbolic (and parabolic) functions detect that c is a function of u (or ux, uy) by passing a vector, u=NaN, to your routine for calculating c. They expect that the computed c-matrix contains at least one NaN. Usually this requirement is satisfied automatically because computations involving NaN produce a NaN as a result. But occasionally you need to treat this NaN call specially.
Finally, I suggest you place a debugger break point in your routine that calculates the c-coefficient and then run your example problem. Obviously, the hyperbolic solver should be calling this routine many times during the analysis and each time it will stop at your break point.
If none of these pointers help you resolve the issue, please post your code and I'll take a look.
Oh, sorry, I thought you were using the documented hyperbolic function. I'm afraid I have no idea why your code isn't doing what you expect.
I think that referring to gradient components in the form e.g. ux(1,:) is fine. But entering complicated nonlinear coefficients for the system (N=2) case in the pdetool GUI can be tricky. I suggest solving the problem using command-line functions instead of the GUI.
I have attached code for a simple linear cantilever beam transient response example that uses the hyperbolic function. If you already have expressions for the geometrically nonlinear c-coeffients, you have already done most of the work. It should be straightforward to code those in the calcCmat function in the example.