A composite material modeling method based on continuous phase evolution and dual-domain independent damage

Abstract

The application discloses a composite material modeling method based on continuous phase evolution and double-domain independent damage, relates to the technical field of composite material modeling, and comprises the following steps: constructing a crystallization kinetics model under the double mechanism of cold crystallization and melt crystallization, and solving a transient crystallinity field controlled by a forming temperature history; introducing an effective phase dimension reduction strategy, fusing the continuously generated new-born crystal phase and the historical crystal phase into a single effective phase elastic strain tensor, combining an amorphous phase viscoelastic-plastic constitutive, and establishing a multiphase parallel mechanical model; establishing a double-domain independent damage evolution model with crystallinity as a weight, and independently evaluating the amorphous phase and the crystal phase; and integrating the model into a finite element framework, realizing continuous simulation prediction of the whole process from hot-press forming to service failure, and solving the technical bottleneck of low multiphase evolution calculation efficiency and lack of physical correlation between crystallinity and damage failure in the prior art, thereby providing an efficient and reliable numerical tool for composite material process-performance integrated design.

Description

Technical Field

[0001] This invention relates to the field of composite material modeling technology, and in particular to a composite material modeling method based on continuous phase evolution and dual-domain independent damage. Background Technology

[0002] Semi-crystalline thermoplastic polymers such as polyetheretherketone (PEEK) have become ideal matrix materials for carbon fiber reinforced composites in the aerospace field due to their excellent specific stiffness, high temperature resistance, and chemical corrosion resistance. Hot pressing is the mainstream manufacturing process for these composite components, in which the material is melted at high temperature and then cooled and solidified under pressure in a mold. Due to the influence of actual process parameter adjustments or equipment cooling environment, the cooling rate experienced by the material during the molding stage varies significantly. This directly determines the crystallization kinetics evolution of the macromolecular chains within the polymer matrix, resulting in different overall crystallinity levels in the molded samples. Differences in crystallinity not only induce varying degrees of volume shrinkage and residual internal stress, but also lead to significant differences in the material's macroscopic mechanical properties such as yield strength, stiffness, and damage tolerance. Therefore, accurately predicting the phase transformation and residual stress evolution during hot pressing at different cooling rates, and assessing the impact of overall crystallinity differences on microcrack initiation and macroscopic tensile failure during service life, is a core challenge in achieving integrated "material-process-performance" design of composite materials.

[0003] However, existing polymer multiphase evolution theories and molding simulation models still face the following technical bottlenecks when dealing with semi-crystalline thermoplastic materials: On the one hand, during the hot pressing process, both melting and cold crystallization occur simultaneously, and existing models are unable to describe the above-mentioned dual phase transition mechanism and its differentiated mapping to mechanical response within a unified framework; On the other hand, accurately tracking the initial stress-free state of the newly formed crystalline phase and its subsequent elastoplastic evolution within each incremental step will cause the number of tensors within the finite element integration point to increase linearly with time step, which will consume a lot of computing power and even cause memory overflow. In addition, conventional multiphase models of composite materials usually adopt the assumption of macroscopic homogenization to directly drive damage evolution, ignoring the essential physical differences between the amorphous phase and the crystalline phase in the polymer matrix in terms of damage mechanism, and cannot accurately reflect the competitive failure relationship between the two phases under different crystallinity. Summary of the Invention

[0004] The purpose of this invention is to provide a composite material modeling method based on continuous phase evolution and dual-domain independent damage, which solves the problems in the above-mentioned prior art of semi-crystalline thermoplastic composite material molding simulation, such as imperfect coupling of dual phase transformation mechanism, computational disaster caused by multiphase constitutive tracking, and lack of cross-scale physical correlation mechanism between crystallinity and damage failure.

[0005] To achieve the above objectives, this invention provides a composite material modeling method based on continuous phase evolution and dual-domain independent damage, comprising the following steps: Step 100: Construct a dual crystallization kinetic model; based on the temperature history of each node in the finite element analysis, adaptively identify the thermal history stage of the current temperature, and activate the corresponding cold crystallization mechanism or melt crystallization mechanism to solve for the total crystallinity; Step 200: Establish a multiphase continuous evolution mechanical constitutive model; use an effective phase dimensionality reduction strategy to merge the newly formed crystalline phase and the historical crystalline phase in each time increment step into a single effective phase elastic strain tensor; establish a hyperelastic-plastic constitutive model for the crystalline phase, and at the same time establish a viscoelastic-plastic constitutive model for the amorphous phase; obtain the total Cauchy stress of the composite material by weighting the amorphous phase and the crystalline phase according to the volume fraction. Step 300: Establish a dual-domain independent damage evolution model; for the amorphous phase, a temperature-sensitive brittle-ductile dual-mechanism damage model is used to calculate the damage variables of the amorphous phase; for the crystalline phase, an independent damage model based on partial equivalent stress is used to calculate... The damage variables of the amorphous phase and the crystalline phase are weighted according to crystallinity to obtain the total macroscopic stress after damage: Step 400: Compile the dual crystallization kinetic model, the multiphase continuous evolution mechanical constitutive model, and the dual-domain independent damage evolution model into a finite element user-defined material subroutine, construct an automated script to perform continuous simulation prediction of the entire process from hot pressing to service failure, and obtain a numerical model of composite materials based on continuous phase evolution and dual-domain independent damage.

