Nonlinear material simulation in pdetool
5 views (last 30 days)
Show older comments
I tried to solve the similar problem explained here with a nonlinear material. And I wrote this code (with a simpler geometry): (And I defined the c coefficient with an interpolation function with respect to the magnitude of the gradient of the solution|grad(solution)|)
u0=4e-7*pi;
model = createpde(1);
gd = [ 3,4,+2.0,-2.0,-2.0,+2.0,1.5,1.5,-1.5,-1.5;...
3,4,+0.8,-0.8,-0.8,+0.8,0.8,0.8,-0.8,-0.8;...
3,4,+0.4,-0.4,-0.4,+0.4,0.4,0.4,-0.4,-0.4;...
3,4,+1.2,+0.8,+0.8,+1.2,0.4,0.4,-0.4,-0.4;...
3,4,-0.8,-1.2,-1.2,-0.8,0.4,0.4,-0.4,-0.4;]'; % geometry gd
ns = char('bound','c1','c2','c3','c4')';
sf = 'bound+c1+c2+c3+c4';
[ dl, bt] = decsg(gd,sf,ns);
model.SolverOptions.ReportStatistics = 'on';
geometryFromEdges(model,dl); % include geometry
generateMesh(model,'Hmax',0.04,'GeometricOrder','Linear');
applyBoundaryCondition(model,'dirichlet','Edge',[1,2,13,14],'u',0);
specifyCoefficients(model,'m',0,'d',0,'c',1,'a',0,'f',0,'Face',1);
cCoef = @(~,state) interp1([0 1-4 2e-4 2],[1/(1000) 1/(500) 1/(100) 1/(10)],sqrt((state.ux).^2 + (state.uy).^2));
specifyCoefficients(model,'m',0,'d',0,'c',cCoef,'a',0,'f',0,'Face',2);
specifyCoefficients(model,'m',0,'d',0,'c',1,'a',0,'f',u0,'Face',3);
specifyCoefficients(model,'m',0,'d',0,'c',1,'a',0,'f',-u0,'Face',4);
specifyCoefficients(model,'m',0,'d',0,'c',1,'a',0,'f',-u0,'Face',5);
results = solvepde(model);
B=sqrt((results.XGradients.^2+results.YGradients.^2)); %B=curl(A);
pdeplot(model.Mesh.Nodes,model.Mesh.Elements,'XYData',B,'Mesh','on');
figure; quiver(results.Mesh.Nodes(1,:),results.Mesh.Nodes(2,:),results.YGradients',-results.XGradients');
I supposed that the solver will handle the problem iteratively, however, I got this instead:
Iteration Residual Step size Jacobian: Full
0 1.0842e-19
And the solution is exactly the same as the time I replace the c coefficient with 1/1000. I expected to see some sort of saturation in the material.
3 Comments
Viktor Szuhai
on 27 Sep 2023
Nice job. But I believe I found a small typo: an "e" is missing in the cCoef definition. Correct:
cCoef = @(~,state) interp1([0 1e-4 2e-4 2],[1/(1000) 1/(500) 1/(100) 1/(10)],sqrt((state.ux).^2 + (state.uy).^2));
Answers (0)
See Also
Categories
Find more on Geometry and Mesh in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!