- Index
- » Usage and Applications
- » OpenModelica Usage and Applications
- » Finite Volume Code - Structurally...
Finite Volume Code - Structurally Singular
Finite Volume Code - Structurally Singular
The code I'm trying to run (shown below) outputs an error of being structurally singular when I run it in equation mode. But when I run it in algorithm mode, it solves and generates incorrect values.
*************************************************************************************************************************
model FVM_Base_v2
import Modelica.SIunits;
//Configuration parameters
parameter SIunits.Length Length = 8 * 0.0254; //Length of tube in m
parameter SIunits.Length ThicknessF = 0.097 * 0.0254; //Thickness of HTF in m
parameter SIunits.Length ThicknessT = 0.028*0.0254; //Thickness of Tube in m
//Nodal definition
parameter Integer zf = 10;//number of nodes in axial direction for the HTF
parameter Integer rf = 5;//number of nodes in radial direction for the HTF
parameter Integer zt=10; //number of nodes in axial direction for the tube
parameter Integer rt=5; //number of nodes in radial direction for the tube
//Fluid properties (water)
parameter SIunits.Density rho_f = 996.5; //Fluid density
parameter SIunits.SpecificHeatCapacity cp_f = 4181; //specific heat
parameter SIunits.ThermalConductivity k_f = 0.6095; //conductivity
parameter SIunits.ThermalDiffusivity alp_f = k_f / (rho_f * cp_f); //Thermal diffusivity
//Metal properties
parameter SIunits.Density rho_t = 8027; //Fluid density
parameter SIunits.SpecificHeatCapacity cp_t = 502.1; //specific heat
parameter SIunits.ThermalConductivity k_t = 16.3; //conductivity
parameter SIunits.ThermalDiffusivity alp_t = k_t/(rho_t*cp_t); //Thermal diffusivity
//Inlet velocity
parameter SIunits.Velocity u_in = 0.0319; //Inlet velocity in m/s
//Computed parameters
parameter SIunits.Length dzf = Length/zf; //Nodal spacing in the axial for HTF
parameter SIunits.Length drf = ThicknessF/rf; //Nodal spacing in the radial for HTF
parameter SIunits.Length dzt = Length/zt; //Nodal spacing in the axial for Tube
parameter SIunits.Length drt = ThicknessT/rt; //Nodal spacing in the radial for Tube
//Interface thermal conductivities
//parameter SIunits.ThermalConductivity k_ft_int = (k_f*k_t)/((k_f*drt/2+k_t*drf/2)/(drf/2+drt/2));
//parameter SIunits.ThermalConductivity k_tf_int = (k_t*k_f)/((k_t*drf/2+k_f*drt/2)/(drt/2+drf/2));
parameter SIunits.Power Q = 100;
//Define counters
Integer i,j;
//Array and Vector Definition
SIunits.Temp_K Th[zf, rf](fixed=false);//(start = fill(293, zf, rf)); //Nodal temperatures
SIunits.Temp_K Tt[zt, rt](fixed=false);//(start=fill(293,zt, rt)); //Nodal temperatures
SIunits.Temp_K Tin = 323; //Fluid inlet temperature
SIunits.Length rn[rf];
SIunits.Length rs[rf];
//Common multipliers
parameter SIunits.Frequency adz = alp_f/dzf ^ 2; //Axial common multiplier
parameter SIunits.Frequency adr = alp_f /drf ^ 2; //Radial common multiplier
parameter SIunits.Velocity a2r = alp_f/(2*drf);
parameter SIunits.Frequency ae = adz;
//Defining vectors
SIunits.Length rof[rf], rot[rt]; //vector of radial positions: Length based on number of radial nodes
SIunits.Velocity uj[rf]; //vector of radial velocities
SIunits.Frequency udz[rf];
SIunits.Frequency ao[rf];
SIunits.Frequency aw[rf];
initial equation
Th[:,:] = fill(298,zf,rf);
Tt[:,:] = fill(298,zt,rt);
equation
//*************************************************************************************************************************************************
//*************************************************************************************************************************************************
// This section of the code is for the heat transfer fluid
//*************************************************************************************************************************************************
//*************************************************************************************************************************************************
for j in 1:rf loop
rof[j] = (2*j - 1)*drf/2;
uj[j] = 2*u_in*(1 - (rof[j]/ThicknessF)^2);
udz[j] = uj[j]/dzf;
ao[j] = udz[j] + 2*adz + 2*adr;
aw[j] = udz[j] + adz;
end for;
for i in 1:zf loop
for j in 1:rf loop
rn[j] = rof[j] + drf/2;
rs[j] = rof[j] - drf/2;
end for;
end for;
//Interior cells
for j in 2:rf-1 loop
for i in 2:zf-1 loop
der(Th[i,j]) = 1/(rof[j]*drf*dzf)*(-Th[i, j]*(alp_f*(rn[j] + rs[j])*dzf/
drf + 2*alp_f*rof[j]*drf/dzf + uj[j]*rof[j]*drf) + Th[i, j + 1]*(alp_f*
rn[j]*dzf/drf) + Th[i, j - 1]*(alp_f*rs[j]*dzf/drf) + Th[i + 1, j]*(
alp_f*rof[j]*drf/dzf) + Th[i - 1, j]*(alp_f*rof[j]*drf/dzf + uj[j]*rof[
j]*drf));
end for;
end for;
//Inlet cells
i= 1;
for j in 2:rf-1 loop
der(Th[i,j]) = 1/(rof[j]*drf*dzf)*(-Th[i, j]*(alp_f*(rn[j] + rs[j])*dzf/drf +
4*alp_f*rof[j]*drf/dzf + uj[j]*rof[j]*drf) + Th[i, j + 1]*(alp_f*rn[j]*
dzf/drf) + Th[i, j - 1]*(alp_f*rs[j]*dzf/drf) + Th[i + 1, j]*(alp_f*4/3*
rof[j]*drf/dzf) + Tin*(8/3*alp_f/dzf + uj[j])*rof[j]*drf);
end for;
//Outlet cells
i= zf;
for j in 2:rf-1 loop
der(Th[i,j]) = 1/(rof[j]*drf*dzf)*(-Th[i, j]*(alp_f*(rn[j] + rs[j])*dzf/drf +
alp_f*rof[j]*drf/dzf + uj[j]*rof[j]*drf) + Th[i, j + 1]*(alp_f*rn[j]*dzf/
drf) + Th[i, j - 1]*(alp_f*rs[j]*dzf/drf) + Th[i - 1, j]*(alp_f*rof[j]*
drf/dzf + uj[j]*rof[j]*drf));
end for;
//Centerline cells
j= 1;
for i in 2:zf-1 loop
der(Th[i,j]) = 1/(rof[j]*drf*dzf)*(-Th[i, j]*(alp_f*rn[j]*dzf/drf + 2*alp_f*
rof[j]*drf/dzf + uj[j]*rof[j]*drf) + Th[i, j + 1]*(alp_f*rn[j]*dzf/drf) +
Th[i + 1, j]*(alp_f*rof[j]*drf/dzf) + Th[i - 1, j]*(alp_f*rof[j]*drf/dzf +
uj[j]*rof[j]*drf));
end for;
//Inlet bottom corner
i = 1; j = 1;
der(Th[i,j]) = 1/(rof[j]*drf*dzf)*(-Th[i, j]*(alp_f*rn[j]*dzf/drf + 4*alp_f*
rof[j]*drf/dzf + uj[j]*rof[j]*drf) + Th[i, j + 1]*(alp_f*rn[j]*dzf/drf) +
Th[i + 1, j]*(4/3*alp_f*rof[j]*drf/dzf) + Tin*(8/3*alp_f*rof[j]*drf/dzf +
uj[j]*rof[j]*drf));
//Outlet bottom corner
i = zf; j = 1;
der(Th[i,j]) = 1/(rof[j]*drf*dzf)*(-Th[i, j]*(alp_f*rn[j]*dzf/drf + alp_f*
rof[j]*drf/dzf + uj[j]*rof[j]*drf) + Th[i, j + 1]*(alp_f*rn[j]*dzf/drf) +
Th[i - 1, j]*(alp_f*rof[j]*drf/dzf + uj[j]*rof[j]*drf));
//HTF-Tube Wall interface
j = rf;
for i in 2:zf-1 loop
der(Th[i,j]) = 1/(rof[j]*drf*dzf)*(-Th[i, j]*((k_t/k_f*alp_f*rn[j]*dzf/(drf/
2 + drt/2) + alp_f*rs[j]*dzf/drf) + 2*alp_f*rof[j]*drf/dzf + uj[j]*rof[j]*
drf) + Th[i, j - 1]*(alp_f*rs[j]*dzf/drf) + Th[i + 1, j]*(alp_f*rof[j]*
drf/dzf) + Th[i - 1, j]*(alp_f*rof[j]*drf/dzf + uj[j]*rof[j]*drf) + Tt[i,
1]*(alp_f*rn[j]*k_t/k_f*dzf/(drf/2 + drt/2)));
end for;
//Inlet top corner
i = 1; j = rf;
der(Th[i,j]) = 1/(rof[j]*drf*dzf)*(-Th[i, j]*(alp_f*rn[j]*k_t/k_f*dzf/(drf/2
+ drt/2) + alp_f*rs[j]*dzf/drf + 4*alp_f*rof[j]*drf/dzf + uj[j]*rof[j]*drf)
+ Th[i, j - 1]*(alp_f*rs[j]*dzf/drf) + Th[i + 1, j]*(4/3*alp_f*rof[j]*drf/
dzf) + Tin*(8/3*alp_f/dzf + uj[j])*rof[j]*drf + Tt[i, 1]*(alp_f*rn[j]*k_t/
k_f*dzf/(drf/2 + drt/2)));
//Outlet top corner
i = zf; j = rf;
der(Th[i,j]) = 1/(rof[j]*drf*dzf)*(-Th[i, j]*(alp_f*rn[j]*k_t/k_f*dzf/(drf/2 +
drt/2) + alp_f*rs[j]*dzf/drf + alp_f*rof[j]*drf/dzf + uj[j]*rof[j]*drf) +
Th[i, j - 1]*(alp_f*rs[j]*dzf/drf) + Th[i - 1, j]*(alp_f*rof[j]*drf/dzf +
uj[j]*rof[j]*drf) + Tt[i, 1]*(alp_f*rn[j]*k_t/k_f)*dzf/(drf/2 + drt/2));
//*************************************************************************************************************************************************
//*************************************************************************************************************************************************
// This section of the code is for the tube wall
//*************************************************************************************************************************************************
//*************************************************************************************************************************************************
for j in 1:rt loop
rot[j] = ThicknessF + (2*j - 1)*drt/2;
end for;
for j in 1:rt loop
for i in 1:zt loop
rn[j] = rot[j] + drt;
rs[j] = rot[j] - drt;
end for;
end for;
//Interior cells
for j in 2:rt-1 loop
for i in 2:zt-1 loop
der(Tt[i,j]) = 1/(rot[j]*drt*dzt)*(-Tt[i, j]*(alp_t*(rn[j] + rs[j])*dzt/
drt + 2*alp_t*rot[j]*drt/dzt) + Tt[i, j + 1]*(alp_t*rn[j]*dzt/drt) + Tt[
i, j - 1]*(alp_t*rs[j]*dzt/drt) + Tt[i + 1, j]*(alp_t*rot[j]*drt/dzt) +
Tt[i - 1, j]*(alp_t*rot[j]*drt/dzt));
end for;
end for;
//Left boundary
i = 1;
for j in 2:rt-1 loop
der(Tt[i,j]) = 1/(rot[j]*drt*dzt)*(-Tt[i, j]*(alp_t*(rn[j] + rs[j])*dzt/drt +
alp_t*rot[j]*drt/dzt) + Tt[i, j + 1]*(alp_t*rn[j]*dzt/drt) + Tt[i, j - 1]*
(alp_t*rs[j]*dzt/drt) + Tt[i + 1, j]*(alp_t*rot[j]*drt/dzt));
end for;
//Right boundary
i = zt;
for j in 2:rt-1 loop
der(Tt[i,j]) = 1/(rot[j]*drt*dzt)*(-Tt[i, j]*(alp_t*(rn[j] + rs[j])*dzt/drt +
alp_t*rot[j]*drt/dzt) + Tt[i, j + 1]*(alp_t*rn[j]*dzt/drt) + Tt[i, j - 1]*
(alp_t*rs[j]*dzt/drt) + Tt[i - 1, j]*(alp_t*rot[j]*drt/dzt));
end for;
//Tube-HTF Interface
j = 1;
for i in 2:zt-1 loop
der(Tt[i,j]) = 1/(rot[j]*drt*dzt)*(-Tt[i, j]*(alp_t*rn[j]*dzt/drt + alp_t*
rs[j]*k_f/k_t*dzt/(drt/2 + drf/2) + 2*alp_t*rot[j]*drt/dzt) + Tt[i, j + 1]
*(alp_t*rn[j]*dzt/drt) + Tt[i + 1, j]*(alp_t*rot[j]*drt/dzt) + Tt[i - 1,
j]*(alp_t*rot[j]*drt/dzt) + Th[i, rf]*(alp_t*rs[j]*k_f/k_t*dzt/(drt/2 +
drf/2)));
end for;
//Tube bottom inlet corner
i = 1; j = 1;
der(Tt[i,j]) = 1/(rot[j]*drt*dzt)*(-Tt[i, j]*(alp_t*rn[j] + dzt/drt + alp_t*
rs[j]*k_f/k_t*dzt/(drt/2 + drf/2) + alp_t*rof[j]*drt/dzt) + Tt[i, j + 1]*(
alp_t*rn[j]*dzt/drt) + Tt[i + 1, j]*(alp_t*rot[j]*drt/dzt) + Th[i, rf]*(
alp_t*rs[j]*k_f/k_t*dzt/(drf/2 + drt/2)));
//Tube bottom right corner
i = zt; j = 1;
der(Tt[i,j]) = 1/(rot[j]*drt*dzt)*(-Tt[i, j]*(alp_t*rn[j]*dzt/drt + alp_t*rs[
j]*k_f/k_t*dzt/(drf/2 + drt/2)) + Tt[i, j + 1]*(alp_t*rn[j]*dzt/drt) + Tt[i -
1, j]*(alp_t*rot[j]*drt/dzt) + Th[i, rf]*(alp_t*rs[j]*k_f/k_t*dzt/(drf/2 +
drt/2)));
//Tube top wall
j = rt;
for i in 2:zt-1 loop
der(Tt[i,j]) = 1/(rot[j]*drt*dzt)*(alp_t*rn[j]*dzt*Q/k_t - Tt[i, j]*(alp_t*
rs[j]*dzt/drt + 2*alp_t*rot[j]*drt/dzt) + Tt[i, j - 1]*(alp_t*rs[j]*dzt/
drt) + Tt[i + 1, j]*(alp_t*rot[j]*drt/dzt) + Tt[i - 1, j]*(alp_t*rot[j]*
drt/dzt));
end for;
//Inlet top corner
i = 1; j = rt;
der(Tt[i,j]) = 1/(rot[j]*drt*dzt)*(alp_t*rn[j]*dzt*Q/k_t - Tt[i, j]*(alp_t*
rs[j]*dzt/drt + alp_t*rot[j]*drt/dzt) + Tt[i, j - 1]*(alp_t*rs[j]*dzt/drt) +
Tt[i + 1, j]*(alp_t*rot[j]*drt/dzt));
//Outlet top corner
i = zt; j = rt;
der(Tt[i,j]) = 1/(rot[j]*drt*dzt)*(alp_t*rn[j]*dzt*Q/k_t - Tt[i, j]*(alp_t*
rs[j]*dzt/drt + alp_t*rot[j]*drt/dzt) + Tt[i, j - 1]*(alp_t*rs[j]*dzt/drt) +
Tt[i - 1, j]*(alp_t*rot[j]*drt/dzt));
end FVM_Base_v2;
- Index
- » Usage and Applications
- » OpenModelica Usage and Applications
- » Finite Volume Code - Structurally...