[0006] Therefore, the composite material modeling method based on continuous phase evolution and dual-domain independent damage described above has the following beneficial effects: 1. By constructing a continuous temperature-displacement coupled analysis step in a single finite element model, seamless transfer of process state variables such as crystallinity field and effective phase elastic strain tensor is achieved between the molding stage and the service stage. This completely avoids the interpolation error and data loss problems caused by offline mapping of multiple software and multiple models. It can directly and quantitatively evaluate the influence of process parameters such as cooling rate on the service failure behavior of composite materials, and provides an accurate numerical tool for the integrated design of composite materials' "process-performance". 2. An effective phase dimensionality reduction strategy is proposed, which compresses the massive historical tensor set that increases over time within each Gaussian integration point into a single core feature tensor. While fully preserving the memory of the initial stress-free state of the newly formed crystalline phase, the computational overhead and storage requirements for multiphase evolution are reduced by one to two orders of magnitude. This fundamentally solves the computational curse problem in the full life cycle simulation of complex composite material components, making efficient full-size simulation of complex components possible. 3. Breaking the barrier of the traditional macroscopic homogenization assumption, a dual-domain independent damage evolution model was established. This model can directly reflect the control effect of crystallinity determined by the molding process on the final failure path. It clearly reveals the competition mechanism between the failure law dominated by viscoplastic damage under low crystallinity and the failure law dominated by crystal phase yielding and tearing under high crystallinity. It provides a highly reliable core prediction tool for damage tolerance design of composite material structures based on molding process.

[0007] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0008] Figure 1 This is a flowchart of a composite material modeling method based on continuous phase evolution and dual-domain independent damage according to the present invention; Figure 2 Here are the boundary conditions and mesh settings diagrams for Abaqus according to the present invention, wherein (a) is the boundary condition setting diagram and (b) is the mesh setting diagram; Figure 3 The crystallinity field distribution diagrams for two working conditions of the present invention are shown, wherein (a) is the crystallinity field distribution diagram for working condition A and (b) is the crystallinity field distribution diagram for working condition B. Figure 4 This is a comparison diagram of stress-strain curves under two working conditions of the present invention; Figure 5 The diagrams show the Mises stress field distribution under two working conditions of the present invention, where (a) is the Mises stress field distribution under working condition A and (b) is the Mises stress field distribution under working condition B. Figure 6 These are amorphous phase damage distribution diagrams for two working conditions of the present invention, wherein (a) is the amorphous phase damage distribution diagram for working condition A, and (b) is the amorphous phase damage distribution diagram for working condition B. Figure 7 The following are crystal phase damage distribution diagrams for two working conditions of the present invention, wherein (a) is the crystal phase damage distribution diagram for working condition A and (b) is the crystal phase damage distribution diagram for working condition B. Figure 8 The diagrams show the comprehensive damage distribution under two working conditions of the present invention, where (a) is the comprehensive damage distribution under working condition A and (b) is the comprehensive damage distribution under working condition B. Detailed Implementation

[0009] The following detailed description of embodiments of the invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely illustrates selected embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the invention without inventive effort are within the scope of protection of the invention.

[0010] To make the technical solution of this invention clearer, the key terms appearing in the text are defined as follows: Please see Figure 1 A composite material modeling method based on continuous phase evolution and dual-domain independent damage includes the following steps: Step 100: Construct a crystallization kinetic model under a dual mechanism; Based on finite element analysis software, the thermal history stage of the current actual temperature is adaptively identified according to the temperature history of each node, and the corresponding crystallization kinetic mechanism is activated to solve the transient crystallinity field controlled by the molding temperature history. The glass transition temperature of the polymer matrix is ​​set as... Melting point is : Step 101: Modeling the cold crystallization mechanism; When the polymer matrix is ​​undergoing a rapid cooling and heating process, and the current actual temperature of the polymer matrix... satisfy At this time, the model activates the cold crystallization mechanism. Relative volume fraction of cold crystals. The dynamic competition process of quantitatively coupling "nucleation" and "growth" is calculated using the following formula: ; in, This represents the relative volume fraction of cold crystallization. The current time; and All are time-integral variables. The nucleation time of the crystal nucleus. This is an intermediate time variable used for integration. Nucleation density as it evolves with supercooling; The crystal growth rate is jointly controlled by the molecular chain segment diffusion term and the nucleation barrier term.

[0011] Crystal growth rate Calculated using the following formula: ; in, This refers to the crystal growth rate; Basic growth rate parameter; It is the activation energy for molecular chain segment diffusion; It is the ideal gas constant; The limiting temperature at which molecular chain segments cease to move; It is the nucleation kinetic constant; The melting point of the polymer matrix; , Supercooling; , This is the temperature correction factor.

[0012] Nucleation density The formula for calculating the evolution of supercooling is: ; in, Nucleation density as it evolves with supercooling; Based on nucleation density; is the thermal nucleation coefficient.

[0013] Furthermore, under the cold crystallization mechanism, the relative volume fraction of cold crystals calculated at the end of the current time increment step... The absolute crystallinity of cold crystallization equivalent to the current time increment step Therefore, the absolute crystallinity of cold crystallization at the end of the current increment step is updated as follows: = .

[0014] Step 102: Modeling the melting and crystallization mechanism; When the polymer matrix is ​​in a molten state and cooling process, and the temperature satisfy At this point, the model activates the melting and crystallization mechanism. The relative rate of change of crystallinity is characterized by the differential form of the Nakamura equation, and the calculation formula is as follows: ; in, This represents the relative rate of change in crystallinity. Relative crystallinity; The Avrami index characterizes the nucleation and growth dimensional features of crystals. is the crystallization rate constant that varies with temperature.

[0015] The relative increase in crystallinity can reach a maximum crystallization enthalpy. The physical increment of the absolute enthalpy of crystallization is calculated using the following formula: ; in, This is the physical increment of the absolute enthalpy of crystallization; The time increment step; It represents the maximum enthalpy of crystallization that the polymer matrix can achieve under fully crystalline conditions.

[0016] Furthermore, under the melt crystallization mechanism, the relative crystallinity change rate is first analyzed. Integrating over the time increment step Δt yields the relative melt crystallinity increment, calculated using the following formula: = ·Δt; in Indicates the relative increase in melt crystallinity; Subsequently, to map the relative melt crystallinity increment to an absolute crystallinity increment that can directly drive the subsequent multiphase mechanical constitutive evolution, the model introduces the ratio of the maximum achievable crystallization enthalpy at the current temperature to the latent heat of complete crystallization of the material as a conversion coefficient, calculated as follows: ; in, This represents the absolute crystallinity increment within the current time increment step. This represents the relative increase in melt crystallinity. This represents the maximum enthalpy of crystallization at the current temperature. This is the latent heat of crystallization of the polymer matrix under fully crystalline conditions.

