stability improvements for v0.0.3 postfix 001

VERSION

@@ -1 +1 @@
-__version__ = "0.0.3 master"
+__version__ = "0.0.3 dev"

doc/conf.py

@@ -63,9 +63,9 @@
 # built documents.
 #
 # The short X.Y version.
-version = '0.0.3'
+version = '0.0.3-001'
 # The full version, including alpha/beta/rc tags.
-release = '0.0.3 master'
+release = '0.0.3-001 master'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

examples/chp.py

@@ -2,13 +2,14 @@
 
 from tespy import con, cmp, nwk
 import matplotlib.pyplot as plt
+import pandas as pd
 
 # %% network
 
 fluids = ['water']
 
 nw = nwk.network(fluids=fluids, p_unit='bar', T_unit='C', h_unit='kJ / kg',
-                 p_range=[0.02, 110], T_range=[20, 800], h_range=[15, 5000])
+                 p_range=[0.01, 110], T_range=[20, 800], h_range=[15, 5000])
 
 # %% components
 
@@ -83,8 +84,8 @@
 turbine_hp.set_attr(eta_s=0.9)
 turbine_lp.set_attr(eta_s=0.9)
 
-condenser.set_attr(pr1=0.95, pr2=0.95, ttd_u=12)
-preheater.set_attr(pr1=0.95, pr2=0.99, ttd_u=7)
+condenser.set_attr(pr1=0.99, pr2=0.99, ttd_u=12)
+preheater.set_attr(pr1=0.99, pr2=0.99, ttd_u=5)
 vessel.set_attr(mode='man')
 
 pump.set_attr(eta_s=0.8)
@@ -124,55 +125,61 @@
 
 nw.solve(init_file=None, mode=mode)
 nw.print_results()
-nw.save('chp_' + mode)
+nw.save('chp')
 
-file = 'chp_' + mode + '_results.csv'
+file = 'chp/results.csv'
 mode = 'offdesign'
 
 # representation of part loads
 P_range = [5.25e6, 5e6, 4.5e6, 4e6, 3.5e6, 3e6]
 
 # temperatures for the heating system
-T_fl = [80, 90, 100, 110, 120]
+T_range = [120, 110, 100, 95, 90, 85, 80, 77.5, 75, 72.5, 70]
 
-P = {}
-Q = {}
+df = pd.DataFrame(columns=P_range)
 
 # iterate over temperatures
-for i in T_fl:
-    cw_out.set_attr(T=i)
-    P[i] = []
-    Q[i] = []
+for T in T_range:
+    cw_out.set_attr(T=T)
+    Q = []
     # iterate over power plant power output
-    for j in P_range:
-        power_bus.set_attr(P=j)
+    for P in P_range:
+        print(T, P)
+        power_bus.set_attr(P=P)
 
         # use an initialisation file with parameters similar to next
         # calculation
-        if j == P_range[0]:
-            init_file = file
+        if T == T_range[0]:
+            nw.solve(init_file=file, design_file=file, mode=mode)
         else:
-            init_file = 'chp_' + mode + '_results.csv'
+            nw.solve(init_file='chp_' + str(P/1e6) + '/results.csv',
+                     design_file=file, mode=mode)
 
-        nw.solve(init_file=init_file, design_file=file, mode=mode)
-        nw.save('chp_' + mode)
-        P[i] += [power_bus.P]
-        Q[i] += [heat_bus.P]
+        nw.save('chp_' + str(P/1e6))
+        Q += [heat_bus.P]
 
+    df.loc[T] = Q
+
+df.to_csv('chp.csv')
 # plotting
-colors = ['#00395b', '#74adc1', '#bfbfbf', '#b54036', '#ec6707']
+df = pd.read_csv('chp.csv', index_col=0)
+
+colors = ['#00395b', '#74adc1', '#b54036', '#ec6707',
+          '#bfbfbf', '#999999', '#010101']
 
 fig, ax = plt.subplots()
-j = 0
-for i in T_fl:
-    plt.plot(Q[i], P[i], '.-', Color=colors[j],
-             label='$T_{VL}$ = '+str(i)+' °C', markersize=15, linewidth=2)
-    j += 1
+
+i = 0
+for T in T_range:
+    if T % 10 == 0:
+        plt.plot(df.loc[T], P_range, '-x', Color=colors[i],
+                 label='$T_{VL}$ = ' + str(T) + ' °C', markersize=7, linewidth=2)
+        i += 1
+
 ax.set_ylabel('$P$ in W')
 ax.set_xlabel('$\dot{Q}$ in W')
 plt.title('P-Q diagram for CHP with backpressure steam turbine')
 plt.legend(loc='lower left')
