Rewriting the quench protection by temperature rise routine #1048

Open jonmaddock opened 4 years ago

In GitLab by @jlion on May 26, 2020, 17:28

Quench Protection `protect`

As the current routine protect (quench protection) is missing a documentation (?) I would like to propose a new more transparent implementation of basically the same formalism as before, but with clearer documentation and clearer assumptions.

Assumptions

Equilibrated heat in the cable space during quench
Heat loss during quench only goes into helium and copper fractions
No external heat loss during quench (not even transported away by the helium fraction)
Exponential decay of $J(t)$ with a decay time $\tau$ and a delay time $t_d$ (this is strictly speaking inconsistent with the rest)
A helium pressure of 6bar

Equations

In thermal equilibrium and without losses the heat produced by the copper resistivity during a quench is equal to the heat needed to rise the temperature in the material by $dT$:

dQ_{heat} = dQ_{temp}

Assuming the materials in the winding pack are equilibrated this is

P(t)dt = \sum_i c_i\rho_i V_i \,dT,

where $P$ is the power deposited by the copper current in time $t$, $i$ runs over all winding pack materials and $V_i$ stands for the volume of those materials. With $p=J^2 \eta$, where $\eta$ is the resistivity of copper and $J$ the current in the copper, it is:

J(t)^2dt = \sum_i \frac{c_i \rho_i}{\eta_{Cu}(T)} \frac{V_i}{V_{Cu}} \,dT,

Now, the quench restriction is to impose

\int J(t)^2dt \overset{!}{<} \int_{T_{op}}^{T_{max}}\sum_i \frac{c_i \rho_i}{\eta_{Cu}(T)} f_i \,dT.

The integral on the left hand side runs over the whole quench time while the integral on the right hand side goes from the operation temperature $T_{op}$ to a maximal $T_{max}$. The difference $T_{max}-T_{op}$ is usually chosen in the order of 150 K. The sum on the right hand side goes over all materials with fractions $f_i$.

The current density decay shall be written in the form:

J(t) = \begin{cases}
    J_0, & \text{if $t<t_d$}.\\
    J_0 e^{-(t-t_{det})\over \tau_{dump}}, & \text{otherwise}.
  \end{cases}

(Note again that this is inconsistent and an assumption.)

Given this form the integral is $\int J(t)^2 dt = J^2 (\frac{1}{2} \tau_{dump} + t_d)$.

Then:

J^2 \left(\frac{1}{2} \tau_{dump} +  t_d\right)   <  q_{Cu} +  \frac{V_{He}}{V_{Cu}} q_{He} +  \frac{V_{scu}}{V_{Cu}}q_{scu}

with

q_{Cu} \equiv \int_{T_0}^{T_{max}} \frac{\rho_{Cu} c_{Cu}(T)}{\eta_{Cu}(T)} \,dT\\
q_{He} \equiv \int_{T_0}^{T_{max}} \frac{\rho_{He}(T) c_{He}(T)}{\eta_{Cu}(T)} \,dT\\
q_{scu} \equiv \int_{T_0}^{T_{max}} \frac{\rho_{scu} c_{scu}}{\eta_{Cu}(T)} \,dT

and using

V_{Cu} = V_{conduit} \,(1-f_{He}) \, f_{Cu}\\
V_{He} = V_{conduit} \,f_{He}\\
V_{scu} = V_{conduit} \,(1-f_{He}) \, (1-f_{Cu})

one ends up with:

J_{cu}    < \sqrt{\frac{1}{(\frac{1}{2} \tau_{dump} +  t_d)} \Big(q_{cu} +  \frac{f_{He}}{(1-f_{He})f_{Cu}} q_{He} +  \frac{1-f_{Cu}}{f_{Cu}}q_{scu}\Big)}

To rewrite this as the winding pack current density and with $1-f_{He}=f_{cond}$ and $J_{WP} = J_{Cu} f_{Cu} f_{cond} (1-f_{case})$:

J_{WP}    < (1-f_{case}) \sqrt{\frac{1}{(\frac{1}{2} \tau_{dump} +  t_d)} \Big( f_{Cu}^2 f_{cond}^2 q_{cu} +  f_{Cu} f_{cond} (1-f_{cond}) q_{He} +   f_{Cu} f_{cond}^2(1-f_{Cu})q_{scu}\Big)}

This is also the form of the previous implementation, apart from the added quench detection term $t_d$.

Prerequisites

Needed are:

Copper Resistivity $\eta_{Cu}(T)$
Copper isobaric heat capacity $c_{Cu}(T)$
Copper density $\rho_{Cu}$
Helium isobaric heat capacity $c_{He}(T)$
Helium density $\rho_{He}(T)$
SC isobaric heat capacity $c_{scu}(T)$
SC density $\rho_{scu}(T)$

I suggest using a Bloch-Grüneisen parametrization for the copper resistivity, along the lines of:

\eta_{Cu}(T[K],RRR)={1.687\cdot 10^{-8}-2.213\cdot10^{-15} RRR \over -1+RRR} m \Omega + 2.8526\cdot10^{-5} m \Omega K  {T^5 \over (343.5 K)^6}\cdot \int_0^{343.5/T}dx \left[{x^5\over (e^x-1)(1-e^{-x})}\right]

References, e.g.

https://www.copper.org/resources/properties/cryogenic/, https://en.wikipedia.org/wiki/Bloch%E2%80%93Gr%C3%BCneisen_temperature

This produces:

CopperResistivity

data from: https://nvlpubs.nist.gov/nistpubs/Legacy/MONO/nistmonograph177.pdf

And store it as a list for a certain RRR value.

I give some values for different temperatures below (helium at 6 bar):

$`T [K]`$	$`\eta_{Cu}^{RRR=100} [n\Omega\,m]`$	$`\eta_{Cu}^{RRR=1000} [n\Omega\,m]`$	$`\rho_{He}[kg/m^3]`$	$`c_{He}[J/(K kg)]`$	$`c_{Cu}[J/(K kg)]`$
3	0.170402	0.016885	152.856	2101.29	0.053032
3.28943	0.170403	0.016886	151.067	2327.67	0.0626
3.60679	0.170403	0.016886	148.769	2608.21	0.07452
3.95477	0.170404	0.016887	145.818	2951.17	0.089486
4.33632	0.170405	0.016888	141.987	3385.88	0.108421
4.75468	0.170407	0.01689	136.908	3975.82	0.132515
5.2134	0.17041	0.016893	129.934	4853.12	0.16335
5.71638	0.170415	0.016898	119.928	6271.55	0.202972
6.26789	0.170423	0.016906	104.918	8961.91	0.254002
6.8726	0.170435	0.016918	80.8147	12511.6	0.319897
7.53566	0.170454	0.016937	61.3544	10143.3	0.404955
8.26269	0.170485	0.016968	48.8461	8607.63	0.516382
9.05986	0.170534	0.017017	40.6341	7513.53	0.663479
9.93393	0.170611	0.017094	34.7607	6884.86	0.857247
10.8923	0.170733	0.017216	30.2513	6480.71	1.11221
11.9432	0.170927	0.01741	26.6352	6196.73	1.45254
13.0955	0.171234	0.017717	23.6489	5987.61	1.91343
14.3589	0.171721	0.018204	21.1289	5829.35	2.5326
15.7442	0.172492	0.018975	18.967	5707.26	3.36621
17.2632	0.173715	0.020197	17.0883	5611.62	4.51107
18.9287	0.175651	0.022134	15.4393	5535.76	6.07632
20.7549	0.178716	0.025199	13.9803	5474.94	8.20812
22.7573	0.183558	0.030041	12.6813	5425.74	11.1119
24.9529	0.191175	0.037658	11.5187	5385.65	15.0603
27.3603	0.203074	0.049557	10.4741	5352.78	20.1398
30	0.221461	0.067944	9.53248	5325.69	26.6421
32.8943	0.249439	0.095922	8.68137	5303.29	34.9171
36.0679	0.291177	0.13766	7.91047	5284.69	45.073
39.5477	0.351986	0.198468	7.21103	5269.22	57.2157
43.3632	0.438251	0.284733	6.57555	5256.33	71.2622
47.5468	0.557186	0.403669	5.99756	5245.56	87.2477
52.134	0.716421	0.562904	5.47138	5236.58	105.06
57.1638	0.923482	0.769965	4.99205	5229.07	124.485
62.6789	1.18526	1.03174	4.55515	5222.79	145.145
68.726	1.50757	1.35405	4.15677	5217.55	166.568
75.3566	1.89492	1.7414	3.79337	5213.18	188.218
82.6269	2.35043	2.19691	3.46182	5209.53	209.698
90.5986	2.87607	2.72256	3.15927	5206.5	230.602
99.3393	3.47297	3.31945	2.88313	5203.97	250.553
108.923	4.14182	3.9883	2.6311	5201.88	269.191
119.432	4.88336	4.72984	2.40104	5200.14	286.35
130.955	5.69873	5.54521	2.19104	5198.71	302.185
143.589	6.58981	6.43629	1.99935	5197.53	316.461
157.442	7.55946	7.40594	1.82437	5196.56	329.185
172.632	8.6117	8.45818	1.66466	5195.76	340.543
189.287	9.75175	9.59823	1.51887	5195.12	350.648
207.549	10.9861	10.8326	1.38581	5194.6	359.571
227.573	12.3227	12.1692	1.26437	5194.18	367.367
249.529	13.7707	13.6172	1.15353	5193.84	374.081
273.603	15.3405	15.187	1.05238	5193.57	379.825
300	17.0442	16.8907	0.960071	5193.36	384.761

Required data (q Integrals)

Only the q integrals will be useful in the function so they can be precalculated in advance and then included in the function later. Assume a maximal temperature rise of 150 K, no matter from which operating temperature. The $q$ integrals are then only dependent on the operating temperature and can be calculated by the above given formula to (RRR=100):

$`T_{op}[K]`$	$`q_{He} [s A^2/m^4]`$	$`q_{Cu} [s A^2 / m^4]`$
4	3.44562*10^16	1.08514*10^17
14	9.92398*10^15	1.12043*10^17
24	4.90462*10^15	1.12406*10^17
34	2.41524*10^15	1.0594*10^17
44	1.26368*10^15	9.49741*10^16
54	7.51617*10^14	8.43757*10^16
64	5.01632*10^14	7.56346*10^16
74	3.63641*10^14	6.85924*10^16
84	2.79164*10^14	6.28575*10^16
94	2.23193*10^14	5.81004*10^16
104	1.83832*10^14	5.40838*10^16
114	1.54863*10^14	5.06414*10^16
124	1.32773*10^14	4.76531*10^16

These values will enter the function.

Implementation

I suggest implementing these two functions in the module superconductors:

real(dp) function j_max_protect_Am2(tau_quench,t_detect,fcu,fcond,temp,acs,aturn)
!! Finds the current density limited by the protection limit
!! acs : input real : Cable space - inside area (m2)
!! aturn : input real : Area per turn (i.e.  entire cable) (m2)
!! tau_quench : input real : Dump time (sec)
!! tau_detect : input real : Quench detection time (sec)
!! fcond : input real : Fraction of cable space containing conductor
!! fcu : input real : Fraction of conductor that is copper
!! temp : input real : Operating He temperature (K)
!! j_max_protect_Am2 : output real :  Winding pack current density from temperature
!! rise protection (A/m2)
!!
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
use maths_library, only: find_y_nonuniform_x

implicit none

