Measuring nonuniform web tension for roll-to-roll manufacturing of flexible and printed electronics

High-quality roll-to-roll (R2R) manufacturing for flexible and printed electronics often requires uniform web tension. Nonuniformity in web tension can lead to nonuniform performance of printed electronic devices across the width of the web, and excessive nonuniformity in web tension can lead to web wrinkling. Here we develop and test a noncontact resonance (NCR) method and a gentle contact stiffness mapping (GCSM) method for measuring the average web tension and its linear variation across the width of the web. The NCR method uses the lowest symmetric and anti-symmetric frequencies of a web with a closed-form expression to obtain its linearly varied tension. The closed-form expression includes the significant effects of air loading on web vibrations through accurate hydrodynamic functions. While the GCSM method is based on nonlinear regression of the contact stiffness on multiple locations of the web. Both methods are accurate, reliable, and inexpensive, and are compatible for a wide range of web properties, web path, web tension, measurement configurations, and environmental conditions. We cross-validate the two methods on a stationary test stand and in-line test the NCR method in two spans of a moving commercial R2R system. We measure up to 35.58% cross-span tension variation in that system, and both average tension and its linear variation can vary in different spans of the same R2R system. We expect the results presented in this article can improve quality control of R2R processes for flexible and printed electronics and maximize device yields.

There are different basic types of web tension measuring methods in literature. Instrumented rollers are the most common tension measurement tools in R2R industry. Koç et al and Jeong et al used tension sensors (load cells and dancer, respectively) with speed sensors to measure and control web tension with velocity in multiple spans under a uniform tension assumption [34,35]. Schultheis invented a method with a pressure sensor wound on a roller to measure continuous web tension [36]. However, instrumented rollers methods need investment in customized rollers and need re-calibration when the web path changes. Moreover, they assume that the tension in the web wrapped around the roller is unchanged by the friction. In addition, these sensors tend to drift and are sensitive to the environmental temperature and vibrations [37]. Several authors have used fundamental vibration frequency measurements to infer web tension; however, they did not consider nonuniform web tension and the effect of air loading [38][39][40][41][42]. Our previous works on vibrations of aircoupled web systems showed air effect significantly changes the frequencies and mode shapes of webs with typical tensions used by commercial R2R systems [43][44][45]. Scientists in VTT Information Technology developed a system to measure the nonuniform web tension using air film pressure [37]. This system only works for high web speed in air R2R processes and needs re-calibration when the web path changes. In addition, they assumed the cross-width web contact stiffness is uniform which, as we will demonstrate in this article, cannot occur in finite width systems even the tension is uniform. Jin et al experimentally fitted web tension and the contact force under a fixed web deflection on the web [46]; however, this approach is specific to the exact web properties, web geometry, and roller configuration of their experiment, and needs re-calibration for each specific web measured. Overall, comprehensive, accurate, reliable, and inexpensive methods for measuring nonuniform web tension in R2R processes are an unsolved problem.
In this work, we develop and test a noncontact resonance (NCR) method and a gentle contact stiffness mapping (GCSM) method for measuring average web tension and its linear variation. These methods can supplement the in-line metrology process of existing in air or in vacuo R2R systems, without the need for expensive instrumented rollers. Both methods can accurately measure web tension for a wide range of web properties, web path, web tension, measurement configurations, and environmental conditions. Both methods are based on first principles mechanics models of a tensioned plate. The NCR method includes the plate's interactions with a surrounding fluid. The methods are cross validated on a stationary test stand and the NCR method is used to study the web tension distribution within a commercial R2R system.

Modeling
We choose a 2D isotropic, linearly elastic, uniaxially tensioned, rectangular Kirchhoff plate model to simulate a single span in R2R systems as shown in figure 1(a). To identify average web tension and the linear variation of tension, we present an NCR method which can distinguish the lowest symmetric and antisymmetric resonance frequencies and a GCSM method with accurate web contact stiffness profile. There is no tension applied in the cross-span direction, thus we use a tensioned plate model instead of a linear membrane model to define the spatial dependence of eigenmodes in that direction, which includes a very small but finite web bending stiffness. In addition, the effect of the web line speed to frequencies is very small in the operation range of R2R systems in flexible electronics manufacturing [13,40,47,48], we choose the partial differential equation of motion for out-of-plane vibrations of a stationary tensioned plate [49]: As shown in figure 1(b), x 1 , x 2 , and x 3 are the coordinates along longitudinal (in the plane of the web and along the direction of web tension), transverse (normal to the web surface), and lateral (in the plane of the web but normal to the direction of web tension) directions, respectively; L is the in-span length of the web in the longitudinal direction; b is the cross-span width of the web in the lateral direction; τ is the time;w(x 1 , x 3 , τ ) is the web deflection in the transverse direction; ρ web is the web areal mass density; is the web bending stiffness, with E, h, and υ denoting Young's modulus, thickness, and Poisson's ratio of the web, respectively; ∇ 4 is the biharmonic operator; N 11 (x 3 ) is the web's uniaxial tension per unit width; P(x 1 , x 3 , τ ) is the web surface pressure, such as air pressure in air-coupled web vibrations and pressure by contact force. The nonuniform tension can be described as a polynomial function of x 3 . In this manuscript we present methods to determine from experiments the 1st two terms of the polynomial as: where N ave 11 is the average web tension, σ is a dimensionless ratio for the linear variation of tension describing the discrepancy between maximum or minimum tension to N ave 11 . Especially, N 11 (±b/2) = 2N ave 11 or 0 when |σ| = 1, which indicates there is no tension applied on one of the free edge. So that, |σ| = 1 is a critical tension variation for local web wrinkling close to one of the free edges. We substitute equation (2) into equation (1): Turnbull et al [50] has proved that simply supported boundary conditions can accurately predict linear vibrations of pre-tensioned webs across finite radius rollers. So that we use the following boundary conditions for both NCR and GCSM analyses: (a) the web is simply supported on the upstream and downstream rollers: (b) there is no shear force or bending moment for the free edges:

NCR method
We develop an NCR method to measure N ave 11 and σ using the lowest symmetric transverse resonance frequency f 11 and the lowest antisymmetric transverse resonance frequency f 12 of a web. The eigenmodes of f 11 and f 12 are coupled by the nonuniform tension.
To describe the method in a consistent way for both in air and in vacuo R2R processes, we derive the key results in the presence of air effect. Our previous work predicted f 11 and f 12 for an air-coupled web system semi-analytically [45] but used a uniform tension assumption. In this paper, we improve on the previous work and assume that the web tension varies linearly. We derive the results and perform the convergence study as shown in appendix A and obtain the lowest two web frequencies as: where the value ofM 11 andM 12 are shown in equations (A32) and (A33) and are described as functions of the hydrodynamic functions in (A25), the web dimensions L and b, the web areal mass density ρ web , and the air density ρ air . We choose ρ air = 0 for the in vacuo web systems, so thatM 11 =M 12 = ρ web and: The linear variation of tension σ changes significantly the frequencies and their corresponding mode shapes compared to the uniform tension case. Figure 2 shows the effect of σ to f 11  To measure N ave 11 and σ by f 11 and f 12 , we rewrite equations (6) and (7) inversely and obtain a closedform expression as: Note that the NCR method can solve the linear variation of tension but not its direction, since both positive and negative value of σ gives the same frequencies as shown in equations (6)- (9). There is no air loading in vacuo R2R systems, and we simplify equations (10) and (11) as: Equations (10)-(13) are applicable for a wide range of web properties, web thickness, web aspect ratio, web path, web tension, measurement configurations, and environmental conditions without the need for recalibration.

GCSM method
The 2nd method we developed is a GCSM method.
In this method, we gently apply multiple contact forces in at least two locations along the web width and measure the deflections. The forces are applied gently to make sure the web deforms far below its plastic deformation region. The web tension distribution (N ave 11 and σ) is obtained by nonlinear regression of local contact stiffness in different locations. In each contact location, multiple local contact forces are fitted with their deflections cubically to extract the local contact stiffness for a linear elastic model. In addition, plastic deformation locally damages the web and reduces the accuracy of measurement. To avoid the local plastic deformation, we apply the contact force gently and with a large contact sphere as a Ping-Pong ball to make sure that all the work done by  [51]. c this property is provided by [52]. d this property is solved by [53].
contact force is absorbed by the web strain energy in the elastic region.
To understand the relation between web tension distribution, web deflections, and web stiffness, we assume the web deflection determined by equation (1) in the quasi static loading case (i.e.ẅ = 0) is a linear combination of admissible basis functions as: where M and N are the number of functions along x 1 and x 3 directions, respectively. A mn and B mn are the amplitude of antisymmetric and symmetric components, respectively. In appendix B, we show the details of derivations from equations (3)-(5) and (14) and obtain equation (15) to describe the relation between A mn , B mn , N ave 11 , σ, and the local contact force F as: where the descriptions of matrices C 1 to C 4 and vectors A, B, C 5 , and C 6 are shown in equations (B8)-(B15). X 1 is the coordinate of the contact location in x 1 direction. Once the deflection of a web under the contact force is obtained, we can solve the local contact stiffness as: where X 3 is the coordinate of the contact location in From table B1 we can see that applying the same contact force on different locations gives different local deflections. To understand the deflection shape of the entire web with respect to the contact location, we plot 3D deflection profiles . of a PET web with a same contact force F = 0.1 N on four different locations with (X 1 , figure 3. The web properties are the same as used in figure 2 but with σ = 0. We choose the number of basis functions as M = 51 and N = 18. The results show under a same contact force, along x 1 direction the deflection is larger on the center than close to the rollers, which agrees with the restriction of simply supported boundary conditions in equation (4); while along x 3 direction, the web deflects more with a contact force on the free edge than on the center, which means there is an edge effect of web stiffness.
To understand the spatial variations of the web's contact stiffness, we plot the stiffness profiles of a web under a same N ave 11 but different σ ′ s in figure 4. The contact stiffness is solved by equations (15) and (16). We choose all other web properties the same as for figures 2 and 3. (a) Along the longitudinal direction, the contact stiffness close to the simply supported boundaries is higher than the center region. Eventually, the contact stiffness is infinity on these boundaries. (b) There are edge effects close to the free edges.
Along the lateral direction, the local contact stiffness close to the edges can drop down to half of it in the central region. (c) When the web is under a uniform distributed tension, its contact stiffness profile is symmetric to both x 1 = L/2 and x 3 = 0. (d) Nonuniform tension changes the contact stiffness profile to asymmetric to x 3 = 0.
With the linear elastic Kirchhoff plate model, we need to extract the local contact stiffness in a linear elastic deformation region from contact force and web deflection before solving N ave 11 and σ. The web deflection by a gentle contact force can easily beyond the web thickness, so that the Kirchhoff hypothesis is not suitable [54]. Taking the example in table B1 and figure 3, a 0.1 N contact force on 127 µm thick web can deflect the web in millimeters. Von Kármán theory is used for large deformation of plates [55]. Based on the strain-displacement relations in von Kármán theory and the antisymmetric contact force and web deflection relation on two sides of the web in transverse direction, we write the following equation: where k ′ 1 and k ′ 3 are the linear and cubic coefficient to the deflection to web thickness ratio. So that, we measure the web deflection under multiple local contact forces on the same location of the web to obtain k ′ 1 and k ′ 3 with least squares fitting. The local contact stiffness in linear elastic deformation region yields: After measuring the contact stiffnesses in at least two locations with different X 3 ′ s, we can solve N ave 11 and σ of the web. Solving contact stiffness from known N ave 11 and σ is straight forward by using equations (15) and (16), however a closed-form solution to solve N ave 11 and σ from the contact stiffness is not available. We use nonlinear regression with the trustregion-reflective algorithm to solve N ave 11 and σ. In the optimization procedure, we minimize the root mean square error between the measured stiffnesses and the model as: where r is the total number of locations being measured. The gradients of e to N ave 11 and σ in each iteration step is chosen to be those between ±10% of current estimated N ave 11 and σ, respectively.

Experimental procedure
We perform two experimental tests of the developed methods. First, we sequentially test both NCR and GCSM on a stationary test stand. This enables cross validation of the two methods. Next, we use the NCR method to do in-line monitoring of a commercial R2R system. This demonstrates the performance of the technique in a realistic manufacturing environment.

Cross validation experiment
Schematics and pictures of the stationary test scan used for cross validation is shown in figure 5.  found from the radio of amplitudes and the differences of phases of a measurement on the web and a measurement directly on the speaker. The resonance frequencies are obtained by a single degree-offreedom (SDOF) fitting of the transfer function to its half-power bandwidth as shown in appendix C. Figures 5(d) and (e) are the system schematic and a picture of our implementation of the GCSM method. A force gauge (VTSYIQI HF-5 Digital Push Pull Force Gauge, Vetus Electronic Technology Co., Ltd, Hefei, China) with a Ping-Pong ball contact head is used to both deform the web and measure the resulting force. The same laser sensor and data acquisition system as the NCR method is used to measure and record the web deformation at the loading location. Use of the Ping-Pong ball contact head with its relatively large contact area minimizes local plastic deformation of the web. Forces from 0.06 N to 0.13 N in increments of between 0.004 N and 0.01 N are applied to the web with corresponding web deformations varying between 0.55 mm and 1.10 mm.

In-line monitoring experimental setup
To demonstrate the applicability of our method on a commercial R2R system, we use the NCR method to in-line measure the web tension distribution in two spans of a DICEweb digital inkjet printer R2R system (Prototype and Production Systems, Inc., Minneapolis, MN). A picture of this setup is shown in figure 6. In the DICEweb system, the web tension is applied by the torque from the unwind roller servo motor and the web line speed is adjusted by the torque from the rewind roller servo motor. Both rollers have sensors to measure the radii of rolls to control the applied torques. The DICEweb system has no feedback control for the tension and speed dynamics. For these measurements we use an opaque PET web, whose properties are shown in table 1. The in-span lengths of span 1 and span 2 are 292.10 mm and 107.95 mm, respectively. As shown in figure 6(a), span 1 is upstream, where a module was originally set up for cutting the web to replace a running out unwind roll with a new one. Span 2 is in the functional area of inkjet printers. Figures 6(b) and (c) show the setup of the laser sensor and speaker in the two spans. We place the laser head with its laser sensor orientation orthogonal to the web moving direction.    (19)) e = 2.4594. Since the NCR method could not identify the sign of σ, we choose σ from the NCR method to have the same sign as that from the GCSM method to plot in figure 7(c) and solve the discrepancies of N ave 11 , maximum tension N ave 11 (1 + |σ|), and minimum tension N ave 11 (1 − |σ|) between two methods to be 3.43%, 2.62%, and 4.78%, respectively. The ratio of cross-span tension variation to average tension (2σ) is 54.48%. This amount of variation can significantly affect the performance of printed electronics devices and was present despite our best effects at web alignment. In addition, observations show the edge effect of contact stiffness agrees with theory. Figure 7(d) shows the comparison of contact stiffness and web tension using the NCR and GCSM methods with the 2.85 kg hanging mass. The environmental conditions for this measurement were 101.84 kPa air pressure, 21 • C air temperature, and 76% relative humidity with a calculated air density ρ air = 1.198 kg m −3 , as shown in We observe high contact stiffness at X 3 = ±45.72 mm from the GCSM method on this web. The average stress measured by the NCR method is 139.0 MPa, which is larger than the yield strength of the PET web [51]. The web likely undergoes local plastic elongation along x 1 direction. Local plastic elongation reduces the local web areal mass density ρ web and the local web thickness h. By examining equations (6)-(9), (A32), and (A33) it can be shown that the resonance frequencies slightly increase when the global web areal mass density slightly reduces due to the local plastic elongation. Hence, the NCR method slightly over predicts N ave 11 by equations (10) and (12). The GCSM method uses local contact stiffnesses extracted by von Kármán theory as shown in equations (18) and (19). The reduced local web thickness increases the predicted linear contact stiffness and over predicts the web tension. Since we measured the contact stiffness in limited locations, the local plastic elongation affects the GCSM method more than the NCR method in areas where local plastic elongation has occurred. Without including the two local anomalously high contact stiffness in figure 7(d), we obtain N ave 11 = 175.33 N m −1 and σ = 0.2894 with e = 2.9621. Compared with the results by the NCR method, the discrepancies of N ave 11 ,   maximum tension, and minimum tension between two methods reduce to 0.72%, 5.88%, and 10.27% when excluding the two anomalously high contact stiffness measurements. The two methods NCR and GCSM can also be used for moving webs so long as the transport speeds are much smaller than the critical speed of the web. As speed of the web approaches critical transport speed, the resonance frequencies drop down to zero [48].

Experimental results and discussion
To demonstrate the performance of the NCR method in a realistic manufacturing environment, we perform single spot location measurements at three different line speeds, and over two different spans on a commercial R2R system. Table 2 summarizes these results. The environmental conditions during the measurements are 102.71 kPa air pressure, 21 • C air temperature, 43% humidity, and an air density ρ air = 1.212 kg m −3 , as shown in table 1. The frequency is swept between 40 Hz and 100 Hz to measure resonances in span 1 and between 140 Hz and 200 Hz to measure resonances in span 2. The cross span linear variation in tension is negligible in span 2, but varies by up to 35.58% in span 1, showing that different spans of a R2R system can have different variations of tension. The average tension is between 21.46% and 28.29% smaller in span 2 compared to span 1.
When the linear variation of tension is small, we observe cases where σ 2 takes on a small negative value (i.e. imaginary σ, see appendix D). This is caused by neglecting higher order variation of tension and errors in fitting and measurements; however, this does not significantly affect the use of the NCR method as it is reasonable to approximate cases with small negative σ 2 as having uniform tension.
We performed the in-line measurements over a range of speed typical for flexible printed electronics manufacturing [13,47]. As described in appendix E, we use an axially moving string model to estimate the web's critical transport speed. We define the normalized line speed as the ratio of the line speed to the critical speed. For our measurements, the maximum tested line speed of 2 m min −1 is much less than the critical transport speeds, which are between 53.57 m s −1 and 55.89 m s −1 for span 1 and between 47.32 m s −1 and 47.48 m s −1 for span 2, indicating it is safe to neglect the effect of web speed in our methods. However, the measured tension of the stationary case in span 1 is significantly smaller than in the moving cases. We hypothesize that this is the effect of the motion of the unwind and rewind motors on the web tension. When the web start to move, the resistance by frictions between web and rollers changes and lead to tension variation in some spans. Figures 8(a)-(c) show the transfer function and the SDOF fittings used to extract resonance frequencies at line speeds set to 0, 1 m min −1 , and 2 m min −1 , respectively. The actual line speed fluctuates by up to 0.26 m min −1 when the web moves. The transfer functions of the moving webs are noisier than those of the stationary web. There are several probable causes for this noise. The 1st is the effect of nonuniformity and roughness of the web surface on the output of the laser sensor. As the web moves, imperfections in the web surface could be falsely detected as web deflection. The 2nd possible cause is noise from the motor driving the web motion being transmitted through the web causing it to vibrate. Such spurious vibrations can be detected by the laser sensor and captured in the data analysis. A final possible source of this noise is variation in web tension induced by the unwind and rewind motors as they move. Ultimately, this noise does not stop us from correctly identifying the web's resonance frequencies and therefore is not an impediment to implementing the NCR method.

