MATLAB: Is the Thermal Model giving me an error

errorPartial Differential Equation Toolboxthermalbcthermalconductivitythermalicthermalmodelwarning

I am trying to solve for a thermal model using the "solve" function. When I try calling solve on my model, I am getting a warning and an error. I use "thermalBC" to set my boundary conditions and "thermalProperties" to set my 'ThermalConductivity' property with an anonymous function. Then I call solve on my model.

Warning: Matrix is singular to working precision.

> In pde.EquationModel/solveStationaryNonlinear (line 21)

In pde.ThermalModel/solve (line 141)

In Thermal (line 28)

 Error using pde.EquationModel/solveStationaryNonlinear (line 27)
 Unsuitable initial guess U0  (default: U0=0).
 Error in pde.ThermalModel/solve (line 141)
             u = self.solveStationaryNonlinear(coefstruct, u0);
 Error in Thermal (line 28)
 results = solve(model);

Best Answer

This warning and error are due to the face that initial conditions are not set and therefore default to zero. The 'ThermalConductivity' function has some part of the equation that is infinite, for example if you are calling log10, then log10(0) is -Inf.

Make sure to set your initial conditions using the "thermalIC" function and to set it to the appropriate value so that the function does not turn out to be infinite and error. See the link below for more information:

https://www.mathworks.com/help/pde/ug/pde.thermalmodel.thermalic.html

Related Solutions

MATLAB: Using a function handle for ThermalConductivity is not changing resulting temperature plot

The Solver is misclassifying thermal conductivity as global assignment with k1=@(location,~) 1*location.x;. This is happening because at location.x = 0, which is the query point used to detect difference in conductivity, both k1 and k2 yield the same value of zero.

In the meantime a workaround is to add a small noise, like k2=@(location,~) 20*location.x + eps, this should give the correct results.

The code with the workaround looks like the following:

%% set geometry up and solve 
model = createpde('thermal','steadystate'); 
R1=[3 4 0 1 1 0 0 0 1 1]'; 
R2=[3 4 0.5 1 1 0.5 0.5 0.5 1 1]'; 
gm=[R1,R2]; 
sf='R1+R2'; 
ns=char('R1','R2') 
ns=ns'; 
g=decsg(gm,sf,ns); 
figure(1) 
pdegplot(g,'FaceLabels','on','EdgeLabels','on')

%% create mesh 
geometryFromEdges(model,g) 
generateMesh(model,'Hmax',0.025); 
figure(2) 
pdemesh(model);

%% set boundary conditions 
temperatureTop=@(location,~) 1*location.x 
thermalBC(model,'Edge',[1],'Temperature',0); 
thermalBC(model,'Edge',[7:8],'Temperature',temperatureTop); 
thermalBC(model,'Edge',[2],'HeatFlux',0); 
thermalBC(model,'Edge',[5:6],'HeatFlux',0);

%Here you can find the function handles for the ThermalConductivity
k1=@(location,~) 1*location.x; 
k2=@(location,~) 20*location.x+eps; %here you see the fix with 'eps'

thermalProperties(model,'Face',2,'ThermalConductivity',k2); 
thermalProperties(model,'Face',1,'ThermalConductivity',k1);

%% solve model 
result = solve(model); 
temp = result.Temperature;

[qx,qy] = result.evaluateHeatFlux;

figure(3) 
pdeplot(model,'XYData',temp,'Contour','on'); 
set(gca,'box','on')

This is the resulting plot:

Our development team is currently working on a fix for a future release. It seems to be that this will be fixed in one of the next releases after R2019a.

MATLAB: Does the Partial Differential Equation Toolbox provide the functionality of “assembleFEMatrices” for the new “modal-solid” workflow introduced in R2018a

For the vertical workflow using "createpde('structural','modal-solid')", the Partial Differential Equation (PDE) Toolbox does not provide the functionality of "assembleFEMatrices". As an alternative, you could solve the 'modal-solid' problem by setting up a PDE model with "createpde()", where you can get access to the matrices with all the checks and internal validatation in place. You can think of PDE model workflow as a super workflow which can also solve ‘modal-solid’ systems. However, the setup is different because of the generality of the approach.

This example shows you how to compute the modes of a cantilever beam using the PDEModel workflow. In the end, you get access to the matrices and can you use them for further calculations, e.g. finding the Frequency Response Function.

% Create PDEModel to solve for 3 PDEs of equilibrium for a solid body
PDEmodel = createpde(3);
% Create a cantilever beam geometry
PDEmodel.Geometry = multicuboid(0.45, 0.02, 0.003);
% Plot geometry
figure(1)
    pdegplot(PDEmodel, 'FaceLabels', 'on'); 
    title('Geometry');
    % NOTE: the beam will be fixed on the left face ([5])
% Generate mesh
PDEmodel.generateMesh;
% Plot mesh
figure(2)
    pdemesh(PDEmodel);
    title('Mesh');
% Set material properties (steel)
E = 210E9;
nu = 0.3;
rho = 7800;
% Compute c-coefficients in the required form using material properties
c = elasticityC3D(E,nu);
PDEmodel.specifyCoefficients('m', rho, 'd', 0, 'c', c, 'a', 0, 'f', [0;0;0]);
% Apply boundary conditions:
% % B.C. 1: fix beam on left face
PDEmodel.applyBoundaryCondition('dirichlet', 'Face', 5, 'u', [0;0;0]);
% % B.C. 2: one face is applied with a traction
%applyBoundaryCondition(model,'neumann','Face',3,'g',[0;0;1000]); 
% NOTE: line is commented out, as this BC is ineffective as loads are not 
%       accounted for in modal analysis.
% Solve for modes and frequency
modalSol = PDEmodel.solvepdeeig([0,3E5]);
% Plot mode shapes
figure(3)
subplot(2,1,1)
    pdeplot3D(PDEmodel, 'ColorMapData', modalSol.Eigenvectors(:,3,1))
    title(['Z-Disp, Mode = 1, frequency = ' num2str(sqrt(modalSol.Eigenvalues(1))/2/pi) ' Hz'])
subplot(2,1,2)
    pdeplot3D(PDEmodel, 'ColorMapData', modalSol.Eigenvectors(:,3,2))
    title(['Z-Disp, Mode = 2, frequency = ' num2str(sqrt(modalSol.Eigenvalues(2))/2/pi) ' Hz'])
% Get access to the matrices used in modal solution
FEmodel = assembleFEMatrices(PDEmodel)

Best Answer

Related Solutions

MATLAB: Using a function handle for ThermalConductivity is not changing resulting temperature plot

MATLAB: Does the Partial Differential Equation Toolbox provide the functionality of “assembleFEMatrices” for the new “modal-solid” workflow introduced in R2018a

Related Question