real(dp), intent(in) :: tau_quench, t_detect, fcu, fcond, temp, acs, aturn
real(dp), dimension(13) :: temp_k, q_he_array_sA2m4, q_cu_array_sA2m4
real(dp) :: q_cu, q_he

temp_k = (/4,14,24,34,44,54,64,74,84,94,104,114,124/)

q_cu_array_sA2m4 = (/ 1.08514d17, 1.12043d17, 1.12406d17, 1.05940d17, &
                      9.49741d16, 8.43757d16, 7.56346d16, 6.85924d16, &
                      6.28575d16, 5.81004d16, 5.40838d16, 5.06414d16, &
                      4.76531d16/)

q_he_array_sA2m4 = (/ 3.44562d16, 9.92398d15, 4.90462d15, 2.41524d15, &
                      1.26368d15, 7.51617d14, 5.01632d14, 3.63641d14, &
                      2.79164d14, 2.23193d14, 1.83832d14, 1.54863d14, &
                      1.32773d14 /)

! Interpolate to find the correct value for the temperature
q_he = find_y_nonuniform_x(temp,temp_k,q_he_array_sA2m4,13)
q_cu = find_y_nonuniform_x(temp,temp_k,q_cu_array_sA2m4,13)

! This leaves out the contribution from the superconductor for now
j_max_protect_Am2 = (acs/aturn) * sqrt(1.0d0/(0.5d0 *tau_quench + t_detect)* &
                    (fcu**2*fcond**2 * q_cu + fcu* fcond* (1.0d0-fcond)* q_he))

end function

The maximum voltage is independent of this and should be treated in a seperate function, e.g.:

real(dp) function u_max_protect_V(tfes, tdump,aio)

!! tfes : input real : Energy stored in one TF coil (J)
!! tdump : input real : Dump time (sec)
!! aio : input real : Operating current (A)
!! u_max_protect_V : output real: Discharge voltage imposed on a TF coil (V)
!!
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none

real, intent(in) :: tfes, tdump,aio

!  Dump voltage
u_max_protect_V = 2.0D0 * tfes/(tdump*aio)

end function

[ ] This assumes, as we do currently, that the dump resistor has a resistance much higher than that of the coil during the quench. I am not certain if I have ever verified this assumption. It probably could be removed with a little work.

Todo

I have only tested this for a few values but it seems to be a bit more optimistic than the previous implementation which might be due to the used copper RRR value.

Todo:

[ ] Decide whether to use
[ ] Test
[ ] Implement
[ ] Run Test Suite

In GitLab by @jlion on May 26, 2020, 17:31

@jmorris-uk, @mkovari please have a look.

In GitLab by @mkovari on May 27, 2020, 16:44

RRR is an issue. Transmutation will gradually introduce impurities in the copper, and radiation damage will gradually cause lattice defects. Both factors will reduce the RRR during the lifetime of the machine.

For normal superconducting wire, 100 seems common. This should certainly be our highest assumption, but maybe we need to lower it for the reasons above.

Magnetoresistance will also be significant, increasing the resistance of the copper further.

Properties_of_Copper_and_Copper_Alloys_nistmonograph177.pdf

I have added a few queries in-line - in this format:

[ ] Comment

In GitLab by @jlion on May 28, 2020, 09:04

[x] Why do you need to approximate rather than integrate exactly?

The integral is exact when taken the form as given and integrated to "infinity"

[x] Is $\tau_{dump}$ the same as $\tau$?

yes

[x] $1.687\cdot 10^{-8}-2.213\cdot10^{-15} RRR$ doesn't look right. The second term is always negligible.

No, its the term that takes the RRR terms into account at low temperatures, I added a plot.

In GitLab by @jlion on May 28, 2020, 09:16

[ ] This assumes, as we do currently, that the dump resistor has a resistance much higher than that of the coil during the quench. I am not certain if I have ever verified this assumption. It probably could be removed with a little work.

Is there any documentation for this voltage formula by the way? I am not so sure how to arrive at this form.

ukaea / PROCESS