-
 ax.set_ylim([0, 6e6])
 ax.set_xlim([0, 14e6])
 

examples/simple_clausius_rankine.py

@@ -76,12 +76,12 @@
 
 mode = 'design'
 
-file = mode + '_results.csv'
+file = 'cr'
 
 # solve the network, print the results to prompt and save
 nw.solve(mode=mode)
 nw.print_results()
-nw.save('design')
+nw.save(file)
 
 # change to offdesign mode
 mode = 'offdesign'
@@ -91,6 +91,7 @@
 
 # the design file holds the information on the design case
 # initialisation from previously design process
-nw.solve(mode=mode, design_file=file, init_file=file)
+nw.solve(mode=mode, design_file=file + '/results.csv',
+         init_file=file + '/results.csv')
 nw.print_results()
-nw.save('offdesign')
+nw.save('cr_OD')

examples/tutorial/step_1.py

@@ -67,5 +67,6 @@
 cons.set_attr(Q=-200e3)
 
 nw.solve('offdesign',
-         init_file='condenser_results.csv', design_file='condenser_results.csv')
+         init_file='condenser/results.csv',
+         design_file='condenser/results.csv')
 nw.print_results()

examples/tutorial/step_3.py

@@ -156,6 +156,6 @@
 
 cons.set_attr(Q=-200e3)
 
-nw.solve('offdesign', init_file='heat_pump_results.csv',
-         design_file='heat_pump_results.csv')
+nw.solve('offdesign', init_file='heat_pump/results.csv',
+         design_file='heat_pump/results.csv')
 nw.print_results()

setup.py