[0017] Therefore, the absolute crystallinity of the melt crystallization at the end of the current time increment step is updated to... = + ,in, This represents the absolute crystallinity contributed by melt crystallization at the end of the previous time increment step. This indicates the absolute crystallinity of the melt at the end of the current time increment step.

[0018] Step 103: Calculate the total crystallinity.

[0019] The absolute crystallinity of melt crystallization and the absolute crystallinity of cold crystallization together constitute the total crystallinity at the end of the current increment step. The calculation formula is: = + .

[0020] The crystallinity field obtained in step 100 will directly drive the continuous evolution of the initial tensor of the effective phase in step 200. These physical quantities, which are highly related to the specific forming process, will reside in the solver as internal state variables to drive the multiphase mechanical constitutive evolution in the subsequent service stage.

[0021] Step 200: Multiphase continuous evolution and effective phase mechanics constitutive modeling; The macroscopic mechanical response within the composite matrix is ​​jointly borne by the parallel interaction of the amorphous and crystalline phases. To address the tensor storage bottleneck caused by multiphase continuous evolution, this method introduces an effective phase dimensionality reduction strategy, specifically including the following steps: Step 201: Volume shrinkage correction; The transition of polymers from an amorphous to a crystalline state is accompanied by an increase in density and a macroscopic volume shrinkage. The model incorporates the density of the amorphous phase. With perfect crystal density of crystal phase The ratio is calculated using the following formula: ; in, The ratio of the physical density of the amorphous phase to that of the crystalline phase; The density of the amorphous phase; Let be the perfect crystal density of the crystalline phase. The corrected volume Jacobian component of the crystalline phase is then calculated using the following formula: ; in, The corrected Jacobian component of the crystal phase volume; This is the total deformation gradient tensor at the end of the current increment step; This is a matrix determinant operator. This formula directly couples the real volume shrinkage effect caused by the crystallization phase transition to the framework of large mechanical deformation.

[0022] Step 202: Evolution of continuous crystal phases and simplification of effective phases; Based on the principle of equivalent mapping, the elastic isochoric left Cauchy-Green tensor of the old crystal phase recorded in the previous increment step is extrapolated to the current configuration through the isochoric relative deformation gradient tensor. The calculation formula is as follows: in, To extrapolate the elastic isocompassive left Cauchy-Green tensor of the old crystal phase to the current configuration; The elastic isocompensative left Cauchy-Green tensor of the old crystal phase recorded in the previous increment step; For isochoric relative deformation gradient tensor; for The transpose of .

[0023] Assuming that the newly formed crystalline phase increment is in a purely stress-free initial state at the moment of formation, i.e., its elastic isochoric left Cauchy-Green tensor is a second-order unity tensor. Based on the volume-weighted rule, the old and new crystals are fused into a single macroscopically effective phase elastic isochoric left-Cauchy-Green tensor, and the calculation formula is as follows: in, The effective phase elasticity isocompensatory left Cauchy-Green tensor after fusion; This represents the total crystallinity at the end of the previous increment step; This represents the total crystallinity at the end of the current increment step; , This represents the crystallinity increment of the newly formed crystalline phase within the current time increment step; It is a second-order unit tensor, representing the stress-free configuration at the instant of formation of the new crystalline phase.

[0024] Step 203: Establish the hyperelastic-plastic constitutive model of the crystalline phase; The elastic deformation response of the crystalline effective phase is described using a compressible Neo-Hookean hyperelastic potential function. The Cauchy stress of the crystalline effective phase is calculated using the following formula: ; in, Cauchy stress is the effective phase of the crystal. The shear modulus of the crystalline phase; To account for the phase transformation volume shrinkage correction of the Jacobian component of the crystal phase; For a tensor, the trace operator is the sum of the elements on the main diagonal; This is the bulk modulus of the crystalline phase.

[0025] It should be noted that, considering that the crystalline phases in highly crystalline regions may still undergo crystal plane slip, lamellar rotation, and local shear yielding under large macroscopic deformation, a crystalline phase plastic yield criterion is introduced here. Ignoring this plastic evolution would lead to irreversible deformation being misjudged as a purely elastic response, thereby overestimating the material's peak load-bearing capacity and underestimating the stress softening effect. Therefore, plastic yield criterion and evolution calculations are performed based on the effective stress of the crystalline phase. Equivalent deviatoric stress is used. The yield condition, which drives plastic evolution, is determined by the following formula: ; ; ; in, The yield function of the crystal phase; The equivalent deviatoric stress of the effective phase of the crystalline phase; This represents the yield stress of the crystalline phase. Calculate the deviatoric stress for the crystal phase. If Then the increment is an elastic step, and the Cauchy stress of the effective phase of the crystal is directly taken as the real stress of the crystal phase. And maintain the current effective elastic configuration of the crystal phase without plastic modification; if Then the effective phase crystal enters a plastic flow state, which has an effect on... The corrected phase deviatoric stress is obtained by projecting radial regression mapping onto the current yield surface. Thus, the true stress of the crystal phase is obtained. The calculation formula is as follows: ; ; in, This is the corrected phase deviatoric stress. This represents the true stress of the crystal phase.

