A full-range formulation for dynamic loss of high-temperature superconductor coated conductors

Dynamic loss has prompted widespread concern recently in the domains related to second generation (2G) high-temperature superconductor (HTS) coated conductors (CC) [1–5]. It is generated when HTS CCs carry a direct current (DC) in an alternating current (AC) magnetic field, which is tightly relevant for machine windings and other power applications. The widely adopted analytical expression for dynamic loss per unit time, Q_dyn,l, can be written as [4]

$$\begin{eqnarray}&&{Q}_{\mathrm{dyn},l}=4wfL{I}_{{\rm{t}}}i\left({B}_{{\rm{ext}}}-{B}_{{\rm{th}}}\right),\end{eqnarray} \tag{ 1 }$$

where w is the half-width of the HTS CC, f is the frequency of the AC magnetic field, L is the length of the CC, I_t represents the transport current, and i denotes the load ratio between I_t and the self-field critical current, I_c0. B_ext is the amplitude of the external AC magnetic field, and B_th is the threshold field, defined by [6]

$$\begin{eqnarray}&&{B}_{{\rm{th}}}=\displaystyle \frac{{\mu }_{0}h{J}_{{\rm{c0}}}}{2\pi }\left[\displaystyle \frac{1}{i}\,\mathrm{ln}\left(\displaystyle \frac{1+i}{1-i}\right)+\,\mathrm{ln}\left(\displaystyle \frac{1-{i}^{2}}{4{i}^{2}}\right)\right],\end{eqnarray} \tag{ 2 }$$

where μ₀ is the free space permeability, h is the thickness of the HTS thin film, and J_c0 is the critical current density in self-field, determined by I_c0/(2wh).

The linearity of dynamic loss with respect to B_ext, as shown in (1), has been well verified by many studies [7–9]. However, recent experimental measurements have shown that, when the HTS CC with a high i is exposed to a high external magnetic field, the dynamic loss exhibits a non-linearity with the augmentation of field intensity and the CC is under the risk of quench [10–12]. A critical-state model for dynamic resistance has been proposed to describe this non-linearity in [8]. When B_ext is far greater than B_th, dynamic loss can be written as [8]

$$\begin{eqnarray}&&{Q}_{\mathrm{dyn},\mathrm{nl}}={{I}_{{\rm{t}}}}^{2}\cdot \displaystyle \frac{4wfL}{{I}_{{\rm{c0}}}}\cdot \left({B}_{{\rm{ext}}}+\displaystyle \frac{{{B}_{{\rm{ext}}}}^{2}}{{B}_{0}}\right),\end{eqnarray} \tag{ 3 }$$

where B₀ is a constant depending on the material property.

However, it appears that the measured dynamic resistance rises much faster than the analytical results obtained by (3) [8]. Therefore, for now, the existing analytical expressions are inapplicable to characterize the non-linearity of dynamic loss, and the magnetic field and current load ratio causing the sudden rise in loss are still unclear. In [9], a time-averaged DC flux flow resistance has been introduced to explain the non-linear fast rise of dynamic resistance, which provides a good description of the deviation from linearity at high currents. On the basis of this work, we have studied the current density and magnetic flux density distributions inside the HTS CC in detail, and propose a new formulation for dynamic loss, which considers a full-range of magnetic fields and load ratios. This paper helps better understand the variation properties of dynamic loss from the perspective of theoretical fundamentals and determine its influencing factors, which is significant for avoiding quench in HTS machine windings, etc.

When a HTS CC carries a DC transport current under an AC magnetic field, the transport current occupies the superconducting layer with width 2iw in the center of the CC (named dynamic region), leaving the rest with width 2(1−i)w free on both sides [10]. Therefore, the dynamic loss per unit time can be formulated by

$$\begin{eqnarray}&&{Q}_{{\rm{dyn}}}=\displaystyle \frac{hL}{T}\displaystyle {\int }_{0}^{T}\displaystyle {\int }_{(1-i)w}^{(1+i)w}E\cdot Jdydt,\end{eqnarray} \tag{ 4 }$$

where E represents the electric field, and J is the current density along the width of the coated conductor.

Considering the field dependence of the critical current density J_c (B) [8]

$$\begin{eqnarray}&&{J}_{{\rm{c}}}\left(B\right)=\displaystyle \frac{{J}_{{\rm{c0}}}}{1+\left|{B}_{\perp }\right|/{B}_{0}},\end{eqnarray} \tag{ 5 }$$

and the E-J power law

$$\begin{eqnarray}&&E={E}_{0}\cdot {\left(\displaystyle \frac{J}{{J}_{{\rm{c}}}\left(B\right)}\right)}^{n},\end{eqnarray} \tag{ 6 }$$

(4) can be transformed to

$$\begin{eqnarray}\begin{array}{lll}{Q}_{{\rm{dyn}}} & = & \displaystyle \frac{hL}{T}\displaystyle {\int }_{0}^{T}\displaystyle {\int }_{(1-i)w}^{(1+i)w}{E}_{0}{\left[\displaystyle \frac{J\cdot \left(1+\left|{B}_{\perp }\right|/{B}_{0}\right)}{{J}_{{\rm{c0}}}}\right]}^{n}Jdydt\\ & = & \displaystyle \frac{{E}_{0}hLf}{{{J}_{{\rm{c0}}}}^{n}}\displaystyle {\int }_{0}^{T}\displaystyle {\int }_{(1-i)w}^{(1+i)w}{J}^{n+1}\cdot f\left(B\right)dydt,\end{array}\end{eqnarray} \tag{ 7 }$$

with

$$\begin{eqnarray}&&f\left(B\right)={\left[1+\displaystyle \frac{\left|{B}_{{\rm{ext}}}\sin \left(2\pi ft\right)-{B}_{{\rm{s}}}\right|}{{B}_{0}}\right]}^{n},\end{eqnarray} \tag{ 8 }$$

where E₀ = 10⁻⁴ V m⁻¹, B_⊥ is the local magnetic field perpendicular to the wide surface of the CC. B_⊥ is decided by both B_ext and the self-field B_s. The negative sign indicates the direction of the field vector.

Although the current density and magnetic distribution along the width of the CC has been studied by Brandt, it is not feasible to calculate the integration in (7) directly by analytical methods [13]. However, the non-linear fast rise of dynamic loss occurs when the real load ratio I_t/I_c(B) is approaching 1, i.e., when the CC is in a critical state. In this case, we can simplify (7) by investigating the current density and magnetic flux density characteristics in the dynamic region.

Firstly, the J properties need to be studied. According to the Bean Model, when the superconductor is in a critical state, the current density distribution along the width of the HTS CC should be characterized by the critical current ±J_c. Here, concerning dynamic loss, we assume that the dynamic region is fully occupied by I_t. To verify this assumption, the T-formulation based numerical modeling method has been adopted here [14, 15]. We define the average current density during one cycle T in the dynamic region, J_avg, as

$$\begin{eqnarray}&&{J}_{{\rm{avg}}}=\displaystyle \frac{1}{T}\cdot \displaystyle \frac{1}{2iw}\cdot \displaystyle {\int }_{0}^{T}\displaystyle {\int }_{(1-i)w}^{(1+i)w}Jdydt.\end{eqnarray} \tag{ 9 }$$

During the simulation, the transport current I_t has been fixed as 60 A. When I_c0 decreases from 160 A to 80 A (i.e. i increases from 0.375 to 0.75), the relationship between i and J_avg/J_t under different B_ext is presented in figure 1. Here, J_t is the average current density determined by the transport current, with J_t = I_t/(2wh). It can be seen that with the increase of i, J_avg will get close to J_t. In other words, the dynamic region with the width of 2iw is the effective region to carry transport current and in the critical state, the dynamic region has been filled with I_t. Therefore, the current density distribution in the dynamic region should be decided by J_t.

Zoom In Zoom Out Reset image size

Next, the magnetic flux density properties in the dynamic region need to be explored. In (8), B_s at any position inside the dynamic region, y₀, is decided by

$$\begin{eqnarray}&&\begin{array}{l}{B}_{{\rm{s}}}\left({y}_{0}\right)=\displaystyle \frac{{\mu }_{0}h}{2\pi }\left[\displaystyle {\int }_{0}^{(1-i)w}\displaystyle \frac{J(y)}{({y}_{0}-y)}dy-\displaystyle {\int }_{(1-i)w}^{{y}_{0}}\displaystyle \frac{J(y)}{({y}_{0}-y)}dy\right.\\ \left.\,\,\,\,\,+\,\displaystyle {\int }_{{y}_{0}}^{2w}\displaystyle \frac{J(y)}{(y-{y}_{0})}dy\right],\end{array}\end{eqnarray} \tag{ 10 }$$

with (1−i)w ≤ y₀ ≤ (1 + i)w.

Based on the above analysis, J(y) is determined by J_t, so (10) can be simplified as

$$\begin{eqnarray}&&{B}_{{\rm{s}}}\left({y}_{0}\right)=\displaystyle \frac{{\mu }_{0}h{J}_{{\rm{t}}}}{2\pi }\cdot \,\mathrm{ln}\left(\displaystyle \frac{\left(2w-{y}_{0}\right)\cdot {y}_{0}}{{\left[{y}_{0}-(1-i)w\right]}^{2}}\right).\end{eqnarray} \tag{ 11 }$$

Then the average self-field flux density in the dynamic region can be calculated by

$$\begin{eqnarray}&&\begin{array}{l}\overline{{B}_{{\rm{s}}}}=\displaystyle \frac{1}{2iw}\displaystyle {\int }_{(1-i)w}^{(1+i)w}{B}_{{\rm{s}}}\left({y}_{0}\right)d{y}_{0}\\ \,\,=\,\displaystyle \frac{{\mu }_{0}h{J}_{{\rm{t}}}}{2\pi }\left[\displaystyle \frac{1}{i}\,\mathrm{ln}\left(\displaystyle \frac{1+i}{1-i}\right)+\,\mathrm{ln}\left(\displaystyle \frac{1-{i}^{2}}{4{i}^{2}}\right)\right]=i\cdot {B}_{{\rm{th}}}.\end{array}\end{eqnarray} \tag{ 12 }$$

According to (2), it is evident that B_th decreases with i. Therefore, at a high load rate, B_s becomes much smaller compared to B_ext, i.e. B_⊥ will be dominated by the external field. In this case, we have

$$\begin{eqnarray}&&f\left(B\right)\approx {\left[1+\displaystyle \frac{\left|{B}_{{\rm{ext}}}\sin \left(2\pi ft\right)\right|}{{B}_{0}}\right]}^{n}.\end{eqnarray} \tag{ 13 }$$

Based on the above formulas, (7) can be transformed to

$$\begin{eqnarray}\begin{array}{lll}{Q}_{\mathrm{dyn},\mathrm{nl}} & = & \displaystyle \frac{{E}_{0}hLf}{{{J}_{{\rm{c0}}}}^{n}}\displaystyle {\int }_{0}^{T}\displaystyle {\int }_{(1-i)w}^{(1+i)w}{\left(\displaystyle \frac{{I}_{{\rm{t}}}}{2wh}\right)}^{n+1}\cdot f\left(B\right)dydt\\ & = & \displaystyle \frac{{E}_{0}Lf{I}_{{\rm{t}}}}{2w}{i}^{n}\displaystyle {\int }_{0}^{T}\displaystyle {\int }_{(1-i)w}^{(1+i)w}f\left(B\right)dydt\\ & = & {E}_{0}L{I}_{{\rm{t}}}{i}^{n+1}\cdot {f}_{{\rm{avg}}}\left(B\right),\end{array}\end{eqnarray} \tag{ 14 }$$

with the average of f (B) defined as

$$\begin{eqnarray}&&{f}_{{\rm{avg}}}\left(B\right)=\displaystyle \frac{1}{T}\cdot \displaystyle \frac{1}{2iw}\cdot \displaystyle {\int }_{0}^{T}\displaystyle {\int }_{(1-i)w}^{(1+i)w}f\left(B\right)dydt,\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{Q}_{\mathrm{dyn},\mathrm{nl}}={E}_{0}L{I}_{{\rm{t}}}{i}^{n+1}\cdot \left\{\begin{array}{l}1+\displaystyle {\sum }_{p=0}^{n/2-1}\displaystyle \frac{n!}{\left(2p+1\right)!\left[n-\left(2p+1\right)\right]!}{\left(\displaystyle \frac{{B}_{{\rm{ext}}}}{{B}_{0}}\right)}^{2p+1}{\left(\displaystyle \frac{1}{2}\right)}^{2p+1}\cdot \displaystyle \frac{{2}^{3p+2}\cdot p!}{\pi \displaystyle {\prod }_{q=0}^{2p+1}\left(2q+1\right)}\\ +\,\displaystyle {\sum }_{p=0}^{n/2-1}\displaystyle \frac{n!}{\left(2p+2\right)!\left[n-\left(2p+2\right)\right]!}{\left(\displaystyle \frac{{B}_{{\rm{ext}}}}{{B}_{0}}\right)}^{2p+2}{\left(\displaystyle \frac{1}{2}\right)}^{2p+2}\cdot \displaystyle \frac{\left(2p+2\right)!}{{\left[\left(p+1\right)!\right]}^{2}}\end{array}\right\},\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\begin{array}{l}{Q}_{{\rm{dyn}}}={Q}_{\mathrm{dyn},l}+{Q}_{\mathrm{dyn},\mathrm{nl}}=4wfL{I}_{{\rm{t}}}i\left({B}_{{\rm{ext}}}-{B}_{{\rm{th}}}\right)+{E}_{0}L{I}_{{\rm{t}}}{i}^{n+1}\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\cdot \,\left\{\begin{array}{l}1+\displaystyle {\sum }_{p=0}^{n/2-1}\displaystyle \frac{n!}{\left(2p+1\right)!\left[n-\left(2p+1\right)\right]!}{\left(\displaystyle \frac{{B}_{{\rm{ext}}}}{{B}_{0}}\right)}^{2p+1}{\left(\displaystyle \frac{1}{2}\right)}^{2p+1}\cdot \displaystyle \frac{{2}^{3p+2}\cdot p!}{\pi \displaystyle {\prod }_{q=0}^{2p+1}\left(2q+1\right)}\\ +\displaystyle {\sum }_{p=0}^{n/2-1}\displaystyle \frac{n!}{\left(2p+2\right)!\left[n-\left(2p+2\right)\right]!}{\left(\displaystyle \frac{{B}_{{\rm{ext}}}}{{B}_{0}}\right)}^{2p+2}{\left(\displaystyle \frac{1}{2}\right)}^{2p+2}\cdot \displaystyle \frac{\left(2p+2\right)!}{{\left[\left(p+1\right)!\right]}^{2}}\end{array}\right\}.\end{array}\end{eqnarray} \tag{ 17 }$$

Based on binomial theorem and Euler's formula, through a series of derivations, the details of which can be found in appendix, (14) can be written as (16). It should be emphasized that Q_dyn,nl in (16) is derived on the basis that the HTS CC is at relatively high i and simultaneous high B_ext, which can only be used to characterize the non-linearly fast-rising part of the Q_dyn(B_ext) curve. In (1), Q_dyn,l is used to describe the linear correlation between dynamic loss and B_ext. There will be a cross point between the Q_dyn,l(B_ext) and Q_dyn,nl(B_ext) curves as at low external magnetic field, B_ext is much smaller than B₀, thus the n-value does not play an important role and Q_dyn,nl will approach 0. However, at high external field, B_ext becomes comparable to B₀, then under the influence of the n-value, Q_dyn will increase rapidly in the form of a power function. In short, before the cross point, Q_dyn,l is far greater than Q_dyn,nl, and vice versa. On the basis of the above analysis, we have proposed a full-range analytical equation for calculating dynamic loss, as (17). It should be noted that, in (17) n is even. When n is odd, we need to modify the upper bound of summation accordingly. To verify the correctness of (17), both numerical modeling and experimental methods have been adopted. In addition, based on the time-averaged DC flux-flow resistance term mentioned in [9], Q_dyn can also be fitted by

$$\begin{eqnarray}&&{Q}_{{\rm{dyn}}}=4wfL{I}_{{\rm{t}}}i\left({B}_{{\rm{ext}}}-{B}_{{\rm{th}}}\right)+Lf{{I}_{{\rm{t}}}}^{n+1}\displaystyle {\int }_{0}^{1/f}\displaystyle \frac{{E}_{0}}{{I}_{{\rm{c}}}{\left(\left|{B}_{\perp }\right|\right)}^{n}}dt.\end{eqnarray} \tag{ 18 }$$

(18) has also been compared with (17) in section 4.

3.1. Numerical modelling

The numerical model of the HTS CC was developed by use of finite element method based on the T-formulation [14, 15]. The diagram of the modelled HTS CC is shown as figure 2. The governing equation can be written as [4]

$$\begin{eqnarray}&&\displaystyle \frac{1}{\sigma }{{\rm{\nabla }}}^{2}T-\displaystyle \frac{{\mu }_{0}h}{2\pi }\cdot \displaystyle \frac{\partial }{\partial t}\displaystyle \sum _{en=1}^{e\,\max }\displaystyle \frac{{l}_{{\rm{en}}}{\rm{\nabla }}\times ({T}_{{\rm{en}}}^{{\rm{e}}}{\bf{n}}^{\prime} )\cdot \mathop{{\bf{x}}}\limits^{\frown {}}}{{r}_{{\rm{en}}}}-\displaystyle \frac{\partial {{\bf{B}}}_{{\rm{ext}}}}{\partial t}\cdot {\bf{n}}=0,\end{eqnarray} \tag{ 19 }$$

where l_en is the width of each element, T_en^e is the current vector potential in each element, r_en is the distance between the current source element and the calculation point, nʹ is the normal vector at the current source element, n is the normal vector at the calculation point, and σ is the equivalent conductivity determined by the E-J power law, σ = J/E.

Zoom In Zoom Out Reset image size

3.2. Experiment

Two measurement systems were used for the experimental results [3, 9]. One of the experimental setups for measuring Q_dyn is shown in figure 3. It is composed of a custom-built AC magnet, which can produce a flux density up to 100 mT, a DC power supply which provides I_t varying between 0–240 A, a nano-voltage meter for averaged voltage measuremnt across the voltage taps attached on the middle of the sample conductor, and a cryostat containing liquid nitrogen. Q_dyn is obtained by multiplying the measured voltage across the sample with the transport current.

Zoom In Zoom Out Reset image size

Two different YBCO CCs have been tested, of which the specifications are shown in table 1.

Table 1.  Specifications for tested HTS CCS.

		Value
Symbols	Parameters	Sample 1	Sample 2
w	half width of the CC	2 mm	2 mm
hHTS	YBCO film thickness	1 μm	1 μm
Ic0	critical current in self-field	105.3 A	87.67 A
It	transport current	95 A	79 A
n	n-value	23	22
f	frequency of the AC field	26.62 Hz	67.89 Hz
B0	magnetic field constant	0.135 T	0.17 T
To	operating temperature	77 K
E0	characteristic E-field	10−4 V m−1
μ0	free space permeability	4π × 10−7 H m−1

For the HTS CCs with different I_c0 varying from 80 to 160 A while carrying I_t = 60 A, they are exposed to AC external magnetic fields with increasing amplitude between 10–500 mT. The simulated Q_dyn (solid symbols) and the analytical results calculated by (1) (short dash lines), (16) (dash-dot lines), (17) (solid lines) as well as (18) (dash lines) are depicted together in figure 4.

Zoom In Zoom Out Reset image size

It shows that there exists a cross-field between the Q_dyn,l(B_ext) and Q_dyn,nl(B_ext) curves, as mentioned above. Before the cross-field, Q_dyn(B_ext) agrees well with (1), which shows a linear correlation. After the cross-field, Q_dyn(B_ext) is in good accordance with (16), which is in the form of a summation of power functions. In general, the simulated Q_dyn(B_ext) agrees well with (17) in the full magnetic field range. Although (18) also describes the non-linearity of Q_dyn, it appears to be much larger in terms of value. (17) has been derived through a rigorous process, providing a more accurate result. Furthermore, (18) is given based on the estimation of the time-average DC flux flow resistance, which cannot intuitively account for the non-linearity of Q_dyn in the form of a summation of power functions. With regard to the cross-field, it increases inversely with i. In fact, for a lower i the CC has a larger capacity to withstand the induced current by the external magnetic field, thus it will be harder for I_t to go beyond the real critical current I_c(B). Therefore, the cross-field increases with I_c0 for a fixed I_t, which complies with the J_c(B) dependence and the conclusions drawn in [4].

The simulated Q_dyn based on T-formulation, the measured experimental data, and the analytical results obtained by (1), (3), (16) and (17) for the two tested CC samples are presented together in figure 5. It can be found that the experimental data is in good agreement with the simulation results and the proposed equation (17). The Q_dyn(B_ext) curve exhibits an evident non-linearity with the increase of the external magnetic field, which is significantly different from the previously widely-adopted analytical equations (1) and (3). It should be pointed out that the relatively poor agreement between the simulation and experiment at low magnetic fields for the CC sample with I_c0 = 87.67 A is due to its effective width, which is different from its physical width [3]. However, it does not affect the validity of (17), in that it describes well the non-linearity of Q_dyn under high magnetic fields.

Zoom In Zoom Out Reset image size

With respect to computational complexity, analytical methods generally require far less computing time and resources than the finite element method (FEM) based numerical modelling. In our case, (17) can be simply calculated in a few seconds using MATLAB. However, it requires more time and effort to build and run a FEM model. For example, the computational time to obtain Q_dyn of a single HTS layer with the T-formulation based numerical model is tens of seconds or longer, depending on the solving accuracy and number of mesh elements. Therefore, the proposed formulation (17) is a more efficient and convenient way to calculate dynamic loss compared with numerical models.

Based on the J_c(B) dependence and E-J power law, through investigation of the characteristics of current density and magnetic flux density distributions in the dynamic region of HTS CCs, a new analytical expression for dynamic loss has been proposed. This expression considers a full-range of magnetic fields and load ratios. At low external fields and simultaneous low load ratios, dynamic loss is in a linear correlation with the external magnetic field. At both high magnetic fields and high load ratios, it is in the form of a summation of power functions. This new formulation can be used to characterize the non-linear variation of dynamic loss, which cannot be accurately predicted and explained by the existing analytical methods. The proposed formulas have been well verified by the T-formulation based numerical modeling method and experiments, which can be achieved in MATLAB and dramatically save computational effort compared with numerical modelling methods. The cross-field has also been studied, of which the variation complies well with the J_c(B) dependence. This paper comprehensively demonstrates the variation of dynamic loss at different load ratios and magnetic fields. It can serve as a guide for loss controlling and avoiding quench in the application of HTS CCs, especially with regards to machine windings.

This work was supported by the joint scholarship from the University of Edinburgh and China Scholarship Council (CSC) under Grant [2018] 3101.

According to binomial theorem

$$\begin{eqnarray}&&{\left(a+b\right)}^{n}=\displaystyle \sum _{k=0}^{n}\displaystyle \frac{n!}{k!(n-k)!}\cdot {a}^{n-k}\cdot {b}^{k}\end{eqnarray} \tag{ A.1 }$$

(8) can be expanded into

$$\begin{eqnarray}\begin{array}{lll}f\left(B\right) & = & {\left[1+\displaystyle \frac{\left|{B}_{{\rm{ext}}}\sin \left(2\pi ft\right)\right|}{{B}_{0}}\right]}^{n}=1+n\displaystyle \frac{\left|{B}_{{\rm{ext}}}\,\sin \left(2\pi ft\right)\right|}{{B}_{0}}\\ & & +\,\displaystyle \frac{n\left(n-1\right)}{2!}{\left(\displaystyle \frac{\left|{B}_{{\rm{ext}}}\sin \left(2\pi ft\right)\right|}{{B}_{0}}\right)}^{2}\\ & & +\,\displaystyle \frac{n\left(n-1\right)\left(n-2\right)}{3!}{\left(\displaystyle \frac{\left|{B}_{{\rm{ext}}}\sin \left(2\pi ft\right)\right|}{{B}_{0}}\right)}^{3}+\cdots \end{array}\end{eqnarray} \tag{ A.2 }$$

According to Euler's formula

