A method for predicting electrical failure of a gilding conductive layer bending

Abstract

The application discloses a gilding conductive layer bending electrical failure prediction method, and belongs to the technical field of flexible electronic reliability simulation and gilding process. The method comprises the following steps: grid modeling, spatial discretization of the gilding conductive layer on the flexible substrate on a D2Q9 discrete velocity lattice; electric conduction lattice Boltzmann solving, iterative solving of the electric potential field distribution through a collision-migration evolution equation; damage evolution, updating of damage variables of each grid point based on the bending mechanics of laminated plates and the cumulative damage criterion; crack dynamic introduction, changing the grid point reaching the critical damage into an insulating boundary and applying a rebound boundary condition; resistance extraction and life prediction, extracting the equivalent resistance between two electrodes, and when the equivalent resistance exceeds the failure threshold, the corresponding bending cycle number is the bending fatigue life. The application can quantitatively evaluate the bending reliability of the gilding conductive layer in the design stage without additional experiments.

Description

Technical Field

[0001] This invention relates to the field of flexible electronics reliability simulation and hot stamping process technology, and in particular to a method for predicting electrical failures caused by bending of the hot stamping conductive layer. Background Technology

[0002] Hot stamping, also known as thermal foil stamping, is a process that transfers a layer of metal foil from a carrier film to the surface of a substrate using heat and pressure. In recent years, this technology has expanded from traditional decorative printing to functional electrical manufacturing. Conductive hot stamping, which uses conductive aluminum foil or metal foil instead of decorative foil for conductive pattern transfer, is widely used in low-cost manufacturing of RFID antennas, flexible circuits, and other fields. Compared to traditional etching and inkjet printing methods, hot stamping offers significant advantages such as a shorter process flow, no wet waste liquid, and continuous roll-to-roll production.

[0003] In the field of flexible electronics manufacturing, transferring a conductive layer of metal foil onto the surface of a flexible substrate using hot stamping to form a bendable conductive circuit pattern has become an important technological approach. This technology boasts a resolution exceeding 300 dpi, eliminates the need for photolithography masks, and offers low manufacturing costs. It has already been used for prototype fabrication in fields such as flexible sensors and wearable electronics. Digital hot stamping technology further eliminates the need for a stamping stencil. Stamping accuracy and positioning are uniformly controlled by a digital motion controller, increasing the utilization rate of electroplated aluminum by 30% to 50% compared to traditional hot stamping. It is suitable for personalized and small-batch precision conductive pattern manufacturing on various electronic-grade substrates, including PVC cards and flexible circuit substrates.

[0004] However, when applying hot stamping conductive layers to flexible electronic applications, there is a technical problem that has not yet been fully solved: the electrical continuity failure of the hot stamping conductive layer on the flexible substrate under repeated bending conditions.

[0005] Specifically, the aluminum foil layer is only about 0.02 micrometers thick, resulting in a significant mismatch in elastic modulus between it and PET or paper-based substrates. Aluminum has a Young's modulus of approximately 70 GPa, while PET's is only about 3 to 5 GPa, a difference of more than ten times. During repeated bending, the strain concentration effect within the aluminum layer causes microcracks to initiate at grain boundaries and propagate along them, eventually leading to the breakage of the conductive path. Because the aluminum layer is extremely thin and covered by a protective layer, the cracks are difficult to detect visually in the early stages, but by then the resistance may have already increased several times to tens of times, severely impacting circuit function.

[0006] Existing bending life testing systems in the field of flexible electronics are mainly established for copper-based conductive inks and liquid metals. Their damage mechanisms, crack morphologies, and failure modes differ fundamentally from those of hot-stamped aluminum layers. The failure of copper-based conductive inks typically stems from the breakage of interparticle connections, while liquid metals maintain conductivity through flow-based self-healing. Hot-stamped aluminum layers, as continuous ultrathin metal films, exhibit a failure mechanism dominated by the initiation, propagation, and penetration of grain boundary microcracks, displaying a completely different damage evolution pattern from the aforementioned two materials. Currently, there are no dedicated failure criteria and life prediction models for hot-stamped aluminum layers, making it impossible to accurately assess their electrical reliability during the engineering design phase. Product design relies heavily on extensive trial-and-error experiments, resulting in long cycles and high costs.

[0007] Therefore, developing a simulation method that can accurately predict the bending fatigue electrical life of the hot stamping conductive layer based on the micro-crack evolution mechanism is of great significance for promoting the engineering application of conductive hot stamping technology in the field of flexible electronics. Summary of the Invention

[0008] The purpose of this invention is to provide a method for predicting electrical failures of hot-stamped conductive layers during bending. This method, by transferring the lattice Boltzmann method to the field of electrical conduction simulation of hot-stamped conductive layers and combining a bending mechanical damage evolution model and dynamic crack boundary treatment technology, can accurately predict the resistance degradation curve and bending fatigue life of the hot-stamped aluminum layer under repeated bending conditions during the design stage. This not only solves the problem of the lack of dedicated failure criteria and life prediction models for hot-stamped aluminum layers in the prior art, but also reveals the nonlinear coupling effect between crack network evolution and current redistribution, providing a quantitative tool for the reliability design of flexible electronic products.

[0009] The above-mentioned technical objective of the present invention is achieved through the following technical solution: a method for predicting electrical failure of hot stamping conductive layer bending, comprising grid-based modeling, simplifying the hot stamping conductive layer on a flexible substrate into a two-dimensional planar conductive film, spatially discretizing the conductive pattern on a D2Q9 discrete velocity grid, assigning a conductor or insulating boundary attribute label to each grid point; Boltzmann solution of the electrical conduction lattice, defining distribution functions in each direction on the lattice and adopting a diffusion-type equilibrium distribution function, iterating to a steady state through a collision-migration evolution equation to solve for the potential field distribution within the conductive layer; Damage evolution is calculated based on the bending mechanics of laminates, which calculates the local strain of each grid point. The damage variables of each grid point are updated cyclically with each bend, combined with the cumulative damage criterion. Crack dynamics are introduced by changing the grid point property to an insulating boundary when the cumulative damage reaches the critical value, and applying a rebound boundary condition to achieve a crack effect with zero current flux. Resistance extraction and lifetime prediction are performed by rerunning the Boltzmann solution of the electrical conduction lattice after each round of crack update to extract the equivalent resistance between the two electrodes. When the equivalent resistance exceeds the preset failure threshold, the corresponding number of bending cycles is the bending fatigue lifetime of the conductive layer.

[0010] The present invention is further configured such that: the collision-migration evolution equation in the electrically conductive lattice Boltzmann solution is Among them, g i (x,t) is the distribution function of grid point x at time t along the i-th discrete velocity direction, c i Let δt be the discrete velocity vector in the i-th direction, and τ be the time step. e The electrical conduction relaxation time is greater than 0.5, g i eq (x,t) is the local equilibrium distribution function; the macroscopic potential is obtained by summing the distribution functions over all directions, i.e. .

[0011] The present invention is further configured such that the equilibrium distribution function adopts a diffusion form. ; where ω i In the D2Q9 model, φ(x,t) represents the grid weight coefficient in the i-th direction, with 4 / 9 for the center direction ω0, 1 / 9 for each of the coordinate axis directions ω1 to ω4, and 1 / 36 for each of the diagonal directions ω5 to ω8. φ(x,t) represents the potential value at grid point x. The electrical conduction relaxation time τ... e With the lattice conductivity σ of the conductive layer LB The relationship between them is One-third of which is the squared value of the lattice sound velocity c in the D2Q9 model. s 2 .

[0012] The present invention is further configured such that the local strain at each grid point in the damage evolution is obtained in the following manner: ; where ε Al (x) represents the local strain of the aluminum layer at grid point x, K sc (x) is the strain concentration factor at grid point x, t PET R is the thickness of the flexible substrate, in meters (m); R is the bending radius, in meters (m); the strain concentration factor K... sc (x) reflects the local strain concentration effect at the grain boundary and is randomly assigned a value at each grid point according to a log-normal distribution.

[0013] The present invention is further configured such that: the cumulative damage criterion is established based on the Coffin-Manson fatigue relationship and Miner's linear cumulative damage theory, and the local bending fatigue life of each grid point is... ; where N f (x) represents the number of bending cycles required to break at grid point x, ε Al For the macroscopic strain of the aluminum layer, ε f c is the fatigue ductility coefficient of the aluminum film. fThe fatigue ductility index is negative; the damage increment at grid point x per bending cycle is Δd(x) = 1 / N. f When the cumulative damage d(x) = ΣΔd(x) reaches 1, the grid point is transformed into a crack grid point.

[0014] The present invention is further configured such that: the damage evolution also includes a crack propagation acceleration mechanism, wherein when a lattice point fractures, the effective strain concentration factor of its adjacent unfractured lattice points is updated according to the following formula. Among them, K sc,eff (x) is the updated effective strain concentration factor, K sc,0 (x) is the initial strain concentration factor at this grid point, α is the stress redistribution amplification factor, and n crack (x) represents the number of broken grid points among the nearest neighbors of this grid point, n. total This represents the total number of nearest neighbors of the given grid point.

[0015] The present invention is further configured such that, in the dynamic introduction of the crack, for partially penetrating cracks that have not yet completely penetrated the thickness direction of the aluminum layer, a partial bounce pattern is adopted. Among them, g i ' is the distribution function in the opposite direction, x w Let t* be the coordinates of the crack grid points, t* be the time after the collision, and β be the crack insulation parameter with a value ranging from 0 to 1; the relationship between β and the crack penetration depth is as follows: , where h crack t represents the crack penetration depth. Al n is the total thickness of the aluminum layer. β This represents the nonlinear index of crack morphology.

[0016] The present invention is further configured such that, in the resistance extraction and lifetime prediction, the equivalent resistance is extracted by applying a fixed potential φ1 to one electrode of the conductive pattern and a fixed potential φ2 to the other end, and after steady state, calculating the total current through any cross-section S and obtaining the equivalent resistance according to the following formula. , Among them, R eq Equivalent resistance, unit is Ω, I total σ is the total current passing through the cross-section, expressed in A. e The surface conductivity of the thin film, δx represents the gradient component of the potential along the direction of current flow, and δx is the grid spacing in meters.

[0017] The present invention is further configured such that: the preset failure threshold is determined according to the application scenario; for RFID antenna application scenarios, the preset failure threshold is 3 to 5 times the initial equivalent resistance; for flexible circuit wire application scenarios, the preset failure threshold is 2 times the initial equivalent resistance; the convergence criterion for iterating to a steady state (i.e., the steady-state convergence criterion) is that the relative L2 norm of the full-field potential change in two adjacent steps is less than the convergence threshold ε. conv , where ε conv The value is 10 -6 Up to 10 -8 .

[0018] The present invention is further configured such that: in the gridded modeling, the grid spacing δx ranges from 50 to 500 nanometers, so that a single grid point approximately corresponds to a grain unit of the hot-stamped aluminum layer; the strain concentration factor K sc The log-normal distribution of (x) is expressed as: ; where μ K is the location parameter of the log-normal distribution, and s is the scale parameter. Both are determined by the grain size distribution and texture characteristics of the hot stamping aluminum layer.

[0019] In summary, the present invention has the following beneficial effects: This invention establishes a complete system for predicting electrical failures due to bending fatigue by transferring the lattice Boltzmann method from the field of fluid mechanics to the field of electrical reliability simulation of hot stamping conductive layers. This fills the gap in the existing technology for the lack of specific failure criteria and life prediction models for hot stamping aluminum layers, enabling designers to quantitatively assess the bending reliability of hot stamping conductive layers during the product development stage.

[0020] The rebound boundary conditions of the lattice Boltzmann method of this invention can handle arbitrarily complex and dynamically evolving crack geometry with extremely low computational cost, without having to regenerate the computational mesh after each crack propagation. In contrast, the traditional finite element method faces huge computational overhead and numerical stability problems caused by frequent remeshing in such dynamic fracture problems.

[0021] Furthermore, this invention utilizes the full-field potential solving capability of the lattice Boltzmann method to reveal the nonlinear resistance growth effect caused by current redistribution in crack networks, particularly the current density concentration phenomenon in the inter-crack bridging region and the resulting damage acceleration positive feedback mechanism. This coupling effect surpasses existing prediction methods based on a simple linear relationship between crack density and resistance, and can more accurately capture the accelerated rise phase in the resistance degradation curve, thus providing more reliable lifetime prediction results. Attached Figure Description

[0022] Figure 1This is a schematic diagram of the overall process of the method in an embodiment of the present invention, illustrating the cyclical coupling relationship of five steps: grid-based modeling, electrical conduction lattice Boltzmann solution, damage evolution, dynamic crack introduction, and resistance extraction and lifetime prediction. Detailed Implementation

[0023] The present invention will be further described in detail below with reference to the accompanying drawings.

[0024] I. Holistic Approach This invention provides a method for predicting electrical failure of a hot-stamped conductive layer during bending. The method comprises five steps: gridded modeling, electrical conduction lattice Boltzmann solution, damage evolution, dynamic crack introduction, and resistance extraction and lifetime prediction. The latter four steps—electrical conduction lattice Boltzmann solution, damage evolution, dynamic crack introduction, and resistance extraction and lifetime prediction—form a closed loop, progressively advancing according to the number of bending cycles until the equivalent resistance exceeds a preset failure threshold.

[0025] Grid-based modeling is the preprocessing step of the entire method. Its task is to transform the actual physical structure of the hot-stamped conductive layer into a discrete computational domain that can be processed by the lattice Boltzmann method. Since the thickness of the hot-stamped aluminum layer is only about 0.02 micrometers, which is negligible compared to its in-plane dimensions, simplifying the three-dimensional conductive layer into a two-dimensional planar conductive film is physically reasonable. When spatially discretizing the conductive pattern on the D2Q9 discrete velocity lattice (a computational model that divides a two-dimensional plane into a grid and stipulates that the information on each grid point can only propagate along nine fixed directions, namely stationary, up, down, left, right and four diagonals; used in this invention to simulate the potential distribution and current path of the conductive layer under the influence of cracks; a standard model in the lattice Boltzmann method), it is necessary to determine the range and boundary of the computational domain according to the actual geometry of the conductive pattern, select an appropriate grid spacing to ensure that the influence of microcracks on the current path can be distinguished, and assign an initial attribute label to each grid point, namely a conductor grid point or an insulating boundary grid point.

[0026] The Boltzmann solution for the electrical conduction lattice is the computational engine of the method. This step utilizes the mathematical isomorphism between the steady-state electrical conduction equation and the steady-state diffusion equation, i.e., both satisfy the Laplace equation. This method directly applies the lattice Boltzmann framework, originally used to solve diffusion problems, to the solution of potential fields. Distribution functions are defined in each direction on the lattice, employing a diffusion-type equilibrium distribution function. The collision-migration evolution equation is iterated repeatedly until the potential field converges to a steady state. The steady-state potential field distribution reflects the potential values ​​at various points within the conductive layer under the current crack distribution, from which the current density distribution and equivalent resistance can be further calculated.

[0027] Damage evolution is responsible for simulating the accumulation of microscopic damage in the aluminum layer during bending. This step first calculates the local strain values ​​at various locations in the aluminum layer at a given bending radius based on the bending mechanics theory of laminates. Then, it determines the damage increment at each grid point for each bending cycle using a cumulative damage criterion, and adds the damage increment to the cumulative damage variable at each grid point. The damage evolution process fully considers the statistical randomness of strain concentration at the grain boundaries of the aluminum layer and the stress redistribution effect during crack propagation.

[0028] The dynamic introduction of cracks serves as a bridge between damage evolution and electrical conduction solutions. When the cumulative damage at a lattice point reaches a critical value of 1, this step changes its property from conductor to insulating boundary and applies a bounce boundary condition. This operation is extremely simple within the lattice Boltzmann framework, requiring only modification of the lattice label and the corresponding distribution function update rules, without any alteration to the computational mesh. For partially penetrating cracks, a partial bounce scheme is used to partially block the current flux.

[0029] Resistance extraction and lifetime prediction are the output components of the method. After each crack update cycle, the electrical conduction lattice Boltzmann solver is rerun to a new steady state. Then, the equivalent resistance between the two electrodes is extracted by applying a fixed potential boundary condition and calculating the steady-state current. The equivalent resistance values ​​corresponding to each bending cycle are recorded to form a degradation curve of resistance as a function of the number of bending cycles. When the equivalent resistance first exceeds the preset failure threshold, the corresponding number of bending cycles is determined as the bending fatigue life of the conductive layer.

[0030] Additionally, regarding the restart strategy for the lattice Boltzmann solver after crack update, this invention adopts a hot-start approach, that is, using the potential field distribution obtained from the previous steady-state solution as the initial condition for the new iteration. Specifically, after the crack lattice point update is completed, the distribution function values ​​on all unaffected lattice points are kept unchanged, and only the distribution functions of the newly generated crack lattice points and their direct neighboring lattice points are locally reset. That is, the distribution functions of all directions of the crack lattice points are set to zero, and the distribution functions of its neighboring lattice points are reset to the equilibrium distribution functions corresponding to the current potential value. This allows the solver to iterate from an initial field close to the true steady-state solution, reducing the number of iterations required for convergence by approximately 2.5 to 5 times compared to a complete re-initialization.

[0031] Regarding the statistical reliability guarantee of the method output results, due to the strain concentration factor K in this method... sc (x) Values ​​are assigned randomly at each grid point; different random seeds will produce different simulation results. To obtain lifetime predictions with engineering reliability, this invention employs a Monte Carlo repeated simulation strategy: running 20 to 50 complete simulations with different random seeds, each simulation generating an independent R-value. eq(N) Degradation curves and corresponding failure lifetime values. The final output is a statistical summary, including the mean, standard deviation, and 95% confidence interval of the failure lifetime. For conservative engineering designs, the mean minus one standard deviation is taken as the design lifetime value.

[0032] II. Collision-Migration Evolution Equation This invention provides specific definitions for the collision-migration evolution equation in the Boltzmann solution of the electrical conduction lattice. This evolution equation is the mathematical model for solving the entire potential field, and its specific expression is as follows:

[0033] ; In computational physics, this equation serves as the fundamental evolution operator of the lattice Boltzmann method. Its function is to drive the distribution function to gradually approach a steady state on a discrete lattice through two sub-processes: collision relaxation and migration propagation. This implicitly solves the macroscopic transport equation at the mesoscopic scale. In this invention, it is specifically used to solve the two-dimensional steady-state potential field with complex crack boundaries.

[0034] This equation describes the temporal evolution of the distribution function on the lattice, and is composed of collision steps and migration steps. The first term g i (x,t) represents the distribution function value at grid point x at the current time along the i-th direction. The second term... This is the BGK collision operator, describing the relaxation process of the distribution function approaching a local equilibrium state under collision. (Left side of the equation) This indicates that after collisions and migrations, the distribution function along c i The direction is shifted to the new value of the next grid point.

[0035] The detailed meanings of each parameter are as follows. g i (x,t) is the distribution function of lattice point x along the i-th discrete velocity direction at time t, and is dimensionless. i Let be the discrete velocity vector in the i-th direction. There are nine directions in the D2Q9 model. δt is the time step of a single iteration, taking the value 1 in a lattice unit. τ e is the electrical conduction relaxation time, a dimensionless parameter whose value must be strictly greater than 0.5 to ensure a positive physical conductivity. Its typical value ranges from 0.7 to 1.2. i eq (x,t) is the local equilibrium distribution function, which is determined by the macroscopic potential value at this specific grid point.

[0036] The recovery of the macroscopic potential is achieved by performing zero-order moment operations on the distribution function in all directions, i.e.: ; This summation relation is used as a standard means of extracting macroscopic quantities in the lattice Boltzmann method. Its function is to aggregate the distribution function information of multiple discrete directions at the mesoscopic level into a single macroscopic scalar field value. In this invention, it is specifically used to extract the potential value at each grid point in real time from the distribution function in 9 directions, providing basic data for subsequent current density calculation and resistance extraction.

[0037] From an algorithm implementation perspective, each iteration time step contains two sub-steps. The first is the collision step, where the distribution function is updated in-situ at each grid point:

[0038] ; The collision step formula serves as a local relaxation operator at the lattice points in the lattice Boltzmann method. Its function is to relax the non-equilibrium distribution function in each direction with a relaxation time τ. e The controlled rate returns to the equilibrium state, which physically corresponds to the process of microscopic particles colliding inside the medium to reach local thermodynamic equilibrium. In this invention, it corresponds to the charge tending to be evenly distributed in a local area due to the effect of conductivity.

[0039] Then comes the migration step, which moves the post-collision distribution functions along their respective directions to adjacent grid points: ; The migration step formula serves as an information propagation operator in the lattice Boltzmann method. Its function is to accurately shift the distribution function after the collision along the discrete velocity direction by a lattice point distance, thereby realizing the information exchange between adjacent lattice points. In this invention, it corresponds to the propagation and diffusion process of potential information in the conductive thin film along various directions.

[0040] This collision-transfer separation algorithm structure is suitable for parallel computing because the collision step is entirely a local operation at the grid point, while the transfer step is a regular data moving operation.

[0041] The specific implementation of the fixed potential boundary condition in the lattice Boltzmann framework adopts the equilibrium state forced assignment method: for an applied fixed potential φ boundary For the electrode grid points, after the collision step is completed and before the migration step is executed at each iteration time step, the distribution function of all directions on the grid point is forcibly reset to: ; This forced assignment formula serves as a standard implementation of the first type of boundary condition (Dirichlet condition) in the lattice Boltzmann method. Its function is to lock the distribution function of the boundary lattice points to the equilibrium value corresponding to the specified macroscopic quantity at each time step, thereby keeping the macroscopic field quantity at the boundary constant during the iteration process. In this invention, it is specifically used to simulate the physical boundary condition of applying a constant voltage at the electrode.

[0042] This method ensures that the potential at the electrode remains at a constant applied value. In subsequent migration steps, the reset distribution function will migrate normally to adjacent grid points, thus establishing the correct potential gradient distribution near the electrode. This method is simple to implement and has good accuracy, making it particularly suitable for computational domains with many dynamically changing internal boundaries.

[0043] III. Relationship between Equilibrium Distribution Function and Conductivity This invention defines the specific form of the equilibrium distribution function and the quantitative correlation between relaxation time and conductivity. The equilibrium distribution function adopts a diffusion-type form:

[0044] ; This equilibrium distribution function serves as a simplified equilibrium form without convection terms in the field of lattice Boltzmann diffusion solutions. Its function is to remove higher-order terms related to macroscopic velocities, retaining only the zeroth-order terms proportional to the scalar field, thus enabling the exact recovery of the Laplace equations from the lattice Boltzmann evolution equations under the Chapman-Enskog expansion. Specifically, in this invention, it is used to ensure that the macroscopic equations solved by the solver are consistent with the steady-state electrical conduction control equations.

[0045] This form is a simplified version of the hydrodynamic equilibrium distribution function in the standard lattice Boltzmann method, removing the nonlinear terms related to macroscopic velocities and retaining only the zeroth-order term proportional to the scalar field value. This simplification is reasonable because the electrical conduction problem is a pure diffusion problem and does not involve convective transport.

[0046] Where ω i In the D2Q9 model, the grid weight coefficients for the i-th direction are a set of predetermined constants. The weight ω0 for the center direction (i=0) is 4 / 9, the weights for the four directions along the coordinate axes (i=1 to 4) are each 1 / 9, and the weights for the four diagonal directions (i=5 to 8) are each 1 / 36. This set of weight coefficients satisfies the normalization condition Σ. i ω i =1, while satisfying the symmetry requirements of each moment, ensuring that the macroscopic diffusion equation can be correctly recovered through the Chapman-Enskog expansion. φ(x,t) is the potential value at the lattice point x.

[0047] Relaxation time τ e With the lattice conductivity σ of the conductive layer LB There is a resolution relationship between them: ; In the lattice Boltzmann method, this correlation formula serves as an analytical mapping between the transport coefficient and the relaxation parameter, and its function is to establish a controllable numerical parameter τ at the lattice scale. eThe quantitative correspondence between the macroscopic physical transport coefficient and the target material transport characteristics can be accurately simulated by selecting an appropriate relaxation time. In this invention, it is specifically used to map the actual conductivity of the hot stamping aluminum layer to the calculation parameters in the lattice Boltzmann solver, ensuring that the simulation results have the correct physical dimensions and values.

[0048] One-third of which is the square of the lattice sound velocity c in the D2Q9 model. s 2 This relationship is a fixed mathematical constant derived through Chapman-Enskog multiscale expansion. When τ e When the value is 1, the lattice conductivity σ LB Equals 1 / 6; when τ e When the value is 0.8, the lattice conductivity σ LB Equals 0.1. Choose τ. e It should be noted that the closer the value is to 0.5, the worse the numerical accuracy; when the value is too large, the iteration convergence speed decreases.

[0049] The calibration relationship between lattice conductivity and physical conductivity is as follows: Given the physical sheet resistance R... sheet The physical surface conductivity is σ e =1 / R sheet Lattice current I LB With physical current I p The conversion relationship between them is I p =I LB ×σ e / σ LB Grid unit resistance R LB With physical resistance R p The transformation relationship is R p =R LB ×σ LB / σ e Verification: Let τ e =1.0, i.e., σ LB =1 / 6, sheet resistance R sheet =2.0Ω / square, the physical resistance R of a rectangular conductor with a grid of 200×100 points. p =R sheet ×L / W=2.0×200 / 100=4.0Ω, the result obtained through lattice unit conversion is consistent, and the calibration relationship is correct.

[0050] IV. Local Strain Calculation This invention defines a specific method for calculating the local strain at each grid point during damage evolution. The formula for calculating local strain is:

[0051] ; In the field of thin film bending mechanics, this formula serves as a simplified application of the neutral plane theory of laminates. Its function is to transform macroscopically measurable bending parameters (substrate thickness and bending radius) into local tensile strain values ​​at various points on the thin film surface. It also introduces strain non-uniformity at microscopic grain boundaries through strain concentration factors. In this invention, it is specifically used to provide local mechanical load inputs that drive fatigue damage accumulation for each grid point, connecting macroscopic bending conditions with microscopic damage evolution.

[0052] The physical basis of this formula is the classical theory of laminate bending mechanics. For a multilayer structure consisting of a hot-stamped aluminum layer, a hot-melt adhesive layer, and a PET substrate, under a given bending radius R, the strain of each layer during pure bending deformation is determined by its distance from the neutral plane. Since the thickness of the hot-stamped aluminum layer is approximately 0.02 micrometers, which is negligible compared to the thickness of the PET substrate, and the neutral plane is approximately located in the middle of the PET substrate, the macroscopic bending strain of the aluminum layer is simplified to t. PET / (2×R).

[0053] The detailed meanings of each parameter are as follows. ε Al (x) represents the local strain of the aluminum layer at grid point x, which is dimensionless and indicates the actual degree of tensile deformation of the aluminum film material at that location. K sc (x) is the strain concentration factor at lattice point x, which is dimensionless and reflects the local strain amplification effect at grain boundaries caused by the difference in crystal orientation between adjacent grains. PET This refers to the thickness of the flexible substrate, in meters (m). A typical value for commonly used PET films is 25 × 10⁻⁶. -6 Up to 125×10 -6 m. R is the bending radius, in meters, defined as the radius of curvature of the neutral surface of the substrate during bending deformation.

[0054] Strain Concentration Factor K sc (x) is one of the physical parameters of this model. In actual hot stamping aluminum layers, the film is composed of a large number of closely packed columnar grains, with grain boundaries between adjacent grains. Due to the different crystal orientations of each grain, strain concentration occurs at the grain boundaries under macroscopic uniform strain, and the strain concentration at the grain boundary junctions is particularly severe. K sc (x) is randomly assigned a value according to a log-normal distribution at each grid point. The assignment process is completed in the gridded modeling. After the assignment, the value remains unchanged as an inherent material property throughout the simulation process, unless affected by the stress redistribution effect in the crack neighborhood.

[0055] V. Cumulative Damage Criterion This invention defines the specific mathematical form of the cumulative damage criterion. This criterion is based on the Coffin-Manson low-cycle fatigue relationship and Miner's linear cumulative damage theory. The formula for the local bending fatigue life at each grid point is:

[0056] ; In the field of metal fatigue mechanics, this formula serves as the inverse function of the Coffin-Manson low-cycle fatigue life relationship. Its function is to predict the number of cycles required for a material to fracture under cyclic loading based on the plastic strain amplitude it experiences. In this invention, it is specifically used to convert the local effective strain level of each grid point into the expected fatigue fracture life of that point, thereby determining the damage accumulation rate and fracture sequence of each grid point in bending cycles.

[0057] This formula is derived from the standard Coffin-Manson relation ε. a =ε f ×(2N f ) cf The algebraic transformation, where ε a N represents the strain amplitude. f This represents the fatigue life cycle count. The meanings of each parameter are as follows: N f (x) represents the number of bending cycles required to break at lattice point x. ε Al For the macroscopic strain of the aluminum layer, by t PET / (2×R) is determined. ε f is the fatigue ductility coefficient of aluminum foil, representing the ultimate ductility of the material under a single reverse loading. For ultra-thin hot stamping aluminum foil, due to the film size effect, dislocation movement is restricted, and its value is significantly lower than that of bulk aluminum material, with typical values ​​between 0.01 and 0.05. f The fatigue ductility index is a negative number, typically between -0.5 and -0.7.

[0058] Because of c f When the effective strain K at the lattice point is negative, it indicates that the strain is negative. sc (x)×ε Al Less than ε f When, the ratio K sc (x)×ε Al / ε f Less than 1, while 1 / c f Since it is a negative number, the calculated result N is... f If (x) is greater than 1, it indicates that multiple cycles are required for fracture. The smaller the effective strain, the better. f The larger (x) is, the longer the lifespan. Based on Miner's linear cumulative damage theory, the damage increment for each bending cycle is:

[0059] Δd(x)=1 / N f (x); This damage increment formula serves as the fundamental expression of Miner's linear cumulative damage theory in the field of structural fatigue reliability. Its function is to quantify the irreversible damage caused to the material by each cycle of load as the reciprocal of the fatigue life at that point, so that the damage contributions under different strain levels can be linearly superimposed. In this invention, it is specifically used to track the damage accumulation process of each grid point in a bending cycle and determine when the critical fracture state is reached.

[0060] The grid point breaks when the cumulative damage d(x) = ΣΔd(x) reaches 1.

[0061] Regarding ε f and c f The experimental method for obtaining fatigue data in this invention employs an indirect film bending fatigue testing method based on substrate support: A series of hot-stamped aluminum layer samples with different PET substrate thicknesses are prepared to correspond to different strain levels. Constant-amplitude bending is performed at a fixed bending radius, and resistance changes are monitored online using a four-probe method. The number of cycles when the resistance increases to twice its initial value is recorded as the fatigue life at that strain level. After obtaining 3 to 5 sets of fatigue life data at different strain amplitudes, the fatigue life is measured using ln(ε... a ) and ln(2×N f Linear fitting in a coordinate system with a slope of c f The intercept corresponds to ln(ε) f This method does not require direct manipulation of the ultrathin aluminum film; the strain level is controlled by changing the substrate thickness, and it can be implemented on a conventional bending tester.

[0062] VI. Crack Propagation Acceleration Mechanism This invention specifies the concrete implementation of the crack propagation acceleration mechanism in damage evolution. When a lattice point fractures, the effective strain concentration factor of its adjacent unfractured lattice points is updated according to the following formula:

[0063] ; In the field of fracture mechanics and damage localization, this formula serves as a discretized approximation of the stress redistribution effect at the crack tip. Its function is to quantify the degree of additional stress concentration caused by an existing crack to its neighboring unbroken material, making crack propagation exhibit spatial correlation and directionality rather than purely random behavior. In this invention, it is specifically used to simulate the experimental observation phenomenon of cracks in hot stamping aluminum layers continuously propagating along the direction perpendicular to bending, and to reproduce the evolution characteristics of cracks from dispersed initiation to banded penetration.

[0064] The meanings of each parameter are as follows. K sc,eff (x) is the updated effective strain concentration factor after considering the influence of neighboring cracks. K sc,0(x) represents the initial strain concentration factor at this grid point, determined by random assignment from a log-normal distribution during gridded modeling. α is the stress redistribution amplification factor, which physically represents the scaling parameter of the strain amplification caused by a single crack neighbor; its typical value ranges from 2 to 4 and can be calibrated using a refined finite element model or experimental data. crack (x) represents the number of broken grid points among the nearest neighbors of grid point x, and its value ranges from 0 to n. total n total This represents the total number of nearest neighbors of grid point x in the D2Q9 model, with a fixed value of 8.

[0065] This acceleration mechanism makes cracks prone to propagate along their initial direction. Once a crack initiates, the effective strain concentration factor of unbroken neighboring lattice points at the crack tip increases significantly, leading to a shortened fatigue life and accelerated damage accumulation at these lattice points, causing the crack to continue propagating along the initial direction. This is consistent with the morphology of bending cracks in aluminum films observed in experiments, where cracks typically appear as elongated, through-cracks perpendicular to the bending direction. In actual calculations, whenever a new crack lattice point is generated, it is necessary to traverse all its nearest neighbor lattice points and update K for neighboring lattice points that are still in a conductive state. sc,eff value.

[0066] VII. Partial Rebound Format This invention defines the treatment method for partially penetrating cracks during dynamic crack introduction. In actual aluminum layer bending damage processes, microcracks may undergo a gradual propagation process from the surface inwards. Before the crack completely penetrates the thickness direction of the aluminum layer, some conductivity is still retained at that point. To describe this partial fracture state, a partial bounce scheme is used:

[0067] ; This formula serves as a standard implementation of the partial bounce boundary condition in the lattice Boltzmann method. Its function is to simulate an internal interface with adjustable permeability by using a linear weighted combination of the bounce component and the transmission component. This allows the flux blocking degree at the boundary to be continuously adjusted between complete blocking and complete transmission. In this invention, it is specifically used to describe the transition state of gradual loss of conductivity during the process of crack initiation to complete penetration of the entire thickness of the aluminum layer, avoiding the distortion caused by the sudden change in resistance due to using only a binary model of complete or unbroken cracks.

[0068] The meanings of each parameter are as follows. g i ' is the distribution function in the opposite direction to the i-th direction. x w t* represents the coordinates of the crack grid points. t* represents the time after the collision. β is the crack insulation parameter, ranging from 0 to 1. When β equals 1, it corresponds to a completely insulating crack; when β equals 0, it corresponds to an intact conductor; when β takes a value between 0 and 1, it represents a partially penetrated state.

[0069] The relationship between β and the physical state of the crack is established through the following: ; In the field of thin film damage mechanics, this formula serves as a phenomenological mapping relationship between crack geometry and macroscopic conductivity. Its function is to transform the physical penetration progress of the crack in the thickness direction into a numerical parameter that can be directly used in the lattice Boltzmann boundary condition. It also reflects the nonlinear influence of crack cross-sectional shape (V-shaped, U-shaped, etc.) on the remaining conductive cross-sectional area through the nonlinear exponent n_β. In this invention, it is specifically used to achieve quantitative coupling between crack penetration depth and local resistance increment.

[0070] Where h crack t represents the depth of the crack penetration along the thickness of the aluminum layer, measured in meters (m). Al This represents the total thickness of the aluminum layer, which is approximately 0.02 × 10⁻⁶ for a typical hot stamping aluminum layer. -6 m. n β It is a nonlinear index of crack morphology, and its value is usually between 1 and 2.

[0071] Crack penetration depth h crack The determination follows these rules: when the cumulative damage d(x) at the lattice points reaches 1, a crack initiates from the outer surface of the aluminum layer. Whether the crack instantaneously penetrates the entire thickness depends on the aluminum layer thickness t. Al The radius r of the plastic zone at the crack tip pz The comparison relationship. The formula for estimating the radius of the plastic zone is:

[0072] ; This formula serves as an approximate expression of plane stress in the field of linear elastic fracture mechanics, representing the modified Irwin plastic zone model. Its function is to estimate the characteristic size of the plastic deformation region at the crack tip caused by stress singularity. In this invention, it is specifically used to determine whether a crack can instantly penetrate the entire thickness of the aluminum layer after initiation. That is, when the thickness of the aluminum layer is less than the radius of the plastic zone, the crack will unstablely propagate and immediately penetrate.

[0073] Where K IC For the fracture toughness of aluminum film, σ y R is the yield strength. For nanocrystalline aluminum thin films, r pz It is approximately 50 to 200 nanometers.

[0074] When t Al Less than or equal to r pz At that time, after crack initiation, the stress concentration at the crack tip covers the entire remaining thickness, and the crack will instantly penetrate the entire thickness, h crack Directly equal to t Al β equals 1. For a hot stamping aluminum layer with a thickness of approximately 20 nanometers in a typical application scenario of this invention, this thickness is much smaller than r. pzTherefore, all crack lattice points are directly subject to the fully rebound boundary condition.

[0075] When t Al Greater than r pz At this point, when the aluminum layer thickness is greater than approximately 0.1 micrometers (a thicker conductive hot stamping layer), additional cycles are required to penetrate the full thickness after crack initiation. At this time, h... crack Initialize to r pz After that, each time a tortuous cycle occurs, h crack The increment is:

[0076] ; This incremental formula draws on the basic idea of ​​strain-driven crack propagation rate in the field of fatigue crack propagation. Its function is to describe the distance that the crack front advances into the material in each cycle under continuous cyclic loading. In this invention, it is specifically used to track the process of cracks gradually penetrating from the surface to the interior in a thicker aluminum layer cycle by cycle. Combined with a partial bounce scheme, it achieves a fine simulation of the gradual increase of resistance with crack depth.

[0077] When h crack Accumulated to t Al When β becomes 1. During this transition period, the lattice Boltzmann solver executes a partial bounce scheme with the updated β value, which can accurately capture the transition behavior of gradually increasing resistance during the crack penetration process.

[0078] VIII. Extraction of Equivalent Resistance This invention specifies the concrete method for extracting equivalent resistance in resistance extraction and lifetime prediction. The extraction of equivalent resistance is based on Ohm's law:

[0079] ; In the field of circuit analysis, this formula is a direct application of Ohm's law. Its function is to calculate the equivalent resistance between two ports using a known applied voltage difference and a measured steady-state current value. In this invention, it is specifically used to convert the full-field potential distribution output by the lattice Boltzmann solver into a single engineering-usable scalar resistance value, which serves as an indicator for determining the electrical failure state of the conductive layer.

[0080] ; In the field of continuous medium electromagnetics, this current summation formula serves as a discrete numerical approximation of the surface integral of current density (i.e., the total current passing through a cross section). Its function is to transform the gridded potential gradient field information into the macroscopic total current value passing through a specified cross section. In this invention, it is specifically used to extract the physical current from the steady-state potential field of the LBM and to complete the quantitative calculation of the equivalent resistance by combining Ohm's law.

[0081] The specific implementation process is as follows: Apply a fixed potential φ1 to all grid points at one electrode of the conductive pattern, and apply a fixed potential φ2 to the other end. Typically, φ1 = 1V and φ2 = 0V. The fixed potential boundary conditions are achieved using the equilibrium-forced assignment method. After running the electrical conduction lattice Boltzmann solver until steady-state convergence, select any complete cross-section S perpendicular to the main current flow direction within the computational domain, and calculate the total current I passing through this cross-section. total For each grid point x on section S, the gradient of the potential in the normal direction is calculated using finite difference. Multiply by the thin film surface conductivity σ e The current contribution represented by a grid point is obtained by summing the contributions of all grid points.

[0082] The meanings of each parameter are as follows. R eq φ1 represents the equivalent resistance between the two electrodes, measured in Ω. φ1 and φ2 are the fixed potential values ​​of the two electrodes, measured in V. total σ represents the total current flowing through cross-section S, expressed in amperes (A). e The sheet conductivity is equal to the sheet resistance R. sheet The reciprocal of the value, expressed in units of S / square. δx represents the gradient component of the potential along the main current flow direction, in V / m. δx is the grid spacing, in meters.

[0083] IX. Failure Threshold and Convergence Criterion This invention defines the method for determining the preset failure threshold and the steady-state convergence criterion. The setting of the failure threshold needs to be determined according to the specific application scenario, because different applications have different tolerances for resistance changes.

[0084] For RFID antenna applications, the effective operation of the antenna depends on a certain quality factor (Q), which is inversely proportional to the antenna resistance. When the antenna resistance increases to a certain level, the Q value decreases to the point where the chip cannot obtain enough energy to be activated, and the tag is considered to have failed. Based on the design margin of a typical RFID system, the failure threshold is set to 3 to 5 times the initial equivalent resistance.

[0085] For flexible circuit conductors, normal operation requires that the voltage drop between nodes does not exceed the system design tolerance. When the conductor resistance increases to twice its initial value, the voltage drop also increases to twice its initial value at the same current, typically exceeding the noise margin of digital circuits or the accuracy requirements of analog circuits. Therefore, the failure threshold for flexible circuit conductors is set to twice the initial equivalent resistance.

[0086] The steady-state convergence criterion uses the relative L2 norm of the total potential change in two adjacent steps: ; This convergence criterion serves as a standard measure of the convergence of steady-state solutions in the field of iterative numerical solutions. Its function is to determine whether the iterative process has sufficiently approached the true steady-state solution by calculating the relative change of the full-field solution between adjacent iteration steps. In this invention, it is specifically used to determine when the lattice Boltzmann solver can stop iterating and output reliable potential field results, thus balancing the relationship between computational accuracy and computation time.

[0087] Where φ n and φ n+1 These are vectors composed of the total field potential values ​​after the nth and (n+1)th iterations, respectively. Represents the L2 norm operation. Convergence threshold ε conv The value is 10 -6 Up to 10 -8 Smaller ε conv The value ensures higher solution accuracy but requires more iterations.

[0088] 10. Grid Spacing and Strain Concentration Factor Distribution This invention limits the selection range of grid spacing and the statistical distribution form of strain concentration factor in gridded modeling. The value of grid spacing δx ranges from 50 to 500 nanometers.

[0089] The selection of grid spacing requires a balance between computational accuracy and computational cost. The basic principle is to ensure that a single grid point approximately corresponds to a grain unit of the hot-stamped aluminum layer. Hot-stamped aluminum layers are typically prepared using a vacuum evaporation process, and their microstructure consists of columnar grains growing perpendicular to the substrate surface. The grain diameter is typically in the range of 50 to 200 nanometers, depending on the evaporation rate and substrate temperature. A lower limit of 50 nanometers corresponds to a fine-grained aluminum film, while an upper limit of 500 nanometers corresponds to a coarser grain structure or engineering estimation scenarios requiring reduced computational load.

[0090] Strain Concentration Factor K sc The statistical distribution of (x) at each grid point follows a log-normal distribution, and its probability density function is: ; This distribution function serves as a standard probabilistic model in the fields of materials micromechanics and reliability engineering, describing the statistical characteristics of grain boundary strain concentration factors in polycrystalline metal thin films. Its function is to provide a basis for random assignment of values ​​for large-scale grid arrays, ensuring that the spatial randomness of crack initiation locations in simulations is consistent with the statistical non-uniformity of the actual aluminum film microstructure. In this invention, it is specifically used to generate K0 on the order of 2 million grid points. sc The initial value field is the foundation upon which the entire method can output statistically significant lifetime prediction results (mean, standard deviation, and confidence interval).

[0091] The physical basis for choosing this distribution is that the strain concentration factor is the product effect of multiple random factors, including the orientation difference angle between adjacent grains, grain boundary type, and the geometric configuration of the grain boundary trihedral point. According to the central limit theorem, the product of multiple independent random variables tends to a normal distribution on a logarithmic scale, that is, the original variables follow a log-normal distribution.

[0092] Where μ K Let be the location parameter of the log-normal distribution, and s be the scale parameter. The mean of the log-normal distribution is exp(μ). K +s 2 / 2). For typical hot stamping aluminum foil, K sc The average value is approximately 1.5 to 3.0.

[0093] μ K There are two methods for obtaining the two parameters, μ and s. The first is the transmission electron microscopy (TEM) cross-sectional observation method: A TEM sample of the aluminum film cross-section is prepared, and the grain orientation distribution is observed through selected area electron diffraction (SID) and dark-field imaging. The orientation difference distribution between adjacent grains is statistically calculated, and then converted into a strain concentration factor distribution using crystal plasticity. The second method is the indirect calibration method: An empirical value for s is selected (s is usually taken as 0.2 to 0.4 for vapor-deposited aluminum films), and then μ is used as the coefficient of variation. K As the only parameter to be calibrated, R is measured in a set of experiments. eq (N) The degradation curve is fitted using least squares, and μ is adjusted. K This method minimizes the deviation between simulation predictions and experimental data. It requires minimal testing equipment, needing only a standard bending tester and a four-probe resistance meter.

[0094] Algorithm conversion process analysis This chapter details the entire process of the cross-domain application of the lattice Boltzmann method from the fields of computational fluid dynamics and statistical physics to the field of hot stamping, including technical modifications and an analysis of the objective technical difficulties.

[0095] Original Algorithm Description: The lattice Boltzmann method, in its original domain—computational fluid dynamics and statistical physics—is a mesoscopic numerical method for solving the Navier-Stokes equations. Its initial inputs are the geometry of the fluid domain, boundary conditions (including inlet velocity distribution or pressure, outlet conditions, wall conditions, etc.), and fluid properties (density and viscosity). Its mathematical model is the lattice Boltzmann equations with the BGK collision operator, which realizes the time evolution of the particle distribution function on a discrete lattice through two steps of collision and migration. Its output is the steady-state or transient flow field distribution, including velocity and pressure fields. Applicable data types include the spatiotemporal distribution of scalar and vector fields on a structured lattice. Constraints include a Mach number much less than 1 to ensure negligible compressibility error, a lattice resolution sufficient to resolve the smallest characteristic scale in the flow, and a relaxation time greater than 0.5 to ensure positive viscosity.

[0096] The specific process for conversion and modification is as follows: Step 1: Replace the complete Navier-Stokes type equilibrium distribution function used in the original algorithm to describe the fluid velocity field and pressure field with a simplified equilibrium distribution function that is only applicable to pure diffusion problems.

[0097] The original form is It contains first- and second-order nonlinear terms related to the macroscopic velocity u. The modified form is g. i eq =ω i ×φ, retaining only the zeroth-order term. This modification allows the evolution equations to revert to the diffusion equations rather than the Navier-Stokes equations at the macroscopic level, corresponding to the Laplace equations for the steady-state electrical conduction problem.

[0098] Step 2: Redefining the recovery methods for macroscopic physical quantities. In the original algorithm, the zeroth moment of the distribution function recovers the fluid density ρ, and the first moment recovers the momentum density ρ×u. After modification, the zeroth moment is redefined to recover the potential value φ, while the non-equilibrium component of the first moment is redefined as information related to the direction of the current density. Specifically, the potential... This is mathematically identical to the original density recovery, but the physical meaning changes from fluid density to electrical potential.

[0099] Step 3: A physical reinterpretation of the relationship between relaxation time and transport coefficient. In the original algorithm, relaxation time τ is determined by... This relationship was previously associated with the kinematic viscosity of the fluid. After modification, this relationship was reinterpreted as a correlation between relaxation time and conductivity. The mathematical form remains the same, but the physical meaning changes from viscous dissipation characteristics to electrical conduction characteristics, and the corresponding calibration experiment changes from rheological measurement to four-probe resistance measurement.

[0100] Step 4: Transformation of the physical meaning of boundary conditions. The impermeable wall bounce boundary condition in the original algorithm corresponds to the physical condition that the fluid velocity is zero on a no-slip wall. After modification, the same mathematical operation is given a new physical meaning: the normal current at the electrical insulation crack is zero. The fixed density (pressure) boundary condition in the original algorithm is converted into a fixed potential electrode boundary condition, implemented using the equilibrium-forced assignment method.

[0101] Step 5: Damage Accumulation-Driven Boundary Dynamic Evolution Mechanism. In traditional fluid dynamics applications, geometric boundaries are typically fixed or move according to a predetermined pattern. However, in the hot stamping application of this invention, crack boundaries are dynamically generated with each bending cycle, and their generation location and timing depend on the distribution of the accumulated damage field. This requires wrapping a damage evolution cycle around the lattice Boltzmann solver and iterating the solver again via a hot-start method after each crack update. This mechanism is absent in the original algorithm and is a novel construct of this invention for the bending failure problem of the hot stamping conductive layer.

[0102] Step 6: Repurposing the Partial Rebound Scheme. In traditional fluid dynamics, the partial rebound scheme is used to simulate semi-permeable walls through which fluid can partially pass in porous media, with the permeability parameter corresponding to porosity. In this invention, this scheme is repurposed to simulate the conductive residual effect of partially penetrating cracks, and the β parameter is redefined from porosity as a function of the crack penetration depth ratio. This repurposing is not merely a simple replacement of the parameter's meaning; it also requires establishing a quantitative correlation model between the β parameter and the crack's physical state, i.e. , and h crack The dynamic rules of cyclical evolution with bending.

[0103] The first problem between the original algorithm and the field of this invention is the difference in physical dimensions and properties. The original algorithm deals with fluid motion problems exhibiting inertia, viscosity, and compressibility, involving the coupled evolution of vector and scalar fields such as density, velocity, and pressure. The electrical conduction problem in the target field is a pure diffusion problem without inertia or convection, involving the steady-state distribution of the potential scalar field. To solve this problem, the equilibrium distribution function must be substantially simplified, all nonlinear terms related to convection and transport must be removed, and the simplified evolution equation must be re-derived to determine which partial differential equation is recovered at the macroscopic level, verifying that it is indeed the Laplace equation required for the target field.

[0104] The second issue is the fundamental difference in time scales. In the original algorithm, there is a direct linear mapping between the iterative time step and physical time; the algorithm's time progression itself simulates the temporal evolution of the physical process. However, in this invention, there are two distinct time scales: one is the extremely short time (on the picosecond scale) for electrical conduction to reach steady state, corresponding to the iterative convergence process within the algorithm; the other is the macroscopic time scale of bending fatigue (with tens to tens of thousands of iterations), corresponding to the external damage-driven loop. To address this issue, a nested structure must be designed between the internal iterative time step and the external bending loop. Within each bending loop, the algorithm runs to steady state, and then damage is advanced and geometry is updated in the outer loop.

[0105] The third problem is the lack of dynamic boundary evolution capability. The original algorithm assumes that the geometric boundary of the computational domain remains fixed or changes according to known rules during the simulation. However, this invention requires the boundary (crack) to be dynamically generated based on the damage information accumulated during the simulation, forming a feedback mechanism where the computational result drives the change of the computational domain. To achieve this function, a completely new damage-crack coupling module must be constructed, and a triggering and execution mechanism for dynamic switching of grid point attributes, as well as a hot-start re-initialization strategy for the distribution function after switching, must be designed.

[0106] The technical verification method for the modified algorithm is mainly carried out in the following ways: First, for a complete conductive layer without cracks, verify whether the equivalent resistance calculated by the lattice Boltzmann solver is consistent with the analytical solution (the resistance of a rectangular conductor is equal to R). sheet ×L / W); secondly, for simple configurations containing cracks of known shape, the lattice Boltzmann results are compared and verified with the finite element method results; finally, the R output of the complete bending fatigue simulation is... eq The (N) curve was compared and calibrated with the experimentally measured bending life data.

[0107] Calculation Example Project Background: An RFID logistics tag uses conductive thermal stamping to form a helical antenna on a PET film. A straight conductive trace with a length L = 200 μm and a width W = 100 μm is selected as the analysis object. Parameters: PET substrate thickness t PET =50μm, hot stamping aluminum layer thickness t Al =20nm, sheet resistance R sheet =2.0Ω / block, bending radius R=10mm, ε f =0.03, c f =-0.6, μ K =0.4, s=0.3, α=3.0, n β =1.5, the failure threshold is defined as 4 times the initial resistance, and the convergence criterion ε conv =10-7 .

[0108] Step 1: Grid-based modeling (one-time preprocessing).

[0109] A grid spacing of δx = 100 nm was selected, with 2000 grid points in the length direction and 1000 grid points in the width direction, for a total of 2 million grid points.

[0110] Boundary conditions are defined as follows: column x=1 represents electrode 1 with φ=1V applied, column x=2000 represents electrode 2 with φ=0V applied, and columns y=0 and y=1000 represent insulating boundaries. A strain concentration factor K is randomly assigned to each conductor grid point according to a log-normal distribution. sc Taking three representative personality points as an example: Grid point A obtains K. sc =2.51, grid point B obtains K sc =1.78, grid point C yields K sc =1.02.

[0111] Calculate the macroscopic strain ε of the aluminum layer Al =t PET / (2R)=2.5×10 -3 Local fatigue life at each grid point N of lattice point A f =6.8 times, N of lattice point B f =12.0 times, N of grid point C f =30.4 times.

[0112] Calculate the radius r of the plastic zone pz =(1 / (2π))×(K IC / σ y ) 2 ≈490nm, due to t Al =20nm is much smaller than r pz After the crack initiates, it instantly penetrates the entire thickness, and β is always taken as 1.

[0113] Step 2: Solve the initial electrical conduction lattice Boltzmann equation. Use the D2Q9 scheme with a relaxation time τ. e =1.0, corresponding to lattice conductivity =(1 / 3)×0.5=1 / 6.

[0114] Initialize the distribution function of each grid point: For any grid point x, first determine the initial potential φ by linear interpolation. init (x) = 1 - (x-1) / 1999, then calculate its equilibrium distribution function g. i (x,0)=ω i ×φ init (x), where ω0=4 / 9, ω1 to ω4=1 / 9, and ω5 to ω8=1 / 36.

[0115] Taking the grid point (1000, 500) as an example, φ init =1-999 / 1999=0.5002, then g0=4 / 9×0.5002=0.2223, g1 to g4=1 / 9×0.5002=0.0556, g5 to g8=1 / 36×0.5002=0.0139.

[0116] Perform a one-step collision-migration operation: collision step g i ×(x,t)=g i (x,t)-(1 / τ e )×[g i (x,t)-ω i [×φ(x,t)], due to the initial state g i It is already in equilibrium, after the collision g i Unchanged; migration step g i (x+c i ×δt,t+δt)=g i (x,t) pushes the distribution functions in each direction along the corresponding grid velocity direction to adjacent grid points. After each collision, g is forcibly reset at the electrode grid points. i =ω i ×φ boundary The insulating boundary grid points employ normal zero-flux bounce. After approximately 5000 steps, the convergence criterion ||φ^(n+1)-φ^(n)|| / ||φ^(n)|| < 10 is met in steady state. -7 After steady state, the potential at each grid point φ(x) = Σg i (x), the potential field is linearly distributed.

[0117] Extracting the initial resistance: Calculate the current density j(x) = Σc at the cross-section x = 1000. i ×g i (x) / (1 / 3), for a uniform linear potential field, the component of j in the x-direction is j. x =σ LB ×(Δφ / δx)=(1 / 6)×(1 / 1999)=8.34×10 -5 (Grid unit).

[0118] Total current Equivalent resistance R eq =ΔV / I total Converted to a physical quantity: R eq =R sheet ×L / W = 2.0 × 2000 / 1000 = 4.0 Ω, determine the failure threshold R. fail =4×4.0=16.0Ω.

[0119] After steps one and two are completed, the remaining four steps (damage evolution, crack introduction, re-solution, and resistance extraction) enter a closed loop. Each bending loop N is executed sequentially: damage evolution → dynamic crack introduction → re-solution of electrical conduction → resistance extraction and failure determination, until R. eq ≥R fail The process will end at that time.

[0120] Cycles 1 through 6 (latency period). Damage evolution is performed in each round, with the cumulative damage at each grid point d(x) += 1 / N. f (x). Perform crack dynamic introduction inspection on all grid points, only a few extreme K points. sc Grid point (K) sc >2.87) Fracture was first triggered when d≥1 at N=5, accounting for approximately 1% of the total grid points, and was dispersed. Each time a new crack occurred, a hot-start electrical conduction calculation was performed to resolve and extract the resistance. Due to the extremely low crack density, the current path was almost unaffected, R eq It slowly increases to about 4.2Ω, which is far below the failure threshold, and the cycle continues.

[0121] The 7th cycle (closed loop and acceleration mechanism).

[0122] Damage evolution: The cumulative damage of grid point A is d = 7 × (1 / 6.8) = 1.03 > 1.

[0123] Crack dynamic introduction: Mark lattice point A as the crack lattice point, and apply a fully rebound boundary condition, that is, all incident distribution functions on this lattice point bounce back in the opposite direction g. i '(x A )=g ī (x A This makes the net current flux zero. The effective strain concentration factor is updated for the first-layer nearest neighbors of grid point A that are not yet fractured (such as grid point B'): if one of B''s eight nearest neighbors is fractured (i.e., grid point A), then K... sc_eff (B')=K sc_0 (B')×[1+3.0×1 / 8]=K sc_0 With a crack density of (B')×1.375, the fatigue life is correspondingly shortened, and subsequent damage accumulation is accelerated. Approximately 40,000 new crack points were added across the entire domain in this round, resulting in a cumulative crack density of approximately 5%.

[0124] The Boltzmann re-solution steps for the electrical conduction lattice are as follows: The distribution functions of all unaffected lattice points from the previous steady-state solution are retained unchanged. Only for the new crack lattice point, the distribution functions in all nine directions are set to zero. For the first-layer nearest neighbor lattice points of the new crack lattice point (i.e., the eight lattice points directly adjacent to the crack lattice point), the distribution functions are reset to the equilibrium value g corresponding to the current potential. i =ω i ×φ(x), while the remaining grid points remain unchanged.

[0125] Starting from these initial conditions, the collision-migration evolution equations are run, and the convergence criterion is satisfied in about 1500 steps, since only the local field is disturbed.

[0126] Resistance Extraction and Lifetime Prediction: Under the new steady state, the equivalent resistance is extracted, and the total current is calculated by summing the values ​​across all 1000 grid points in the width direction at the x=1000 section. R eq =ΔV / I total Because the cracks are dispersed and have a density of only 5%, the current re-converges after bypassing the fracture lattice points, R eq (N=7)≈4.5Ω<16.0Ω, the cycle continues.

[0127] Cycles 8 through 14 (acceleration phase). Due to the positive feedback effect of the crack acceleration mechanism, the number of new cracks increases with each cycle, and crack bands begin to form along the width direction. Each cycle fully executes the four steps of a closed loop. The resistance increases rapidly with each cycle: R eq (N=10)≈5.6Ω, R eq (N=12)≈7.2Ω, R eq (N=14)≈9.5Ω. None of them reached the failure threshold.

[0128] Cycles 15 through 19 (rapid failure period). Cracks extend across the width of the conductive trace, forcing current to bypass the residual bridging region, causing a sharp increase in resistance: R eq (N=15)≈10.8Ω, R eq (N=17)≈13.8Ω, R eq (N=18)≈15.2Ω. R after the 19th cycle eq ≈18.5Ω>R fail =16.0Ω, exceeding the failure threshold for the first time. The failure lifetime N is calculated by linear interpolation between N=18 and N=19. fail Approximately 18.2 times, then the cycle terminates.

[0129] Monte Carlo statistics were performed. The entire procedure was repeated with 30 different random seeds, yielding a mean failure life of 18.6 cycles and a standard deviation of 2.4 cycles. The mean minus one standard deviation was taken as the conservative design value for 16 cycles. Bending tests on samples of the same specification were conducted to verify the method. The mean failure life of 5 samples was 17.8 cycles, and the relative deviation of the simulation prediction was 4.5%, thus verifying the effectiveness of the method.

[0130] Additional scenario: Partial rebound pattern. When the aluminum layer thickness is large (e.g., t...). Al =5μm) when t Al >r pz The crack is approximately 102 nm thick and will not penetrate instantly. At a certain lattice point, when d reaches 1, the initial crack depth h... crack =rpz =102nm, β=(102 / 5000) 1.5 =0.003, which hardly blocks the current. A partial bounce scheme g is applied to this grid point. i '(x w )=β×g ī (x w )+(1-β)×g i eq (x w The current flows through almost unimpeded. The crack depth increases by approximately 2020 nm per cycle, and after 3 cycles, h... crack More than t Al When β reaches 1, it transitions to a complete bounce. During this process, a warm start is performed after each β update to resolve the problem, resulting in a smooth increase in resistance rather than a step jump.

[0131] This specific embodiment is merely an explanation of the present invention and is not intended to limit the invention. After reading this specification, those skilled in the art can make modifications to this embodiment without contributing any inventive step, but such modifications are protected by patent law as long as they are within the scope of the claims of the present invention.

Claims

1. A method for predicting electrical failure of a hot-stamped conductive layer due to bending, characterized in that, include: Grid-based modeling simplifies the hot stamping conductive layer on the flexible substrate into a two-dimensional planar conductive film. The conductive pattern is spatially discretized on the D2Q9 discrete velocity grid, and each grid point is assigned a property label of conductor or insulating boundary. The Boltzmann solution for the electrical conduction lattice defines distribution functions in each direction on the lattice and adopts a diffusion-type equilibrium distribution function. The potential field distribution within the conductive layer is solved by iterating to a steady state through the collision-migration evolution equation. Damage evolution is calculated based on the bending mechanics of laminated plates to determine the local strain of each grid point, and the damage variables of each grid point are updated cyclically by bending, combined with the cumulative damage criterion. The crack dynamics are introduced, and the grid properties that have accumulated damage to the critical value are changed to insulating boundaries. Rebound boundary conditions are applied to them to achieve the crack effect with zero current flux. Resistance extraction and lifetime prediction: After each crack update, the electrical conduction lattice Boltzmann solution is rerun to extract the equivalent resistance between the two electrodes. When the equivalent resistance exceeds the preset failure threshold, the corresponding number of bending cycles is the bending fatigue lifetime of the conductive layer.

2. The method for predicting electrical failure of hot stamping conductive layer by bending according to claim 1, characterized in that: The collision-migration evolution equation in the electrically conductive lattice Boltzmann solution is: ； Among them, g i (x,t) is the distribution function of grid point x at time t along the i-th discrete velocity direction, c i Let δt be the discrete velocity vector in the i-th direction, and τ be the time step. e The electrical conduction relaxation time is greater than 0.5, g i eq (x,t) is the local equilibrium distribution function; the macroscopic potential is obtained by summing the distribution functions over all directions, i.e. .

3. The method for predicting electrical failure of hot stamping conductive layer by bending according to claim 2, characterized in that: The equilibrium distribution function adopts a diffusion form: ； Where, ω i In the D2Q9 model, φ(x,t) represents the grid weight coefficient in the i-th direction, with 4 / 9 for the center direction ω0, 1 / 9 for each of the coordinate axis directions ω1 to ω4, and 1 / 36 for each of the diagonal directions ω5 to ω8. φ(x,t) represents the potential value at grid point x. The electrical conduction relaxation time τ... e With the lattice conductivity σ of the conductive layer LB The relationship between them is as follows: ； Where 1 / 3 is the squared value of the lattice sound velocity c in the D2Q9 model. s 2 .

4. The method for predicting electrical failure of hot stamping conductive layer by bending according to claim 1, characterized in that: The local strain at each grid point during the damage evolution is obtained in the following way: ； Where, ε Al (x) represents the local strain of the aluminum layer at grid point x, K sc (x) is the strain concentration factor at grid point x, t PET R is the thickness of the flexible substrate, in meters (m); R is the bending radius, in meters (m); the strain concentration factor K... sc (x) reflects the local strain concentration effect at the grain boundary and is randomly assigned a value at each grid point according to a log-normal distribution.

5. The method for predicting electrical failure of hot stamping conductive layer by bending according to claim 1, characterized in that: The cumulative damage criterion is established based on the Coffin-Manson fatigue relationship and Miner's linear cumulative damage theory. The local bending fatigue life at each grid point is: ； Where, N f (x) represents the number of bending cycles required to break at grid point x, ε Al For the macroscopic strain of the aluminum layer, ε f c is the fatigue ductility coefficient of the aluminum film. f The fatigue ductility index is negative; the damage increment at grid point x per bending cycle is Δd(x) = 1 / N. f When the cumulative damage d(x) = ΣΔd(x) reaches 1, the grid point is transformed into a crack grid point.

6. The method for predicting electrical failure of hot stamping conductive layer by bending according to claim 1, characterized in that: The damage evolution also includes a crack propagation acceleration mechanism. When a lattice point fractures, the effective strain concentration factor of its adjacent unfractured lattice points is updated according to the following formula: ； Among them, K sc,eff (x) is the updated effective strain concentration factor, K sc,0 (x) is the initial strain concentration factor at this grid point, α is the stress redistribution amplification factor, and n crack (x) represents the number of broken grid points among the nearest neighbors of this grid point, n. total This represents the total number of nearest neighbors of the given grid point.

7. The method for predicting electrical failure of hot stamping conductive layer by bending according to claim 1, characterized in that: In the dynamic introduction of cracks, for partially penetrating cracks that have not yet completely penetrated the thickness direction of the aluminum layer, a partial bounce pattern is adopted: ； Among them, g i ' is the distribution function in the opposite direction, x w Let t* be the coordinates of the crack grid points, t* be the time after the collision, and β be the crack insulation parameter with a value ranging from 0 to 1; the relationship between β and the crack penetration depth is as follows: , where h crack t represents the crack penetration depth. Al n is the total thickness of the aluminum layer. β This represents the nonlinear index of crack morphology.

8. The method for predicting electrical failure of hot stamping conductive layer by bending according to claim 1, characterized in that: In the resistance extraction and lifetime prediction, the equivalent resistance is extracted by applying a fixed potential φ1 to one end of the conductive pattern electrode and a fixed potential φ2 to the other end. After steady state, the total current passing through any cross-section S is calculated, and the equivalent resistance is obtained according to the following formula: ； Among them, R eq Equivalent resistance, unit is Ω, I total σ is the total current passing through the cross-section, expressed in A. e The surface conductivity of the thin film, δx represents the gradient component of the potential along the direction of current flow, and δx is the grid spacing in meters.

9. The method for predicting electrical failure of hot stamping conductive layer by bending according to claim 1, characterized in that: The preset failure threshold is determined according to the application scenario. For RFID antenna application scenarios, the preset failure threshold is 3 to 5 times the initial equivalent resistance. For flexible circuit wire applications, the preset failure threshold is twice the initial equivalent resistance; the convergence criterion for iterating to a steady state is that the relative L2 norm of the full-field potential change in two adjacent steps is less than the convergence threshold ε. conv , where ε conv The value is 10 -6 Up to 10 -8 .

10. The method for predicting electrical failure of hot stamping conductive layer by bending according to claim 4, characterized in that: In the gridded modeling, the grid spacing δx ranges from 50 to 500 nanometers, so that a single grid point approximately corresponds to a grain unit of the hot-stamped aluminum layer; the strain concentration factor K sc The log-normal distribution of (x) is expressed as: ； Where, μ K is the location parameter of the log-normal distribution, and s is the scale parameter. Both are determined by the grain size distribution and texture characteristics of the hot stamping aluminum layer.

Images