[0026] Step 204: Amorphous phase viscoelastic-plastic constitutive modeling; To address the slip phenomenon at molecular chain crosslinking points in amorphous phases under large deformations, this method constructs a constitutive model consisting of a generalized Maxwell viscoelastic branch (characterizing the viscoelastic response) and an elastoplastic slip branch in parallel, as follows: (a) Generalized Maxwell viscoelastic branch; The formula for calculating the instantaneous Neo-Hookean pseudo-stress of amorphous phase viscoelastic deformation using the Neo-Hookean hyperelastic model is as follows: ; in, For the instantaneous pseudo-stress tensor of viscoelastic deformation in amorphous phase; This represents the total initial shear modulus of the amorphous phase. For tensors, trace operators; It is a second-order unit tensor; The bulk modulus of the amorphous phase; The Jacobian component of the total volume deformation is calculated using the following formula: , This represents the total deformation gradient tensor at the end of the current time increment step. The total elastic isotropic left Cauchy-Green tensor is calculated using the following formula: , , The isochoric deformation gradient. Total deformation gradient tensor. The deformation gradient at the integration point is obtained by differentiating the deformation mapping of the current configuration relative to the initial configuration at the integration point based on the nodal displacement field of each element at the end of the current time increment step by the finite element solver, and is passed as the input variable to the user material subroutine UMAT; the total deformation gradient at the end of the previous time increment step is denoted as... .

[0027] The Cauchy stress contributed by the viscoelastic branching of the amorphous phase is solved using the convolution integral, and the calculation formula is as follows: ; in, Cauchy stress contributes to the viscoelastic branching of amorphous phases; The current time; This is the normalized relaxation modulus function; Let be the time integral variable. The normalized relaxation modulus represents the complete relaxation spectrum from rapid relaxation to long-term equilibrium. The calculation formula, obtained through Prony series expansion, is as follows: ; in, To normalize the relaxation modulus; The current time; The equilibrium weighting factor represents the proportion of the long-term equilibrium modulus to the initial modulus. This refers to the Maxwell branch number; For the first The weight coefficients of each Maxwell branch; For the first Relaxation time of each Maxwell branch; This represents the total number of branches in Maxwell.

[0028] To characterize the nonlinear effect of the temperature field on the relaxation rate of macromolecular chain segments, the model uses the WLF (Williams-Landel-Ferry) equation to calculate the temperature shift factor. The calculation formula is as follows: ; in, It is the temperature shift factor; , All are WLF model constants; This is the actual temperature of the polymer matrix; This is a reference temperature.

[0029] The crystal structure resulting from the crystallization phase transition acts as additional physical cross-linking points, hindering the relaxation motion of molecular chain segments. Therefore, a crystallization shift factor is introduced into the model. The calculation formula is as follows: ; in, It is the crystallization shift factor; The crystal barrier parameter characterizes the strength of the restriction on the relaxation motion of molecular chain segments by the crystal structure; The transient total crystallinity; Boltzmann's constant; This represents the actual temperature of the polymer matrix.

[0030] The temperature shift factor and the crystallization shift factor are used together to correct the equivalent time scale of the viscoelastic relaxation process in amorphous phases. Specifically, the model calculates the temperature shift factor based on the current actual temperature. Calculate the crystallization shift factor based on the current total crystallinity. And construct the total shift factor The calculation formula is: Based on this, the real time increment will be... Convert to equivalent time increment The calculation formula is: The equivalent time increment is further used to update the recursive integrals of each branch of the generalized Maxwell model, thereby determining the viscoelastic stress decay, stress increment response, and tangential stiffness contribution of each branch, and ultimately transferring them to the homogenization stress of the amorphous phase, the total stress assembly, and subsequent damage evolution analysis.

[0031] (b) Elastic-plastic slip branch; Construct an elastoplastic constitutive model of slip at physical crosslinking points in the amorphous phase, and solve for the true Cauchy stress contributed by the slip network: in, The actual Cauchy stress tensor contributed to the slip network; The shear modulus of the slip network; The bulk modulus of the slip network; For the elasticity of the sliding network, the left-inclusive Cauchy-Green tensor; The volume Jacobian component of the elastic deformation of the slip network is calculated using the following formula: , This is the elastic deformation gradient tensor of the slip network. The finite element solver calculates the incremental deformation information of the current configuration relative to the configuration of the previous increment step at the integration point based on the displacement field of each element node at the end of the current time increment step. It then inputs the total deformation gradient and its incremental information as input variables to the user material subroutine UMAT. UMAT then updates the elastic deformation of the slip network by combining it with the elastic state of the slip network saved in the previous increment step, thus obtaining the elastic deformation gradient tensor of the slip network corresponding to the current increment step.

[0032] Since amorphous phase physical crosslinking network mainly exhibits reversible segment stretching and recovery at low stress levels, but triggers irreversible slip dissipation at high eccentric stress levels, it is necessary to distinguish between elastic and plastic slip responses based on the slip network yield condition, and update the stress state and accumulated plastic strain accordingly. This is based on equivalent eccentric stress. Constructing yield functions for amorphous phase slip networks to drive plastic evolution: ; ; ; in, The yield function of the amorphous phase slip network; This represents the equivalent deviatoric stress of the slip network; The yield stress of the amorphous phase slip network; Calculate the deviatoric stress for an amorphous phase slip network. If... If the increment is the elastic response, then the actual stress of the slip network is directly taken as the actual stress contributed by the elastoplastic slip branch of the amorphous phase. The calculation formula is: ;like Then, irreversible slippage occurs at the physical cross-linking points, affecting... The corrected deviatoric stress is obtained by projecting it onto the current yield surface using radial regression mapping. This allows us to obtain the true stress contributed by the elastoplastic slip branch of the amorphous phase. The calculation formula is as follows: ; ; in, For the corrected amorphous phase deviatoric stress, The actual stress contributed to the elastoplastic slip branch of the amorphous phase.

[0033] Step 205: Total stress assembly.