$$\begin{eqnarray}&&{e}^{ix}=\,\cos \,x+i\,\sin \,x\end{eqnarray} \tag{ A.3 }$$

we have

$$\begin{eqnarray}\begin{array}{lll}{\sin }^{N}x & = & {\left[\displaystyle \frac{1}{2i}\left({e}^{ix}-{e}^{-ix}\right)\right]}^{N}\\ & = & {\left(\displaystyle \frac{1}{2i}\right)}^{N}\cdot \displaystyle \sum _{k=0}^{N}\displaystyle \frac{N!}{k!\left(N-k\right)!}{\left({e}^{ix}\right)}^{N-k}{\left(-{e}^{-ix}\right)}^{k}\\ & = & {\left(\displaystyle \frac{1}{2i}\right)}^{N}\cdot \displaystyle \sum _{k=0}^{N}\displaystyle \frac{N!}{k!\left(N-k\right)!}{\left(-1\right)}^{k}\\ & \times & \left\{\begin{array}{l}\cos \left[\left(N-2k\right)x\right]\\ +i\,\sin \left[\left(N-2k\right)x\right]\end{array}\right\}.\end{array}\end{eqnarray} \tag{ A.4 }$$

When N is odd, we have

$$\begin{eqnarray}&&\left|{\sin }^{N}x\right|={\left(\displaystyle \frac{1}{2}\right)}^{N}\left|\displaystyle \sum _{k=0}^{N}\displaystyle \frac{N!}{k!\left(N-k\right)!}{\left(-1\right)}^{k}\,\sin \left[\left(N-2k\right)x\right]\right|.\end{eqnarray} \tag{ A.5 }$$

On the contrary, when N is even we get

$$\begin{eqnarray}&&\left|{\sin }^{N}x\right|={\left(\displaystyle \frac{1}{2}\right)}^{N}\left|\displaystyle \sum _{k=0}^{N}\displaystyle \frac{N!}{k!\left(N-k\right)!}{\left(-1\right)}^{k}\,\cos \left[\left(N-2k\right)x\right]\right|.\end{eqnarray} \tag{ A.6 }$$

Based on (A.1)–(A.6), (15) can be expanded to

$$\begin{eqnarray}\begin{array}{lll}{f}_{{\rm{avg}}}\left(B\right) & = & \displaystyle \frac{1}{T}\cdot \displaystyle \frac{1}{2iw}\cdot \displaystyle {\int }_{0}^{T}\displaystyle {\int }_{(1-i)w}^{(1+i)w}f\left(B\right)dydt\\ & = & \displaystyle \frac{1}{T}\cdot \displaystyle \frac{1}{2iw}\cdot \displaystyle {\int }_{0}^{T}\displaystyle {\int }_{(1-i)w}^{(1+i)w}{\left[1+\displaystyle \frac{\left|{B}_{{\rm{ext}}}\sin \left(2\pi ft\right)\right|}{{B}_{0}}\right]}^{n}dydt\\ & = & \displaystyle \sum _{k=0}^{n}\displaystyle \frac{n!}{k!\left(n-k\right)!}{\left(\displaystyle \frac{{B}_{{\rm{ext}}}}{{B}_{0}}\right)}^{k}\cdot \displaystyle \frac{1}{T}\displaystyle {\int }_{0}^{T}\left|{\sin }^{k}\left(2\pi ft\right)\right|dt\\ & = & \displaystyle \sum _{k=0}^{n}\displaystyle \frac{n!}{k!\left(n-k\right)!}{\left(\displaystyle \frac{{B}_{{\rm{ext}}}}{{B}_{0}}\right)}^{k}\cdot \displaystyle \frac{1}{T}\displaystyle {\int }_{0}^{T}\left|{\left(\displaystyle \frac{1}{2i}\right)}^{k}\cdot \displaystyle \sum _{m=0}^{k}\displaystyle \frac{k!}{m!\left(k-m\right)!}{\left(-1\right)}^{m}\left\{\begin{array}{l}\cos \left[2\pi \left(k-2m\right)ft\right]\\ +i\,\sin \left[2\pi \left(k-2m\right)ft\right]\end{array}\right\}\right|dt\\ & = & \displaystyle \sum _{p=0}^{n/2-1}\displaystyle \frac{n!}{\left(2p+1\right)!\left[n-\left(2p+1\right)\right]!}{\left(\displaystyle \frac{{B}_{{\rm{ext}}}}{{B}_{0}}\right)}^{2p+1}{\left(\displaystyle \frac{1}{2}\right)}^{2p+1}\cdot \displaystyle \frac{1}{T}\displaystyle {\int }_{0}^{T}\left|\displaystyle \sum _{m=0}^{2p+1}\displaystyle \frac{\left(2p+1\right)!}{m!\left(2p+1-m\right)!}{\left(-1\right)}^{m}\left\{\sin \left[2\pi \left(2p+1-2m\right)ft\right]\right\}\right|dt\\ & & +\,\displaystyle \sum _{p=0}^{n/2-1}\displaystyle \frac{n!}{\left(2p+2\right)!\left[n-\left(2p+2\right)\right]!}{\left(\displaystyle \frac{{B}_{{\rm{ext}}}}{{B}_{0}}\right)}^{2p+2}{\left(\displaystyle \frac{1}{2}\right)}^{2p+2}\cdot \displaystyle \frac{1}{T}\displaystyle {\int }_{0}^{T}\left|\displaystyle \sum _{m=0}^{2p+2}\displaystyle \frac{\left(2p+2\right)!}{m!\left(2p+2-m\right)!}{\left(-1\right)}^{m}\left\{\cos \left[2\pi \left(2p+2-2m\right)ft\right]\right\}\right|dt\\ & = & 1+\displaystyle \sum _{p=0}^{n/2-1}\left\{\displaystyle \frac{n!}{\left(2p+1\right)!\left[n-\left(2p+1\right)\right]!}{\left(\displaystyle \frac{{B}_{{\rm{ext}}}}{{B}_{0}}\right)}^{2p+1}\displaystyle \frac{{2}^{p+1}\cdot p!}{\pi \displaystyle \prod _{q=0}^{2p+1}\left(2q+1\right)}+\displaystyle \frac{n!}{\left(2p+2\right)!\left[n-\left(2p+2\right)\right]!}{\left(\displaystyle \frac{{B}_{{\rm{ext}}}}{{B}_{0}}\right)}^{2p+2}\displaystyle \frac{\left(2p+2\right)!}{{2}^{2p+2}\cdot {\left[\left(p+1\right)!\right]}^{2}}\right\}.\end{array}\end{eqnarray} \tag{ A.7 }$$

Therefore, (14) can be written as

$$\begin{eqnarray}&&\begin{array}{l}{Q}_{\mathrm{dyn},\mathrm{nl}}={E}_{0}L{I}_{{\rm{t}}}{i}^{n+1}\cdot {f}_{{\rm{avg}}}\left(B\right)={E}_{0}L{I}_{{\rm{t}}}{i}^{n+1}\\ \,\,\,\,\,\,\,\,\,\,\cdot \,\left\{\begin{array}{l}1+\displaystyle {\sum }_{p=0}^{n/2-1}\displaystyle \frac{n!}{\left(2p+1\right)!\left[n-\left(2p+1\right)\right]!}{\left(\displaystyle \frac{{B}_{{\rm{ext}}}}{{B}_{0}}\right)}^{2p+1}{\left(\displaystyle \frac{1}{2}\right)}^{2p+1}\cdot \displaystyle \frac{{2}^{3p+2}\cdot p!}{\pi \displaystyle {\prod }_{q=0}^{2p+1}\left(2q+1\right)}\\ +\displaystyle {\sum }_{p=0}^{n/2-1}\displaystyle \frac{n!}{\left(2p+2\right)!\left[n-\left(2p+2\right)\right]!}{\left(\displaystyle \frac{{B}_{{\rm{ext}}}}{{B}_{0}}\right)}^{2p+2}{\left(\displaystyle \frac{1}{2}\right)}^{2p+2}\cdot \displaystyle \frac{\left(2p+2\right)!}{{\left[\left(p+1\right)!\right]}^{2}}\end{array}\right\}.\end{array}\end{eqnarray} \tag{ A.8 }$$