@@ -7,7 +7,7 @@ def read(fname):
 
 
 setup(name='TESPy',
-      version='0.0.3',
+      version='0.0.3-001',
       description='Thermal Engineering Systems in Python (TESPy)',
       url='http://github.com/oemof/tespy',
       author='Francesco Witte',

tespy/components/components.py

@@ -2249,14 +2249,17 @@ def convergence_check(self, nw):
         """
         i, o = nw.comps.loc[self].i, nw.comps.loc[self].o
 
-        if i[0].p.val_SI < o[0].p.val_SI and not o[0].p.val_set:
+        if i[0].p.val_SI <= 1e5 and not i[0].p.val_set:
+            i[0].p.val_SI = 1e5
+
+        if i[0].p.val_SI <= o[0].p.val_SI and not o[0].p.val_set:
             o[0].p.val_SI = i[0].p.val_SI / 2
 
         if i[0].h.val_SI < 10e5 and not i[0].h.val_set:
             i[0].h.val_SI = 10e5
 
-        if o[0].h.val_SI < 8e5 and not o[0].h.val_set:
-            o[0].h.val_SI = 8e5
+        if o[0].h.val_SI < 5e5 and not o[0].h.val_set:
+            o[0].h.val_SI = 5e5
 
         if i[0].h.val_SI <= o[0].h.val_SI and not o[0].h.val_set:
             o[0].h.val_SI = i[0].h.val_SI * 0.75
@@ -3515,7 +3518,7 @@ def convergence_check(self, nw):
         for o in nw.comps.loc[self].o:
             fluids = [f for f in o.fluid.val.keys() if not o.fluid.val_set[f]]
             for f in fluids:
-                if f not in [self.o2, self.co2, self.h2o, self.fuel]:
+                if f not in [self.o2, self.co2, self.h2o, self.fuel.val]:
                     m_f = 0
                     for i in nw.comps.loc[self].i:
                         m_f += i.fluid.val[f] * i.m.val_SI
@@ -3541,7 +3544,7 @@ def convergence_check(self, nw):
                     if o.fluid.val[f] < 0.02:
                         o.fluid.val[f] = 0.02
 
-                elif f == self.fuel:
+                elif f == self.fuel.val:
                     if o.fluid.val[f] > 0:
                         o.fluid.val[f] = 0
 
@@ -3728,13 +3731,14 @@ def outlets(self):
         return ['out1']
 
     def attr(self):
-        return ['fuel', 'fuel_alias', 'air', 'air_alias', 'path', 'lamb', 'ti']
+        return ['fuel', 'fuel_alias', 'air', 'air_alias', 'path',
+                'lamb', 'ti', 'S']
 
     def attr_prop(self):
         return {'fuel': dc_cp(), 'fuel_alias': dc_cp(),
                 'air': dc_cp(), 'air_alias': dc_cp(),
                 'path': dc_cp(),
-                'lamb': dc_cp(), 'ti': dc_cp()}
+                'lamb': dc_cp(), 'ti': dc_cp(), 'S': dc_cp()}
 
     def fuels(self):
         return ['methane', 'ethane', 'propane', 'butane',
@@ -4594,7 +4598,8 @@ class heat_exchanger_simple(component):
       factor is used
     - D: diameter of the pipes, :math:`[D]=\text{m}`
     - L: length of the pipes, :math:`[L]=\text{m}`
-    - ks: pipes roughness
+    - ks: pipes roughness , :math:`[ks]=\text{m}` in case of darcy friction,
+      :math:`[ks]=\text{1}` in case of hazen williams
     - kA: area independent heat transition coefficient,
       :math:`[kA]=\frac{\text{W}}{\text{K}}`
     - t_a: ambient temperature, provide parameter in network's temperature unit
@@ -4909,7 +4914,7 @@ def hw_func(self, inl, outl):
 
             0 = p_{in} - p_{out} - \frac{10.67 \cdot \dot{m}_{in} ^ {1.852}
             \cdot L}{ks^{1.852} \cdot D^{4.871}} \cdot g \cdot
-            \frac{v_{in} + v_{out}}{2}^{0.852}
+            \left(\frac{v_{in} + v_{out}}{2}\right)^{0.852}
 
             \text{note: g is set to } 9.81 \frac{m}{s^2}
         """

tespy/networks.py

@@ -1041,7 +1041,6 @@ def solve_loop(self):
 
         **Improvememts**
         """
-
         print('iter\t| residual')
         for self.iter in range(self.max_iter):
 
@@ -1099,6 +1098,9 @@ def solve_loop(self):
                                                len(self.conns),))
                     break
 
+            if self.lin_dep:
+                break
+
     def solve_control(self):
         """
         calculation step of newton algorithm
@@ -1121,34 +1123,15 @@ def solve_control(self):
         self.solve_connections()
         self.solve_busses()
 
-        lin_dep = True
+        self.lin_dep = True
         try:
             self.vec_z = inv(self.mat_deriv).dot(-np.asarray(self.vec_res))
-            lin_dep = False
+            self.lin_dep = False
         except:
             pass
 
         # check for linear dependency
-        if lin_dep:
-            if self.iter > 0 and self.num_restart < 3:
-                for c in self.conns.index:
-                    c.m.val_SI = c.m.val0
-                    c.p.val_SI = c.p.val0
-                    c.h.val_SI = c.h.val0
-                    c.fluid.val = c.fluid.val0.copy()
-
-                self.relax *= 0.8
-                self.num_restart += 1
-                self.iter = 0
-                self.convergence[0] = np.zeros((len(self.conns), 0))
-                self.convergence[1] = np.zeros((len(self.conns), 0))
-                self.convergence[2] = np.zeros((len(self.conns), 0))
-                self.res = np.array([])
-                self.vec_res = []
-                print('Adjusting relaxation factors and'
-                      'restarting calculation.')
-                self.solve_loop()
-            else:
+        if self.lin_dep:
                 msg = ('error calculating the network:\n'
                        'singularity in jacobian matrix, possible reasons are\n'
                        '-> given Temperature with given pressure in two phase '
@@ -1160,7 +1143,8 @@ def solve_control(self):
                        '-> bad starting value for fuel mass flow of '
                        'combustion chamber, provide small (near to zero, '
                        'but not zero) starting value.')
-                raise hlp.MyNetworkError(msg)
+                print(msg)
+                return
 
         # add increment
         i = 0
@@ -1189,11 +1173,13 @@ def solve_control(self):
             i += 1
 
         # check properties for consistency
-        if self.iter < 3:
+        if self.iter < 5:
+            for c in self.conns.index:
+                self.solve_check_properties(c)
+
             for cp in self.comps.index:
                 cp.convergence_check(self)
 
-        if self.iter < 3:
             for c in self.conns.index:
                 self.solve_check_properties(c)
 
@@ -1228,7 +1214,8 @@ def solve_check_properties(self, c):
         # acutal value
         # for pure fluids:
         # obtain maximum temperature from fluid properties directly
-        if c.T.val_set and not c.h.val_set:
+        ############ this part seems to cause bad convergence ############
+        if c.T.val_set and not c.h.val_set and not c.p.val_set:
             self.solve_check_temperature(c, 'min')
             self.solve_check_temperature(c, 'max')