[0034] Total Cauchy stress in composite materials The amorphous phase and the crystalline phase share the burden in parallel according to their volume fractions: ; in, The total Cauchy stress of the composite material; This represents the total crystallinity at the end of the current increment step; Cauchy stress contributes to the viscoelastic branching of amorphous phases; The actual stress contributed to the elastoplastic slip branch of the amorphous phase; This represents the true stress of the crystalline phase. The macroscopically effective total Cauchy stress... The ideal theoretical bearing limit of the composite material without microscopic damage was characterized and used as the evolution driving source and reduction reference state for the subsequent dual-domain independent damage model, thereby completing the calculation of the macroscopic nonlinear mechanical degradation behavior of the composite material.

[0035] Step 300: Construct a dual-domain independent damage evolution model; Abandoning the traditional assumption of macroscopic homogeneous damage, and taking crystallinity as the core weight for damage allocation, the formula for calculating the total stress after macroscopic damage is as follows: ; in, This represents the total macroscopic stress after damage. and These are respectively the amorphous phase damage variables and the crystalline phase damage variables; , The homogenization stress component contributed by the amorphous phase at the macroscopic integration point; , This represents the homogenization stress component contributed by the crystalline phase at the macroscopic integration point. The transient total crystallinity is used to characterize the volume fraction of crystalline phases and to perform bi-phase force weighting.

[0036] The macroscopic total damage state variable is D, defined as follows: The macroscopic total damage state variable can effectively characterize the comprehensive stiffness degradation of composite materials at the macroscopic scale. This variable not only serves as a core output index for evaluating the macroscopic failure region and damage degree in finite element post-processing, and provides a comprehensive judgment threshold for macroscopic crack initiation and element failure, but also establishes a cross-scale physical mapping relationship between microscopic dual-domain damage evolution and macroscopic mechanical property degradation, thereby ensuring the physical consistency between numerical simulation and macroscopic mechanical tensile test results.

[0037] Step 301: Damage evolution of amorphous phase under dual thermosensitive mechanisms; The failure behavior of amorphous phases is significantly driven by temperature, so a temperature-related weighting factor is introduced. Achieving a smooth transition from low-temperature brittleness to high-temperature ductility. Amorphous phase damage. Brittle mechanism based on maximum principal stress Ductility mechanism based on cumulative plastic strain The fusion is calculated using the following formula: ; in, For damage variables in amorphous phases; As a temperature-related weighting factor, it approaches zero at low temperatures to activate the brittle mechanism, and approaches one at high temperatures to activate the ductile mechanism. This represents the brittle damage component of the amorphous phase. This represents the ductile damage component of the amorphous phase.

[0038] The calculation formula for the brittle damage mechanism is as follows: ; in, For amorphous phase brittle damage variables; This is the brittle damage evolution rate constant; The maximum principal stress of the true stress tensor of the amorphous phase; This is the threshold value for brittle fracture strength.

[0039] The calculation formula for the ductile damage mechanism is as follows: ; in, For amorphous phase ductility damage variables; This is the ductile damage evolution rate constant; This represents the accumulated equivalent plastic strain of the amorphous phase. This is the strain threshold for the initiation of ductile failure.

[0040] Step 302: Crystal phase-independent damage evolution.

[0041] Crystallographic damage is mainly driven by crystal plane slip or tearing caused by high stress. Based on the equivalent stress evolution of the true phase bias, the calculation formula is as follows: ; in, ; The crystal phase damage evolution rate constant; The deviatoric equivalent stress is the actual stress tensor of the crystal phase. This represents the real-time yield strength of the crystalline phase.

[0042] Step 400: Finite element full-process mapping and tensile failure prediction simulation.

[0043] The constitutive relations and damage evolution information from steps 100 to 300 are written into a custom UMAT (User Material Subroutine) built using the Fortran language, and an automated script is built using the Python language to achieve closed-loop simulation.

[0044] (a) Process simulation field reconstruction: During the cooling stage of process forming, the finite element solver processes the heat transfer field information, and the UMAT subroutine updates the crystallinity field distribution in real time according to the temperature history of each integration point, and calculates the volume change caused by the macroscopic volume shrinkage due to the crystallization phase transformation.

[0045] (b) Inheritance of state variables: After the hot pressing analysis step is completed, the physical state variables of the entire sample reside in the solver memory and are directly inherited as initial conditions into the subsequent service loading analysis step without the need for importing external data or variable mapping operations.

[0046] (c) Non-uniform loading and failure assessment: After applying external loads in the service loading analysis step, the mechanical properties and damage modes at each integration point are strictly controlled by their local crystallinity level. When local stress or strain triggers the damage threshold, the stiffness reduction at that integration point is updated synchronously, and a complete macroscopic stress-strain curve and microscopic failure location prediction are output, thus obtaining a numerical model of composite materials based on continuous phase evolution and dual-domain independent damage.

[0047] The technical solution of the present invention will be described in detail below through a specific embodiment: This embodiment demonstrates the simulation of the entire process of polyether ether ketone (PEEK) composite material specimens from hot pressing to tensile failure through a complete numerical simulation process, in order to verify the accuracy and robustness of the method of the present invention in predicting the nonlinear mechanical response and associated damage evolution behavior of materials under different process crystallization states.

[0048] I. Finite Element Model Construction and Preprocessing: This embodiment combines the finite element analysis software Abaqus with an automated script built using Python to construct a thermo-mechanical coupled three-dimensional finite element model of a tensile specimen conforming to the ISO 527-2 standard. The gauge length is set to 60mm × 10mm × 4mm. To accurately capture the local deformation gradient during stress softening, an eight-node thermo-coupled hexahedral element (C3D8T) is used to mesh the entire model, resulting in a total of 2208 finite element elements. The geometric boundary conditions and mesh settings of the specimen are as follows. Figure 2 As shown.

[0049] At the bottom layer of the solver, the core algorithms proposed in this invention, including the dual crystallization dynamic evolution model, the amorphous phase viscoelastic-plastic constitutive model, the crystalline effective phase simplified constitutive model, and the dual-domain independent damage evolution model, are compiled in a unified manner and linked to the finite element solver in the form of a custom material subroutine (UMAT), thus completing the deep coupling of theoretical constitutive relation information to the physical finite element model.