Conclusions
Nonuniform web tension in R2R processes can lead to nonuniform device performance crossing the width of web. We develop and cross validate an NCR method and a GCSM method to measure the average web tension and its linear variation. These methods are in-expensive and friendly to adapt to different spans with different web path flexibly without the needs of new calibration. In addition, we in-line measure the tension distribution in two spans of a commercial R2R system. The main results are as follows: (a) Increasing linear variation of tension reduces the lowest symmetric frequency and increases the lowest antisymmetric frequency, who are tightly clustered in vacuo with uniform tension. (b) The linear varied tension coupled symmetric and antisymmetric basis functions. (c) We developed an NCR method with a closedform expression to measure average tension and its linear variation by the lowest symmetric and antisymmetric frequencies. (d) The local contact stiffness close to the free edges of a web is smaller than its central region. On the other word, the web deforms more with a contact force on its edge than with a same contact force on its central region. (e) Nonuniform tension changes the contact stiffness profile to asymmetric to its cross-span center. (f) We developed a GCSM method to measure the average web tension and its linear variation from local contact stiffness by nonlinear regression. (g) The NCR method and GCSM method are experimentally cross validated with 3.43% discrepancy in the average web tension and 4.12% discrepancy in the linear variation of tension when there is no local plastically elongation. (h) We use the NCR method to in-line measure the tension in two spans of a commercial R2R system under three different line speeds, respectively. (i) Different spans in a same R2R system can have different average tension and different linear variation of tension. (j) In the commercial R2R system we tested, the cross-span linear variation can be up to 35.58%.
The limits of the proposed methods are as follows: (a) Both methods strictly speaking only work for webs with uniform layers. However, if the printed devices create local nonuniformities in stress or web mass, we expect that the method will continue to apply so long as these local perturbations are small in some sense. Overall, we expect these methods provide accurate, reliable, and inexpensive manners for measuring nonuniform web tension to improve quality control of R2R processes for flexible and printed electronics and maximize device yields.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).

