Does hyperbolic-solver really use a nonlinear solver?!

6 次查看（过去 30 天）
Lena
Lena 2014-1-27
评论： Hassan Diab 2014-3-14
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
Bill Greene 2014-1-29
Hi,
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:
http://www.mathworks.com/help/pde/ug/nonlinear-equations.html?searchHighlight=nonlinear+equations
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.
Regards,
Bill
  4 个评论
显示 2更早的评论
Lena
Lena 2014-1-30
I have just found the error!!! Instead of
u = u(:);
We of course have to write
u = [u(1,:)'; u(2,:)'];
Now we get
Max internal force = 3.05176e-05
I'm now trying to solve my original Problem, and inform you about my results! Thank you Bill!
Hassan Diab
Hassan Diab 2014-3-14
I am using matlab 2011, I have the pdenonlin command but i didn't find the line to change,
I have the same problem, pdenonlin try to divide H and get a memory error when there is not any Dirichlet condition.
do you think this issue is solved in the new release of Matlab.
Thanks

请先登录，再进行评论。

更多回答（4 个）

Bill Greene
Bill Greene 2014-1-27
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.
Bill
  1 个评论
Lena
Lena 2014-1-27
编辑：Lena 2014-1-27
Thanks for your quick answer! I've checked the things you recommended, but the difficulties continue to exist. As a model problem I considered a rotation around one fixed point, so we consider large deformations. In the attached example I left this fixed point out for simplicity, since otherwise it has to be given manually to the odesolver. But the main effect is also visible from this non-fixed case.
I don't use the hyperbolic solver itself, because I want to implement my own time integrator. Only for this example I used the ode15s, which I have replaced in my work by another solver. But the problem of giving over the right-hand-side, mass and jacobi is the same. I have just proceeded as in the code of hyperbolic.m and it's subroutines.
If you run the code with both c = @cCoeff_rotation; and c = @cCoeff_rotation_linear; you get the same results: The block moves down and it's size increases!
This increasing of the size under rotation should be avoided by use of the nonlinear coefficient c (nonlinear strain measure), so I imply there was never used a nonlinear solver.
function [] = rotation()
clear all; close all;
global p; global t; global e; global c; global b; global a; global f;
global E1; global nu1; global Q0; global K0;
Tspan = [0.0, 6];
E1 = 7.0e10;
rho1 = 2700;
nu1 = 0.34;
%Coefficients
a = 0.0;
d = rho1;
f = [0.0; 0.0];
c = @cCoeff_rotation;
%c = @cCoeff_rotation_linear;
b = @wbound_rotation;
%Geometry
gd = [2.0000;5.0000;-1.0000;-0.5000;1.0000;1.0000;-1.0000;-0.5000;...
-0.5000;-0.5000;0.5000;0.5000];
ns =[ 80; 49];
sf = 'P1';
gd = decsg(gd, sf, ns);
%Mesh
p =[1.0000 1.0000 -1.0000 -1.0000 -0.5000;
-0.5000 0.5000 0.5000 -0.5000 -0.5000];
e = [1 2 3 4 5; 2 3 4 5 1; 0 0 0 0 0; 1 1 1 1 1; 1 2 3 4 5; 1 1 1 1 1;...
0 0 0 0 0];
t = [1 2 3; 2 3 4; 5 5 5; 1 1 1];
np = size(p,2);
ne = size(e, 2);
u0 = zeros(2*np, 1);
ut0 = zeros(2*np, 1);
%Space-discretization
t0 = min(Tspan);
[~, MM0] = assema(p, t, 0, d, f, t0);
[K0, M0, F0] = assema(p, t, c, a, f, u0, t0);
[Q0, G0, ~, ~] = assemb(b, p, e, u0, t0);
nu = size(MM0,1);
mass = [speye(nu,nu) sparse(nu,nu); sparse(nu,nu) MM0];
%Odesovler
uu0 = [u0; ut0];
options = odeset('Jacobian', @Jac_rotation);
options = odeset(options, 'Mass', mass);
reltol = 1e-3;
options = odeset(options, 'RelTol', reltol, 'Stats', 'on', 'MaxOrder', 2);
[tt, uu] = ode15s(@ODEfct_rotation, Tspan, uu0, options);
uu = uu';
uSol = uu(1:2*np,:);
%Plots
result = uSol(:,end);
pdeplot(p,e,t)
hold on
pdeplot(p+[result(1:np) result(np+1:np+np)]', e, t, 'mesh', 'on', 'title', 'deformed mesh at last timestep');
xlabel('x-axis')
ylabel('y-axis')
axis equal
end
function [c] = cCoeff_rotation(p, t, u, time)
global E1;
global nu1;
[UX, UY] = pdegrad(p,t,u);
u1x = UX(1,:);
u1y = UY(1,:);
u2x = UX(2,:);
u2y = UY(2,:);
g1 = (1-nu1)/2;
e1 = E1/(1-nu1^2);
tn = size(t, 2);
c = zeros(16, tn);
for i=1:tn
if(isnan(time) || any(isnan(u)))
c(:, i) = NaN*ones(16, 1);
else
c(:, i) = e1*[1+0.5*u1x(i);
0.5*nu1*u1y(i);
g1*u1y(i);
g1;
0.5*u2x(i);
nu1+0.5*nu1*u2y(i);
g1+g1*u2y(i);
0;
g1*u1y(i);
g1;
nu1+0.5*nu1*u1x(i);
0.5*u1y(i);
g1+g1*u2y(i);
0;
0.5*nu1*u2x(i);
1+0.5*u2y(i)];
end
end
end
function [c] = cCoeff_rotation_linear(p, t, u, time)
global E1;
global nu1;
Gm1 = E1/(2*(1+nu1));
mu1 = 2*Gm1*(nu1/(1-nu1));
tn = size(t, 2);
c = zeros(16, tn);
for i=1:tn
c(:,i) = [2*Gm1+mu1 0 0 Gm1, 0 mu1 Gm1 0, 0 Gm1 mu1 0, Gm1 0 0 2*Gm1+mu1]';
end
end
function [df] = Jac_rotation(time,u) % according to pdehypdf.m
global Q0; global p; global c; global a; global t;
nu = size(Q0,1);
uu = u(1:nu);
[K0, M0, ~] = assema(p, t, c, a, [0;0], uu, time);
KK = K0+M0+Q0; %only G0 varies due to time-dependent NBC
df = - [sparse(nu,nu) -speye(nu,nu); KK sparse(nu,nu)];
end
function [fct]= ODEfct_rotation(time,u) % according to pdehypf.m
global K0; global b; global p; global e; global c; global a;
global f; global t;
nu = size(K0,1);
uu = u(1:nu);
[K0, M0, F0] = assema(p, t, c, a, f, uu, time);
[Q0, G0, ~, ~] = assemb(b, p, e, uu, time);
KK = K0 + M0 + Q0;
K = [sparse(nu,nu) -speye(nu,nu); KK sparse(nu,nu)];
FF =F0+G0;
F = [sparse(nu,1); FF];
fct = -K*u+F;
end
function [q,g,h,r]=wbound_rotation(p,e,u,time)
%WBOUND_ROTATION Boundary condition data.
bl=[
2 2 2 2 2
0 0 0 0 0
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 20 1
1 1 1 20 1
48 48 48 48 48
48 48 48 48 48
48 48 48 48 48
48 48 48 48 48
48 48 48 40 48
48 48 48 53 48
48 48 48 48 48
48 48 48 48 48
48 48 48 42 48
48 48 48 115 48
48 48 48 105 48
48 48 48 110 48
49 49 49 40 49
48 48 48 48 48
48 48 48 46 48
49 49 49 53 49
48 48 48 42 48
48 48 48 116 48
0 0 0 41 0
0 0 0 41 0
0 0 0 46 0
0 0 0 42 0
0 0 0 110 0
0 0 0 120 0
0 0 0 40 0
0 0 0 53 0
0 0 0 48 0
0 0 0 48 0
0 0 0 42 0
0 0 0 115 0
0 0 0 105 0
0 0 0 110 0
0 0 0 40 0
0 0 0 48 0
0 0 0 46 0
0 0 0 53 0
0 0 0 42 0
0 0 0 116 0
0 0 0 41 0
0 0 0 41 0
0 0 0 46 0
0 0 0 42 0
0 0 0 110 0
0 0 0 121 0
];
if any(size(u))
[q,g,h,r]=pdeexpd(p,e,u,time,bl);
else
[q,g,h,r]=pdeexpd(p,e,time,bl);
end
end

请先登录，再进行评论。

Bill Greene
Bill Greene 2014-1-27
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.
Bill
  1 个评论
Lena
Lena 2014-1-27
I have just tried to write my problem in the pdetool-GUI (=> use hyperbolic function). So I have to enter the coefficiets c11 , c12 ,... I tried for example to refer to the first component du1/dx of the gradient with ux(1,:) e.g.
c11 = 7.9150e10.*(1+0.5.*ux(1,:)) 7.9150e10.*(0.5*0.34.*uy(1,:)) 7.9150e10.*(0.33.*uy(1,:)) 7.9150e10*0.33
But I got "PDE Toolbox Error: Matrix dimensions must agree". What is the correct handling for the gradient in the system case? I think that, this is also one problem in my code.
Another remark to my code: My implementation follows the code hyperbolic.m . In the code of hyperbolic.m and it's subroutines I can nowhere find any call of a nonlinear solver such as pdenonlin.m , which uses Gauss-Newton or something like that to treat the nonlinear equation. Where exactly in hyperbolic.m (or parabolic.m ) is treated a nonlinearity like K(u)u=F like in the elliptic case? Thank you very much, Lena

请先登录，再进行评论。

Bill Greene
Bill Greene 2014-1-28
Hi,
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.
Bill
Bill Greene
Bill Greene 2014-1-28
Sorry about that! Apparently this GUI is too sophisticated for me ;-) I think I got it this time.
Bill
  1 个评论
Lena
Lena 2014-1-28
Hi Bill, thank you for your code, I have adapted my code to the hyperbolic-solver now (but keeping in mind that I want to replace its ode-solver). I simply used my implementation of the c -Coefficient from the above code:
Line 15: c = @cCoeff_rotation;
Line 16: %c = @cCoeff_rotation_linear;
If you run it with Line 15 (the nonlinear Version) I get the error message
Error using +
Matrix dimensions must agree.
With Line 16 it runs, but with the block's size increasing due to the use of a linear strain measure.
What is wrong in the nonlinear c-Coefficient-Function? Thank you very much for your help! Lena

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Structural Mechanics 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by