cane/main/main_FEMThermalConductionAnalysis/main_steadyStateThermalConductionAnalysis.m at master · andreas-apostolatos/cane · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
%% Licensing
%
% License:         BSD License
%                  cane Multiphysics default license: cane/license.txt
%
% Main authors:    Andreas Apostolatos
%                  Marko Leskovar
%
%% Script documentation
%
% Task : Steady-state thermal conduction analysis
%
% Date : 14.04.2020
%
%% Preamble
clear;
clc;
close all;

%% Includes

% Add general math functions
addpath('../../generalMath/');

% Add all functions related to parsing
addpath('../../parsers/');

% Add all functions related to the low order basis functions
addpath('../../basisFunctions/');

% Add the equation system solvers
addpath('../../equationSystemSolvers/');

% Add all the efficient computation functions
addpath('../../efficientComputation/');

% Add all functions related to body forces
addpath('../../FEMPlateInMembraneActionAnalysis/loads/');

% Add all functions related to heat transfer analysis
addpath('../../FEMThermalConductionAnalysis/solvers/',...
        '../../FEMThermalConductionAnalysis/solutionMatricesAndVectors/',...
        '../../FEMThermalConductionAnalysis/loads/',...
        '../../FEMThermalConductionAnalysis/graphics/',...
        '../../FEMThermalConductionAnalysis/output/',...
        '../../FEMThermalConductionAnalysis/postprocessing/');

%% Parse data from GiD input file

% Define the path to the case
pathToCase = '../../inputGiD/FEMThermalConductionAnalysis/';
caseName = 'steadyStateSquareCavity';
% caseName = 'steadyStateWallConduction';
% caseName = 'rectangularPlateWithTwoHoles';
% caseName = 'rectangularPlateWithCenterHole';
% caseName = 'trapezoidalPlateHeatFlux';
% caseName = 'rectangularPlateHeatFlux';

% Parse the data from the GiD input file
[strMsh, homDOFs, inhomDOFs, valuesInhomDOFs, propNBC, propAnalysis, ...
    parameters, propNLinearAnalysis, ~, propGaussInt] = ...
    parse_ThermalModelFromGid(pathToCase, caseName, 'outputEnabled');

%% GUI

% On the body force computation
computeBodyForces = 'undefined';

% Choose solver for the linear equation system
solve_LinearSystem = @solve_LinearSystemMatlabBackslashSolver;
% solve_LinearSystem = @solve_LinearSystemGMResWithIncompleteLUPreconditioning;

% On the writing the output function
propVTK.isOutput = true;
propVTK.writeOutputToFile = @writeOutputFEMThermalConductionAnalysisToVTK;
propVTK.VTKResultFile = 'undefined';

% Choose computation of the stiffness matrix
computeStiffMtxLoadVct = @computeStiffMtxAndLoadVctFEMThermalConductionAnalysisCST;

% Linear analysis
propHeatDynamics = 'undefined';

% Initialize graphics index
graph.index = 1;

%% Assign heat flux (load)
% computeConstantFlux
if strcmp(caseName, 'rectangularPlateWithTwoHoles')
    propNBC.flux = 5e4;
elseif strcmp(caseName,'trapezoidalPlateHeatFlux') || strcmp(caseName,'rectangularPlateHeatFlux')
    propNBC.flux = 100;
end

%% Output data to a VTK format
pathToOutput = '../../outputVTK/FEMThermalConductionAnalysis/';

%% Initialize solution
numNodes = length(strMsh.nodes(:,1));
numDOFs = numNodes;
dHat = zeros(numDOFs,1);

%% Solve the plate in membrane action problem
[dHat, WComplete, minElSize] = solve_FEMThermalConductionSteadyState ...
    (propAnalysis, strMsh, dHat, homDOFs, inhomDOFs, valuesInhomDOFs, ...
    propNBC, computeBodyForces, parameters, computeStiffMtxLoadVct, ...
    solve_LinearSystem, propNLinearAnalysis, propGaussInt, propVTK, ...
    caseName, pathToOutput, 'outputEnabled');

%% Define a function to compute analytical results
if strcmp(caseName, 'steadyStateSquareCavity')
    propPostproc.computeAnalytical = @(x,y,t,propPostproc) ...
        propPostproc.T1+(propPostproc.T2-propPostproc.T1) * (2/pi) * ...
        sum( ((( (-1).^( (1:propPostproc.k) +1) )+1)./ (1:propPostproc.k) ) .* ...
        sin(( (1:propPostproc.k) *pi*x)/propPostproc.width) .* ...
        ( sinh(( (1:propPostproc.k) *pi*y)/propPostproc.width) ./ ...
        sinh(( (1:propPostproc.k) *pi*propPostproc.height)/propPostproc.width) ) );

    % Assign temperatures
    propPostproc.T1 = min(valuesInhomDOFs);
    propPostproc.T2 = max(valuesInhomDOFs);
end

%% Visualize analytical solution
if strcmp(caseName, 'steadyStateSquareCavity')
    % At steady state time is infinite
    t = inf;
    graph.index = plot_analyticalSolutionAtTimeInstance...
        (strMsh,t,propPostproc,graph,'outputEnabled');
end

%% END OF THE SCRIPT