Appendix A. Derivation of NCR method
We assume the surrounding air is initially quiescent (no bulk velocity or transport), and is incompressible, inviscid, and irrotational to choose a 3D potential flow theory for the aerodynamics of it. The surface pressure differential between the bottom surface and upper surface of the vibrating web [43][44][45]: where the factor '2' comes from the anti-symmetric aerodynamic pressure on the bottom surface and upper surface; ρ air is the air density (we choose ρ air = 0 for the in vacuo R2R processes); ϕ(x 1 , 0 + , x 3 , τ ) is the velocity potential of air on the web. Substituting equation (A1) into equation (3), the governing equation for the air-coupled web system is: The aerodynamics of surrounding air follows the continuity equation: where ∇ 2 is the Laplace operator, ϕ(x 1 , x 2 , x 3 , τ ) is the 3D velocity potential of surrounding air. The air velocity on the web surface and normal to the web is identical to that of the vibrating web: We nondimensionalize the governing equations using the definitions: where κ is the in-span length to width ratio (i.e. aspect ratio) of web, ε is the bending stiffness to tension ratio, Λ is a measure of air density to areal mass density of web, the primes denote the nondimensional quantities. Thus, the nondimensional forms of equations (4), (5), (A2)-(A4) are: The boundary conditions for the 3D air are: (a) the air can exchange outside the area of web, i.e. the web is unbaffled: (b) the air is stationary in the far field: Equations (A8)-(A10) are discretized using AMM with the separable forms: is the nondimensional admissible basis functions for web vibration and normalized as´1 is the generalized coordinate, m and n are their notations along longitudinal and lateral directions, respectively. The odd-numbered m and n values describe symmetric basis functions with respect to x ′ 1 = 0.5 and x ′ 3 = 0, respectively; while the even-numbered ones correspond to anti-symmetric basis functions with respect to x ′ 1 = 0.5 and x ′ 3 = 0, respectively. φ ′ mn is the velocity potential with respect to web vibration in the corresponding basis function W ′ mn (x ′ 1 , x ′ 3 ). We do not include initial cross flow along longitudinal or lateral direction in this model [44,45]. Substituting equations (A13) and (A14) into equation (A10) yields: We combine equations (A2), (A13)-(A16) and apply AMM with inner products with W ′ mn : where: where M ′ air is the added air mass matrix [45], which is zero for the in vacuo R2R processes. q ′ is the vector of the generalized coordinates q ′ mn (τ ′ ). K ′ web and K ′ σ are the stiffness matrices corresponding to the average web tension and σ, respectively. Based on the convergence study we presented in our previous works [45], we start from 4 × 4 AMM with polynomial admissible function basis to find the analytical formulas that relate f 11 and f 12 to N ave 11 and σ. We choose an orthogonal set of admissible basis functions: We substitute equation (A21) into equation (A17) and obtain: where F ijmn (κ) ′ s are the hydrodynamic functions fitted following the derivation we presented in our previous work with κ from 0.1 to 10 [45]. Their formulas are: Analytical solutions to the eigenvalue problem with 4 × 4 AMM are not easily available, however we develop a higher order mass approximation method to obtain an analytical solution with the same accuracy but reduced size. In R2R processes for flexible and printed electronics, the web bending stiffness to tension ratio ε is very small, typically from 10 −6 to 10 −4 , so that we assume ε ≈ 0 and start from the 2 × 2 AMM with W ′ 11 and W ′ 12 to solve the eigenvalue problem of equation (A17), and obtain: Solving the eigenvalue problem of equation (A17), we obtain: where: Similarly, an equivalent valueM 12  . (A33) Updating M 11 and M 12 byM 11 andM 12 respectively in equations (A26) and (A27), we obtain the new frequencies as shown in equations (6) and (7). Table A1 shows the convergence study of different AMM orders from equation (A17) and the higher order mass approximation from equations (6) and (7). We choose a PET web with the same system properties as used in figure 2. Recalling our previous works [45] show that 4 × 4 AMM with polynomial admissible function basis can accurately predict the lowest two frequencies, we compare the higher order approximation results with 4 × 4 AMM. In the range of small linear variation of tension with σ from 0 to 0.2, higher order mass approximation estimates f 11 with less than 1.7% error and f 12 with less than 2.6% error in comparison to those by 4 × 4 AMM. When there is a non-zero σ, higher order mass approximation under predicts f 11 and over predicts f 12 .