[0050] II. Implementation of Continuous Coupled Simulation of the Entire "Forming-Service" Process: To achieve the "integrated prediction" goal of this invention, two temperature-displacement coupled analysis steps are sequentially established using a Python script in a single simulation task. At the starting point of the analysis step (InitialStep), the global initial temperature field is set to 400°C (i.e., fully molten state), and the transient total crystallinity is initialized. and the state of no stress throughout the entire domain.

[0051] Phase 1: Simulation of hot pressing molding process; This stage aims to simulate the cooling and solidification process of the sample in a hot press mold. By changing the convective heat transfer boundary conditions, samples with different degrees of crystallinity are prepared: Operating Condition A (Slow Cooling): Set the convective heat transfer coefficient W / (m 2 •K), simulating natural cooling conditions, with a cooling time of 300s.

[0052] Operating Condition B (Quick Cooling): Set the convective heat transfer coefficient W / (m 2 •K), simulating forced air cooling intervention conditions, with a cooling time of 60s.

[0053] At the end of the first analysis step, the physical state variables residing in the finite element solver include at least: the transient total crystallinity c, the crystallinity increment Δc of the newly formed crystalline phase, and the effective phase elastic isocompensatory Cauchy-Green tensor. Temperature shift factor and crystallization shift factor The aforementioned physical state variables will serve as the initial conditions for the continuous physical field, seamlessly inherited into the analysis step of stage two, and used to drive the constitutive stress recursive update, viscoelastic relaxation time spectrum evolution, and dual-domain independent damage threshold determination during the subsequent service loading process.

[0054] The cooling and solidification process in stage one mainly executes steps 100, 200, and (a) and (b) of step 400. During this process, the system inputs the initial melting temperature field, the mold heat transfer boundary conditions, and the material thermodynamic parameters, and outputs the physical state variables at each integration point.

[0055] Phase Two: Service Tensile Loading Simulation; This analysis step builds upon the calculations in Stage 1, requiring no import of external data or variable inheritance. While maintaining the room temperature boundary conditions across the entire specimen domain, all degrees of freedom at one end of the specimen are constrained, and a constant tensile load controlled by displacement is applied to the other end at a tensile rate of 3 mm / min until macroscopic failure occurs.

[0056] The stretching service process in stage two takes the output of stage one as the initial condition and executes steps 202 to 205 and steps 301 to 302 continuously to achieve a continuous mapping of "forming-service".

[0057] III. Results Analysis and Method Validation: Based on visualization post-processing technology, the macroscopic stress-strain curves and state variable cloud maps output from continuous simulation are visualized and compared and analyzed to fully verify the effectiveness of the "process-service" integrated modeling method proposed in this invention in predicting the nonlinear mechanical behavior of semi-crystalline composite materials.

[0058] 3.1 The decisive role of process cooling rate on the crystallinity field: like Figure 3 As shown, the distribution of crystallinity field after the molding stage is completed is illustrated. For condition A, the slower cooling rate provides ample time for polymer molecular chain rearrangement, resulting in a uniform high crystallinity field across the entire domain, with a final total crystallinity of 34.4%. In condition B, the extremely fast cooling rate effectively suppresses crystal growth, and the final total crystallinity is uniformly fixed at 27%. This prediction result is highly consistent with the actual phase transition kinetics of PEEK materials, laying a physical foundation for subsequent mechanical property difference analysis.

[0059] 3.2 The effect of crystallinity differences on macroscopic load-bearing capacity and local stress: like Figure 4As shown, the stress-strain curves of specimens with different crystallinities are compared. The highly crystallinity specimen (condition A, red line) exhibits a higher peak stress (approximately 112 MPa) and initial stiffness. Conversely, the low crystallinity specimen (condition B, blue line) shows a peak stress that drops to approximately 110 MPa, and exhibits a steeper stress softening slope (strain range of 0.1–0.2) after exceeding the peak. This difference in load-bearing capacity stems from the load-bearing sharing mechanism of the internal microstructure. The Mises stress field at the moment of tensile failure is shown in the figure. Figure 5 As shown, the maximum local Mises stress in specimen A under condition A reached 128 MPa, significantly higher than the 51.3 MPa in condition B. This result strongly demonstrates that the well-developed crystalline phase network under high crystallinity is the main framework resisting external loads; while under low crystallinity, due to the lack of crystalline phase support, the overall stress distribution level is lower, and the material is more prone to yielding.

[0060] 3.3 Dual-domain independent damage and macroscopic composite damage: To explain Figure 4 The more severe macroscopic stress softening phenomenon under medium working condition B was analyzed in depth using the dual-domain independent damage output of the present invention to investigate the microscopic failure mechanism.

[0061] Amorphous phase damage evolution: such as Figure 6 As shown, the volume of amorphous phase is extremely high in condition B, which undergoes rapid cooling. Under large deformation and stretching, a large number of amorphous domains are forced to bear intense viscoplastic energy dissipation, leading to rapid deterioration and large-area saturation of amorphous phase damage in the core region of the gauge length. This directly explains the physical reason why the macroscopic curve of condition B falls into severe softening earlier; while the amorphous phase damage in condition A is relatively slight and dispersed.

[0062] Crystal phase damage evolution: such as Figure 7 As shown, in condition A, the crystal network is well-developed, capable of continuously distributing and transferring loads, resulting in widespread but controllable crystal phase damage. However, in condition B, due to the overall low crystallinity, a very small number of crystal structures encounter extremely severe real stress concentration at the moment of failure, leading to a sudden burst of local crystal phase damage peaks. The model of this invention accurately overcomes the problem of tangential stiffness degradation under low crystallinity and robustly captures this stress anomaly.

