Archived OpenModelica forums. Posting is disabled.

Alternative forums include GitHub discussions or StackOverflow (make sure to read the Stack Overflow rules; you need to have well-formed questions)


Forgot password? | Forgot username? | Register

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;

Edited by: shig1108 - Jul-30-20 01:38:38
There are 0 guests and 0 other users also viewing this topic
You are here: