(self, test_data_path)
| 914 | _mech = "sofc.yaml" |
| 915 | |
| 916 | def test_sofc(self, test_data_path): |
| 917 | mech = self._mech |
| 918 | T = 1073.15 # T in K |
| 919 | P = ct.one_atm |
| 920 | TPB_length_per_area = 1.0e7 # TPB length per unit area [1/m] |
| 921 | |
| 922 | def newton_solve(f, xstart, C=0.0): |
| 923 | """ Solve f(x) = C by Newton iteration. """ |
| 924 | x0 = xstart |
| 925 | dx = 1.0e-6 |
| 926 | |
| 927 | n = 0 |
| 928 | while True: |
| 929 | n += 1 |
| 930 | f0 = f(x0) - C |
| 931 | x0 -= f0/(f(x0 + dx) - C - f0)*dx |
| 932 | if n > 1000: |
| 933 | raise Exception('No convergence in Newton solve') |
| 934 | if abs(f0) < 0.00001: |
| 935 | return x0 |
| 936 | |
| 937 | # Anode-side phases |
| 938 | tpb_a = ct.Interface(mech, "tpb") |
| 939 | anode_surf = tpb_a.adjacent["metal_surface"] |
| 940 | gas_a = anode_surf.adjacent["gas"] |
| 941 | oxide_surf_a = tpb_a.adjacent["oxide_surface"] |
| 942 | oxide_a = oxide_surf_a.adjacent["oxide_bulk"] |
| 943 | anode_bulk = tpb_a.adjacent["metal"] |
| 944 | |
| 945 | # Cathode-side phases |
| 946 | tpb_c = ct.Interface(mech, "tpb") |
| 947 | cathode_surf = tpb_c.adjacent["metal_surface"] |
| 948 | gas_c = cathode_surf.adjacent["gas"] |
| 949 | oxide_surf_c = tpb_c.adjacent["oxide_surface"] |
| 950 | oxide_c = oxide_surf_c.adjacent["oxide_bulk"] |
| 951 | cathode_bulk = tpb_c.adjacent["metal"] |
| 952 | |
| 953 | kElectron_a = tpb_a.kinetics_species_index("electron") |
| 954 | def anode_curr(E): |
| 955 | anode_bulk.electric_potential = E |
| 956 | w = tpb_a.net_production_rates |
| 957 | return ct.faraday * w[kElectron_a] * TPB_length_per_area |
| 958 | |
| 959 | kElectron_c = tpb_c.kinetics_species_index("electron") |
| 960 | def cathode_curr(E): |
| 961 | cathode_bulk.electric_potential = E + oxide_c.electric_potential |
| 962 | w = tpb_c.net_production_rates |
| 963 | return -ct.faraday * w[kElectron_c] * TPB_length_per_area |
| 964 | |
| 965 | # initialization |
| 966 | gas_a.TPX = T, P, 'H2:0.97, H2O:0.03' |
| 967 | gas_a.equilibrate('TP') |
| 968 | gas_c.TPX = T, P, 'O2:1.0, H2O:0.001' |
| 969 | gas_c.equilibrate('TP') |
| 970 | |
| 971 | for p in [anode_bulk, anode_surf, oxide_surf_a, oxide_a, cathode_bulk, |
| 972 | cathode_surf, oxide_surf_c, oxide_c, tpb_a, tpb_c]: |
| 973 | p.TP = T, P |
nothing calls this directly
no test coverage detected