[0063] Macroscopic Comprehensive Damage Verification: To intuitively assess the impact of independent dual-domain damage on the overall material failure behavior, this embodiment outputs a macroscopic total damage state variable based on volume fraction weighting. For example... Figure 8 As shown, the peak value of the total macroscopic damage in the core area of ​​the gauge length section under condition B is significantly higher than that under condition A, and the failure zone exhibits a highly concentrated trend. This numerical field distribution of macroscopic damage is completely consistent with the source tracing conclusions of the microscopic dual-domain damage, and even more so with… Figure 3 The steep macroscopic stress softening trajectory under medium working condition B forms a perfect closed-loop confirmation.

[0064] This embodiment fully demonstrates that the modeling method and finite element framework proposed in this invention completely break down the physical barriers between complex process thermal histories and macroscopic nonlinear failures of multiphase composite materials. This model can not only accurately characterize macroscopic mechanical responses, but also delve into the material's interior to reveal the competing damage mechanism between amorphous phase viscoplastic dissipation and crystalline phase local stress concentration caused by differences in cooling rates. Without requiring expensive trial-and-error experiments, it provides highly accurate numerical support for process optimization and damage-tolerant design of high-performance composite materials.

[0065] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit them. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the technical solutions of the present invention, and these modifications or equivalent substitutions cannot cause the modified technical solutions to deviate from the spirit and scope of the technical solutions of the present invention.

Claims

1. A composite material modeling method based on continuous phase evolution and dual-domain independent damage, characterized in that, Includes the following steps: Step 100: Construct a dual crystallization kinetic model; based on the temperature history of each node in the finite element analysis, adaptively identify the thermal history stage of the current temperature, and activate the corresponding cold crystallization mechanism or melt crystallization mechanism to solve for the total crystallinity; Step 200: Establish a multiphase continuous evolution mechanical constitutive model; use an effective phase dimensionality reduction strategy to merge the newly formed crystalline phase and the historical crystalline phase in each time increment step into a single effective phase elastic strain tensor; establish a hyperelastic-plastic constitutive model for the crystalline phase, and at the same time establish a viscoelastic-plastic constitutive model for the amorphous phase; obtain the total Cauchy stress of the composite material by weighting the amorphous phase and the crystalline phase according to the volume fraction. Step 300: Establish a dual-domain independent damage evolution model; for the amorphous phase, a temperature-sensitive brittle-ductile dual-mechanism damage model is used to calculate the damage variables of the amorphous phase; for the crystalline phase, an independent damage model based on partial equivalent stress is used to calculate... ; By weighting the damage variables of the amorphous phase and the crystalline phase according to crystallinity, the total macroscopic stress after damage is obtained: Step 400: Compile the dual crystallization kinetic model, the multiphase continuous evolution mechanical constitutive model, and the dual-domain independent damage evolution model into a finite element user-defined material subroutine, construct an automated script to perform continuous simulation prediction of the entire process from hot pressing to service failure, and obtain a numerical model of composite materials based on continuous phase evolution and dual-domain independent damage.

2. The composite material modeling method based on continuous phase evolution and dual-domain independent damage according to claim 1, characterized in that, In step 100, based on the temperature history of each node in the finite element analysis, the current thermal history stage of the temperature is adaptively identified, and the corresponding cold crystallization mechanism or melt crystallization mechanism is activated. Solving the transient crystallinity field controlled by the molding temperature history includes: Step 101: When the polymer matrix is ​​in a rapid cooling and heating process, and the current actual temperature of the polymer matrix... satisfy The cold crystallization mechanism is activated, and the absolute crystallinity of cold crystals is calculated using the following formula: = ； ； in, This represents the relative volume fraction of cold crystallization. The current time; and All are time-integral variables. The nucleation time of the crystal nucleus. This is an intermediate time variable used for integration. Nucleation density as it evolves with supercooling; The crystal growth rate is jointly controlled by the molecular chain segment diffusion term and the nucleation barrier term. This represents the absolute crystallinity of cold crystallization at the end of the current increment step. and These are the glass transition temperature and melting point of the polymer matrix, respectively. Step 102: When the polymer matrix is ​​in a molten state and cooling process, and the temperature... satisfy At this time, the model activates the melting and crystallization mechanism, using the differential form of the Nakamura equation to characterize the rate of change in relative crystallinity. The rate of change in relative crystallinity is integrated over the time increment step to obtain the relative melt crystallinity increment, thereby calculating the absolute crystallinity of the melt crystallization. The calculation formula is as follows: ； = + ； in, This represents the absolute crystallinity increment within the current time increment step. This represents the relative increase in melt crystallinity. This represents the maximum enthalpy of crystallization at the current temperature. The latent heat of crystallization of the polymer matrix under fully crystalline conditions; This represents the absolute crystallinity contributed by melt crystallization at the end of the previous time increment step. This indicates the absolute crystallinity of the melt at the end of the current time increment step; Step 103: Calculate the total crystallinity based on the absolute crystallinity of melt crystallization and the absolute crystallinity of cold crystallization. The calculation formula is as follows: = + ;in The total crystallinity is denoted as .

3. The composite material modeling method based on continuous phase evolution and dual-domain independent damage according to claim 2, characterized in that, In step 200, the effective phase dimensionality reduction strategy is used to fuse the newly formed crystal phase and the historical crystal phase into a single effective phase elastic strain tensor within each time increment step, which includes the following steps: Step 201: Consider the phase transformation volume shrinkage and calculate the corrected Jacobian component of the crystal phase volume. The calculation formula is as follows: ； ； in, The corrected Jacobian component of the crystal phase volume; This is the total deformation gradient tensor at the end of the current increment step; It is a matrix determinant operator; The ratio of the physical density of the amorphous phase to that of the crystalline phase; The density of the amorphous phase; The perfect crystal density of the crystalline phase; Step 202: Based on the principle of equivalent mapping, the elastic isochoric left Cauchy-Green tensor of the old crystal phase recorded in the previous incremental step is extrapolated to the current configuration through the isochoric relative deformation gradient tensor. Then, based on the volume weighting rule, the old and new crystals are merged into a single macroscopic effective phase elastic isochoric left Cauchy-Green tensor. The calculation formula is as follows: ； ； in, To extrapolate the elastic isocompassive left Cauchy-Green tensor of the old crystal phase to the current configuration; The elastic isocompensative left Cauchy-Green tensor of the old crystal phase recorded in the previous increment step; For isochoric relative deformation gradient tensor; for Transpose of; The effective phase elasticity isocompensatory left Cauchy-Green tensor after fusion; It is a second-order unit tensor.

4. The composite material modeling method based on continuous phase evolution and dual-domain independent damage according to claim 3, characterized in that, In step 200, establishing the hyperelastic-plastic constitutive model of the crystalline phase includes: The elastic deformation response of the effective phase of the crystal is described by a compressible Neo-Hookean hyperelastic potential function. The formula for calculating the Cauchy stress of the effective phase of the crystal is as follows: ； in, Cauchy stress is the effective phase of the crystal. This represents the shear modulus of the crystalline phase. For tensors, trace operators; The bulk modulus of the crystalline phase; The true stress of a crystal phase is determined based on its yield function, specifically as follows: If the yield function of the crystal phase is less than or equal to 0, then the Cauchy stress of the effective phase of the crystal is directly taken as the real stress of the crystal phase. If the yield function of the crystal phase is greater than 0, the Cauchy stress of the effective phase of the crystal is projected onto the current yield surface using radial regression mapping to obtain the corrected deviatoric stress of the crystal phase as the true stress of the crystal phase.

5. The composite material modeling method based on continuous phase evolution and dual-domain independent damage according to claim 4, characterized in that, In step 200, the established viscoelastic-plastic constitutive model of the amorphous phase adopts a parallel architecture of generalized Maxwell viscoelastic branches and elastoplastic slip branches. The specific construction process is as follows: For the generalized Maxwell viscoelastic branch, the Neo-Hookean hyperelastic model is used to calculate the instantaneous pseudo-stress, and the viscoelastic stress is solved by using the convolution integral and the normalized relaxation modulus of the Prony series expansion to solve for the Cauchy stress contributed by the viscoelastic branch of the amorphous phase. Calculate the total shift factor, convert the real time increment into an equivalent time increment, and complete the recursive update of the viscoelastic stress for each Maxwell branch. The calculation formula is as follows: ； ； ； ； in, It is the temperature shift factor; , All are WLF model constants; This is the actual temperature of the polymer matrix; For reference temperature; It is the crystallization shift factor; The energy barrier parameter of the crystal; The transient total crystallinity; Boltzmann's constant; Total shift factor; Increment to actual time; This is the equivalent time increment; For the elastoplastic slip branch, the true Cauchy stress contributed by the slip network is solved. A yield function for the amorphous phase slip network is constructed based on the equivalent deviatoric stress. The true stress contributed by the amorphous phase elastoplastic slip branch is determined according to the yield function. The formula for solving the true Cauchy stress contributed by the slip network is as follows: ； in, The actual Cauchy stress tensor contributed to the slip network; The shear modulus of the slip network; The bulk modulus of the slip network; For the elasticity of the sliding network, use the left-inclusive Cauchy-Green tensor; For the volumetric Jacobian component of the elastic deformation of the slip network; If the yield function of the amorphous phase slip network is less than or equal to 0, then the actual stress of the slip network is directly taken as the actual stress contributed by the amorphous phase elastoplastic slip branch. If the yield function of the amorphous phase slip network is greater than 0, the true Cauchy stress tensor contributed by the slip network is projected onto the current yield surface using radial regression mapping to obtain the corrected deviatoric stress as the true stress contributed by the amorphous phase elastoplastic slip branch.

6. The composite material modeling method based on continuous phase evolution and dual-domain independent damage according to claim 5, characterized in that, In step 200, the formula for calculating the total Cauchy stress of the composite material obtained by weighting the amorphous phase and crystalline phase according to their volume fractions is as follows: ； in, The total Cauchy stress of the composite material; Cauchy stress contributes to the viscoelastic branching of amorphous phases; The actual stress contributed to the elastoplastic slip branch of the amorphous phase; This represents the true stress of the crystal phase.

7. The composite material modeling method based on continuous phase evolution and dual-domain independent damage according to claim 6, characterized in that, For amorphous phases, a temperature-sensitive brittle-ductile dual-mechanism damage model is used to calculate the damage variables. The calculation method is as follows: ； in, For damage variables in amorphous phases; Temperature-related weighting factors; This represents the brittle damage component of the amorphous phase. This represents the ductile damage component of the amorphous phase. The calculation formula for the brittle damage mechanism is as follows: ； in, This is the brittle damage evolution rate constant; The maximum principal stress of the true stress tensor of the amorphous phase; This is the threshold value for brittle fracture strength. The calculation formula for the ductile damage mechanism is as follows: ； in, This is the ductile damage evolution rate constant; This represents the accumulated equivalent plastic strain of the amorphous phase. This is the strain threshold for the initiation of ductile failure.

8. The composite material modeling method based on continuous phase evolution and dual-domain independent damage according to claim 7, characterized in that, The crystal phase damage variable was calculated using an independent damage model based on biased equivalent stress. The calculation method is as follows: ； in, For crystal phase damage variables; The crystal phase damage evolution rate constant; The deviatoric equivalent stress is the actual stress tensor of the crystal phase. This represents the real-time yield strength of the crystalline phase.

9. The composite material modeling method based on continuous phase evolution and dual-domain independent damage according to claim 8, characterized in that, The formula for calculating the total stress after macroscopic damage is as follows: ； in, This represents the total macroscopic stress after damage. The homogenization stress component contributed by the amorphous phase at the macroscopic integration point; This represents the homogenization stress component contributed by the crystalline phase at the macroscopic integration point. This represents the transient total crystallinity.

Images