Appendix B. Derivation of GCSM method
To mathematically find the relation between local contact stiffness and web tension distribution, we choose a point load to describe the contact pressure in the linearly elastic model as: where F is the contact force, δ is the Dirac delta function, X 1 and X 3 are the coordinates of the contact location in x 1 and x 3 , respectively. Substituting equation (B1) into equation (3) and evaluating the time derivative to zero give: Applying inner products of equation (B2) with a trail function φ i and the boundary conditions in equations (4) and (5) to its weak form, we obtain: We assume the web deflection as shown in equation (14) and choose two trail functions as: where p and q are the integers from 1 to M and 1 to N, respectively. Substituting equations (14) and (B4) into equation (B3) gives: Similarly, substituting equations (14) and (B5) into equation (B3) gives: We combine and rewrite equations (B6) and (B7) in a matrix form as equation (15), where: N) . . . . . . . . . . . . . . . C 2 (q, 1) · · · C 2 (q, n) · · · C 2 (q, N) . . . . . . . . . . . . . . .
C 3 (q, n) = C 2 (n, q), C 6 = C 6 (1) · · · C 6 (q) · · · C 6 (N) T , Once the deflection of a web under the contact force is obtained, we can solve the local contact stiffness as in equation (16).
To identify the number of M and N, we compare the convergence of web deflection under a known contact force and with known web tension. A same PET web as used in figure 2 is chosen here with the web tension as N ave 11 = 150.47 N m −1 and σ = 0. We choose a contact force F = 0.1 N applied on the center of the web (X 1 , X 3 ) = L 2 , 0 and the center of a free edge (X 1 , X 3 ) = L 2 , − b 2 , which are the two critical locations on the web where the maximum number of basis functions may be needed to ensure convergence. We solve A mn and B mn from equations (15) and (B8)-(B15). Substituting the values of A mn and B mn into equation (14), we obtain the web deflection. The convergence results with σ ̸ = 0 is similar to those with σ = 0, so we choose to only present those with σ = 0 in table B1. The deflections in table B1 are in the contact locations, which is the maximum local deflections of the web under that contact force. The discrepancies are from the comparison of deflections with the chosen numbers of M and N to those with M = 1E4 and N = 1E4. We choose the results with discrepancies just under 20%, 15%, 10%, 5%, and 1% and with minimum total number of admissible basis functions (M × N) to present. To make sure the entire is under a 1% discrepancy, we choose M = 51 and N = 18 in our investigation. However, for the convenience of a R2R process users M = 3 and N = 5 can also be chosen with a 20% error in their calculations.

Appendix C. SDOF fitting
We obtain the frequencies from frequency response functions by a SDOF fitting. The equation of motion of a damped SDOF system is: where m, c, and k are the mass, viscous damping coefficient, and stiffness, respectively. A is the magnitude of harmonic excitation with frequency ω in rad s −1 . By assuming the particular solution as: where and X is the magnitude of response, θ is the phase between excitation and response, and substituting equation (C2) into equation (C1), we obtain the transfer function G as: where ω n = k m is the nature frequency of the system with a unit rad s −1 , Q = √ km c is the quality factor. So that we can rewrite equation (C3) as: where a 0 = k 2 , a 2 = 1 Q 2 − 2 k 2 ω 2 n , and a 4 = k 2 ω 4 n . We obtain a 0 , a 2 , and a 4 by least squares regression in the half-power bandwidth of a peak magnitude in the frequency response functions. So that we can solve the nature frequency and quality factor as:

Appendix E. Critical transport speed
We consider the web to be stationary in our analysis since the R2R processes for flexible and printed electronics are operated at extremely low speeds [13,47] compared with their critical transport speeds. Since the web bending stiffness to tension ratio ε is very small, typically from 10 −6 to 10 −4 , so that we assume ε ≈ 0 and use the critical transport speed of an axially moving string by Wickert and Mote [48] to estimate the effect of web line speed to web vibrations: The ratio of frequency of moving web f (U11) to frequency of stationary web f (0) is [48]: where U 11 is the web line speed, is the normalized line speed. Thus, when the web line speed increases to its critical transport speed, the frequency of moving web reduces to zero.